#include <MeshGeneration.hpp>
Classes | |
struct | BrickOpts |
Public Member Functions | |
MeshGeneration (Interface *mbi, EntityHandle rset=0) | |
virtual | ~MeshGeneration () |
ErrorCode | BrickInstance (BrickOpts &opts) |
Private Attributes | |
Interface * | mb |
EntityHandle | cset |
Definition at line 18 of file MeshGeneration.hpp.
moab::MeshGeneration::MeshGeneration | ( | Interface * | mbi, |
EntityHandle | rset = 0 |
||
) |
Definition at line 30 of file MeshGeneration.cpp.
35 : mb( impl ),
36 #ifdef MOAB_HAVE_MPI
37 pc( comm ),
38 #endif
39 cset( rset )
40 {
41 // ErrorCode error;
42
43 #ifdef MOAB_HAVE_MPI
44 // Get the Parallel Comm instance to prepare all new sets to work in parallel
45 // in case the user did not provide any arguments
46 if( !comm ) pc = moab::ParallelComm::get_pcomm( mb, 0 );
47 #endif
48 }
References moab::ParallelComm::get_pcomm(), and mb.
|
virtual |
Definition at line 50 of file MeshGeneration.cpp.
50 {}
ErrorCode moab::MeshGeneration::BrickInstance | ( | MeshGeneration::BrickOpts & | opts | ) |
Definition at line 52 of file MeshGeneration.cpp.
53 {
54 int A = opts.A, B = opts.B, C = opts.C, M = opts.M, N = opts.N, K = opts.K;
55 int blockSize = opts.blockSize;
56 double xsize = opts.xsize, ysize = opts.ysize, zsize = opts.zsize; // The size of the region
57 bool newMergeMethod = opts.newMergeMethod;
58 bool quadratic = opts.quadratic;
59 bool tetra = opts.tetra;
60 bool adjEnts = opts.adjEnts;
61 bool parmerge = opts.parmerge;
62
63 #ifdef MOAB_HAVE_MPI
64 int GL = opts.GL; // number of ghost layers
65 bool keep_skins = opts.keep_skins;
66 #endif
67
68 int rank = 0, size = 1;
69 #ifndef _MSC_VER /* windows */
70 clock_t tt = clock();
71 #endif
72 #ifdef MOAB_HAVE_MPI
73 rank = pc->rank();
74 size = pc->size();
75 #endif
76
77 if( M * N * K != size )
78 {
79 if( 0 == rank ) std::cout << "M*N*K = " << M * N * K << " != size = " << size << "\n";
80
81 return MB_FAILURE;
82 }
83 // Determine m, n, k for processor rank
84 int m, n, k;
85 k = rank / ( M * N );
86 int leftover = rank % ( M * N );
87 n = leftover / M;
88 m = leftover % M;
89
90 // Used for nodes increments
91 int q = ( quadratic ) ? 2 : 1;
92 // Used for element increments
93 int factor = ( tetra ) ? 6 : 1;
94
95 double dx = xsize / ( A * M * blockSize * q ); // Distance between 2 nodes in x direction
96 double dy = ysize / ( B * N * blockSize * q ); // Distance between 2 nodes in y direction
97 double dz = zsize / ( C * K * blockSize * q ); // Distance between 2 nodes in z direction
98
99 int NX = ( q * M * A * blockSize + 1 );
100 int NY = ( q * N * B * blockSize + 1 );
101 int nex = M * A * blockSize; // Number of elements in x direction, used for global id on element
102 int ney = N * B * blockSize; // Number of elements in y direction ...
103 // int NZ = (K * C * blockSize + 1); // Not used
104 int blockSize1 = q * blockSize + 1; // Used for vertices
105 long num_total_verts = (long)NX * NY * ( K * C * blockSize + 1 );
106 if( 0 == rank )
107 {
108 std::cout << "Generate mesh on " << size << " processors \n";
109 std::cout << "Total number of vertices: " << num_total_verts << "\n";
110 }
111 // int xstride = 1;
112 int ystride = blockSize1;
113
114 int zstride = blockSize1 * blockSize1;
115 // Generate the block at (a, b, c); it will represent a partition, it will get a partition tag
116
117 ReadUtilIface* iface;
118 ErrorCode rval = mb->query_interface( iface );MB_CHK_SET_ERR( rval, "Can't get reader interface" );
119
120 Tag global_id_tag;
121 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" );
122
123 // set global ids
124 Tag new_id_tag;
125 if( !parmerge )
126 {
127 rval =
128 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" );
129 }
130 Tag part_tag;
131 int dum_id = -1;
132 rval =
133 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" );
134
135 Range wsets; // write only part sets
136 Range localVerts;
137 Range all3dcells;
138 for( int a = 0; a < A; a++ )
139 {
140 for( int b = 0; b < B; b++ )
141 {
142 for( int c = 0; c < C; c++ )
143 {
144 // We will generate (q*block + 1)^3 vertices, and block^3 hexas; q is 1 for linear,
145 // 2 for quadratic the global id of the vertices will come from m, n, k, a, b, c x
146 // will vary from m*A*q*block + a*q*block to m*A*q*block + (a+1)*q*block etc;
147 int num_nodes = blockSize1 * blockSize1 * blockSize1;
148
149 vector< double* > arrays;
150 EntityHandle startv;
151 rval = iface->get_node_coords( 3, num_nodes, 0, startv, arrays );MB_CHK_SET_ERR( rval, "Can't get node coords" );
152
153 // Will start with the lower corner:
154 int x = m * A * q * blockSize + a * q * blockSize;
155 int y = n * B * q * blockSize + b * q * blockSize;
156 int z = k * C * q * blockSize + c * q * blockSize;
157 int ix = 0;
158 vector< int > gids( num_nodes );
159 vector< long > lgids( num_nodes );
160 Range verts( startv, startv + num_nodes - 1 );
161 for( int kk = 0; kk < blockSize1; kk++ )
162 {
163 for( int jj = 0; jj < blockSize1; jj++ )
164 {
165 for( int ii = 0; ii < blockSize1; ii++ )
166 {
167 arrays[0][ix] = ( x + ii ) * dx;
168 arrays[1][ix] = ( y + jj ) * dy;
169 arrays[2][ix] = ( z + kk ) * dz;
170 gids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + ( z + kk ) * ( NX * NY );
171 if( !parmerge )
172 lgids[ix] = 1 + ( x + ii ) + ( y + jj ) * NX + (long)( z + kk ) * ( NX * NY );
173 // Set int tags, some nice values?
174
175 ix++;
176 }
177 }
178 }
179
180 rval = mb->tag_set_data( global_id_tag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids to vertices" );
181 if( !parmerge )
182 {
183 rval = mb->tag_set_data( new_id_tag, verts, &lgids[0] );MB_CHK_SET_ERR( rval, "Can't set the new handle id tags" );
184 }
185 localVerts.merge( verts );
186 int num_hexas = blockSize * blockSize * blockSize;
187 int num_el = num_hexas * factor;
188
189 EntityHandle starte; // Connectivity
190 EntityHandle* conn;
191 int num_v_per_elem = 8;
192 if( quadratic )
193 {
194 num_v_per_elem = 27;
195 rval = iface->get_element_connect( num_el, 27, MBHEX, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
196 }
197 else if( tetra )
198 {
199 num_v_per_elem = 4;
200 rval = iface->get_element_connect( num_el, 4, MBTET, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
201 }
202 else
203 {
204 rval = iface->get_element_connect( num_el, 8, MBHEX, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
205 }
206
207 Range cells( starte, starte + num_el - 1 ); // Should be elements
208 // Fill cells
209 ix = 0;
210 // Identify the elements at the lower corner, for their global ids
211 int xe = m * A * blockSize + a * blockSize;
212 int ye = n * B * blockSize + b * blockSize;
213 int ze = k * C * blockSize + c * blockSize;
214 gids.resize( num_el );
215 lgids.resize( num_el );
216 int ie = 0; // Index now in the elements, for global ids
217 for( int kk = 0; kk < blockSize; kk++ )
218 {
219 for( int jj = 0; jj < blockSize; jj++ )
220 {
221 for( int ii = 0; ii < blockSize; ii++ )
222 {
223 EntityHandle corner = startv + q * ii + q * jj * ystride + q * kk * zstride;
224 // These could overflow for large numbers
225 gids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + ( ze + kk ) * ( nex * ney ) ) *
226 factor; // 6 more for tetra
227 lgids[ie] = 1 + ( ( xe + ii ) + ( ye + jj ) * nex + (long)( ze + kk ) * ( nex * ney ) ) *
228 factor; // 6 more for tetra
229 // EntityHandle eh = starte + ie;
230
231 ie++;
232 if( quadratic )
233 {
234 // 4 ----- 19 ----- 7
235 // . | . |
236 // 16 25 18 |
237 // . | . |
238 // 5 ----- 17 ----- 6 |
239 // | 12 | 23 15
240 // | | |
241 // | 20 | 26 | 22 |
242 // | | |
243 // 13 21 | 14 |
244 // | 0 ----- 11 ----- 3
245 // | . | .
246 // | 8 24 | 10
247 // | . | .
248 // 1 ----- 9 ----- 2
249 //
250 conn[ix] = corner;
251 conn[ix + 1] = corner + 2;
252 conn[ix + 2] = corner + 2 + 2 * ystride;
253 conn[ix + 3] = corner + 2 * ystride;
254 conn[ix + 4] = corner + 2 * zstride;
255 conn[ix + 5] = corner + 2 + 2 * zstride;
256 conn[ix + 6] = corner + 2 + 2 * ystride + 2 * zstride;
257 conn[ix + 7] = corner + 2 * ystride + 2 * zstride;
258 conn[ix + 8] = corner + 1; // 0-1
259 conn[ix + 9] = corner + 2 + ystride; // 1-2
260 conn[ix + 10] = corner + 1 + 2 * ystride; // 2-3
261 conn[ix + 11] = corner + ystride; // 3-0
262 conn[ix + 12] = corner + zstride; // 0-4
263 conn[ix + 13] = corner + 2 + zstride; // 1-5
264 conn[ix + 14] = corner + 2 + 2 * ystride + zstride; // 2-6
265 conn[ix + 15] = corner + 2 * ystride + zstride; // 3-7
266 conn[ix + 16] = corner + 1 + 2 * zstride; // 4-5
267 conn[ix + 17] = corner + 2 + ystride + 2 * zstride; // 5-6
268 conn[ix + 18] = corner + 1 + 2 * ystride + 2 * zstride; // 6-7
269 conn[ix + 19] = corner + ystride + 2 * zstride; // 4-7
270 conn[ix + 20] = corner + 1 + zstride; // 0154
271 conn[ix + 21] = corner + 2 + ystride + zstride; // 1265
272 conn[ix + 22] = corner + 1 + 2 * ystride + zstride; // 2376
273 conn[ix + 23] = corner + ystride + zstride; // 0374
274 conn[ix + 24] = corner + 1 + ystride; // 0123
275 conn[ix + 25] = corner + 1 + ystride + 2 * zstride; // 4567
276 conn[ix + 26] = corner + 1 + ystride + zstride; // center
277 ix += 27;
278 }
279 else if( tetra )
280 {
281 // E H
282 // F G
283 //
284 // A D
285 // B C
286 EntityHandle AA = corner;
287 EntityHandle BB = corner + 1;
288 EntityHandle CC = corner + 1 + ystride;
289 EntityHandle D = corner + ystride;
290 EntityHandle E = corner + zstride;
291 EntityHandle F = corner + 1 + zstride;
292 EntityHandle G = corner + 1 + ystride + zstride;
293 EntityHandle H = corner + ystride + zstride;
294
295 // tet EDHG
296 conn[ix] = E;
297 conn[ix + 1] = D;
298 conn[ix + 2] = H;
299 conn[ix + 3] = G;
300
301 // tet ABCF
302 conn[ix + 4] = AA;
303 conn[ix + 5] = BB;
304 conn[ix + 6] = CC;
305 conn[ix + 7] = F;
306
307 // tet ADEF
308 conn[ix + 8] = AA;
309 conn[ix + 9] = D;
310 conn[ix + 10] = E;
311 conn[ix + 11] = F;
312
313 // tet CGDF
314 conn[ix + 12] = CC;
315 conn[ix + 13] = G;
316 conn[ix + 14] = D;
317 conn[ix + 15] = F;
318
319 // tet ACDF
320 conn[ix + 16] = AA;
321 conn[ix + 17] = CC;
322 conn[ix + 18] = D;
323 conn[ix + 19] = F;
324
325 // tet DGEF
326 conn[ix + 20] = D;
327 conn[ix + 21] = G;
328 conn[ix + 22] = E;
329 conn[ix + 23] = F;
330 ix += 24;
331 for( int ff = 0; ff < factor - 1; ff++ )
332 {
333 gids[ie] = gids[ie - 1] + 1; // 6 more for tetra
334
335 // eh = starte + ie;
336
337 ie++;
338 }
339 }
340 else
341 { // Linear hex
342 conn[ix] = corner;
343 conn[ix + 1] = corner + 1;
344 conn[ix + 2] = corner + 1 + ystride;
345 conn[ix + 3] = corner + ystride;
346 conn[ix + 4] = corner + zstride;
347 conn[ix + 5] = corner + 1 + zstride;
348 conn[ix + 6] = corner + 1 + ystride + zstride;
349 conn[ix + 7] = corner + ystride + zstride;
350 ix += 8;
351 }
352 }
353 }
354 }
355
356 EntityHandle part_set;
357 rval = mb->create_meshset( MESHSET_SET, part_set );MB_CHK_SET_ERR( rval, "Can't create mesh set" );
358 rval = mb->add_entities( part_set, cells );MB_CHK_SET_ERR( rval, "Can't add entities to set" );
359 all3dcells.merge( cells );
360 // update adjacencies now, because some elements are new;
361 rval = iface->update_adjacencies( starte, num_el, num_v_per_elem, conn );MB_CHK_SET_ERR( rval, "Can't update adjacencies" );
362 // If needed, add all edges and faces
363 if( adjEnts )
364 {
365 // Generate all adj entities dimension 1 and 2 (edges and faces/ tri or qua)
366 Range edges, faces;
367 rval = mb->get_adjacencies( cells, 1, true, edges, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get edges" );
368 rval = mb->get_adjacencies( cells, 2, true, faces, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get faces" );
369 // rval = mb->add_entities(part_set, edges);MB_CHK_SET_ERR(rval, "Can't add
370 // edges to partition set"); rval = mb->add_entities(part_set,
371 // faces);MB_CHK_SET_ERR(rval, "Can't add faces to partition set");
372 }
373
374 rval = mb->tag_set_data( global_id_tag, cells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids to elements" );
375 if( !parmerge )
376 {
377 rval = mb->tag_set_data( new_id_tag, cells, &lgids[0] );MB_CHK_SET_ERR( rval, "Can't set new ids to elements" );
378 }
379 int part_num = a + m * A + ( b + n * B ) * ( M * A ) + ( c + k * C ) * ( M * A * N * B );
380 rval = mb->tag_set_data( part_tag, &part_set, 1, &part_num );MB_CHK_SET_ERR( rval, "Can't set part tag on set" );
381 wsets.insert( part_set );
382 }
383 }
384 }
385
386 mb->add_entities( cset, all3dcells );
387 rval = mb->add_entities( cset, wsets );MB_CHK_SET_ERR( rval, "Can't add entity sets" );
388 #ifdef MOAB_HAVE_MPI
389 pc->partition_sets() = wsets;
390 #endif
391
392 /*
393 // Before merge locally
394 rval = mb->write_file("test0.h5m", 0, ";;PARALLEL=WRITE_PART");MB_CHK_SET_ERR(rval, "Can't write
395 in parallel, before merging");
396 */
397 // After the mesh is generated on each proc, merge the vertices
398 MergeMesh mm( mb );
399
400 // rval = mb->get_entities_by_dimension(0, 3, all3dcells);MB_CHK_SET_ERR(rval, "Can't get all 3d
401 // cells elements");
402
403 if( 0 == rank )
404 {
405 #ifndef _MSC_VER /* windows */
406 std::cout << "generate local mesh: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
407 tt = clock();
408 #endif
409
410 std::cout << "number of elements on rank 0: " << all3dcells.size() << endl;
411 std::cout << "Total number of elements " << all3dcells.size() * size << endl;
412 std::cout << "Element type: " << ( tetra ? "MBTET" : "MBHEX" )
413 << " order:" << ( quadratic ? "quadratic" : "linear" ) << endl;
414 }
415
416 if( A * B * C != 1 )
417 { // Merge needed
418 if( newMergeMethod )
419 {
420 rval = mm.merge_using_integer_tag( localVerts, global_id_tag );MB_CHK_SET_ERR( rval, "Can't merge" );
421 }
422 else
423 {
424 rval = mm.merge_entities( all3dcells, 0.0001 );MB_CHK_SET_ERR( rval, "Can't merge" );
425 }
426 #ifndef _MSC_VER /* windows */
427 if( 0 == rank )
428 {
429 std::cout << "merge locally: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
430 tt = clock();
431 }
432 #endif
433 }
434 // if adjEnts, add now to each set
435 if( adjEnts )
436 {
437 for( Range::iterator wsit = wsets.begin(); wsit != wsets.end(); ++wsit )
438 {
439 EntityHandle ws = *wsit; // write set
440 Range cells, edges, faces;
441 rval = mb->get_entities_by_dimension( ws, 3, cells );MB_CHK_SET_ERR( rval, "Can't get cells" );
442 rval = mb->get_adjacencies( cells, 1, false, edges, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get edges" );
443 rval = mb->get_adjacencies( cells, 2, false, faces, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get faces" );
444 rval = mb->add_entities( ws, edges );MB_CHK_SET_ERR( rval, "Can't add edges to partition set" );
445 rval = mb->add_entities( ws, faces );MB_CHK_SET_ERR( rval, "Can't add faces to partition set" );
446 }
447 }
448 #ifdef MOAB_HAVE_MPI
449 if( size > 1 )
450 {
451
452 // rval = mb->create_meshset(MESHSET_SET, mesh_set);MB_CHK_SET_ERR(rval, "Can't create new
453 // set");
454
455 if( parmerge )
456 {
457 ParallelMergeMesh pm( pc, 0.00001 );
458 rval = pm.merge();MB_CHK_SET_ERR( rval, "Can't resolve shared ents" );
459 #ifndef _MSC_VER /* windows */
460 if( 0 == rank )
461 {
462 std::cout << "parallel mesh merge: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
463 tt = clock();
464 }
465 #endif
466 }
467 else
468 {
469 rval = pc->resolve_shared_ents( cset, -1, -1, &new_id_tag );MB_CHK_SET_ERR( rval, "Can't resolve shared ents" );
470 #ifndef _MSC_VER /* windows */
471 if( 0 == rank )
472 {
473 std::cout << "resolve shared entities: " << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds"
474 << endl;
475 tt = clock();
476 }
477 #endif
478 }
479 if( !keep_skins )
480 { // Default is to delete the 1- and 2-dimensional entities
481 // Delete all quads and edges
482 Range toDelete;
483 rval = mb->get_entities_by_dimension( cset, 1, toDelete );MB_CHK_SET_ERR( rval, "Can't get edges" );
484
485 rval = mb->get_entities_by_dimension( cset, 2, toDelete );MB_CHK_SET_ERR( rval, "Can't get faces" );
486
487 rval = pc->delete_entities( toDelete );MB_CHK_SET_ERR( rval, "Can't delete entities" );
488 rval = mb->remove_entities( cset, toDelete );MB_CHK_SET_ERR( rval, "Can't remove entities from base set" );
489 if( 0 == rank )
490 {
491
492 std::cout << "delete edges and faces \n";
493 toDelete.print( std::cout );
494 #ifndef _MSC_VER /* windows */
495 std::cout << ( clock() - tt ) / (double)CLOCKS_PER_SEC << " seconds" << endl;
496 tt = clock();
497 #endif
498 }
499 }
500 // do some ghosting if required
501 if( GL > 0 )
502 {
503 rval = pc->exchange_ghost_cells( 3, // int ghost_dim
504 0, // int bridge_dim
505 GL, // int num_layers
506 0, // int addl_ents
507 true );MB_CHK_ERR( rval ); // bool store_remote_handles
508 #ifndef _MSC_VER /* windows */
509 if( 0 == rank )
510 {
511 std::cout << "exchange " << GL << " ghost layer(s) :" << ( clock() - tt ) / (double)CLOCKS_PER_SEC
512 << " seconds" << endl;
513 tt = clock();
514 }
515 #endif
516 }
517 }
518 #endif
519
520 if( !parmerge )
521 {
522 rval = mb->tag_delete( new_id_tag );MB_CHK_SET_ERR( rval, "Can't delete new ID tag" );
523 }
524
525 return MB_SUCCESS;
526 }
References moab::MeshGeneration::BrickOpts::A, moab::Interface::add_entities(), moab::MeshGeneration::BrickOpts::adjEnts, moab::MeshGeneration::BrickOpts::B, moab::Range::begin(), moab::MeshGeneration::BrickOpts::blockSize, moab::MeshGeneration::BrickOpts::C, CC, moab::Interface::create_meshset(), cset, moab::E, moab::Range::end(), ErrorCode, moab::F, moab::Interface::get_adjacencies(), moab::Interface::get_entities_by_dimension(), moab::MeshGeneration::BrickOpts::GL, iface, moab::Range::insert(), moab::MeshGeneration::BrickOpts::K, moab::MeshGeneration::BrickOpts::keep_skins, moab::MeshGeneration::BrickOpts::M, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TAG_SPARSE, MB_TYPE_INTEGER, MB_TYPE_OPAQUE, MBHEX, MBTET, moab::Range::merge(), moab::ParallelMergeMesh::merge(), moab::MergeMesh::merge_entities(), moab::MergeMesh::merge_using_integer_tag(), MESHSET_SET, moab::MeshGeneration::BrickOpts::N, moab::MeshGeneration::BrickOpts::newMergeMethod, moab::MeshGeneration::BrickOpts::parmerge, moab::Range::print(), moab::MeshGeneration::BrickOpts::quadratic, moab::Interface::query_interface(), moab::Interface::remove_entities(), size, moab::Range::size(), moab::Interface::tag_delete(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), moab::MeshGeneration::BrickOpts::tetra, moab::Interface::UNION, moab::MeshGeneration::BrickOpts::xsize, moab::MeshGeneration::BrickOpts::ysize, and moab::MeshGeneration::BrickOpts::zsize.
Referenced by main().
|
private |
Definition at line 53 of file MeshGeneration.hpp.
Referenced by BrickInstance().
|
private |
Definition at line 49 of file MeshGeneration.hpp.
Referenced by BrickInstance(), and MeshGeneration().