MOAB: Mesh Oriented datABase  (version 5.5.0)
urefine_mesh_test.cpp
Go to the documentation of this file.
1 /*This unit test is for the uniform refinement capability based on AHF datastructures*/
2 #include <iostream>
3 #include <string>
4 #include <sstream>
5 #if defined( __MINGW32__ )
6 #include <sys/time.h>
7 #else
8 #include <ctime>
9 #endif
10 
11 #include <vector>
12 #include <algorithm>
13 #include "moab/Core.hpp"
14 #include "moab/Range.hpp"
15 #include "moab/MeshTopoUtil.hpp"
16 #include "moab/ReadUtilIface.hpp"
17 #include "moab/NestedRefine.hpp"
18 #include "TestUtil.hpp"
19 
20 #ifdef MOAB_HAVE_MPI
21 #include "moab/ParallelComm.hpp"
22 #include "MBParallelConventions.h"
23 #include "ReadParallel.hpp"
24 #include "moab/FileOptions.hpp"
25 #include "MBTagConventions.hpp"
26 #include "moab_mpi.h"
27 #endif
28 
29 using namespace moab;
30 
31 #ifdef MOAB_HAVE_MPI
32 std::string read_options;
33 #endif
34 
35 int number_tests_successful = 0;
36 int number_tests_failed = 0;
37 
38 void handle_error_code( ErrorCode rv, int& number_failed, int& number_successful )
39 {
40  if( rv == MB_SUCCESS )
41  {
42 #ifdef MOAB_HAVE_MPI
43  int rank = 0;
44  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
45  if( rank == 0 ) std::cout << "Success";
46 #else
47  std::cout << "Success";
48 #endif
49  number_successful++;
50  }
51  else
52  {
53  std::cout << "Failure";
54  number_failed++;
55  }
56 }
57 
58 ErrorCode test_adjacencies( Interface* mbImpl, NestedRefine* nr, const Range& all_ents )
59 {
60  MeshTopoUtil mtu( mbImpl );
61  ErrorCode error;
62  Range verts, edges, faces, cells;
63  verts = all_ents.subset_by_dimension( 0 );
64  edges = all_ents.subset_by_dimension( 1 );
65  faces = all_ents.subset_by_dimension( 2 );
66  cells = all_ents.subset_by_dimension( 3 );
67 
68  std::vector< EntityHandle > adjents;
69  Range mbents, ahfents;
70 
71  // Update the moab native data structures
72  ReadUtilIface* read_face;
73  error = mbImpl->query_interface( read_face );CHECK_ERR( error );
74  std::vector< EntityHandle > ents, conn;
75 
76  if( !edges.empty() )
77  {
78  conn.clear();
79  ents.clear();
80 
81  for( Range::iterator it = edges.begin(); it != edges.end(); it++ )
82  ents.push_back( *it );
83  // std::copy(edges.begin(), edges.end(), ents.begin());
84  error = mbImpl->get_connectivity( &ents[0], (int)ents.size(), conn );CHECK_ERR( error );
85  error = read_face->update_adjacencies( *ents.begin(), (int)ents.size(), 2, &conn[0] );CHECK_ERR( error );
86  }
87 
88  if( !faces.empty() )
89  {
90  conn.clear();
91  ents.clear();
92 
93  for( Range::iterator it = faces.begin(); it != faces.end(); it++ )
94  ents.push_back( *it );
95 
96  // std::copy(faces.begin(), faces.end(), ents.begin());
97  error = mbImpl->get_connectivity( &ents[0], 1, conn );CHECK_ERR( error );
98  int nvF = conn.size();
99  conn.clear();
100  error = mbImpl->get_connectivity( &ents[0], (int)ents.size(), conn );CHECK_ERR( error );
101  error = read_face->update_adjacencies( *ents.begin(), (int)ents.size(), nvF, &conn[0] );CHECK_ERR( error );
102  }
103 
104  if( !cells.empty() )
105  {
106  conn.clear();
107  ents.clear();
108 
109  for( Range::iterator it = cells.begin(); it != cells.end(); it++ )
110  ents.push_back( *it );
111  // std::copy(cells.begin(), cells.end(), ents.begin());
112  error = mbImpl->get_connectivity( &ents[0], 1, conn );CHECK_ERR( error );
113  int nvF = conn.size();
114  conn.clear();
115  error = mbImpl->get_connectivity( &ents[0], (int)ents.size(), conn );CHECK_ERR( error );
116  error = read_face->update_adjacencies( *ents.begin(), (int)ents.size(), nvF, &conn[0] );CHECK_ERR( error );
117  }
118 
119  if( !edges.empty() )
120  {
121  // 1D Queries //
122  // IQ1: For every vertex, obtain incident edges
123  for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
124  {
125  adjents.clear();
126  mbents.clear();
127  ahfents.clear();
128  error = nr->get_adjacencies( *i, 1, adjents );CHECK_ERR( error );
129  error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error );
130  std::sort( adjents.begin(), adjents.end() );
131  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
132 
133  if( ahfents.size() != mbents.size() )
134  {
135  std::cout << "ahf results = " << std::endl;
136  ahfents.print();
137  std::cout << "native results = " << std::endl;
138  mbents.print();
139  }
140 
141  CHECK_EQUAL( adjents.size(), mbents.size() );
142  mbents = subtract( mbents, ahfents );
143  if( ahfents.size() != mbents.size() )
144  {
145  // std::cout<<"ahf results = "<<std::endl;
146  // ahfents.print();
147  // std::cout<<"native results = "<<std::endl;
148  // mbents.print();
149  }
150  // CHECK_EQUAL(adjents.size(),mbents.size());
151  CHECK( !mbents.size() );
152  }
153 
154  // NQ1: For every edge, obtain neighbor edges
155  for( Range::iterator i = edges.begin(); i != edges.end(); ++i )
156  {
157  adjents.clear();
158  mbents.clear();
159  ahfents.clear();
160  error = nr->get_adjacencies( *i, 1, adjents );CHECK_ERR( error );
161  error = mtu.get_bridge_adjacencies( *i, 0, 1, mbents );CHECK_ERR( error );
162  CHECK_EQUAL( adjents.size(), mbents.size() );
163  std::sort( adjents.begin(), adjents.end() );
164  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
165  mbents = subtract( mbents, ahfents );
166  CHECK( !mbents.size() );
167  }
168  }
169 
170  if( !faces.empty() )
171  {
172  // IQ21: For every vertex, obtain incident faces
173  for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
174  {
175  adjents.clear();
176  mbents.clear();
177  ahfents.clear();
178  error = nr->get_adjacencies( *i, 2, adjents );CHECK_ERR( error );
179  error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error );
180  CHECK_EQUAL( adjents.size(), mbents.size() );
181  std::sort( adjents.begin(), adjents.end() );
182  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
183  mbents = subtract( mbents, ahfents );
184  CHECK( !mbents.size() );
185  }
186 
187  // IQ22: For every edge, obtain incident faces
188  if( !edges.empty() )
189  {
190  for( Range::iterator i = edges.begin(); i != edges.end(); ++i )
191  {
192  adjents.clear();
193  mbents.clear();
194  ahfents.clear();
195  error = nr->get_adjacencies( *i, 2, adjents );CHECK_ERR( error );
196  error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error );
197  CHECK_EQUAL( adjents.size(), mbents.size() );
198  std::sort( adjents.begin(), adjents.end() );
199  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
200  mbents = subtract( mbents, ahfents );
201  CHECK( !mbents.size() );
202  }
203  }
204 
205  // NQ2: For every face, obtain neighbor faces
206  for( Range::iterator i = faces.begin(); i != faces.end(); ++i )
207  {
208  adjents.clear();
209  mbents.clear();
210  ahfents.clear();
211  error = nr->get_adjacencies( *i, 2, adjents );CHECK_ERR( error );
212  error = mtu.get_bridge_adjacencies( *i, 1, 2, mbents );CHECK_ERR( error );
213  CHECK_EQUAL( adjents.size(), mbents.size() );
214  std::sort( adjents.begin(), adjents.end() );
215  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
216  mbents = subtract( mbents, ahfents );
217  CHECK( !mbents.size() );
218  }
219 
220  if( !edges.empty() )
221  {
222  for( Range::iterator i = faces.begin(); i != faces.end(); ++i )
223  {
224  adjents.clear();
225  mbents.clear();
226  ahfents.clear();
227  error = nr->get_adjacencies( *i, 1, adjents );CHECK_ERR( error );
228  error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error );
229  CHECK_EQUAL( adjents.size(), mbents.size() );
230  std::sort( adjents.begin(), adjents.end() );
231  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
232  mbents = subtract( mbents, ahfents );
233  CHECK( !mbents.size() );
234  }
235  }
236  }
237 
238  if( !cells.empty() )
239  {
240  // IQ 31: For every vertex, obtain incident cells
241  for( Range::iterator i = verts.begin(); i != verts.end(); ++i )
242  {
243  adjents.clear();
244  mbents.clear();
245  ahfents.clear();
246  error = nr->get_adjacencies( *i, 3, adjents );CHECK_ERR( error );
247  error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error );
248 
249  if( adjents.size() != mbents.size() )
250  {
251  // std::cout<<"ahf results = "<<std::endl;
252  // ahfents.print();
253  // std::cout<<"native results = "<<std::endl;
254  // mbents.print();
255  }
256 
257  CHECK_EQUAL( adjents.size(), mbents.size() );
258  std::sort( adjents.begin(), adjents.end() );
259  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
260  mbents = subtract( mbents, ahfents );
261  CHECK( !mbents.size() );
262  }
263 
264  if( !edges.empty() )
265  {
266  for( Range::iterator i = edges.begin(); i != edges.end(); ++i )
267  {
268  adjents.clear();
269  mbents.clear();
270  ahfents.clear();
271  error = nr->get_adjacencies( *i, 3, adjents );CHECK_ERR( error );
272  error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error );
273  CHECK_EQUAL( adjents.size(), mbents.size() );
274  std::sort( adjents.begin(), adjents.end() );
275  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
276  mbents = subtract( mbents, ahfents );
277  CHECK( !mbents.size() );
278  }
279  }
280 
281  if( !faces.empty() )
282  {
283  for( Range::iterator i = faces.begin(); i != faces.end(); ++i )
284  {
285  adjents.clear();
286  mbents.clear();
287  ahfents.clear();
288  error = nr->get_adjacencies( *i, 3, adjents );CHECK_ERR( error );
289  error = mbImpl->get_adjacencies( &*i, 1, 3, false, mbents );CHECK_ERR( error );
290  CHECK_EQUAL( adjents.size(), mbents.size() );
291  std::sort( adjents.begin(), adjents.end() );
292  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
293  mbents = subtract( mbents, ahfents );
294  CHECK( !mbents.size() );
295  }
296  }
297 
298  // NQ3: For every cell, obtain neighbor cells
299  for( Range::iterator i = cells.begin(); i != cells.end(); ++i )
300  {
301  adjents.clear();
302  mbents.clear();
303  ahfents.clear();
304  error = nr->get_adjacencies( *i, 3, adjents );CHECK_ERR( error );
305  error = mtu.get_bridge_adjacencies( *i, 2, 3, mbents );CHECK_ERR( error );
306  CHECK_EQUAL( adjents.size(), mbents.size() );
307  std::sort( adjents.begin(), adjents.end() );
308  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
309  mbents = subtract( mbents, ahfents );
310  CHECK( !mbents.size() );
311  }
312 
313  if( !edges.empty() )
314  {
315  for( Range::iterator i = cells.begin(); i != cells.end(); ++i )
316  {
317  adjents.clear();
318  mbents.clear();
319  ahfents.clear();
320  error = nr->get_adjacencies( *i, 1, adjents );CHECK_ERR( error );
321  error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents );CHECK_ERR( error );
322  CHECK_EQUAL( adjents.size(), mbents.size() );
323  std::sort( adjents.begin(), adjents.end() );
324  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
325  mbents = subtract( mbents, ahfents );
326  CHECK( !mbents.size() );
327  }
328  }
329 
330  if( !faces.empty() )
331  {
332  for( Range::iterator i = cells.begin(); i != cells.end(); ++i )
333  {
334  adjents.clear();
335  mbents.clear();
336  ahfents.clear();
337  error = nr->get_adjacencies( *i, 2, adjents );CHECK_ERR( error );
338  error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents );CHECK_ERR( error );
339  CHECK_EQUAL( adjents.size(), mbents.size() );
340  std::sort( adjents.begin(), adjents.end() );
341  std::copy( adjents.begin(), adjents.end(), range_inserter( ahfents ) );
342  mbents = subtract( mbents, ahfents );
343  CHECK( !mbents.size() );
344  }
345  }
346  }
347 
348  return MB_SUCCESS;
349 }
350 
351 ErrorCode refine_entities( Interface* mb,
352  ParallelComm* pc,
353  EntityHandle fset,
354  int* level_degrees,
355  const int num_levels,
356  bool output )
357 {
358  ErrorCode error;
359 
360  // Get the range of entities in the initial mesh
361  Range init_ents[4];
362 
363  int dim[3] = { 1, 2, 3 };
364 
365  if( output )
366  {
367  // int inents = init_ents.size();
368  std::stringstream file;
369 #ifdef MOAB_HAVE_MPI
370  file << "MESH_LEVEL_0.h5m";
371 #else
372  file << "MESH_LEVEL_0.vtk";
373 #endif
374  std::string str = file.str();
375  const char* output_file = str.c_str();
376 #ifdef MOAB_HAVE_MPI
377  error = mb->write_file( output_file, 0, ";;PARALLEL=WRITE_PART" );CHECK_ERR( error );
378 #else
379  error = mb->write_file( output_file, 0, NULL );CHECK_ERR( error );
380 #endif
381  }
382 
383  // Create an hm object and generate the hierarchy
384  std::cout << "Creating a hm object" << std::endl;
385 
386 #ifdef MOAB_HAVE_MPI
387  Range averts, aedges, afaces, acells;
388  error = mb->get_entities_by_dimension( fset, 0, averts );MB_CHK_ERR( error );
389  error = mb->get_entities_by_dimension( fset, 1, aedges );MB_CHK_ERR( error );
390  error = mb->get_entities_by_dimension( fset, 2, afaces );MB_CHK_ERR( error );
391  error = mb->get_entities_by_dimension( fset, 3, acells );MB_CHK_ERR( error );
392 
393  /* filter based on parallel status */
394  if( pc )
395  {
396  error = pc->filter_pstatus( averts, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[0] );MB_CHK_ERR( error );
397  error = pc->filter_pstatus( aedges, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[1] );MB_CHK_ERR( error );
398  error = pc->filter_pstatus( afaces, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[2] );MB_CHK_ERR( error );
399  error = pc->filter_pstatus( acells, PSTATUS_GHOST, PSTATUS_NOT, -1, &init_ents[3] );MB_CHK_ERR( error );
400  }
401  else
402  {
403  init_ents[0] = averts;
404  init_ents[1] = aedges;
405  init_ents[2] = afaces;
406  init_ents[3] = acells;
407  }
408 #else
409  error = mb->get_entities_by_dimension( fset, 0, init_ents[0] );CHECK_ERR( error );
410  error = mb->get_entities_by_dimension( fset, 1, init_ents[1] );CHECK_ERR( error );
411  error = mb->get_entities_by_dimension( fset, 2, init_ents[2] );CHECK_ERR( error );
412  error = mb->get_entities_by_dimension( fset, 3, init_ents[3] );CHECK_ERR( error );
413 #endif
414 
415  NestedRefine uref( dynamic_cast< Core* >( mb ), pc, fset );
416  std::vector< EntityHandle > set;
417 
418  std::cout << "Starting hierarchy generation" << std::endl;
419  bool opt = true;
420  // bool opt = false;
421  error = uref.generate_mesh_hierarchy( num_levels, level_degrees, set, opt );CHECK_ERR( error );
422  std::cout << "Finished hierarchy generation in " << uref.timeall.tm_total << " secs" << std::endl;
423 #ifdef MOAB_HAVE_MPI
424  if( pc )
425  {
426  std::cout << "Time taken for refinement " << uref.timeall.tm_refine << " secs" << std::endl;
427  std::cout << "Time taken for resolving shared interface " << uref.timeall.tm_resolve << " secs" << std::endl;
428  }
429 #endif
430 
431  // error = uref.exchange_ghosts(set, 1); CHECK_ERR(error);
432 
433  std::cout << std::endl;
434  std::cout << "Mesh size for level 0 :: inverts = " << init_ents[0].size() << ", inedges = " << init_ents[1].size()
435  << ", infaces = " << init_ents[2].size() << ", incells = " << init_ents[3].size() << std::endl;
436 
437  Range prev_ents[4];
438  for( int i = 0; i < 4; i++ )
439  prev_ents[i] = init_ents[i];
440 
441  // Loop over each mesh level and check its topological properties
442  for( int l = 0; l < num_levels; l++ )
443  {
444  Range all_ents;
445  error = mb->get_entities_by_handle( set[l + 1], all_ents );CHECK_ERR( error );
446 
447  Range ents[4];
448  for( int k = 0; k < 4; k++ )
449  ents[k] = all_ents.subset_by_dimension( k );
450 
451  if( ents[0].empty() || all_ents.empty() ) std::cout << "Something is not right" << std::endl;
452 
453  std::cout << std::endl;
454  std::cout << "Mesh size for level " << l + 1 << " :: nverts = " << ents[0].size()
455  << ", nedges = " << ents[1].size() << ", nfaces = " << ents[2].size()
456  << ", ncells = " << ents[3].size() << std::endl;
457 
458  // Check if the number of new entities created are correct.
459 
460  for( int type = 1; type < 3; type++ )
461  {
462  int factor = 1;
463  if( !ents[type + 1].empty() )
464  {
465 
466  for( int p = 0; p <= l; p++ )
467  {
468  for( int d = 0; d < dim[type]; d++ )
469  factor *= level_degrees[p];
470  }
471  int expected_nents = factor * init_ents[type + 1].size();
472  CHECK_EQUAL( expected_nents, (int)ents[type + 1].size() );
473  }
474  }
475 
476  // Check adjacencies
477  error = test_adjacencies( mb, &uref, all_ents );CHECK_ERR( error );
478 
479  // Check interlevel child-parent query between previous and current level
480  for( int type = 1; type < 3; type++ )
481  {
482  if( !prev_ents[type + 1].empty() )
483  {
484  for( Range::iterator e = prev_ents[type + 1].begin(); e != prev_ents[type + 1].end(); e++ )
485  {
486  std::vector< EntityHandle > children;
487  error = uref.parent_to_child( *e, l, l + 1, children );CHECK_ERR( error );
488  for( int i = 0; i < (int)children.size(); i++ )
489  {
490  EntityHandle parent;
491  error = uref.child_to_parent( children[i], l + 1, l, &parent );CHECK_ERR( error );
492  assert( parent == *e );
493  }
494  }
495  }
496  }
497 
498  for( int i = 0; i < 4; i++ )
499  prev_ents[i] = ents[i];
500 
501  // Print out the boundary vertices
502  /* EntityHandle bnd_set;
503  error = mb->create_meshset(MESHSET_SET, bnd_set); MB_CHK_ERR(error);
504  std::vector<EntityHandle> vbnd;
505 
506  for (int k=0; k<3; k++){
507  std::cout<<"Finding boundary of dimension "<<k<<" with size
508  "<<ents[k].size()<<std::endl; if (ents[k].size() != 0){ for (Range::iterator v =
509  ents[k].begin(); v != ents[k].end(); v++)
510  {
511  EntityHandle ent = *v;
512  bool bnd = uref.is_entity_on_boundary(ent);
513  if (bnd)
514  vbnd.push_back(*v);
515  }
516  }
517  }
518 
519  std::cout<<"vbnd.size = "<<vbnd.size()<<std::endl;
520  error = mb->add_entities(bnd_set, &vbnd[0], (int)vbnd.size()); MB_CHK_ERR(error);
521  if (output)
522  {
523  std::stringstream file;
524  file << "VBND_LEVEL_" <<l+1<<".vtk";
525  std::string str = file.str();
526  const char* output_file = str.c_str();
527  char * write_opts = NULL;
528  error = mb->write_file(output_file, 0, write_opts, &bnd_set, 1); CHECK_ERR(error);
529  }*/
530 
531  // Print out the mesh
532  if( output )
533  {
534  std::stringstream file;
535 
536 #ifdef MOAB_HAVE_MPI
537  file << "MESH_LEVEL_" << l + 1 << ".h5m";
538 #else
539  file << "MESH_LEVEL_" << l + 1 << ".vtk";
540 #endif
541  std::string str = file.str();
542  const char* output_file = str.c_str();
543 #ifdef MOAB_HAVE_MPI
544  error = mb->write_file( output_file, 0, ";;PARALLEL=WRITE_PART", &set[l + 1], 1 );CHECK_ERR( error );
545 #else
546  error = mb->write_file( output_file, 0, NULL, &set[l + 1], 1 );CHECK_ERR( error );
547 #endif
548  }
549  }
550 
551  // Check interlevel child-parent query between initial and most refined mesh
552  for( int type = 1; type < 3; type++ )
553  {
554  if( !init_ents[type + 1].empty() )
555  {
556  for( Range::iterator e = init_ents[type + 1].begin(); e != init_ents[type + 1].end(); e++ )
557  {
558  std::vector< EntityHandle > children;
559  error = uref.parent_to_child( *e, 0, num_levels, children );CHECK_ERR( error );
560  for( int i = 0; i < (int)children.size(); i++ )
561  {
562  EntityHandle parent;
563  error = uref.child_to_parent( children[i], num_levels, 0, &parent );CHECK_ERR( error );
564  assert( parent == *e );
565  }
566  }
567  }
568  }
569 
570  // Print out the whole hierarchy into a single file
571  /* if (output)
572  {
573  std::stringstream file;
574  file << "MESH_HIERARCHY.vtk";
575  std::string str = file.str();
576  const char* output_file = str.c_str();
577  error = mb->write_file(output_file); CHECK_ERR(error);
578  }*/
579 
580  return MB_SUCCESS;
581 }
582 
583 ErrorCode create_single_entity( Interface* mbImpl, EntityType type )
584 {
585  ErrorCode error;
586  if( type == MBEDGE )
587  {
588  const double coords[] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 };
589  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
590 
591  const int conn[] = { 0, 1 };
592  const size_t num_elems = sizeof( conn ) / sizeof( conn[0] ) / 2;
593 
594  std::cout << "Specify verts and ents" << std::endl;
595 
596  EntityHandle verts[num_vtx], edges[num_elems];
597  for( size_t i = 0; i < num_vtx; ++i )
598  {
599  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
600  }
601 
602  std::cout << "Created vertices" << std::endl;
603 
604  for( size_t i = 0; i < num_elems; ++i )
605  {
606  EntityHandle c[2];
607  c[0] = verts[conn[0]];
608  c[1] = verts[conn[1]];
609 
610  error = mbImpl->create_element( MBEDGE, c, 2, edges[i] );CHECK_ERR( error );
611  }
612 
613  std::cout << "Created ents" << std::endl;
614  }
615  else if( type == MBTRI )
616  {
617  const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0 };
618  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
619 
620  const int conn[] = { 0, 1, 2 };
621  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3;
622 
623  EntityHandle verts[num_vtx], faces[num_elems];
624  for( size_t i = 0; i < num_vtx; ++i )
625  {
626  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );
627  if( error != MB_SUCCESS ) return error;
628  }
629 
630  for( size_t i = 0; i < num_elems; ++i )
631  {
632  EntityHandle c[3];
633  for( int j = 0; j < 3; j++ )
634  c[j] = verts[conn[3 * i + j]];
635 
636  error = mbImpl->create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error );
637  }
638  }
639  else if( type == MBQUAD )
640  {
641  const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0 };
642  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
643 
644  const int conn[] = { 0, 1, 2, 3 };
645  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3;
646 
647  EntityHandle verts[num_vtx], faces[num_elems];
648  for( size_t i = 0; i < num_vtx; ++i )
649  {
650  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
651  }
652 
653  for( size_t i = 0; i < num_elems; ++i )
654  {
655  EntityHandle c[4];
656  for( int j = 0; j < 4; j++ )
657  c[j] = verts[conn[j]];
658 
659  error = mbImpl->create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error );
660  }
661  }
662  else if( type == MBTET )
663  {
664  const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 };
665 
666  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
667 
668  const int conn[] = { 0, 1, 2, 3 };
669  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
670 
671  EntityHandle verts[num_vtx], cells[num_elems];
672  for( size_t i = 0; i < num_vtx; ++i )
673  {
674  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
675  }
676 
677  for( size_t i = 0; i < num_elems; ++i )
678  {
679  EntityHandle c[4];
680  for( int j = 0; j < 4; j++ )
681  c[j] = verts[conn[j]];
682 
683  error = mbImpl->create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error );
684  }
685  }
686  else if( type == MBHEX )
687  {
688  const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1 };
689  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
690 
691  const int conn[] = { 0, 1, 2, 3, 4, 5, 6, 7 };
692  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8;
693 
694  EntityHandle verts[num_vtx], cells[num_elems];
695  for( size_t i = 0; i < num_vtx; ++i )
696  {
697  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
698  }
699 
700  for( size_t i = 0; i < num_elems; ++i )
701  {
702  EntityHandle c[8];
703  for( int j = 0; j < 8; j++ )
704  c[j] = verts[conn[j]];
705 
706  error = mbImpl->create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error );
707  }
708  }
709  return MB_SUCCESS;
710 }
711 
712 ErrorCode create_mesh( Interface* mbImpl, EntityType type )
713 {
714  ErrorCode error;
715  if( type == MBEDGE )
716  {
717  const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0 };
718  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
719 
720  const int conn[] = { 1, 0, 0, 3, 2, 0, 0, 4 };
721  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 2;
722 
723  EntityHandle verts[num_vtx], edges[num_elems];
724  for( size_t i = 0; i < num_vtx; ++i )
725  {
726  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
727  }
728 
729  for( size_t i = 0; i < num_elems; ++i )
730  {
731  EntityHandle c[2];
732  c[0] = verts[conn[2 * i]];
733  c[1] = verts[conn[2 * i + 1]];
734 
735  error = mbImpl->create_element( MBEDGE, c, 2, edges[i] );CHECK_ERR( error );
736  }
737  }
738  else if( type == MBTRI )
739  {
740  const double coords[] = { 0, 0, 0, 1, -1, 0, 1, 1, 0, -1, 1, 0, -1, -1, 0, 0, 0, 1 };
741 
742  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
743 
744  const int conn[] = { 0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 1, 0, 5, 3, 0, 2, 5, 0, 4, 5, 0, 5, 1 };
745 
746  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3;
747 
748  EntityHandle verts[num_vtx], faces[num_elems];
749  for( size_t i = 0; i < num_vtx; ++i )
750  {
751  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
752  }
753 
754  for( size_t i = 0; i < num_elems; ++i )
755  {
756  EntityHandle c[3];
757  for( int j = 0; j < 3; j++ )
758  c[j] = verts[conn[3 * i + j]];
759 
760  error = mbImpl->create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error );
761  }
762  }
763  else if( type == MBQUAD )
764  {
765  const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, -1, 1, 0, -1, 0, 0,
766  -1, -1, 0, 0, -1, 0, 1, -1, 0, 0, 1, 1, 0, 0, 1, 0, -1, 1 };
767 
768  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
769 
770  const int conn[] = { 0, 1, 2, 3, 0, 3, 4, 5, 7, 8, 1, 0, 6, 7, 0, 5, 0, 3, 9, 10, 0, 10, 11, 7 };
771 
772  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
773 
774  EntityHandle verts[num_vtx], faces[num_elems];
775  for( size_t i = 0; i < num_vtx; ++i )
776  {
777  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
778  }
779 
780  for( size_t i = 0; i < num_elems; ++i )
781  {
782  EntityHandle c[4];
783  for( int j = 0; j < 4; j++ )
784  c[j] = verts[conn[4 * i + j]];
785 
786  error = mbImpl->create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error );
787  }
788  }
789  else if( type == MBTET )
790  {
791  const double coords[] = { 0, -1, 0, 0, 2, 0, 1, 0, 0, -0.5, 0, 0, 0, 0, 1, 0, 0, -2 };
792 
793  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
794 
795  const int conn[] = { 0, 2, 1, 4, 3, 0, 1, 4, 5, 2, 1, 0, 5, 0, 1, 3 };
796 
797  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
798 
799  EntityHandle verts[num_vtx], cells[num_elems];
800  for( size_t i = 0; i < num_vtx; ++i )
801  {
802  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
803  }
804 
805  for( size_t i = 0; i < num_elems; ++i )
806  {
807  EntityHandle c[4];
808  for( int j = 0; j < 4; j++ )
809  c[j] = verts[conn[4 * i + j]];
810 
811  error = mbImpl->create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error );
812  }
813  }
814  else if( type == MBHEX )
815  {
816  const double coords[] = { 0, -1, 0, 1, -1, 0, 1, 1, 0, 0, 1, 0, -1, 1, 0, -1, -1, 0,
817  0, -1, 1, 1, -1, 1, 1, 1, 1, 0, 1, 1, -1, 1, 1, -1, -1, 1,
818  0, -1, -1, 1, -1, -1, 1, 1, -1, 0, 1, -1, -1, 1, -1, -1, -1, -1 };
819  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
820 
821  const int conn[] = { 0, 1, 2, 3, 6, 7, 8, 9, 5, 0, 3, 4, 11, 6, 9, 10,
822  12, 13, 14, 15, 0, 1, 2, 3, 17, 12, 15, 16, 5, 0, 3, 4 };
823  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8;
824 
825  EntityHandle verts[num_vtx], cells[num_elems];
826  for( size_t i = 0; i < num_vtx; ++i )
827  {
828  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
829  }
830 
831  for( size_t i = 0; i < num_elems; ++i )
832  {
833  EntityHandle c[8];
834  for( int j = 0; j < 8; j++ )
835  c[j] = verts[conn[8 * i + j]];
836 
837  error = mbImpl->create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error );
838  }
839  }
840  return MB_SUCCESS;
841 }
842 
843 ErrorCode create_simple_mesh( Interface* mbImpl, EntityType type )
844 {
845  ErrorCode error;
846  if( type == MBEDGE )
847  {
848  const double coords[] = { 0, 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 0, 4, 0, 0,
849  4, 1, 0, 3, 1, 0, 2, 1, 0, 1, 1, 0, 0, 1, 0 };
850  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
851 
852  const int conn[] = { 1, 0, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 0 };
853  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 2;
854 
855  EntityHandle verts[num_vtx], edges[num_elems];
856  for( size_t i = 0; i < num_vtx; ++i )
857  {
858  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
859  }
860 
861  for( size_t i = 0; i < num_elems; ++i )
862  {
863  EntityHandle c[2];
864  c[0] = verts[conn[2 * i]];
865  c[1] = verts[conn[2 * i + 1]];
866 
867  error = mbImpl->create_element( MBEDGE, c, 2, edges[i] );CHECK_ERR( error );
868  }
869  }
870  else if( type == MBTRI )
871  {
872  const double coords[] = { 0, 0, 0, 1, 0, 0, 2, 0, 0, 2.5, 1, 0, 1.5, 1, 0, 0.5, 1,
873  0, -0.5, 1, 0, -0.5, -1, 0, 0.5, -1, 0, 1.5, -1, 0, 2.5, -1, 0 };
874 
875  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
876 
877  const int conn[] = { 0, 5, 6, 0, 1, 5, 1, 4, 5, 1, 2, 4, 2, 3, 4,
878  7, 8, 0, 8, 1, 0, 8, 9, 1, 9, 2, 1, 9, 10, 2 };
879 
880  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 3;
881 
882  EntityHandle verts[num_vtx], faces[num_elems];
883  for( size_t i = 0; i < num_vtx; ++i )
884  {
885  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
886  }
887 
888  for( size_t i = 0; i < num_elems; ++i )
889  {
890  EntityHandle c[3];
891  for( int j = 0; j < 3; j++ )
892  c[j] = verts[conn[3 * i + j]];
893 
894  error = mbImpl->create_element( MBTRI, c, 3, faces[i] );CHECK_ERR( error );
895  }
896  }
897  else if( type == MBQUAD )
898  {
899  const double coords[] = { 0, 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 0, 4, 0, 0, 5, 0, 0, 0, 1, 0, 1, 1, 0, 2, 1, 0,
900  3, 1, 0, 4, 1, 0, 5, 1, 0, 0, 2, 0, 1, 2, 0, 2, 2, 0, 3, 2, 0, 4, 2, 0, 5, 2, 0 };
901 
902  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
903 
904  const int conn[] = { 0, 1, 7, 6, 1, 2, 8, 7, 2, 3, 9, 8, 3, 4, 10, 9, 4, 5, 11, 10,
905  6, 7, 13, 12, 7, 8, 14, 13, 8, 9, 15, 14, 9, 10, 16, 15, 10, 11, 17, 16 };
906 
907  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
908 
909  EntityHandle verts[num_vtx], faces[num_elems];
910  for( size_t i = 0; i < num_vtx; ++i )
911  {
912  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
913  }
914 
915  for( size_t i = 0; i < num_elems; ++i )
916  {
917  EntityHandle c[4];
918  for( int j = 0; j < 4; j++ )
919  c[j] = verts[conn[4 * i + j]];
920 
921  error = mbImpl->create_element( MBQUAD, c, 4, faces[i] );CHECK_ERR( error );
922  }
923  }
924  else if( type == MBTET )
925  {
926  const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1 };
927 
928  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
929 
930  const int conn[] = { 0, 1, 2, 5, 3, 0, 2, 5, 4, 1, 0, 5, 4, 0, 3, 5 };
931 
932  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 4;
933 
934  EntityHandle verts[num_vtx], cells[num_elems];
935  for( size_t i = 0; i < num_vtx; ++i )
936  {
937  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
938  }
939 
940  for( size_t i = 0; i < num_elems; ++i )
941  {
942  EntityHandle c[4];
943  for( int j = 0; j < 4; j++ )
944  c[j] = verts[conn[4 * i + j]];
945 
946  error = mbImpl->create_element( MBTET, c, 4, cells[i] );CHECK_ERR( error );
947  }
948  }
949  else if( type == MBHEX )
950  {
951  const double coords[] = { 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 1, 0, 1, 1, 0, 2, 1, 0, 0, 2, 0, 1, 2, 0, 2, 2, 0,
952  0, 0, 1, 1, 0, 1, 2, 0, 1, 0, 1, 1, 1, 1, 1, 2, 1, 1, 0, 2, 1, 1, 2, 1, 2, 2, 1 };
953  const size_t num_vtx = sizeof( coords ) / sizeof( double ) / 3;
954 
955  const int conn[] = { 0, 1, 4, 3, 9, 10, 13, 12, 1, 2, 5, 4, 10, 11, 14, 13,
956  3, 4, 7, 6, 12, 13, 16, 15, 4, 5, 8, 7, 13, 14, 17, 16 };
957  const size_t num_elems = sizeof( conn ) / sizeof( int ) / 8;
958 
959  EntityHandle verts[num_vtx], cells[num_elems];
960  for( size_t i = 0; i < num_vtx; ++i )
961  {
962  error = mbImpl->create_vertex( coords + 3 * i, verts[i] );CHECK_ERR( error );
963  }
964 
965  for( size_t i = 0; i < num_elems; ++i )
966  {
967  EntityHandle c[8];
968  for( int j = 0; j < 8; j++ )
969  c[j] = verts[conn[8 * i + j]];
970 
971  error = mbImpl->create_element( MBHEX, c, 8, cells[i] );CHECK_ERR( error );
972  }
973  }
974  return MB_SUCCESS;
975 }
976 
977 ErrorCode test_entities( int mesh_type, EntityType type, int* level_degrees, int num_levels, bool output )
978 {
979  ErrorCode error;
980  Core mb;
981  Interface* mbimpl = &mb;
982 
983  // Create entities
984  if( mesh_type == 1 )
985  {
986  error = create_single_entity( mbimpl, type );
987  if( error != MB_SUCCESS ) return error;
988  std::cout << "Entity created successfully" << std::endl;
989  }
990  else if( mesh_type == 2 )
991  {
992  error = create_mesh( mbimpl, type );
993  if( error != MB_SUCCESS ) return error;
994  std::cout << "Small mesh created successfully" << std::endl;
995  }
996  else if( mesh_type == 3 )
997  {
998  error = create_simple_mesh( mbimpl, type );
999  if( error != MB_SUCCESS ) return error;
1000  std::cout << "Small simple mesh created successfully" << std::endl;
1001  }
1002 
1003  // Generate hierarchy
1004  error = refine_entities( mbimpl, 0, 0, level_degrees, num_levels, output );
1005  if( error != MB_SUCCESS ) return error;
1006 
1007  return MB_SUCCESS;
1008 }
1009 ErrorCode test_1D()
1010 {
1011  ErrorCode error;
1012 
1013  std::cout << "Testing EDGE" << std::endl;
1014  EntityType type = MBEDGE;
1015 
1016  std::cout << "Testing single entity" << std::endl;
1017  int deg[3] = { 2, 3, 5 };
1018  int len = sizeof( deg ) / sizeof( int );
1019  error = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
1020 
1021  std::cout << std::endl;
1022  std::cout << "Testing a small mesh" << std::endl;
1023  error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
1024 
1025  std::cout << std::endl;
1026  std::cout << "Testing a small simple mesh" << std::endl;
1027  int degree[4] = { 5, 5, 2, 2 };
1028  len = sizeof( degree ) / sizeof( int );
1029  error = test_entities( 3, type, degree, len, false );CHECK_ERR( error );
1030 
1031  return MB_SUCCESS;
1032 }
1033 
1034 ErrorCode test_2D()
1035 {
1036  ErrorCode error;
1037 
1038  std::cout << "Testing TRI" << std::endl;
1039  EntityType type = MBTRI;
1040 
1041  std::cout << "Testing single entity" << std::endl;
1042  int deg[3] = { 2, 3, 5 };
1043  int len = sizeof( deg ) / sizeof( int );
1044  error = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
1045 
1046  std::cout << std::endl;
1047  std::cout << "Testing a small mesh" << std::endl;
1048  error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
1049 
1050  std::cout << std::endl;
1051  std::cout << "Testing a small simple mesh" << std::endl;
1052  int degree[2] = { 5, 2 };
1053  int length = sizeof( degree ) / sizeof( int );
1054  error = test_entities( 3, type, degree, length, false );CHECK_ERR( error );
1055 
1056  std::cout << std::endl;
1057  std::cout << "Testing QUAD" << std::endl;
1058  type = MBQUAD;
1059 
1060  std::cout << "Testing single entity" << std::endl;
1061  error = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
1062 
1063  std::cout << std::endl;
1064  std::cout << "Testing a small mesh" << std::endl;
1065  error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
1066 
1067  std::cout << std::endl;
1068  std::cout << "Testing a small simple mesh" << std::endl;
1069  error = test_entities( 3, type, degree, length, false );CHECK_ERR( error );
1070 
1071  return MB_SUCCESS;
1072 }
1073 
1074 ErrorCode test_3D()
1075 {
1076  ErrorCode error;
1077 
1078  std::cout << "Testing TET" << std::endl;
1079  EntityType type = MBTET;
1080  int deg[2] = { 2, 3 };
1081  int len = sizeof( deg ) / sizeof( int );
1082 
1083  std::cout << "Testing single entity" << std::endl;
1084  error = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
1085 
1086  std::cout << std::endl;
1087  std::cout << "Testing a small mesh" << std::endl;
1088  error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
1089 
1090  std::cout << std::endl;
1091  std::cout << "Testing a small simple mesh" << std::endl;
1092  int degree[4] = { 2, 2, 2, 2 };
1093  int length = sizeof( degree ) / sizeof( int );
1094  error = test_entities( 3, type, degree, length, false );CHECK_ERR( error );
1095 
1096  std::cout << std::endl;
1097  std::cout << "Testing HEX" << std::endl;
1098  type = MBHEX;
1099 
1100  std::cout << "Testing single entity" << std::endl;
1101  error = test_entities( 1, type, deg, len, false );CHECK_ERR( error );
1102 
1103  std::cout << std::endl;
1104  std::cout << "Testing a small mesh" << std::endl;
1105  error = test_entities( 2, type, deg, len, false );CHECK_ERR( error );
1106 
1107  std::cout << std::endl;
1108  std::cout << "Testing a small simple mesh" << std::endl;
1109  error = test_entities( 3, type, degree, length, false );CHECK_ERR( error );
1110 
1111  return MB_SUCCESS;
1112 }
1113 
1114 ErrorCode test_mesh( const char* filename, int* level_degrees, int num_levels )
1115 {
1116  Core moab;
1117  Interface* mbImpl = &moab;
1118  ParallelComm* pc = NULL;
1119  EntityHandle fileset;
1120 
1121  ErrorCode error;
1122  error = mbImpl->create_meshset( moab::MESHSET_SET, fileset );CHECK_ERR( error );
1123 
1124 #ifdef MOAB_HAVE_MPI
1125  MPI_Comm comm = MPI_COMM_WORLD;
1126  EntityHandle partnset;
1127  error = mbImpl->create_meshset( moab::MESHSET_SET, partnset );MB_CHK_ERR( error );
1128  pc = moab::ParallelComm::get_pcomm( mbImpl, partnset, &comm );
1129 
1130  int procs = 1;
1131  MPI_Comm_size( MPI_COMM_WORLD, &procs );
1132 
1133  if( procs > 1 )
1134  {
1135  read_options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS";
1136 
1137  error = mbImpl->load_file( filename, &fileset, read_options.c_str() );CHECK_ERR( error );
1138 
1139  // DBG
1140  std::set< unsigned int > shprocs;
1141  error = pc->get_comm_procs( shprocs );CHECK_ERR( error );
1142  std::cout << "#sprocs = " << shprocs.size() << std::endl;
1143  // DBG
1144  }
1145  else if( procs == 1 )
1146  {
1147 #endif
1148  error = mbImpl->load_file( filename, &fileset );CHECK_ERR( error );
1149 #ifdef MOAB_HAVE_MPI
1150  }
1151 #endif
1152 
1153  // Generate hierarchy
1154  error = refine_entities( &moab, pc, fileset, level_degrees, num_levels, false );CHECK_ERR( error );
1155 
1156  return MB_SUCCESS;
1157 }
1158 
1159 int main( int argc, char* argv[] )
1160 {
1161 #ifdef MOAB_HAVE_MPI
1162  MPI_Init( &argc, &argv );
1163 
1164  int nprocs, rank;
1165  MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
1166  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
1167 #endif
1168 
1169  ErrorCode result;
1170  if( argc == 1 )
1171  {
1172  result = test_1D();
1173  handle_error_code( result, number_tests_failed, number_tests_successful );
1174  std::cout << "\n";
1175 
1176  result = test_2D();
1177  handle_error_code( result, number_tests_failed, number_tests_successful );
1178  std::cout << "\n";
1179 
1180  result = test_3D();
1181  handle_error_code( result, number_tests_failed, number_tests_successful );
1182  std::cout << "\n";
1183  }
1184  else if( argc == 2 )
1185  {
1186  const char* filename = argv[1];
1187  int deg[1] = { 2 };
1188  int len = sizeof( deg ) / sizeof( int );
1189  result = test_mesh( filename, deg, len );
1190  handle_error_code( result, number_tests_failed, number_tests_successful );
1191  std::cout << "\n";
1192  }
1193 
1194 #ifdef MOAB_HAVE_MPI
1195  MPI_Finalize();
1196 #endif
1197 
1198  return number_tests_failed;
1199 }