 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
00001 /*
00002 * MGen.cpp
00003 *
00004 */
00005
00006 #include "moab/MeshGeneration.hpp"
00007 #include "moab/MergeMesh.hpp"
00008 #include
00009 #include
00010 #include
00011 #ifdef WIN32 /* windows */
00012 #include
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.
Dark theme by Tilen Majerle. All rights reserved.