MOAB: Mesh Oriented datABase  (version 5.3.1)
imoab_test.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines