Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
moab::MeshGeneration Class Reference

#include <MeshGeneration.hpp>

+ Collaboration diagram for moab::MeshGeneration:

Classes

struct  BrickOpts
 

Public Member Functions

 MeshGeneration (Interface *mbi, EntityHandle rset=0)
 
virtual ~MeshGeneration ()
 
ErrorCode BrickInstance (BrickOpts &opts)
 

Private Attributes

Interfacemb
 
EntityHandle cset
 

Detailed Description

Definition at line 18 of file MeshGeneration.hpp.

Constructor & Destructor Documentation

◆ MeshGeneration()

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.

◆ ~MeshGeneration()

moab::MeshGeneration::~MeshGeneration ( )
virtual

Definition at line 50 of file MeshGeneration.cpp.

50 {}

Member Function Documentation

◆ BrickInstance()

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

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(), 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(), 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().

Member Data Documentation

◆ cset

EntityHandle moab::MeshGeneration::cset
private

Definition at line 53 of file MeshGeneration.hpp.

Referenced by BrickInstance().

◆ mb

Interface* moab::MeshGeneration::mb
private

Definition at line 49 of file MeshGeneration.hpp.

Referenced by BrickInstance(), and MeshGeneration().


The documentation for this class was generated from the following files: