MOAB: Mesh Oriented datABase  (version 5.2.1)
MeshGeneration.cpp
Go to the documentation of this file.
00001 /*
00002  * MGen.cpp
00003  *
00004  */
00005 
00006 #include "moab/MeshGeneration.hpp"
00007 #include "moab/MergeMesh.hpp"
00008 #include <iostream>
00009 #include <ctime>
00010 #include <vector>
00011 #ifdef WIN32 /* windows */
00012 #include <time.h>
00013 #endif
00014 
00015 #ifdef MOAB_HAVE_MPI
00016 #include "moab_mpi.h"
00017 #include "moab/ParallelComm.hpp"
00018 #include "MBParallelConventions.h"
00019 #include "moab/ParallelMergeMesh.hpp"
00020 #endif
00021 #include "moab/ReadUtilIface.hpp"
00022 
00023 using std::endl;
00024 using std::string;
00025 using std::vector;
00026 
00027 namespace moab
00028 {
00029 
00030 MeshGeneration::MeshGeneration( Interface* impl, ParallelComm* comm, EntityHandle rset )
00031     : mb( impl ), pc( comm ), cset( rset )
00032 {
00033     // ErrorCode error;
00034 
00035 #ifdef MOAB_HAVE_MPI
00036     // Get the Parallel Comm instance to prepare all new sets to work in parallel
00037     // in case the user did not provide any arguments
00038     if( !comm ) pc = moab::ParallelComm::get_pcomm( mb, 0 );
00039 #endif
00040 }
00041 
00042 MeshGeneration::~MeshGeneration() {}
00043 
00044 ErrorCode MeshGeneration::BrickInstance( MeshGeneration::BrickOpts& opts )
00045 {
00046     int A = opts.A, B = opts.B, C = opts.C, M = opts.M, N = opts.N, K = opts.K;
00047     int blockSize = opts.blockSize;
00048     double xsize = opts.xsize, ysize = opts.ysize, zsize = opts.zsize;  // The size of the region
00049     int GL = opts.GL;                                                   // number of ghost layers
00050 
00051     bool newMergeMethod = opts.newMergeMethod;
00052     bool quadratic      = opts.quadratic;
00053     bool keep_skins     = opts.keep_skins;
00054     bool tetra          = opts.tetra;
00055     bool adjEnts        = opts.adjEnts;
00056     bool parmerge       = opts.parmerge;
00057 
00058     int rank = 0, size = 1;
00059     clock_t tt = clock();
00060 
00061 #ifdef MOAB_HAVE_MPI
00062     rank = pc->rank();
00063     size = pc->size();
00064 #endif
00065 
00066     if( M * N * K != size )
00067     {
00068         if( 0 == rank ) std::cout << "M*N*K = " << M * N * K << " != size = " << size << "\n";
00069 
00070         return MB_FAILURE;
00071     }
00072     // Determine m, n, k for processor rank
00073     int m, n, k;
00074     k            = rank / ( M * N );
00075     int leftover = rank % ( M * N );
00076     n            = leftover / M;
00077     m            = leftover % M;
00078 
00079     // Used for nodes increments
00080     int q = ( quadratic ) ? 2 : 1;
00081     // Used for element increments
00082     int factor = ( tetra ) ? 6 : 1;
00083 
00084     double dx = xsize / ( A * M * blockSize * q );  // Distance between 2 nodes in x direction
00085     double dy = ysize / ( B * N * blockSize * q );  // Distance between 2 nodes in y direction
00086     double dz = zsize / ( C * K * blockSize * q );  // Distance between 2 nodes in z direction
00087 
00088     int NX  = ( q * M * A * blockSize + 1 );
00089     int NY  = ( q * N * B * blockSize + 1 );
00090     int nex = M * A * blockSize;  // Number of elements in x direction, used for global id on element
00091     int ney = N * B * blockSize;  // Number of elements in y direction ...
00092     // int NZ = (K * C * blockSize + 1); // Not used
00093     int blockSize1       = q * blockSize + 1;  // Used for vertices
00094     long num_total_verts = (long)NX * NY * ( K * C * blockSize + 1 );
00095     if( 0 == rank )
00096     {
00097         std::cout << "Generate mesh on " << size << " processors \n";
00098         std::cout << "Total number of vertices: " << num_total_verts << "\n";
00099     }
00100     // int xstride = 1;
00101     int ystride = blockSize1;
00102 
00103     int zstride = blockSize1 * blockSize1;
00104     // Generate the block at (a, b, c); it will represent a partition, it will get a partition tag
00105 
00106     ReadUtilIface* iface;
00107     ErrorCode rval = mb->query_interface( iface );MB_CHK_SET_ERR( rval, "Can't get reader interface" );
00108 
00109     Tag global_id_tag;
00110     rval = mb->tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, global_id_tag );MB_CHK_SET_ERR( rval, "Can't get global id tag" );
00111 
00112     // set global ids
00113     Tag new_id_tag;
00114     if( !parmerge )
00115     {
00116         rval =
00117             mb->tag_get_handle( "HANDLEID", sizeof( long ), MB_TYPE_OPAQUE, new_id_tag, MB_TAG_CREAT | MB_TAG_DENSE );MB_CHK_SET_ERR( rval, "Can't get handle id tag" );
00118     }
00119     Tag part_tag;
00120     int dum_id = -1;
00121     rval =
00122         mb->tag_get_handle( "PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, part_tag, MB_TAG_CREAT | MB_TAG_SPARSE, &dum_id );MB_CHK_SET_ERR( rval, "Can't get parallel partition tag" );
00123 
00124     Range wsets;  // write only part sets
00125     Range localVerts;
00126     Range all3dcells;
00127     for( int a = 0; a < A; a++ )
00128     {
00129         for( int b = 0; b < B; b++ )
00130         {
00131             for( int c = 0; c < C; c++ )
00132             {
00133                 // We will generate (q*block + 1)^3 vertices, and block^3 hexas; q is 1 for linear,
00134                 // 2 for quadratic the global id of the vertices will come from m, n, k, a, b, c x
00135                 // will vary from  m*A*q*block + a*q*block to m*A*q*block + (a+1)*q*block etc;
00136                 int num_nodes = blockSize1 * blockSize1 * blockSize1;
00137 
00138                 vector< double* > arrays;
00139                 EntityHandle startv;
00140                 rval = iface->get_node_coords( 3, num_nodes, 0, startv, arrays );MB_CHK_SET_ERR( rval, "Can't get node coords" );
00141 
00142                 // Will start with the lower corner:
00143                 int x  = m * A * q * blockSize + a * q * blockSize;
00144                 int y  = n * B * q * blockSize + b * q * blockSize;
00145                 int z  = k * C * q * blockSize + c * q * blockSize;
00146                 int ix = 0;
00147                 vector< int > gids( num_nodes );
00148                 vector< long > lgids( num_nodes );
00149                 Range verts( startv, startv + num_nodes - 1 );
00150                 for( int kk = 0; kk < blockSize1; kk++ )
00151                 {
00152                     for( int jj = 0; jj < blockSize1; jj++ )
00153                     {
00154                         for( int ii = 0; ii < blockSize1; ii++ )
00155                         {
00156                             arrays[0][ix] = ( x + ii ) * dx;
00157                             arrays[1][ix] = ( y + jj ) * dy;
00158                             arrays[2][ix] = ( z + kk ) * dz;
00159                             gids[ix]      = 1 + ( x + ii ) + ( y + jj ) * NX + ( z + kk ) * ( NX * NY );
00160                             if( !parmerge )
00161                                 lgids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + (long)( z + kk ) * ( NX * NY );
00162                             // Set int tags, some nice values?
00163 
00164                             ix++;
00165                         }
00166                     }
00167                 }
00168 
00169                 rval = mb->tag_set_data( global_id_tag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids to vertices" );
00170                 if( !parmerge )
00171                 {
00172                     rval = mb->tag_set_data( new_id_tag, verts, &lgids[0] );MB_CHK_SET_ERR( rval, "Can't set the new handle id tags" );
00173                 }
00174                 localVerts.merge( verts );
00175                 int num_hexas = blockSize * blockSize * blockSize;
00176                 int num_el    = num_hexas * factor;
00177 
00178                 EntityHandle starte;  // Connectivity
00179                 EntityHandle* conn;
00180                 int num_v_per_elem = 8;
00181                 if( quadratic )
00182                 {
00183                     num_v_per_elem = 27;
00184                     rval           = iface->get_element_connect( num_el, 27, MBHEX, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
00185                 }
00186                 else if( tetra )
00187                 {
00188                     num_v_per_elem = 4;
00189                     rval           = iface->get_element_connect( num_el, 4, MBTET, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
00190                 }
00191                 else
00192                 {
00193                     rval = iface->get_element_connect( num_el, 8, MBHEX, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
00194                 }
00195 
00196                 Range cells( starte, starte + num_el - 1 );  // Should be elements
00197                 // Fill cells
00198                 ix = 0;
00199                 // Identify the elements at the lower corner, for their global ids
00200                 int xe = m * A * blockSize + a * blockSize;
00201                 int ye = n * B * blockSize + b * blockSize;
00202                 int ze = k * C * blockSize + c * blockSize;
00203                 gids.resize( num_el );
00204                 lgids.resize( num_el );
00205                 int ie = 0;  // Index now in the elements, for global ids
00206                 for( int kk = 0; kk < blockSize; kk++ )
00207                 {
00208                     for( int jj = 0; jj < blockSize; jj++ )
00209                     {
00210                         for( int ii = 0; ii < blockSize; ii++ )
00211                         {
00212                             EntityHandle corner = startv + q * ii + q * jj * ystride + q * kk * zstride;
00213                             // These could overflow for large numbers
00214                             gids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + ( ze + kk ) * ( nex * ney ) ) *
00215                                                factor;  // 6 more for tetra
00216                             lgids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + (long)( ze + kk ) * ( nex * ney ) ) *
00217                                                 factor;  // 6 more for tetra
00218                             // EntityHandle eh = starte + ie;
00219 
00220                             ie++;
00221                             if( quadratic )
00222                             {
00223                                 //                    4   ----- 19   -----  7
00224                                 //                .   |                 .   |
00225                                 //            16         25         18      |
00226                                 //         .          |          .          |
00227                                 //      5   ----- 17   -----  6             |
00228                                 //      |            12       | 23         15
00229                                 //      |                     |             |
00230                                 //      |     20      |  26   |     22      |
00231                                 //      |                     |             |
00232                                 //     13         21  |      14             |
00233                                 //      |             0   ----- 11   -----  3
00234                                 //      |         .           |         .
00235                                 //      |      8         24   |     10
00236                                 //      |  .                  |  .
00237                                 //      1   -----  9   -----  2
00238                                 //
00239                                 conn[ix]      = corner;
00240                                 conn[ix + 1]  = corner + 2;
00241                                 conn[ix + 2]  = corner + 2 + 2 * ystride;
00242                                 conn[ix + 3]  = corner + 2 * ystride;
00243                                 conn[ix + 4]  = corner + 2 * zstride;
00244                                 conn[ix + 5]  = corner + 2 + 2 * zstride;
00245                                 conn[ix + 6]  = corner + 2 + 2 * ystride + 2 * zstride;
00246                                 conn[ix + 7]  = corner + 2 * ystride + 2 * zstride;
00247                                 conn[ix + 8]  = corner + 1;                              // 0-1
00248                                 conn[ix + 9]  = corner + 2 + ystride;                    // 1-2
00249                                 conn[ix + 10] = corner + 1 + 2 * ystride;                // 2-3
00250                                 conn[ix + 11] = corner + ystride;                        // 3-0
00251                                 conn[ix + 12] = corner + zstride;                        // 0-4
00252                                 conn[ix + 13] = corner + 2 + zstride;                    // 1-5
00253                                 conn[ix + 14] = corner + 2 + 2 * ystride + zstride;      // 2-6
00254                                 conn[ix + 15] = corner + 2 * ystride + zstride;          // 3-7
00255                                 conn[ix + 16] = corner + 1 + 2 * zstride;                // 4-5
00256                                 conn[ix + 17] = corner + 2 + ystride + 2 * zstride;      // 5-6
00257                                 conn[ix + 18] = corner + 1 + 2 * ystride + 2 * zstride;  // 6-7
00258                                 conn[ix + 19] = corner + ystride + 2 * zstride;          // 4-7
00259                                 conn[ix + 20] = corner + 1 + zstride;                    // 0154
00260                                 conn[ix + 21] = corner + 2 + ystride + zstride;          // 1265
00261                                 conn[ix + 22] = corner + 1 + 2 * ystride + zstride;      // 2376
00262                                 conn[ix + 23] = corner + ystride + zstride;              // 0374
00263                                 conn[ix + 24] = corner + 1 + ystride;                    // 0123
00264                                 conn[ix + 25] = corner + 1 + ystride + 2 * zstride;      // 4567
00265                                 conn[ix + 26] = corner + 1 + ystride + zstride;          // center
00266                                 ix += 27;
00267                             }
00268                             else if( tetra )
00269                             {
00270                                 //        E      H
00271                                 //     F     G
00272                                 //
00273                                 //        A     D
00274                                 //     B     C
00275                                 EntityHandle AA = corner;
00276                                 EntityHandle BB = corner + 1;
00277                                 EntityHandle CC = corner + 1 + ystride;
00278                                 EntityHandle D  = corner + ystride;
00279                                 EntityHandle E  = corner + zstride;
00280                                 EntityHandle F  = corner + 1 + zstride;
00281                                 EntityHandle G  = corner + 1 + ystride + zstride;
00282                                 EntityHandle H  = corner + ystride + zstride;
00283 
00284                                 // tet EDHG
00285                                 conn[ix]     = E;
00286                                 conn[ix + 1] = D;
00287                                 conn[ix + 2] = H;
00288                                 conn[ix + 3] = G;
00289 
00290                                 // tet ABCF
00291                                 conn[ix + 4] = AA;
00292                                 conn[ix + 5] = BB;
00293                                 conn[ix + 6] = CC;
00294                                 conn[ix + 7] = F;
00295 
00296                                 // tet ADEF
00297                                 conn[ix + 8]  = AA;
00298                                 conn[ix + 9]  = D;
00299                                 conn[ix + 10] = E;
00300                                 conn[ix + 11] = F;
00301 
00302                                 // tet CGDF
00303                                 conn[ix + 12] = CC;
00304                                 conn[ix + 13] = G;
00305                                 conn[ix + 14] = D;
00306                                 conn[ix + 15] = F;
00307 
00308                                 // tet ACDF
00309                                 conn[ix + 16] = AA;
00310                                 conn[ix + 17] = CC;
00311                                 conn[ix + 18] = D;
00312                                 conn[ix + 19] = F;
00313 
00314                                 // tet DGEF
00315                                 conn[ix + 20] = D;
00316                                 conn[ix + 21] = G;
00317                                 conn[ix + 22] = E;
00318                                 conn[ix + 23] = F;
00319                                 ix += 24;
00320                                 for( int ff = 0; ff < factor - 1; ff++ )
00321                                 {
00322                                     gids[ie] = gids[ie - 1] + 1;  // 6 more for tetra
00323 
00324                                     // eh = starte + ie;
00325 
00326                                     ie++;
00327                                 }
00328                             }
00329                             else
00330                             {  // Linear hex
00331                                 conn[ix]     = corner;
00332                                 conn[ix + 1] = corner + 1;
00333                                 conn[ix + 2] = corner + 1 + ystride;
00334                                 conn[ix + 3] = corner + ystride;
00335                                 conn[ix + 4] = corner + zstride;
00336                                 conn[ix + 5] = corner + 1 + zstride;
00337                                 conn[ix + 6] = corner + 1 + ystride + zstride;
00338                                 conn[ix + 7] = corner + ystride + zstride;
00339                                 ix += 8;
00340                             }
00341                         }
00342                     }
00343                 }
00344 
00345                 EntityHandle part_set;
00346                 rval = mb->create_meshset( MESHSET_SET, part_set );MB_CHK_SET_ERR( rval, "Can't create mesh set" );
00347                 rval = mb->add_entities( part_set, cells );MB_CHK_SET_ERR( rval, "Can't add entities to set" );
00348                 all3dcells.merge( cells );
00349                 // update adjacencies now, because some elements are new;
00350                 rval = iface->update_adjacencies( starte, num_el, num_v_per_elem, conn );MB_CHK_SET_ERR( rval, "Can't update adjacencies" );
00351                 // If needed, add all edges and faces
00352                 if( adjEnts )
00353                 {
00354                     // Generate all adj entities dimension 1 and 2 (edges and faces/ tri or qua)
00355                     Range edges, faces;
00356                     rval = mb->get_adjacencies( cells, 1, true, edges, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get edges" );
00357                     rval = mb->get_adjacencies( cells, 2, true, faces, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get faces" );
00358                     // rval = mb->add_entities(part_set, edges);MB_CHK_SET_ERR(rval, "Can't add
00359                     // edges to partition set"); rval = mb->add_entities(part_set,
00360                     // faces);MB_CHK_SET_ERR(rval, "Can't add faces to partition set");
00361                 }
00362 
00363                 rval = mb->tag_set_data( global_id_tag, cells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids to elements" );
00364                 if( !parmerge )
00365                 {
00366                     rval = mb->tag_set_data( new_id_tag, cells, &lgids[0] );MB_CHK_SET_ERR( rval, "Can't set new ids to elements" );
00367                 }
00368                 int part_num = a + m * A + ( b + n * B ) * ( M * A ) + ( c + k * C ) * ( M * A * N * B );
00369                 rval         = mb->tag_set_data( part_tag, &part_set, 1, &part_num );MB_CHK_SET_ERR( rval, "Can't set part tag on set" );
00370                 wsets.insert( part_set );
00371             }
00372         }
00373     }
00374 
00375     mb->add_entities( cset, all3dcells );
00376     rval = mb->add_entities( cset, wsets );MB_CHK_SET_ERR( rval, "Can't add entity sets" );
00377 #ifdef MOAB_HAVE_MPI
00378     pc->partition_sets() = wsets;
00379 #endif
00380 
00381     /*
00382     // Before merge locally
00383     rval = mb->write_file("test0.h5m", 0, ";;PARALLEL=WRITE_PART");MB_CHK_SET_ERR(rval, "Can't write
00384     in parallel, before merging");
00385     */
00386     // After the mesh is generated on each proc, merge the vertices
00387     MergeMesh mm( mb );
00388 
00389     // rval = mb->get_entities_by_dimension(0, 3, all3dcells);MB_CHK_SET_ERR(rval, "Can't get all 3d
00390     // cells elements");
00391 
00392     if( 0 == rank )
00393     {
00394         std::cout << "generate local mesh: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
00395         tt = clock();
00396         std::cout << "number of elements on rank 0: " << all3dcells.size() << endl;
00397         std::cout << "Total number of elements " << all3dcells.size() * size << endl;
00398         std::cout << "Element type: " << ( tetra ? "MBTET" : "MBHEX" )
00399                   << " order:" << ( quadratic ? "quadratic" : "linear" ) << endl;
00400     }
00401 
00402     if( A * B * C != 1 )
00403     {  // Merge needed
00404         if( newMergeMethod )
00405         {
00406             rval = mm.merge_using_integer_tag( localVerts, global_id_tag );MB_CHK_SET_ERR( rval, "Can't merge" );
00407         }
00408         else
00409         {
00410             rval = mm.merge_entities( all3dcells, 0.0001 );MB_CHK_SET_ERR( rval, "Can't merge" );
00411         }
00412 
00413         if( 0 == rank )
00414         {
00415             std::cout << "merge locally: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
00416             tt = clock();
00417         }
00418     }
00419     // if adjEnts, add now to each set
00420     if( adjEnts )
00421     {
00422         for( Range::iterator wsit = wsets.begin(); wsit != wsets.end(); ++wsit )
00423         {
00424             EntityHandle ws = *wsit;  // write set
00425             Range cells, edges, faces;
00426             rval = mb->get_entities_by_dimension( ws, 3, cells );MB_CHK_SET_ERR( rval, "Can't get cells" );
00427             rval = mb->get_adjacencies( cells, 1, false, edges, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get edges" );
00428             rval = mb->get_adjacencies( cells, 2, false, faces, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get faces" );
00429             rval = mb->add_entities( ws, edges );MB_CHK_SET_ERR( rval, "Can't add edges to partition set" );
00430             rval = mb->add_entities( ws, faces );MB_CHK_SET_ERR( rval, "Can't add faces to partition set" );
00431         }
00432     }
00433 #ifdef MOAB_HAVE_MPI
00434     if( size > 1 )
00435     {
00436 
00437         // rval = mb->create_meshset(MESHSET_SET, mesh_set);MB_CHK_SET_ERR(rval, "Can't create new
00438         // set");
00439 
00440         if( parmerge )
00441         {
00442             ParallelMergeMesh pm( pc, 0.00001 );
00443             rval = pm.merge();MB_CHK_SET_ERR( rval, "Can't resolve shared ents" );
00444             if( 0 == rank )
00445             {
00446                 std::cout << "parallel mesh merge: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
00447                 tt = clock();
00448             }
00449         }
00450         else
00451         {
00452             rval = pc->resolve_shared_ents( cset, -1, -1, &new_id_tag );MB_CHK_SET_ERR( rval, "Can't resolve shared ents" );
00453 
00454             if( 0 == rank )
00455             {
00456                 std::cout << "resolve shared entities: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds"
00457                           << endl;
00458                 tt = clock();
00459             }
00460         }
00461         if( !keep_skins )
00462         {  // Default is to delete the 1- and 2-dimensional entities
00463             // Delete all quads and edges
00464             Range toDelete;
00465             rval = mb->get_entities_by_dimension( cset, 1, toDelete );MB_CHK_SET_ERR( rval, "Can't get edges" );
00466 
00467             rval = mb->get_entities_by_dimension( cset, 2, toDelete );MB_CHK_SET_ERR( rval, "Can't get faces" );
00468 
00469             rval = pc->delete_entities( toDelete );MB_CHK_SET_ERR( rval, "Can't delete entities" );
00470             rval = mb->remove_entities( cset, toDelete );MB_CHK_SET_ERR( rval, "Can't remove entities from base set" );
00471             if( 0 == rank )
00472             {
00473 
00474                 std::cout << "delete edges and faces \n";
00475                 toDelete.print( std::cout );
00476 
00477                 std::cout << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
00478                 tt = clock();
00479             }
00480         }
00481         // do some ghosting if required
00482         if( GL > 0 )
00483         {
00484             rval = pc->exchange_ghost_cells( 3,   // int ghost_dim
00485                                              0,   // int bridge_dim
00486                                              GL,  // int num_layers
00487                                              0,   // int addl_ents
00488                                              true );MB_CHK_ERR( rval );  // bool store_remote_handles
00489             if( 0 == rank )
00490             {
00491                 std::cout << "exchange  " << GL << " ghost layer(s) :" << ( clock() - tt ) / (double)CLOCKS_PER_SEC
00492                           << " seconds" << endl;
00493                 tt = clock();
00494             }
00495         }
00496     }
00497 #endif
00498 
00499     if( !parmerge )
00500     {
00501         rval = mb->tag_delete( new_id_tag );MB_CHK_SET_ERR( rval, "Can't delete new ID tag" );
00502     }
00503 
00504     return MB_SUCCESS;
00505 }
00506 
00507 } /* namespace moab */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines