StructuredMeshSimple.cpp File Reference

Example demonstrating creation and query of structured meshes in MOAB. More...

#include "moab/Core.hpp"
#include "moab/ScdInterface.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CN.hpp"
#include <array>
#include <cmath>
#include <iostream>
#include <memory>
#include <stdexcept>
#include <vector>

Include dependency graph for StructuredMeshSimple.cpp:

Go to the source code of this file.

Namespaces
	anonymous_namespace{StructuredMeshSimple.cpp}

Functions
int	main (int argc, char **argv)

Variables
constexpr int	anonymous_namespace{StructuredMeshSimple.cpp}::DEFAULT_DIMENSION = 3
constexpr int	anonymous_namespace{StructuredMeshSimple.cpp}::DEFAULT_ELEMENTS_PER_SIDE = 10

Detailed Description

Example demonstrating creation and query of structured meshes in MOAB.

This example shows how to:

• Create structured meshes in 1D, 2D, or 3D
• Work with structured mesh interfaces (ScdInterface)
• Handle both serial and parallel structured mesh creation
• Query mesh entities and their connectivity
• Access element coordinates and connectivity
• Use parametric indexing for structured meshes

In serial mode, a single N*N*N block of elements is created. In parallel mode, each processor gets an N*N*N block arranged in a 1D column, sharing vertices and faces at interfaces.

Author

MOAB Development Team

Date

2024

Show creation and query of structured mesh, serial or parallel, through MOAB's structured mesh interface. This is an example showing creation and query of a 3D structured mesh. In serial, a single N*N*N block of elements is created; in parallel, each proc gets an N*N*N block, with blocks arranged in a 1d column, sharing vertices and faces at their interfaces (proc 0 has no left neighbor and proc P-1 no right neighbor). Each square block of hex elements is then referenced by its ijk parameterization. 1D and 2D examples could be made simply by changing the dimension parameter passed into the MOAB functions.
This example :

1. Instantiate MOAB and get the structured mesh interface
2. Decide what the local parameters of the mesh will be, based on parallel/serial and rank.
3. Create a N^d structured mesh, which includes (N+1)^d vertices and N^d elements.
4. Get the vertices and elements from moab and check their numbers against (N+1)^d and N^d, resp.
5. Loop over elements in d nested loops over i, j, k; for each (i,j,k):
   1. Get the element corresponding to (i,j,k)
   2. Get the connectivity of the element
   3. Get the coordinates of the vertices comprising that element
6. Release the structured mesh interface and destroy the MOAB instance

To run: ./StructuredMeshSimple [d [N] ]
(default values so can run w/ no user interaction)

Parameters

argc	Number of command line arguments
argv	Command line arguments array

Returns

0 on success, 1 on failure

Definition in file StructuredMeshSimple.cpp.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 81 of file StructuredMeshSimple.cpp.

{
    // Initialize MPI if available
#ifdef MOAB_HAVE_MPI
    if( MPI_Init( &argc, &argv ) != MPI_SUCCESS )
    {
        std::cerr << "MPI_Init failed" << std::endl;
        return 1;
    }
 
    // Ensure MPI is properly finalized when we exit
    struct MPIFinalizer
    {
        void operator()( int* ) const
        {
            MPI_Finalize();
        }
    };
    std::unique_ptr< int, MPIFinalizer > mpi_guard( nullptr );
#endif
 
    // Parse command line options
    int dimension         = DEFAULT_DIMENSION;
    int elements_per_side = DEFAULT_ELEMENTS_PER_SIDE;
 
    try
    {
        ProgOptions opts;
        opts.addOpt< int >( "dim,d", "Dimension of mesh (default=3)", &dimension );
        opts.addOpt< int >( ",n", "Number of elements on a side (default=10)", &elements_per_side );
        opts.parseCommandLine( argc, argv );
 
        // Validate input
        if( dimension < 1 || dimension > 3 )
        {
            throw std::invalid_argument( "Dimension must be 1, 2, or 3" );
        }
        if( elements_per_side < 1 )
        {
            throw std::invalid_argument( "Number of elements must be positive" );
        }
    }
    catch( const std::exception& e )
    {
        std::cerr << "Error parsing command line: " << e.what() << std::endl;
        return 1;
    }
 
    // Initialize MOAB
    std::unique_ptr< Core > moab_instance = std::make_unique< Core >();
    if( !moab_instance )
    {
        std::cerr << "Failed to create MOAB instance" << std::endl;
        return 1;
    }
 
    // Get the structured mesh interface
    ScdInterface* scd_interface = nullptr;
    MB_CHK_SET_ERR( moab_instance->query_interface( scd_interface ), "Failed to get ScdInterface" );
    if( !scd_interface )
    {
        std::cerr << "Failed to get ScdInterface: null pointer returned" << std::endl;
        return 1;
    }
 
    // 1. Determine local parameters based on parallel/serial and rank
    int rank  = 0;
    int ilow  = 0;
    int ihigh = elements_per_side;
 
#ifdef MOAB_HAVE_MPI
    int nprocs = 1;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
 
    // Calculate local range for this process
    const int elements_per_proc = elements_per_side / nprocs;
    const int remainder         = elements_per_side % nprocs;
 
    // Distribute elements among processes
    ilow  = rank * elements_per_proc + std::min( rank, remainder );
    ihigh = ilow + elements_per_proc + ( rank < remainder ? 1 : 0 );
#endif
 
    std::string write_options = "";
    if( nprocs > 1 ) write_options = "PARALLEL=WRITE_PART";
 
    // 2. Create a N^d structured mesh, which includes (N+1)^d vertices and N^d elements.
    ScdBox* box = nullptr;
 
    // Calculate box bounds based on dimension
    const moab::HomCoord low( ilow, ( dimension > 1 ) ? 0 : -1, ( dimension > 2 ) ? 0 : -1 );
 
    const moab::HomCoord high( ihigh, ( dimension > 1 ) ? elements_per_side : -1,
                               ( dimension > 2 ) ? elements_per_side : -1 );
 
    // Create the structured mesh box
    MB_CHK_SET_ERR( scd_interface->construct_box( low, high,
                                                  nullptr,  // No coordinates array
                                                  0,        // No coordinates
                                                  box,      // Output parameter
                                                  nullptr,  // Periodicity
                                                  nullptr,  // Parallel data
                                                  true,     // No coordinates
                                                  0         // Resolve dimensionality
                                                  ),
                    "Failed to construct structured box" );
 
    if( !box )
    {
        std::cerr << "Failed to construct structured box: null box returned" << std::endl;
        return 1;
    }
 
    // 3. Get the vertices and elements from moab and check their numbers against (N+1)^d and N^d
    Range vertices, elements;
 
    // Get all vertices (dimension = 0) and elements (dimension = mesh dimension)
    MB_CHK_SET_ERR( moab_instance->get_entities_by_dimension( 0, 0, vertices ), "Failed to get vertices" );
    MB_CHK_SET_ERR( moab_instance->get_entities_by_dimension( 0, dimension, elements ), "Failed to get elements" );
 
#ifdef MOAB_HAVE_MPI
    std::size_t local_entities_size[2] = { vertices.size(), elements.size() }, global_entities_size[2] = { 0, 0 };
    MPI_Allreduce( local_entities_size, global_entities_size, 2, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD );
#endif
 
    // Calculate expected number of vertices and elements
    const std::size_t expected_vertices = static_cast< std::size_t >( std::pow( elements_per_side + 1, dimension ) );
    const std::size_t expected_elements = static_cast< std::size_t >( std::pow( elements_per_side, dimension ) );
 
    // Only print from rank 0 in parallel
    auto print_message = [rank]( const auto& message ) {
        if( rank == 0 )
        {
            std::cout << message << std::endl;
        }
    };
 
    // Verify the number of created elements and vertices
    if( expected_elements == global_entities_size[1] )
    {
        const auto element_type = moab_instance->type_from_handle( *elements.begin() );
        std::ostringstream msg;
        msg << "Created " << global_entities_size[1] << " " << moab::CN::EntityTypeName( element_type )
            << " elements and " << global_entities_size[0] << " vertices.";
        print_message( msg.str() );
    }
    else
    {
        std::ostringstream err_msg;
        err_msg << "Error: Expected " << expected_elements << " elements and " << expected_vertices
                << " vertices, but got " << global_entities_size[1] << " elements and " << global_entities_size[0]
                << " vertices.";
        print_message( err_msg.str() );
        return 1;
    }
 
    // 4. Loop over elements in nested loops over i, j, k based on dimension
    // const int i_max = elements_per_side - 1;  // 0-based indexing
    const int j_max = ( dimension > 1 ) ? elements_per_side - 1 : 0;
    const int k_max = ( dimension > 2 ) ? elements_per_side - 1 : 0;
 
    // Pre-allocate vectors to avoid reallocation
 
    std::vector< double > coordinates;
    for( int k = 0; k <= k_max; ++k )
    {
        for( int j = 0; j <= j_max; ++j )
        {
            for( int i = ilow; i < ihigh; ++i )
            {
                // 4a. Get the element corresponding to (i,j,k)
                const EntityHandle element = box->get_element( i, j, k );
                if( 0 == element )
                {
                    std::cerr << "Failed to get element at (" << i << ", " << j << ", " << k << ")" << std::endl;
                    // return 1;
                    continue;
                }
 
                // 4b. Get the connectivity of the element
                std::vector< EntityHandle > connectivity;
                MB_CHK_SET_ERR( moab_instance->get_connectivity( &element, 1, connectivity ),
                                "Failed to get connectivity for element at (" << i << ", " << j << ", " << k
                                                                              << ") with handle " << element );
 
                // 4c. Get the coordinates of the vertices comprising that element
                coordinates.resize( 3 * connectivity.size() );
                MB_CHK_SET_ERR( moab_instance->get_coords( connectivity.data(), connectivity.size(),
                                                           coordinates.data() ),
                                "Failed to get coordinates for element at (" << i << ", " << j << ", " << k
                                                                             << ") with handle " << element );
            }
        }
    }
 
    // 5. Write the mesh file to disk
    MB_CHK_SET_ERR( moab_instance->write_file( "structured_mesh.h5m", "h5m", write_options.c_str() ),
                    "Failed to write mesh file" );
 
    // 6. Clean up
    MB_CHK_SET_ERR( moab_instance->release_interface( scd_interface ), "Failed to release interface" );
 
    // MPI_Finalize is handled by the mpi_guard destructor if MPI is enabled
    return 0;
}

References ProgOptions::addOpt(), moab::Range::begin(), moab::ScdInterface::construct_box(), anonymous_namespace{StructuredMeshSimple.cpp}::DEFAULT_DIMENSION, anonymous_namespace{StructuredMeshSimple.cpp}::DEFAULT_ELEMENTS_PER_SIDE, moab::CN::EntityTypeName(), moab::ScdBox::get_element(), MB_CHK_SET_ERR, ProgOptions::parseCommandLine(), and moab::Range::size().