MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 #include "moab/MOABConfig.h" 00002 00003 #ifdef MOAB_HAVE_MPI 00004 #include "moab_mpi.h" 00005 #endif 00006 00007 #include "moab/iMOAB.h" 00008 #include "TestUtil.hpp" 00009 00010 // for malloc, free: 00011 #include <cstdlib> 00012 #include <cstdio> 00013 #include <cstring> 00014 00015 #define CHECKRC( rc, message ) \ 00016 if( 0 != ( rc ) ) \ 00017 { \ 00018 printf( "%s", message ); \ 00019 return 1; \ 00020 } 00021 00022 int main( int argc, char* argv[] ) 00023 { 00024 int nprocs = 1, rank = 0; 00025 std::string filen; 00026 #ifdef MOAB_HAVE_MPI 00027 MPI_Init( &argc, &argv ); 00028 00029 MPI_Comm comm = MPI_COMM_WORLD; 00030 MPI_Comm_size( comm, &nprocs ); 00031 MPI_Comm_rank( comm, &rank ); 00032 #endif 00033 00034 #ifdef MOAB_HAVE_HDF5 00035 filen = TestDir + "/io/p8ex1.h5m"; 00036 #endif 00037 00038 if( argc > 1 ) filen = std::string( argv[1] ); 00039 00040 if( !filen.size() ) return 1; 00041 /* 00042 * MOAB needs to be initialized; A MOAB instance will be created, and will be used by each 00043 * application in this framework. There is no parallel context yet. 00044 */ 00045 ErrCode rc = iMOAB_Initialize( argc, argv ); 00046 00047 CHECKRC( rc, "failed to initialize MOAB" ); 00048 int num_global_vertices = 0, num_global_elements = 0, num_dimension = 0, num_parts = 0; 00049 /* 00050 * Header information is cheap to retrieve from an hdf5 file; it can be called by each process 00051 * or can be called by one master task and broadcasted. The hdf5 file has to be in MOAB native 00052 * format, and has to be partitioned in advance. There should be at least as many partitions in 00053 * the file as number of tasks in the communicator. 00054 */ 00055 rc = iMOAB_ReadHeaderInfo( filen.c_str(), &num_global_vertices, &num_global_elements, &num_dimension, &num_parts, 00056 (int)( filen.size() ) ); 00057 00058 CHECKRC( rc, "failed to read header info" ); 00059 00060 if( 0 == rank ) 00061 { 00062 printf( "file %s has %d vertices, %d elements, %d parts in partition\n", filen.c_str(), num_global_vertices, 00063 num_global_elements, num_parts ); 00064 } 00065 int appID, appDupID; 00066 iMOAB_AppID pid = &appID; 00067 iMOAB_AppID pidDup = &appDupID; 00068 /* 00069 * Each application has to be registered once. A mesh set and a parallel communicator will be 00070 * associated with each application. A unique application id will be returned, and will be used 00071 * for all future mesh operations/queries. 00072 */ 00073 int compid = 11, compdupid = 12; 00074 rc = iMOAB_RegisterApplication( "MBAPP", 00075 #ifdef MOAB_HAVE_MPI 00076 &comm, 00077 #endif 00078 &compid, pid ); 00079 CHECKRC( rc, "failed to register application" ); 00080 00081 rc = iMOAB_RegisterApplication( "MBDUPAPP", 00082 #ifdef MOAB_HAVE_MPI 00083 &comm, 00084 #endif 00085 &compdupid, pidDup ); 00086 CHECKRC( rc, "failed to register duplicate application" ); 00087 00088 #ifdef MOAB_HAVE_MPI 00089 const char* read_opts = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS"; 00090 int num_ghost_layers[1] = { 1 }; 00091 #else 00092 const char* read_opts = ""; 00093 int* num_ghost_layers = nullptr; 00094 #endif 00095 00096 /* 00097 * Loading the mesh is a parallel IO operation. Ghost layers can be exchanged too, and default 00098 * MOAB sets are augmented with ghost elements. By convention, blocks correspond to MATERIAL_SET 00099 * sets, side sets with NEUMANN_SET sets, node sets with DIRICHLET_SET sets. Each element and 00100 * vertex entity should have a GLOBAL ID tag in the file, which will be available for visible 00101 * entities 00102 */ 00103 rc = iMOAB_LoadMesh( pid, filen.c_str(), read_opts, num_ghost_layers, (int)( filen.size() ), strlen( read_opts ) ); 00104 CHECKRC( rc, "failed to load mesh" ); 00105 00106 rc = iMOAB_LoadMesh( pidDup, filen.c_str(), read_opts, num_ghost_layers, (int)( filen.size() ), 00107 strlen( read_opts ) ); 00108 CHECKRC( rc, "failed to load mesh" ); 00109 00110 rc = iMOAB_SetGlobalInfo( pid, &num_global_vertices, &num_global_elements ); 00111 CHECKRC( rc, "failed to set global info" ); 00112 00113 int nverts[3], nelem[3], nblocks[3], nsbc[3], ndbc[3]; 00114 /* 00115 * Each process in the communicator will have access to a local mesh instance, which will 00116 * contain the original cells in the local partition and ghost entities. Number of vertices, 00117 * primary cells, visible blocks, number of sidesets and nodesets boundary conditions will be 00118 * returned in size 3 arrays, for local, ghost and total numbers. 00119 */ 00120 rc = iMOAB_GetMeshInfo( pid, nverts, nelem, nblocks, nsbc, ndbc ); 00121 CHECKRC( rc, "failed to get mesh info" ); 00122 00123 iMOAB_GlobalID* vGlobalID = (iMOAB_GlobalID*)malloc( nverts[2] * sizeof( iMOAB_GlobalID ) ); 00124 /* 00125 * All visible vertices have a unique global ID. All arrays are allocated by client. 00126 * Currently, global ID is a 32 bit integer, total number of visible vertices was returned by 00127 * iMOAB_GetMeshInfo 00128 */ 00129 rc = iMOAB_GetVertexID( pid, &nverts[2], vGlobalID ); 00130 CHECKRC( rc, "failed to get vertex id info" ); 00131 00132 int* vranks = (int*)malloc( nverts[2] * sizeof( int ) ); 00133 /* 00134 * In MOAB terminology, a vertex can be owned, shared or ghost. An owned vertex will be 00135 * owned by the current process. A shared vertex appears at the interface between partitions, it 00136 * can be owned by another process, while a ghost vertex is owned by another process always. All 00137 * vertex arrays have a natural local order, in which the ghost vertices are at the end of the 00138 * arrays. The local index of the vertex will vary from 0 to number of visible vertices (in this 00139 * example nverts[2]-1) 00140 */ 00141 rc = iMOAB_GetVertexOwnership( pid, &nverts[2], vranks ); 00142 CHECKRC( rc, "failed to get vertex ranks" ); 00143 00144 double* coords = (double*)malloc( 3 * nverts[2] * sizeof( double ) ); 00145 int size_coords = 3 * nverts[2]; 00146 /* 00147 * Coordinates are returned interleaved, for all visible vertices. Size is dimension times 00148 * number of visible vertices 00149 */ 00150 rc = iMOAB_GetVisibleVerticesCoordinates( pid, &size_coords, coords ); 00151 CHECKRC( rc, "failed to get coordinates" ); 00152 00153 /* 00154 * Block IDs correspond to MATERIAL_SET tag values in hdf5 file 00155 * They also correspond to block IDs in a Cubit model. 00156 * Local visible blocks might contain ghost cells, owned by other 00157 * processes. 00158 * Inside MOAB, a block will correspond to a meshset with a MATERIAL_SET tag, 00159 * with value the block ID. 00160 * A MATERIAL_SET tag in moab if SPARSE and type INTEGER. 00161 */ 00162 iMOAB_GlobalID* gbIDs = (iMOAB_GlobalID*)malloc( nblocks[2] * sizeof( iMOAB_GlobalID ) ); 00163 rc = iMOAB_GetBlockID( pid, &nblocks[2], gbIDs ); 00164 CHECKRC( rc, "failed to get block info" ); 00165 00166 /* 00167 * The 2 tags used in this example exist in the file, already. 00168 * If a user needs a new tag, it can be defined with the same call as this one 00169 * this method, iMOAB_DefineTagStorage, will return a local index for the tag. 00170 * The name of the tag is case sensitive. 00171 * This method is collective. 00172 */ 00173 int tagIndex[2]; 00174 int entTypes[2] = { 0, 1 }; /* first is on vertex, second is on elements; */ 00175 int tagTypes[2] = { DENSE_INTEGER, DENSE_DOUBLE }; 00176 int num_components = 1; 00177 rc = iMOAB_DefineTagStorage( pid, "INTFIELD", &tagTypes[0], &num_components, &tagIndex[0], strlen( "INTFIELD" ) ); 00178 CHECKRC( rc, "failed to get tag INTFIELD " ); 00179 00180 rc = iMOAB_DefineTagStorage( pid, "DFIELD", &tagTypes[1], &num_components, &tagIndex[1], strlen( "DFIELD" ) ); 00181 CHECKRC( rc, "failed to get tag DFIELD " ); 00182 00183 // synchronize one of the tags only, just to see what happens 00184 int num_tags_to_sync = 1; 00185 /* 00186 * tags are not exchanged for ghost entities by default. The user has to explicitly 00187 * synchronize them with a call like this one. The tags to synchronize are referenced by the 00188 * index returned at their "definition" with iMOAB_DefineTagStorage method. 00189 */ 00190 rc = iMOAB_SynchronizeTags( pid, &num_tags_to_sync, &tagIndex[0], &tagTypes[0] ); 00191 CHECKRC( rc, "failed to sync tag INTFIELD " ); 00192 00193 for( int irank = 0; irank < nprocs; irank++ ) 00194 { 00195 if( irank == rank ) 00196 { 00197 /* printf some of the block info */ 00198 printf( "on rank %d, there are \n" 00199 " %3d visible vertices of which %3d local %3d ghost \n" 00200 " %3d visible elements of which %3d owned %3d ghost \n" 00201 " %3d visible blocks\n" 00202 " %3d visible neumann BCs\n" 00203 " %3d visible dirichlet BCs\n", 00204 rank, nverts[2], nverts[0], nverts[1], nelem[2], nelem[0], nelem[1], nblocks[2], nsbc[2], ndbc[2] ); 00205 /* printf some of the vertex id infos */ 00206 int numToPrint = nverts[2]; 00207 printf( "on rank %d vertex info:\n", rank ); 00208 for( int i = 0; i < numToPrint; i++ ) 00209 printf( " vertex local id: %3d, rank ID:%d global ID: %3d coords: %g, %g, %g\n", i, vranks[i], 00210 vGlobalID[i], coords[3 * i], coords[3 * i + 1], coords[3 * i + 2] ); 00211 00212 iMOAB_GlobalID* element_global_IDs = (iMOAB_GlobalID*)malloc( nelem[2] * sizeof( iMOAB_GlobalID ) ); 00213 iMOAB_GlobalID* block_IDs = (iMOAB_GlobalID*)malloc( nelem[2] * sizeof( iMOAB_GlobalID ) ); 00214 int* ranks = (int*)malloc( nelem[2] * sizeof( int ) ); 00215 /* 00216 * visible elements info is available for all visible elements, with this call. Similar 00217 * information can be retrieved by looping over visible blocks 00218 */ 00219 rc = iMOAB_GetVisibleElementsInfo( pid, &nelem[2], element_global_IDs, ranks, block_IDs ); 00220 CHECKRC( rc, "failed to get all elem info" ); 00221 for( int i = 0; i < nelem[2]; i++ ) 00222 printf( " element local id: %3d, global ID: %3d rank:%d block ID: %2d \n", i, element_global_IDs[i], 00223 ranks[i], block_IDs[i] ); 00224 free( element_global_IDs ); 00225 free( ranks ); 00226 free( block_IDs ); 00227 00228 // get first element connectivity 00229 int conn[27]; 00230 int nv = 27; 00231 int eindex = 0; 00232 /* 00233 * element connectivity can be retrieved for a single element. Local vertex indices are 00234 * returned. 00235 */ 00236 rc = iMOAB_GetElementConnectivity( pid, &eindex, &nv, conn ); 00237 CHECKRC( rc, "failed to get first element connectivity" ); 00238 printf( " conn for first element: \n" ); 00239 for( int i = 0; i < nv; i++ ) 00240 printf( " %3d", conn[i] ); 00241 printf( "\n" ); 00242 00243 // test neighbors 00244 int local_index = 0; // test element with local index 0 00245 int num_adjacent_elements = 10; // we can have maximum 6 actually 00246 iMOAB_LocalID adjacent_element_IDs[10]; 00247 /* 00248 * Neighbor elements that share a face with current element can be determined with this 00249 * call Elements are identified by their local index 00250 */ 00251 rc = iMOAB_GetNeighborElements( pid, &local_index, &num_adjacent_elements, adjacent_element_IDs ); 00252 CHECKRC( rc, "failed to get first element neighbors" ); 00253 printf( " neighbors for first element:\n" ); 00254 for( int i = 0; i < num_adjacent_elements; i++ ) 00255 { 00256 printf( " %4d", adjacent_element_IDs[i] ); 00257 } 00258 printf( "\n" ); 00259 00260 for( int i = 0; i < nblocks[2]; i++ ) 00261 { 00262 printf( " block index: %3d, block ID: %3d \n", i, gbIDs[i] ); 00263 int vertices_per_element, num_elements_in_block; 00264 /* 00265 * blocks should have the same type of primary elements. Number of elements in block 00266 * and number of vertices per element are found with this call. The method refers 00267 * only to the visible elements on the process. The block is identified with its 00268 * block ID. Should it be identified with its local index in the list of visible 00269 * blocks? 00270 */ 00271 rc = iMOAB_GetBlockInfo( pid, &gbIDs[i], &vertices_per_element, &num_elements_in_block ); 00272 CHECKRC( rc, "failed to elem block info" ); 00273 printf( " has %4d elements with %d vertices per element\n", num_elements_in_block, 00274 vertices_per_element ); 00275 int size_conn = num_elements_in_block * vertices_per_element; 00276 iMOAB_LocalID* element_connectivity = (iMOAB_LocalID*)malloc( sizeof( iMOAB_LocalID ) * size_conn ); 00277 /* 00278 * Connectivity for all elements in the block can be determined with this method. 00279 * Vertices are identified by their local index (from 0 to number of visible 00280 * vertices -1) 00281 */ 00282 rc = iMOAB_GetBlockElementConnectivities( pid, &gbIDs[i], &size_conn, element_connectivity ); 00283 CHECKRC( rc, "failed to get block elem connectivity" ); 00284 int* element_ownership = (int*)malloc( sizeof( int ) * num_elements_in_block ); 00285 00286 /* 00287 * Each element is owned by a particular process. Ownership can be returned for all 00288 * elements in a block with this call. 00289 */ 00290 rc = iMOAB_GetElementOwnership( pid, &gbIDs[i], &num_elements_in_block, element_ownership ); 00291 CHECKRC( rc, "failed to get block elem ownership" ); 00292 iMOAB_GlobalID* global_element_ID = 00293 (iMOAB_GlobalID*)malloc( sizeof( iMOAB_GlobalID ) * num_elements_in_block ); 00294 iMOAB_LocalID* local_element_ID = 00295 (iMOAB_LocalID*)malloc( sizeof( iMOAB_LocalID ) * num_elements_in_block ); 00296 00297 /* 00298 * Global element IDs are determined with this call. Local indices within the 00299 * process are returned too. 00300 */ 00301 rc = iMOAB_GetElementID( pid, &gbIDs[i], &num_elements_in_block, global_element_ID, local_element_ID ); 00302 CHECKRC( rc, "failed to get block elem IDs" ); 00303 for( int j = 0; j < num_elements_in_block; j++ ) 00304 { 00305 printf( " elem %3d owned by %d gid: %4d lid: %4d -- ", j, element_ownership[j], 00306 global_element_ID[j], local_element_ID[j] ); 00307 for( int k = 0; k < vertices_per_element; k++ ) 00308 printf( " %5d", element_connectivity[j * vertices_per_element + k] ); 00309 printf( "\n" ); 00310 } 00311 free( global_element_ID ); 00312 free( local_element_ID ); 00313 free( element_connectivity ); 00314 free( element_ownership ); 00315 } 00316 /* 00317 * query integer tag values on vertices 00318 * this tag was synchronized (see iMOAB_SynchronizeTags above) 00319 * Being a dense tag, all visible vertices have a value. 00320 * As usual, the client allocates the returned array, which is filled with data with 00321 * this call (deep copy) 00322 */ 00323 int* int_tag_vals = (int*)malloc( sizeof( int ) * nverts[2] ); // for all visible vertices on the rank 00324 rc = 00325 iMOAB_GetIntTagStorage( pid, "INTFIELD", &nverts[2], &entTypes[0], int_tag_vals, strlen( "INTFIELD" ) ); 00326 CHECKRC( rc, "failed to get INTFIELD tag" ); 00327 printf( "INTFIELD tag values:\n" ); 00328 for( int i = 0; i < nverts[2]; i++ ) 00329 { 00330 printf( " %4d", int_tag_vals[i] ); 00331 if( i % 20 == 19 ) printf( "\n" ); 00332 } 00333 printf( "\n" ); 00334 free( int_tag_vals ); 00335 00336 /* 00337 * query double tag values on elements 00338 * This tag was not synchronized, so ghost elements have a default value of 0. 00339 */ 00340 double* double_tag_vals = 00341 (double*)malloc( sizeof( double ) * nelem[2] ); // for all visible elements on the rank 00342 rc = iMOAB_GetDoubleTagStorage( pid, "DFIELD", &nelem[2], &entTypes[1], double_tag_vals, 00343 strlen( "DFIELD" ) ); 00344 CHECKRC( rc, "failed to get DFIELD tag" ); 00345 printf( "DFIELD tag values: (not exchanged) \n" ); 00346 for( int i = 0; i < nelem[2]; i++ ) 00347 { 00348 printf( " %f", double_tag_vals[i] ); 00349 if( i % 8 == 7 ) printf( "\n" ); 00350 } 00351 printf( "\n" ); 00352 free( double_tag_vals ); 00353 00354 /* query surface BCs */ 00355 iMOAB_LocalID* surfBC_ID = (iMOAB_LocalID*)malloc( sizeof( iMOAB_LocalID ) * nsbc[2] ); 00356 int* ref_surf = (int*)malloc( sizeof( int ) * nsbc[2] ); 00357 int* bc_value = (int*)malloc( sizeof( int ) * nsbc[2] ); 00358 /* 00359 * Surface boundary condition information is returned for all visible Neumann conditions 00360 * the reference surfaces are indexed from 1 to 6 for hexahedrons, 1 to 4 to 00361 * tetrahedrons, in the MOAB canonical ordering convention 00362 */ 00363 rc = iMOAB_GetPointerToSurfaceBC( pid, &nsbc[2], surfBC_ID, ref_surf, bc_value ); 00364 CHECKRC( rc, "failed to get surf boundary conditions" ); 00365 printf( " Surface boundary conditions:\n" ); 00366 for( int i = 0; i < nsbc[2]; i++ ) 00367 { 00368 printf( " elem_localID %4d side:%d BC:%2d\n", surfBC_ID[i], ref_surf[i], bc_value[i] ); 00369 } 00370 free( surfBC_ID ); 00371 free( ref_surf ); 00372 free( bc_value ); 00373 /* query vertex BCs */ 00374 iMOAB_LocalID* vertBC_ID = (iMOAB_LocalID*)malloc( sizeof( iMOAB_LocalID ) * ndbc[2] ); 00375 int* vertBC_value = (int*)malloc( sizeof( int ) * ndbc[2] ); 00376 rc = iMOAB_GetPointerToVertexBC( pid, &ndbc[2], vertBC_ID, vertBC_value ); 00377 CHECKRC( rc, "failed to get vertex boundary conditions" ); 00378 printf( " Vertex boundary conditions:\n" ); 00379 for( int i = 0; i < ndbc[2]; i++ ) 00380 { 00381 printf( " vertex %4d BC:%2d\n", vertBC_ID[i], vertBC_value[i] ); 00382 } 00383 free( vertBC_ID ); 00384 free( vertBC_value ); 00385 } 00386 #ifdef MOAB_HAVE_MPI 00387 MPI_Barrier( comm ); // to avoid printing problems, as we write all procs data 00388 #endif 00389 } 00390 00391 // free allocated data 00392 free( coords ); 00393 free( vGlobalID ); 00394 free( vranks ); 00395 char outputFile[] = "fnew.h5m"; 00396 #ifdef MOAB_HAVE_MPI 00397 char writeOptions[] = "PARALLEL=WRITE_PART"; 00398 #else 00399 char writeOptions[] = ""; 00400 #endif 00401 /* 00402 * the file can be written in parallel, and it will contain additional tags defined by the user 00403 * we may extend the method to write only desired tags to the file 00404 */ 00405 rc = iMOAB_WriteMesh( pid, outputFile, writeOptions, strlen( outputFile ), strlen( writeOptions ) ); 00406 00407 /* 00408 * deregistering application will delete all mesh entities associated with the application and 00409 * will free allocated tag storage. 00410 */ 00411 rc = iMOAB_DeregisterApplication( pid ); 00412 CHECKRC( rc, "failed to de-register application" ); 00413 00414 rc = iMOAB_DeregisterApplication( pidDup ); 00415 CHECKRC( rc, "failed to de-register application" ); 00416 00417 /* 00418 * this method will delete MOAB instance 00419 */ 00420 rc = iMOAB_Finalize(); 00421 CHECKRC( rc, "failed to finalize MOAB" ); 00422 #ifdef MOAB_HAVE_MPI 00423 MPI_Finalize(); 00424 #endif 00425 00426 return 0; 00427 }