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