StructuredMeshSimple.cpp File Reference

#include "moab/Core.hpp"
#include "moab/ScdInterface.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CN.hpp"
#include <iostream>
#include <vector>

Include dependency graph for StructuredMeshSimple.cpp:

Go to the source code of this file.

Macros
#define	MYSTREAM(a)    if( !rank ) cout << a << endl

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

Macro Definition Documentation

MYSTREAM

#define MYSTREAM

(

a

)

if( !rank ) cout << a << endl

Examples

StructuredMeshSimple.cpp.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Examples

StructuredMeshSimple.cpp.

Definition at line 39 of file StructuredMeshSimple.cpp.

   40 {
   41     int N = 10, dim = 3;
   42  
   43 #ifdef MOAB_HAVE_MPI
   44     MPI_Init( &argc, &argv );
   45 #endif
   46  
   47     ProgOptions opts;
   48     opts.addOpt< int >( string( "dim,d" ), string( "Dimension of mesh (default=3)" ), &dim );
   49     opts.addOpt< int >( string( ",n" ), string( "Number of elements on a side (default=10)" ), &N );
   50     opts.parseCommandLine( argc, argv );
   51  
   52     // 0. Instantiate MOAB and get the structured mesh interface
   53     Interface* mb = new( std::nothrow ) Core;
   54     if( NULL == mb ) return 1;
   55     ScdInterface* scdiface;
   56     ErrorCode rval = mb->query_interface( scdiface );MB_CHK_ERR( rval );  // Get a ScdInterface object through moab instance
   57  
   58     // 1. Decide what the local parameters of the mesh will be, based on parallel/serial and rank.
   59 #ifdef MOAB_HAVE_MPI
   60     int rank = 0, nprocs = 1;
   61     MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
   62     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
   63     int ilow = rank * N, ihigh = ilow + N;
   64 #else
   65     int rank = 0;
   66     int ilow = 0, ihigh = N;
   67 #endif
   68  
   69     // 2. Create a N^d structured mesh, which includes (N+1)^d vertices and N^d elements.
   70     ScdBox* box;
   71     rval = scdiface->construct_box(
   72         HomCoord( ilow, ( dim > 1 ? 0 : -1 ),
   73                   ( dim > 2 ? 0 : -1 ) ),  // Use in-line logical tests to handle dimensionality
   74         HomCoord( ihigh, ( dim > 1 ? N : -1 ), ( dim > 2 ? N : -1 ) ), NULL,
   75         0,  // NULL coords vector and 0 coords (don't specify coords for now)
   76         box );MB_CHK_ERR( rval );  // box is the structured box object providing the parametric
   77                          // structured mesh interface for this rectangle of elements
   78  
   79     // 3. Get the vertices and elements from moab and check their numbers against (N+1)^d and N^d,
   80     // resp.
   81     Range verts, elems;
   82     rval = mb->get_entities_by_dimension( 0, 0, verts );MB_CHK_ERR( rval );  // First '0' specifies "root set", or entire MOAB instance, second the
   83     // entity dimension being requested
   84     rval = mb->get_entities_by_dimension( 0, dim, elems );MB_CHK_ERR( rval );
   85  
   86 #define MYSTREAM( a ) \
   87     if( !rank ) cout << a << endl
   88  
   89     if( pow( N, dim ) == (int)elems.size() && pow( N + 1, dim ) == (int)verts.size() )
   90     {  // Expected #e and #v are N^d and (N+1)^d, resp.
   91 #ifdef MOAB_HAVE_MPI
   92         MYSTREAM( "Proc 0: " );
   93 #endif
   94         MYSTREAM( "Created " << elems.size() << " " << CN::EntityTypeName( mb->type_from_handle( *elems.begin() ) )
   95                              << " elements and " << verts.size() << " vertices." << endl );
   96     }
   97     else
   98         cout << "Created the wrong number of vertices or hexes!" << endl;
   99  
  100     // 4. Loop over elements in 3 nested loops over i, j, k; for each (i,j,k):
  101     vector< double > coords( 3 * pow( N + 1, dim ) );
  102     vector< EntityHandle > connect;
  103     for( int k = 0; k < ( dim > 2 ? N : 1 ); k++ )
  104     {
  105         for( int j = 0; j < ( dim > 1 ? N : 1 ); j++ )
  106         {
  107             for( int i = 0; i < N - 1; i++ )
  108             {
  109                 // 4a. Get the element corresponding to (i,j,k)
  110                 EntityHandle ehandle = box->get_element( i, j, k );
  111                 if( 0 == ehandle ) return MB_FAILURE;
  112                 // 4b. Get the connectivity of the element
  113                 rval = mb->get_connectivity( &ehandle, 1, connect );MB_CHK_ERR( rval );  // Get the connectivity, in canonical order
  114                 // 4c. Get the coordinates of the vertices comprising that element
  115                 rval = mb->get_coords( &connect[0], connect.size(), &coords[0] );MB_CHK_ERR( rval );  // Get the coordinates of those vertices
  116             }
  117         }
  118     }
  119  
  120     // 5. Release the structured mesh interface and destroy the MOAB instance
  121     mb->release_interface( scdiface );  // Tell MOAB we're done with the ScdInterface
  122  
  123 #ifdef MOAB_HAVE_MPI
  124     MPI_Finalize();
  125 #endif
  126  
  127     return 0;
  128 }

References ProgOptions::addOpt(), moab::Range::begin(), moab::ScdInterface::construct_box(), dim, moab::CN::EntityTypeName(), ErrorCode, moab::Core::get_connectivity(), moab::Core::get_coords(), moab::ScdBox::get_element(), moab::Core::get_entities_by_dimension(), mb, MB_CHK_ERR, MYSTREAM, ProgOptions::parseCommandLine(), moab::Interface::query_interface(), moab::Interface::release_interface(), moab::Range::size(), and moab::Core::type_from_handle().