 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /** @example StructuredMeshSimple.cpp
2 * \brief Show creation and query of structured mesh, serial or parallel, through MOAB's structured
3 * mesh interface. This is an example showing creation and query of a 3D structured mesh. In
4 * serial, a single N*N*N block of elements is created; in parallel, each proc gets an N*N*N block,
5 * with blocks arranged in a 1d column, sharing vertices and faces at their interfaces (proc 0 has
6 * no left neighbor and proc P-1 no right neighbor). Each square block of hex elements is then
7 * referenced by its ijk parameterization. 1D and 2D examples could be made simply by changing the
8 * dimension parameter passed into the MOAB functions. \n
9 *
10 * <b>This example </b>:
11 * -# Instantiate MOAB and get the structured mesh interface
12 * -# Decide what the local parameters of the mesh will be, based on parallel/serial and rank.
13 * -# Create a N^d structured mesh, which includes (N+1)^d vertices and N^d elements.
14 * -# Get the vertices and elements from moab and check their numbers against (N+1)^d and N^d,
15 * resp.
16 * -# Loop over elements in d nested loops over i, j, k; for each (i,j,k):
17 * -# Get the element corresponding to (i,j,k)
18 * -# Get the connectivity of the element
19 * -# Get the coordinates of the vertices comprising that element
20 * -# Release the structured mesh interface and destroy the MOAB instance
21 *
22 * <b> To run: </b> ./StructuredMeshSimple [d [N] ] \n
23 * (default values so can run w/ no user interaction)
24 */
25
26 #include "moab/Core.hpp"
27 #include "moab/ScdInterface.hpp"
28 #include "moab/ProgOptions.hpp"
29 #include "moab/CN.hpp"
30 #ifdef MOAB_HAVE_MPI
31 #include "moab_mpi.h"
32 #endif
33 #include <iostream>
34 #include <vector>
35
36 using namespace moab;
37 using namespace std;
38
39 int main( int argc, char** argv )
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 }
1.9.1.