|
MOAB: Mesh Oriented datABase
(version 5.2.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 */
1.7.6.1 
