MOAB: Mesh Oriented datABase  (version 5.3.1)
StructuredMeshSimple.cpp
Go to the documentation of this file.
00001 /** @example StructuredMeshSimple.cpp
00002  * \brief Show creation and query of structured mesh, serial or parallel, through MOAB's structured
00003  * mesh interface. This is an example showing creation and query of a 3D structured mesh.  In
00004  * serial, a single N*N*N block of elements is created; in parallel, each proc gets an N*N*N block,
00005  * with blocks arranged in a 1d column, sharing vertices and faces at their interfaces (proc 0 has
00006  * no left neighbor and proc P-1 no right neighbor). Each square block of hex elements is then
00007  * referenced by its ijk parameterization. 1D and 2D examples could be made simply by changing the
00008  * dimension parameter passed into the MOAB functions. \n
00009  *
00010  * <b>This example </b>:
00011  *    -# Instantiate MOAB and get the structured mesh interface
00012  *    -# Decide what the local parameters of the mesh will be, based on parallel/serial and rank.
00013  *    -# Create a N^d structured mesh, which includes (N+1)^d vertices and N^d elements.
00014  *    -# Get the vertices and elements from moab and check their numbers against (N+1)^d and N^d,
00015  * resp.
00016  *    -# Loop over elements in d nested loops over i, j, k; for each (i,j,k):
00017  *      -# Get the element corresponding to (i,j,k)
00018  *      -# Get the connectivity of the element
00019  *      -# Get the coordinates of the vertices comprising that element
00020  *    -# Release the structured mesh interface and destroy the MOAB instance
00021  *
00022  * <b> To run: </b> ./StructuredMeshSimple [d [N] ] \n
00023  * (default values so can run w/ no user interaction)
00024  */
00025 
00026 #include "moab/Core.hpp"
00027 #include "moab/ScdInterface.hpp"
00028 #include "moab/ProgOptions.hpp"
00029 #include "moab/CN.hpp"
00030 #ifdef MOAB_HAVE_MPI
00031 #include "moab_mpi.h"
00032 #endif
00033 #include <iostream>
00034 #include <vector>
00035 
00036 using namespace moab;
00037 using namespace std;
00038 
00039 int main( int argc, char** argv )
00040 {
00041     int N = 10, dim = 3;
00042 
00043 #ifdef MOAB_HAVE_MPI
00044     MPI_Init( &argc, &argv );
00045 #endif
00046 
00047     ProgOptions opts;
00048     opts.addOpt< int >( string( "dim,d" ), string( "Dimension of mesh (default=3)" ), &dim );
00049     opts.addOpt< int >( string( ",n" ), string( "Number of elements on a side (default=10)" ), &N );
00050     opts.parseCommandLine( argc, argv );
00051 
00052     // 0. Instantiate MOAB and get the structured mesh interface
00053     Interface* mb = new( std::nothrow ) Core;
00054     if( NULL == mb ) return 1;
00055     ScdInterface* scdiface;
00056     ErrorCode rval = mb->query_interface( scdiface );MB_CHK_ERR( rval );  // Get a ScdInterface object through moab instance
00057 
00058     // 1. Decide what the local parameters of the mesh will be, based on parallel/serial and rank.
00059 #ifdef MOAB_HAVE_MPI
00060     int rank = 0, nprocs = 1;
00061     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
00062     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00063     int ilow = rank * N, ihigh = ilow + N;
00064 #else
00065     int rank = 0;
00066     int ilow = 0, ihigh = N;
00067 #endif
00068 
00069     // 2. Create a N^d structured mesh, which includes (N+1)^d vertices and N^d elements.
00070     ScdBox* box;
00071     rval = scdiface->construct_box(
00072         HomCoord( ilow, ( dim > 1 ? 0 : -1 ),
00073                   ( dim > 2 ? 0 : -1 ) ),  // Use in-line logical tests to handle dimensionality
00074         HomCoord( ihigh, ( dim > 1 ? N : -1 ), ( dim > 2 ? N : -1 ) ), NULL,
00075         0,  // NULL coords vector and 0 coords (don't specify coords for now)
00076         box );MB_CHK_ERR( rval );  // box is the structured box object providing the parametric
00077                          // structured mesh interface for this rectangle of elements
00078 
00079     // 3. Get the vertices and elements from moab and check their numbers against (N+1)^d and N^d,
00080     // resp.
00081     Range verts, elems;
00082     rval = mb->get_entities_by_dimension( 0, 0, verts );MB_CHK_ERR( rval );  // First '0' specifies "root set", or entire MOAB instance, second the
00083     // entity dimension being requested
00084     rval = mb->get_entities_by_dimension( 0, dim, elems );MB_CHK_ERR( rval );
00085 
00086 #define MYSTREAM( a ) \
00087     if( !rank ) cout << ( a ) << endl
00088 
00089     if( pow( N, dim ) == (int)elems.size() && pow( N + 1, dim ) == (int)verts.size() )
00090     {  // Expected #e and #v are N^d and (N+1)^d, resp.
00091 #ifdef MOAB_HAVE_MPI
00092         MYSTREAM( "Proc 0: " );
00093 #endif
00094         MYSTREAM( "Created " << elems.size() << " " << CN::EntityTypeName( mb->type_from_handle( *elems.begin() ) )
00095                              << " elements and " << verts.size() << " vertices." << endl );
00096     }
00097     else
00098         cout << "Created the wrong number of vertices or hexes!" << endl;
00099 
00100     // 4. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k):
00101     vector< double > coords( 3 * pow( N + 1, dim ) );
00102     vector< EntityHandle > connect;
00103     for( int k = 0; k < ( dim > 2 ? N : 1 ); k++ )
00104     {
00105         for( int j = 0; j < ( dim > 1 ? N : 1 ); j++ )
00106         {
00107             for( int i = 0; i < N - 1; i++ )
00108             {
00109                 // 4a. Get the element corresponding to (i,j,k)
00110                 EntityHandle ehandle = box->get_element( i, j, k );
00111                 if( 0 == ehandle ) return MB_FAILURE;
00112                 // 4b. Get the connectivity of the element
00113                 rval = mb->get_connectivity( &ehandle, 1, connect );MB_CHK_ERR( rval );  // Get the connectivity, in canonical order
00114                 // 4c. Get the coordinates of the vertices comprising that element
00115                 rval = mb->get_coords( &connect[0], connect.size(), &coords[0] );MB_CHK_ERR( rval );  // Get the coordinates of those vertices
00116             }
00117         }
00118     }
00119 
00120     // 5. Release the structured mesh interface and destroy the MOAB instance
00121     mb->release_interface( scdiface );  // Tell MOAB we're done with the ScdInterface
00122 
00123 #ifdef MOAB_HAVE_MPI
00124     MPI_Finalize();
00125 #endif
00126 
00127     return 0;
00128 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines