Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
MeshGeneration.cpp
Go to the documentation of this file.
1 /*
2  * MGen.cpp
3  *
4  */
5 
6 #include "moab/MeshGeneration.hpp"
7 #include "moab/MergeMesh.hpp"
8 #include <iostream>
9 #include <ctime>
10 #include <vector>
11 #ifdef WIN32 /* windows */
12 #include <time.h>
13 #endif
14 
15 #ifdef MOAB_HAVE_MPI
16 #include "moab_mpi.h"
17 #include "moab/ParallelComm.hpp"
18 #include "MBParallelConventions.h"
19 #include "moab/ParallelMergeMesh.hpp"
20 #endif
21 #include "moab/ReadUtilIface.hpp"
22 
23 using std::endl;
24 using std::string;
25 using std::vector;
26 
27 namespace moab
28 {
29 
30 MeshGeneration::MeshGeneration( Interface* impl,
31 #ifdef MOAB_HAVE_MPI
32  ParallelComm* comm,
33 #endif
34  EntityHandle rset )
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 }
49 
50 MeshGeneration::~MeshGeneration() {}
51 
52 ErrorCode MeshGeneration::BrickInstance( MeshGeneration::BrickOpts& opts )
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 }
527 
528 } /* namespace moab */