Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
NestedRefine.cpp
Go to the documentation of this file.
1 #include "moab/NestedRefine.hpp"
3 #include "moab/HalfFacetRep.hpp"
4 #include "moab/CpuTimer.hpp"
5 #include "moab/ReadUtilIface.hpp"
6 #include "Internals.hpp"
7 #include "MBTagConventions.hpp"
8 
9 #ifdef MOAB_HAVE_MPI
10 #include "moab/ParallelComm.hpp"
12 #include "moab/Skinner.hpp"
13 #endif
14 
15 #include <iostream>
16 #include <cassert>
17 #include <vector>
18 #include <limits>
19 #include <cmath>
20 
21 namespace moab
22 {
23 
25  : mbImpl( impl ), pcomm( comm ), _rset( rset )
26 {
28  assert( NULL != impl );
29 
30 #ifdef MOAB_HAVE_MPI
31  // Get the Parallel Comm instance to prepare all new sets to work in parallel
32  // in case the user did not provide any arguments
33  if( !comm ) pcomm = moab::ParallelComm::get_pcomm( mbImpl, 0 );
34 #endif
35  error = initialize();
36  if( error != MB_SUCCESS )
37  {
38  std::cout << "Error initializing NestedRefine\n" << std::endl;
39  exit( 1 );
40  }
41 }
42 
44 {
45 #ifdef MOAB_HAVE_AHF
46  ahf = NULL;
47 #else
48  delete ahf;
49 #endif
50  delete tm;
51 }
52 
54 {
56 
57  tm = new CpuTimer();
58  if( !tm ) return MB_MEMORY_ALLOCATION_FAILED;
59 
60 #ifdef MOAB_HAVE_AHF
61  ahf = mbImpl->a_half_facet_rep();
62 #else
63  ahf = new HalfFacetRep( mbImpl, pcomm, _rset, true );
64  if( !ahf ) return MB_MEMORY_ALLOCATION_FAILED;
65 #endif
66 
67  // Check for mixed entity type
68  bool chk_mixed = ahf->check_mixed_entity_type();
69  if( chk_mixed ) MB_SET_ERR( MB_NOT_IMPLEMENTED, "Encountered a mesh with mixed entity types" );
70 
73 
74  // Check for supported entity type
75  if( !_incells.empty() )
76  {
77  EntityType type = mbImpl->type_from_handle( _incells[0] );
78  if( type != MBTET && type != MBHEX )
79  MB_SET_ERR( MB_FAILURE, "Not supported 3D entity types: MBPRISM, MBPYRAMID, MBKNIFE, MBPOLYHEDRON" );
80 
81  meshdim = 3;
82  elementype = type;
83  }
84  else if( !_infaces.empty() )
85  {
86  EntityType type = mbImpl->type_from_handle( _infaces[0] );
87  if( type == MBPOLYGON ) MB_SET_ERR( MB_FAILURE, "Not supported 2D entity type: POLYGON" );
88 
89  meshdim = 2;
90  elementype = type;
91  }
92  else if( !_inedges.empty() )
93  {
94  meshdim = 1;
96  }
97  else
98  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Encountered a mixed-dimensional or invalid mesh" );
99 
100  // Initialize std::map to get indices of degrees.
101  deg_index[2] = 0;
102  deg_index[3] = 1;
103  deg_index[5] = 2;
104 
105  // Set ghost flag to false
106  hasghost = false;
107  return MB_SUCCESS;
108 }
109 
110 /************************************************************
111  * Interface Functions *
112  ************************************************************/
113 
115  int* level_degrees,
116  std::vector< EntityHandle >& level_sets,
117  bool optimize )
118 {
119  assert( num_level > 0 );
120  nlevels = num_level;
121 
123  std::vector< moab::EntityHandle > hmsets( num_level );
124 
125  if( meshdim <= 2 )
126  {
127  for( int i = 0; i < num_level; i++ )
128  {
129  assert( ( level_degrees[i] == 2 ) || ( level_degrees[i] == 3 ) || ( level_degrees[i] == 5 ) );
130  level_dsequence[i] = level_degrees[i];
131  }
132  }
133  else
134  {
135  for( int i = 0; i < num_level; i++ )
136  {
137  assert( ( level_degrees[i] == 2 ) || ( level_degrees[i] == 3 ) );
138  level_dsequence[i] = level_degrees[i];
139  }
140  }
141 
142  error = generate_hm( level_degrees, num_level, &hmsets[0], optimize );MB_CHK_ERR( error );
143 
144  // copy the entity handles
145  level_sets.resize( num_level + 1 );
146  level_sets[0] = _rset;
147  for( int i = 0; i < num_level; i++ )
148  level_sets[i + 1] = hmsets[i];
149 
150  return MB_SUCCESS;
151 }
152 
153 ErrorCode NestedRefine::get_connectivity( EntityHandle ent, int level, std::vector< EntityHandle >& conn )
154 {
156  EntityType type = mbImpl->type_from_handle( ent );
157  EntityHandle start_ent;
158  if( !conn.empty() ) conn.clear();
159  if( level > 0 )
160  {
161  if( type == MBEDGE )
162  {
163  conn.reserve( 2 );
164  start_ent = level_mesh[level - 1].start_edge;
165  EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
166  conn.push_back( level_mesh[level - 1].edge_conn[2 * offset] );
167  conn.push_back( level_mesh[level - 1].edge_conn[2 * offset + 1] );
168  }
169  else if( type == MBTRI || type == MBQUAD )
170  {
171  int num_corners = ahf->lConnMap2D[type - 2].num_verts_in_face;
172  conn.reserve( num_corners );
173  start_ent = level_mesh[level - 1].start_face;
174  EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
175 
176  for( int i = 0; i < num_corners; i++ )
177  conn.push_back( level_mesh[level - 1].face_conn[num_corners * offset + i] );
178  }
179  else if( type == MBTET || type == MBHEX )
180  {
181  int index = ahf->get_index_in_lmap( *_incells.begin() );
182  int num_corners = ahf->lConnMap3D[index].num_verts_in_cell;
183  conn.reserve( num_corners );
184  start_ent = level_mesh[level - 1].start_cell;
185  EntityID offset = ID_FROM_HANDLE( ent ) - ID_FROM_HANDLE( start_ent );
186  for( int i = 0; i < num_corners; i++ )
187  conn.push_back( level_mesh[level - 1].cell_conn[num_corners * offset + i] );
188  }
189  else
190  MB_SET_ERR( MB_FAILURE, "Requesting connectivity for an unsupported entity type" );
191  }
192  else
193  {
194  error = mbImpl->get_connectivity( &ent, 1, conn );MB_CHK_ERR( error );
195  }
196 
197  return MB_SUCCESS;
198 }
199 
200 ErrorCode NestedRefine::get_coordinates( EntityHandle* verts, int num_verts, int level, double* coords )
201 {
202  if( level > 0 )
203  {
204  EntityID vstart = ID_FROM_HANDLE( level_mesh[level - 1].start_vertex );
205  for( int i = 0; i < num_verts; i++ )
206  {
207  const EntityHandle& vid = verts[i];
208  EntityID offset = ID_FROM_HANDLE( vid ) - vstart;
209  coords[3 * i] = level_mesh[level - 1].coordinates[0][offset];
210  coords[3 * i + 1] = level_mesh[level - 1].coordinates[1][offset];
211  coords[3 * i + 2] = level_mesh[level - 1].coordinates[2][offset];
212  }
213  }
214  else
215  {
217  error = mbImpl->get_coords( verts, num_verts, coords );MB_CHK_ERR( error );
218  }
219 
220  return MB_SUCCESS;
221 }
222 
224  const unsigned int target_dimension,
225  std::vector< EntityHandle >& target_entities )
226 
227 {
229  error = ahf->get_adjacencies( source_entity, target_dimension, target_entities );MB_CHK_ERR( error );
230 
231  return MB_SUCCESS;
232 }
233 
234 ErrorCode NestedRefine::child_to_parent( EntityHandle child, int child_level, int parent_level, EntityHandle* parent )
235 {
236  assert( ( child_level > 0 ) && ( child_level > parent_level ) );
237  EntityType type = mbImpl->type_from_handle( child );
238  assert( type != MBVERTEX );
239 
240  int child_index;
241  if( type == MBEDGE )
242  child_index = child - level_mesh[child_level - 1].start_edge;
243  else if( type == MBTRI || type == MBQUAD )
244  child_index = child - level_mesh[child_level - 1].start_face;
245  else if( type == MBTET || type == MBHEX )
246  child_index = child - level_mesh[child_level - 1].start_cell;
247  else
248  MB_SET_ERR( MB_FAILURE, "Requesting parent for unsupported entity type" );
249 
250  int parent_index;
251  int l = child_level - parent_level;
252  for( int i = 0; i < l; i++ )
253  {
254  int d = get_index_from_degree( level_dsequence[child_level - i - 1] );
255  int nch = refTemplates[type - 1][d].total_new_ents;
256  child_index = child_index / nch;
257  }
258  parent_index = child_index;
259 
260  if( type == MBEDGE )
261  {
262  if( parent_level > 0 )
263  *parent = level_mesh[parent_level - 1].start_edge + parent_index;
264  else
265  *parent = _inedges[parent_index];
266  }
267  else if( type == MBTRI || type == MBQUAD )
268  {
269  if( parent_level > 0 )
270  *parent = level_mesh[parent_level - 1].start_face + parent_index;
271  else
272  *parent = _infaces[parent_index];
273  }
274  else if( type == MBTET || type == MBHEX )
275  {
276  if( parent_level > 0 )
277  *parent = level_mesh[parent_level - 1].start_cell + parent_index;
278  else
279  *parent = _incells[parent_index];
280  }
281 
282  return MB_SUCCESS;
283 }
284 
286  int parent_level,
287  int child_level,
288  std::vector< EntityHandle >& children )
289 {
290  assert( ( child_level > 0 ) && ( child_level > parent_level ) );
291  EntityType type = mbImpl->type_from_handle( parent );
292  assert( type != MBVERTEX );
293 
294  int parent_index;
295  if( type == MBEDGE )
296  {
297  if( parent_level > 0 )
298  parent_index = parent - level_mesh[parent_level - 1].start_edge;
299  else
300  parent_index = _inedges.index( parent );
301  }
302  else if( type == MBTRI || type == MBQUAD )
303  {
304  if( parent_level > 0 )
305  parent_index = parent - level_mesh[parent_level - 1].start_face;
306  else
307  parent_index = _infaces.index( parent );
308  }
309  else if( type == MBTET || type == MBHEX )
310  {
311  if( parent_level > 0 )
312  parent_index = parent - level_mesh[parent_level - 1].start_cell;
313  else
314  parent_index = _incells.index( parent );
315  }
316  else
317  MB_SET_ERR( MB_FAILURE, "Requesting children for unsupported entity type" );
318 
319  int start, end;
320  start = end = parent_index;
321  for( int i = parent_level; i < child_level; i++ )
322  {
324  int nch = refTemplates[type - 1][d].total_new_ents;
325  start = start * nch;
326  end = end * nch + nch - 1;
327  }
328 
329  int num_child = end - start;
330  children.reserve( num_child );
331 
332  for( int i = start; i <= end; i++ )
333  {
335  if( type == MBEDGE )
336  child = level_mesh[child_level - 1].start_edge + i;
337  else if( type == MBTRI || type == MBQUAD )
338  child = level_mesh[child_level - 1].start_face + i;
339  else if( type == MBTET || type == MBHEX )
340  child = level_mesh[child_level - 1].start_cell + i;
341 
342  children.push_back( child );
343  }
344 
345  return MB_SUCCESS;
346 }
347 
349  int vert_level,
350  int parent_level,
351  std::vector< EntityHandle >& incident_entities )
352 {
353  assert( vert_level > parent_level );
355 
356  // Step 1: Get the incident entities at the current level
357  std::vector< EntityHandle > inents;
358  if( meshdim == 1 )
359  {
361  }
362  else if( meshdim == 2 )
363  {
365  }
366  else if( meshdim == 3 )
367  {
369  }
370 
371  // Step 2: Loop over all the incident entities at the current level and gather their parents
372  for( int i = 0; i < (int)inents.size(); i++ )
373  {
374  EntityHandle ent = inents[i];
375  EntityHandle parent;
376  error = child_to_parent( ent, vert_level, parent_level, &parent );MB_CHK_ERR( error );
377  incident_entities.push_back( parent );
378  }
379 
380  // Step 3: Sort and remove duplicates
381  std::sort( incident_entities.begin(), incident_entities.end() );
382  incident_entities.erase( std::unique( incident_entities.begin(), incident_entities.end() ),
383  incident_entities.end() );
384 
385  return MB_SUCCESS;
386 }
387 
389  int vert_level,
390  int child_level,
391  std::vector< EntityHandle >& incident_entities )
392 {
393  assert( vert_level < child_level );
395 
396  // Step 1: Get the incident entities at the current level
397  std::vector< EntityHandle > inents;
398  if( meshdim == 1 )
399  {
401  }
402  else if( meshdim == 2 )
403  {
405  }
406  else if( meshdim == 3 )
407  {
409  }
410 
411  // Step 2: Loop over all the incident entities at the current level and gather their parents
412  std::vector< EntityHandle > childs;
413  for( int i = 0; i < (int)inents.size(); i++ )
414  {
415  childs.clear();
416  EntityHandle ent = inents[i];
417  error = parent_to_child( ent, vert_level, child_level, childs );MB_CHK_ERR( error );
418  for( int j = 0; j < (int)childs.size(); j++ )
419  incident_entities.push_back( childs[j] );
420  }
421 
422  return MB_SUCCESS;
423 }
424 
426 {
427  if( ( vertex - *_inverts.begin() ) > _inverts.size() )
428  MB_SET_ERR( MB_FAILURE, "Requesting duplicates for non-coarse vertices" );
429 
430  dupvertex = level_mesh[level - 1].start_vertex + ( vertex - *_inverts.begin() );
431 
432  return MB_SUCCESS;
433 }
434 
436 {
437  bool is_border = false;
438  EntityType type = mbImpl->type_from_handle( entity );
439 
440  if( type == MBVERTEX )
441  is_border = is_vertex_on_boundary( entity );
442  else if( type == MBEDGE )
443  is_border = is_edge_on_boundary( entity );
444  else if( type == MBTRI || type == MBQUAD )
445  is_border = is_face_on_boundary( entity );
446  else if( type == MBTET || type == MBHEX )
447  is_border = is_cell_on_boundary( entity );
448  else
449  MB_SET_ERR( MB_FAILURE, "Requesting boundary information for unsupported entity type" );
450 
451  return is_border;
452 }
453 
454 ErrorCode NestedRefine::exchange_ghosts( std::vector< EntityHandle >& lsets, int num_glayers )
455 {
457 
458  if( hasghost ) return MB_SUCCESS;
459 
460  hasghost = true;
461 #ifdef MOAB_HAVE_MPI
462  error = pcomm->exchange_ghost_cells( meshdim, 0, num_glayers, 0, true, false );MB_CHK_ERR( error );
463  {
464  Range empty_range;
466  // error = pcomm->assign_global_ids(lsets[i], 0, 1, false, true, false);MB_CHK_ERR(error);
467  }
468 #else
469  MB_SET_ERR( MB_FAILURE, "Requesting ghost layers for a serial mesh" );
470 #endif
471 
472  Range* lverts = new Range[lsets.size()];
473  Range* lents = new Range[lsets.size()];
474  for( size_t i = 0; i < lsets.size(); i++ )
475  {
476  error = mbImpl->get_entities_by_dimension( lsets[i], meshdim, lents[i] );MB_CHK_ERR( error );
477  error = mbImpl->get_connectivity( lents[i], lverts[i] );MB_CHK_ERR( error );
478 
479  for( int gl = 0; gl < num_glayers; gl++ )
480  {
481  error = mbImpl->get_adjacencies( lverts[i], meshdim, false, lents[i], Interface::UNION );MB_CHK_ERR( error );
482  error = mbImpl->get_connectivity( lents[i], lverts[i] );MB_CHK_ERR( error );
483  }
484  }
485  for( size_t i = 0; i < lsets.size(); i++ )
486  {
487  error = mbImpl->add_entities( lsets[i], lverts[i] );MB_CHK_ERR( error );
488  error = mbImpl->add_entities( lsets[i], lents[i] );MB_CHK_ERR( error );
489  }
490 
491  delete[] lverts;
492  delete[] lents;
493  return MB_SUCCESS;
494 }
495 
497 {
498  assert( level > 0 && level < nlevels + 1 );
499 
501  std::vector< Tag > mtags( 3 );
502 
506 
507  for( int i = 0; i < 3; i++ )
508  {
509  // Gather sets of a particular tag
510  Range sets;
511  error = mbImpl->get_entities_by_type_and_tag( _rset, MBENTITYSET, &mtags[i], NULL, 1, sets );MB_CHK_ERR( error );
512 
513  // Loop over all sets, gather entities in each set and add their children at all levels to
514  // the set
515  Range set_ents;
516  Range::iterator set_it;
517  std::vector< EntityHandle > childs;
518 
519  for( set_it = sets.begin(); set_it != sets.end(); ++set_it )
520  {
521  // Get the entities in the set, recursively
522  set_ents.clear();
523  childs.clear();
524  error = mbImpl->get_entities_by_handle( *set_it, set_ents, true );MB_CHK_ERR( error );
525 
526  // Gather child entities at the input level
527  for( Range::iterator sit = set_ents.begin(); sit != set_ents.end(); sit++ )
528  {
529  EntityType type = mbImpl->type_from_handle( *sit );
530  if( type == MBVERTEX )
531  {
532  Range conn;
533  std::vector< EntityHandle > cents;
534  error = vertex_to_entities_down( *sit, 0, level, cents );MB_CHK_ERR( error );
535  error = mbImpl->get_connectivity( &cents[0], (int)cents.size(), conn, true );MB_CHK_ERR( error );
536  childs.insert( childs.end(), cents.begin(), cents.end() );
537  }
538  else
539  {
540  error = parent_to_child( *sit, 0, level, childs );MB_CHK_ERR( error );
541  }
542 
543  std::sort( childs.begin(), childs.end() );
544  childs.erase( std::unique( childs.begin(), childs.end() ), childs.end() );
545 
546  // Add child entities to tagged sets
547  error = mbImpl->add_entities( *set_it, &childs[0], childs.size() );MB_CHK_ERR( error );
548  }
549 
550  // Remove the coarse entities
551  error = mbImpl->remove_entities( *set_it, set_ents );MB_CHK_ERR( error );
552 
553  // Add
554  error = mbImpl->add_entities( lset, &( *set_it ), 1 );MB_CHK_ERR( error );
555  }
556  }
557  return MB_SUCCESS;
558 }
559 
560 /***********************************************
561  * Basic functionalities: generate HM *
562  ***********************************************/
563 
564 ErrorCode NestedRefine::estimate_hm_storage( EntityHandle set, int level_degree, int cur_level, int hmest[4] )
565 {
567 
568  // Obtain the size of input mesh.
569  int nverts_prev, nedges_prev, nfaces_prev, ncells_prev;
570  if( cur_level )
571  {
572  nverts_prev = level_mesh[cur_level - 1].num_verts;
573  nedges_prev = level_mesh[cur_level - 1].num_edges;
574  nfaces_prev = level_mesh[cur_level - 1].num_faces;
575  ncells_prev = level_mesh[cur_level - 1].num_cells;
576  }
577  else
578  {
579  nverts_prev = _inverts.size();
580  nedges_prev = _inedges.size();
581  nfaces_prev = _infaces.size();
582  ncells_prev = _incells.size();
583  }
584 
585  // Estimate mesh size of current level mesh.
586  int nedges = 0, nfaces = 0;
587  error = count_subentities( set, cur_level - 1, &nedges, &nfaces );MB_CHK_ERR( error );
588 
589  int d = get_index_from_degree( level_degree );
590  int nverts = refTemplates[MBEDGE - 1][d].nv_edge * nedges;
591  hmest[0] = nverts_prev + nverts;
592  hmest[1] = nedges_prev * refTemplates[MBEDGE - 1][d].total_new_ents;
593  hmest[2] = 0;
594  hmest[3] = 0;
595 
596  int findex, cindex;
597  if( nfaces_prev != 0 )
598  {
599  EntityHandle start_face;
600  if( cur_level )
601  start_face = level_mesh[cur_level - 1].start_face;
602  else
603  start_face = *_infaces.begin();
604  findex = mbImpl->type_from_handle( start_face ) - 1;
605  hmest[2] = nfaces_prev * refTemplates[findex][d].total_new_ents;
606 
607  if( meshdim == 2 ) hmest[0] += refTemplates[findex][d].nv_face * nfaces_prev;
608 
609  if( meshdim == 3 ) hmest[1] += nfaces_prev * intFacEdg[findex - 1][d].nie;
610  }
611 
612  if( ncells_prev != 0 )
613  {
614  cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
615  hmest[3] = ncells_prev * refTemplates[cindex][d].total_new_ents;
616 
617  hmest[0] += refTemplates[cindex][d].nv_face * nfaces;
618  hmest[0] += refTemplates[cindex][d].nv_cell * ncells_prev;
619  }
620 
621  return MB_SUCCESS;
622 }
623 
625 {
626  // Obtain chunks of memory for the current level. Add them to a particular meshset.
627  EntityHandle set_handle;
628  ErrorCode error = mbImpl->create_meshset( MESHSET_SET, set_handle );MB_CHK_SET_ERR( error, "Cannot create mesh for the current level" );
629  *set = set_handle;
630 
631  ReadUtilIface* read_iface;
632  error = mbImpl->query_interface( read_iface );MB_CHK_ERR( error );
633 
634  // Vertices
635  error = read_iface->get_node_coords( 3, estL[0], 0, level_mesh[cur_level].start_vertex,
636  level_mesh[cur_level].coordinates );MB_CHK_ERR( error );
637  level_mesh[cur_level].num_verts = estL[0];
638 
639  Range newverts( level_mesh[cur_level].start_vertex, level_mesh[cur_level].start_vertex + estL[0] - 1 );
640  error = mbImpl->add_entities( *set, newverts );MB_CHK_ERR( error );
641  level_mesh[cur_level].verts = newverts;
642 
643  Tag gidtag;
645  error = read_iface->assign_ids( gidtag, newverts, level_mesh[cur_level].start_vertex );MB_CHK_ERR( error );
646 
647  // Edges
648  if( estL[1] )
649  {
650  error = read_iface->get_element_connect( estL[1], 2, MBEDGE, 0, level_mesh[cur_level].start_edge,
651  level_mesh[cur_level].edge_conn );MB_CHK_ERR( error );
652  level_mesh[cur_level].num_edges = estL[1];
653 
654  Range newedges( level_mesh[cur_level].start_edge, level_mesh[cur_level].start_edge + estL[1] - 1 );
655  error = mbImpl->add_entities( *set, newedges );MB_CHK_ERR( error );
656  level_mesh[cur_level].edges = newedges;
657  }
658  else
659  level_mesh[cur_level].num_edges = 0;
660 
661  // Faces
662  if( estL[2] )
663  {
664  EntityType type = mbImpl->type_from_handle( *( _infaces.begin() ) );
665  int nvpf = ahf->lConnMap2D[type - 2].num_verts_in_face;
666  error = read_iface->get_element_connect( estL[2], nvpf, type, 0, level_mesh[cur_level].start_face,
667  level_mesh[cur_level].face_conn );MB_CHK_ERR( error );
668  level_mesh[cur_level].num_faces = estL[2];
669 
670  Range newfaces( level_mesh[cur_level].start_face, level_mesh[cur_level].start_face + estL[2] - 1 );
671  error = mbImpl->add_entities( *set, newfaces );MB_CHK_ERR( error );
672  level_mesh[cur_level].faces = newfaces;
673  }
674  else
675  level_mesh[cur_level].num_faces = 0;
676 
677  // Cells
678  if( estL[3] )
679  {
680  EntityType type = mbImpl->type_from_handle( *( _incells.begin() ) );
681  int index = ahf->get_index_in_lmap( *_incells.begin() );
682  int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
683  error = read_iface->get_element_connect( estL[3], nvpc, type, 0, level_mesh[cur_level].start_cell,
684  level_mesh[cur_level].cell_conn );MB_CHK_ERR( error );
685  level_mesh[cur_level].num_cells = estL[3];
686 
687  Range newcells( level_mesh[cur_level].start_cell, level_mesh[cur_level].start_cell + estL[3] - 1 );
688  error = mbImpl->add_entities( *set, newcells );MB_CHK_ERR( error );
689  level_mesh[cur_level].cells = newcells;
690  }
691  else
692  level_mesh[cur_level].num_cells = 0;
693 
694  // Resize the ahf maps
695  error = ahf->resize_hf_maps( level_mesh[cur_level].start_vertex, level_mesh[cur_level].num_verts,
696  level_mesh[cur_level].start_edge, level_mesh[cur_level].num_edges,
697  level_mesh[cur_level].start_face, level_mesh[cur_level].num_faces,
698  level_mesh[cur_level].start_cell, level_mesh[cur_level].num_cells );MB_CHK_ERR( error );
699 
701 
702  // If the mesh type changes, then update the member variable in ahf to use the applicable
703  // adjacency matrix
704  MESHTYPE nwmesh = ahf->get_mesh_type( level_mesh[cur_level].num_verts, level_mesh[cur_level].num_edges,
705  level_mesh[cur_level].num_faces, level_mesh[cur_level].num_cells );MB_CHK_ERR( error );
706  if( ahf->thismeshtype != nwmesh ) ahf->thismeshtype = nwmesh;
707 
708  return MB_SUCCESS;
709 }
710 
711 /**********************************
712  * Hierarchical Mesh Generation *
713  * *********************************/
714 
715 ErrorCode NestedRefine::generate_hm( int* level_degrees, int num_level, EntityHandle* hm_set, bool optimize )
716 {
718 
719  Tag gidtag;
721 
722  nlevels = num_level;
723 
724  timeall.tm_total = 0;
725  timeall.tm_refine = 0;
726  timeall.tm_resolve = 0;
727 
728  for( int l = 0; l < num_level; l++ )
729  {
730  double tstart;
731  tstart = tm->time_elapsed();
732 
733  // Estimate storage
734  int hmest[4] = { 0, 0, 0, 0 };
735  EntityHandle set;
736  if( l )
737  set = hm_set[l - 1];
738  else
739  set = _rset;
740  error = estimate_hm_storage( set, level_degrees[l], l, hmest );MB_CHK_ERR( error );
741 
742  // Create arrays for storing the current level
743  error = create_hm_storage_single_level( &hm_set[l], l, hmest );MB_CHK_ERR( error );
744 
745  // Copy the old vertices along with their coordinates
747 
748  // Create the new entities and new vertices
749  error = construct_hm_entities( l, level_degrees[l] );MB_CHK_ERR( error );
750 
751  timeall.tm_refine += tm->time_elapsed() - tstart;
752 
753  // Go into parallel communication
754  if( !optimize )
755  {
756 #ifdef MOAB_HAVE_MPI
757  if( pcomm && ( pcomm->size() > 1 ) )
758  {
759  double tpstart = tm->time_elapsed();
760  error = resolve_shared_ents_parmerge( l, hm_set[l] );MB_CHK_ERR( error );
761  timeall.tm_resolve += tm->time_elapsed() - tpstart;
762  }
763 #endif
764  }
765  }
766 
767  if( optimize )
768  {
769 #ifdef MOAB_HAVE_MPI
770  if( pcomm && ( pcomm->size() > 1 ) )
771  {
772  double tpstart = tm->time_elapsed();
773  error = resolve_shared_ents_opt( hm_set, nlevels );MB_CHK_ERR( error );
774  timeall.tm_resolve = tm->time_elapsed() - tpstart;
775  }
776 #endif
777  }
779 
780  return MB_SUCCESS;
781 }
782 
784 {
786 
787  // Generate mesh for current level by refining previous level.
788  if( ahf->thismeshtype == CURVE )
789  {
790  error = construct_hm_1D( cur_level, deg );MB_CHK_ERR( error );
791  }
792  else if( ahf->thismeshtype == SURFACE || ahf->thismeshtype == SURFACE_MIXED )
793  {
794  error = construct_hm_2D( cur_level, deg );MB_CHK_ERR( error );
795  }
796  else
797  {
798  error = construct_hm_3D( cur_level, deg );MB_CHK_ERR( error );
799  }
800 
801  return MB_SUCCESS;
802 }
803 
804 ErrorCode NestedRefine::construct_hm_1D( int cur_level, int deg )
805 {
807  int nverts_prev, nents_prev;
808  if( cur_level )
809  {
810  nverts_prev = level_mesh[cur_level - 1].num_verts;
811  nents_prev = level_mesh[cur_level - 1].num_edges;
812  }
813  else
814  {
815  nverts_prev = _inverts.size();
816  nents_prev = _inedges.size();
817  }
818 
819  int d = get_index_from_degree( deg );
820  int vtotal = 2 + refTemplates[0][d].total_new_verts;
821  std::vector< EntityHandle > vbuffer( vtotal );
822 
823  std::vector< EntityHandle > conn;
824  int count_nents = 0;
825  int count_verts = nverts_prev;
826 
827  // Step 1: Create the subentities via refinement of the previous mesh
828  for( int eid = 0; eid < nents_prev; eid++ )
829  {
830  conn.clear();
831 
832  // EntityHandle of the working edge
834  if( cur_level )
835  edge = level_mesh[cur_level - 1].start_edge + eid;
836  else
837  edge = _inedges[eid]; // Makes the assumption initial mesh is contiguous in memory
838 
839  error = get_connectivity( edge, cur_level, conn );MB_CHK_ERR( error );
840 
841  // Add the vertex handles to vbuffer for the current level for the working edge
842 
843  // Since the old vertices are copied first, their local indices do not change as new levels
844  // are added. Clearly the local indices of the new vertices introduced in the current level
845  // is still the same when the old vertices are copied. Thus, there is no need to explicitly
846  // store another map between the old and duplicates in the subsequent levels. The second
847  // part in the following sum is the local index in the previous level.
848 
849  // Add the corners to the vbuffer first.
850 
851  for( int i = 0; i < (int)conn.size(); i++ )
852  {
853  if( cur_level )
854  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
855  else
856  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
857  }
858 
859  // Adding rest of the entityhandles to working buffer for vertices.
860  int num_new_verts = refTemplates[0][d].total_new_verts;
861 
862  for( int i = 0; i < num_new_verts; i++ )
863  {
864  vbuffer[i + 2] = level_mesh[cur_level].start_vertex + count_verts;
865  count_verts += 1;
866  }
867 
868  // Use the template to obtain the subentities
869  int id1, id2;
870  int etotal = refTemplates[0][d].total_new_ents;
871  std::vector< EntityHandle > ent_buffer( etotal );
872 
873  for( int i = 0; i < etotal; i++ )
874  {
875  id1 = refTemplates[0][d].ents_conn[i][0];
876  id2 = refTemplates[0][d].ents_conn[i][1];
877  level_mesh[cur_level].edge_conn[2 * ( count_nents )] = vbuffer[id1];
878  level_mesh[cur_level].edge_conn[2 * ( count_nents ) + 1] = vbuffer[id2];
879  ent_buffer[i] = level_mesh[cur_level].start_edge + count_nents;
880  count_nents += 1;
881  };
882 
883  error = update_local_ahf( deg, MBEDGE, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
884 
885  // Compute the coordinates of the new vertices: Linear interpolation
886  int idx;
887  double xi;
888  for( int i = 0; i < num_new_verts; i++ )
889  {
890  xi = refTemplates[0][d].vert_nat_coord[i][0];
891  idx = vbuffer[i + 2] - level_mesh[cur_level].start_vertex; // index of new vertex in current level
892  if( cur_level )
893  {
894  id1 = conn[0] - level_mesh[cur_level - 1].start_vertex; // index of old end vertices in current level
895  id2 = conn[1] - level_mesh[cur_level - 1].start_vertex;
896  }
897  else
898  {
899  id1 = _inverts.index( conn[0] );
900  id2 = _inverts.index( conn[1] );
901  }
902 
903  level_mesh[cur_level].coordinates[0][idx] =
904  ( 1 - xi ) * level_mesh[cur_level].coordinates[0][id1] + xi * level_mesh[cur_level].coordinates[0][id2];
905  level_mesh[cur_level].coordinates[1][idx] =
906  ( 1 - xi ) * level_mesh[cur_level].coordinates[1][id1] + xi * level_mesh[cur_level].coordinates[1][id2];
907  level_mesh[cur_level].coordinates[2][idx] =
908  ( 1 - xi ) * level_mesh[cur_level].coordinates[2][id1] + xi * level_mesh[cur_level].coordinates[2][id2];
909  }
910  }
911 
912  error = update_global_ahf( MBEDGE, cur_level, deg );MB_CHK_ERR( error );
913 
914  return MB_SUCCESS;
915 }
916 
918  int deg,
919  EntityType type,
920  std::vector< EntityHandle >& trackverts )
921 {
923 
924  int nedges_prev;
925  if( cur_level )
926  nedges_prev = level_mesh[cur_level - 1].num_edges;
927  else
928  nedges_prev = _inedges.size();
929 
930  int d = get_index_from_degree( deg );
931  int nve = refTemplates[0][d].nv_edge;
932  int vtotal = 2 + refTemplates[0][d].total_new_verts;
933  int etotal = refTemplates[0][d].total_new_ents;
934  int ne = 0, dim = 0, index = 0;
935  if( type == MBTRI || type == MBQUAD )
936  {
937  index = type - 2;
938  ne = ahf->lConnMap2D[index].num_verts_in_face;
939  dim = 2;
940  }
941  else if( type == MBTET || type == MBHEX )
942  {
943  index = ahf->get_index_in_lmap( *( _incells.begin() ) );
944  ne = ahf->lConnMap3D[index].num_edges_in_cell;
945  dim = 3;
946  }
947 
948  std::vector< EntityHandle > vbuffer( vtotal );
949  std::vector< EntityHandle > ent_buffer( etotal );
950 
951  std::vector< EntityHandle > adjents, econn, fconn;
952  std::vector< int > leids;
953  int count_nents = 0;
954 
955  // Loop over all the edges and gather the vertices to be used for refinement
956  for( int eid = 0; eid < nedges_prev; eid++ )
957  {
958  adjents.clear();
959  leids.clear();
960  econn.clear();
961  fconn.clear();
962  for( int i = 0; i < vtotal; i++ )
963  vbuffer[i] = 0;
964  for( int i = 0; i < etotal; i++ )
965  ent_buffer[i] = 0;
966 
968  if( cur_level )
969  edge = level_mesh[cur_level - 1].start_edge + eid;
970  else
971  edge = _inedges[eid];
972 
973  error = get_connectivity( edge, cur_level, econn );MB_CHK_ERR( error );
974 
975  for( int i = 0; i < (int)econn.size(); i++ )
976  {
977  if( cur_level )
978  vbuffer[i] = level_mesh[cur_level].start_vertex + ( econn[i] - level_mesh[cur_level - 1].start_vertex );
979  else
980  vbuffer[i] = level_mesh[cur_level].start_vertex + ( econn[i] - *_inverts.begin() );
981  }
982 
983  int fid = -1, lid = -1, idx1 = -1, idx2 = -1;
984 
985  if( dim == 2 )
986  {
987  error = ahf->get_up_adjacencies_2d( edge, adjents, &leids );MB_CHK_ERR( error );
988  if( cur_level )
989  fid = adjents[0] - level_mesh[cur_level - 1].start_face;
990  else
991  fid = _infaces.index( adjents[0] );
992 
993  lid = leids[0];
994  idx1 = lid;
995  idx2 = ahf->lConnMap2D[index].next[lid];
996  }
997  else if( dim == 3 )
998  {
999  error = ahf->get_up_adjacencies_edg_3d( edge, adjents, &leids );MB_CHK_ERR( error );
1000  if( cur_level )
1001  fid = adjents[0] - level_mesh[cur_level - 1].start_cell;
1002  else
1003  fid = _incells.index( adjents[0] );
1004 
1005  lid = leids[0];
1006  idx1 = ahf->lConnMap3D[index].e2v[lid][0];
1007  idx2 = ahf->lConnMap3D[index].e2v[lid][1];
1008  }
1009 
1010  error = get_connectivity( adjents[0], cur_level, fconn );MB_CHK_ERR( error );
1011 
1012  bool orient = false;
1013  if( ( fconn[idx1] == econn[0] ) && ( fconn[idx2] == econn[1] ) ) orient = true;
1014 
1015  if( orient )
1016  {
1017  for( int j = 0; j < nve; j++ )
1018  vbuffer[j + 2] = trackverts[fid * ne * nve + nve * lid + j];
1019  }
1020  else
1021  {
1022  for( int j = 0; j < nve; j++ )
1023  vbuffer[( nve - j - 1 ) + 2] = trackverts[fid * ne * nve + nve * lid + j];
1024  }
1025 
1026  // Use the template to obtain the subentities
1027  int id1, id2;
1028 
1029  for( int i = 0; i < etotal; i++ )
1030  {
1031  id1 = refTemplates[0][d].ents_conn[i][0];
1032  id2 = refTemplates[0][d].ents_conn[i][1];
1033  level_mesh[cur_level].edge_conn[2 * ( count_nents )] = vbuffer[id1];
1034  level_mesh[cur_level].edge_conn[2 * ( count_nents ) + 1] = vbuffer[id2];
1035  ent_buffer[i] = level_mesh[cur_level].start_edge + count_nents;
1036 
1037  count_nents += 1;
1038  };
1039 
1040  error = update_local_ahf( deg, MBEDGE, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1041  }
1042 
1043  error = update_global_ahf_1D_sub( cur_level, deg );MB_CHK_ERR( error );
1044 
1045  return MB_SUCCESS;
1046 }
1047 
1048 ErrorCode NestedRefine::construct_hm_2D( int cur_level, int deg )
1049 {
1050  ErrorCode error;
1051  int nverts_prev, nents_prev;
1052  if( cur_level )
1053  {
1054  nverts_prev = level_mesh[cur_level - 1].num_verts;
1055  nents_prev = level_mesh[cur_level - 1].num_faces;
1056  }
1057  else
1058  {
1059  nverts_prev = _inverts.size();
1060  nents_prev = _infaces.size();
1061  }
1062 
1063  // Create some book-keeping arrays over the old mesh to avoid introducing duplicate vertices and
1064  // calculating vertices more than once.
1065  EntityType ftype = mbImpl->type_from_handle( *_infaces.begin() );
1066  int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
1067  int findex = ftype - 1;
1068 
1069  int d = get_index_from_degree( deg );
1070  int tnv = refTemplates[findex][d].total_new_verts;
1071  int vtotal = nepf + tnv;
1072  std::vector< EntityHandle > vbuffer( vtotal );
1073  int etotal = refTemplates[findex][d].total_new_ents;
1074  std::vector< EntityHandle > ent_buffer( etotal );
1075 
1076  int nve = refTemplates[findex][d].nv_edge;
1077  std::vector< EntityHandle > trackvertsF( nents_prev * nepf * nve, 0 );
1078  int cur_nverts = level_mesh[cur_level].num_verts;
1079  std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
1080 
1081  int count_nverts = nverts_prev;
1082  int count_nents = 0;
1083  std::vector< EntityHandle > conn, cur_conn;
1084 
1085  // Step 1: Create the subentities via refinement of the previous mesh
1086  for( int fid = 0; fid < nents_prev; fid++ )
1087  {
1088  conn.clear();
1089  cur_conn.clear();
1090  for( int i = 0; i < vtotal; i++ )
1091  vbuffer[i] = 0;
1092  for( int i = 0; i < etotal; i++ )
1093  ent_buffer[i] = 0;
1094 
1095  // EntityHandle of the working face
1097  if( cur_level )
1098  face = level_mesh[cur_level - 1].start_face + fid;
1099  else
1100  face = _infaces[fid];
1101 
1102  error = get_connectivity( face, cur_level, conn );MB_CHK_ERR( error );
1103 
1104  // Step 1: Add vertices from the current level for the working face that will be used for
1105  // subdivision.
1106  // Add the corners to vbuffer
1107  for( int i = 0; i < (int)conn.size(); i++ )
1108  {
1109  if( cur_level )
1110  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
1111  else
1112  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
1113 
1114  cur_conn.push_back( vbuffer[i] );
1115  }
1116 
1117  // Gather vertices already added to tracking array due to refinement of the sibling faces
1118 
1119  for( int i = 0; i < nepf; i++ )
1120  {
1121  for( int j = 0; j < nve; j++ )
1122  {
1123  int id = refTemplates[findex][d].vert_on_edges[i][j];
1124  vbuffer[id] = trackvertsF[fid * nve * nepf + nve * i + j];
1125  }
1126  }
1127 
1128  // Add the remaining vertex handles to vbuffer for the current level for the working face
1129  for( int i = 0; i < tnv; i++ )
1130  {
1131  if( !vbuffer[i + nepf] )
1132  {
1133  vbuffer[i + nepf] = level_mesh[cur_level].start_vertex + count_nverts;
1134  count_nverts += 1;
1135  }
1136  }
1137 
1138  // Step 2: Create the subentities using the template and the vbuffer
1139  int idx;
1140  for( int i = 0; i < etotal; i++ )
1141  {
1142  for( int k = 0; k < nepf; k++ )
1143  {
1144  idx = refTemplates[findex][d].ents_conn[i][k];
1145  level_mesh[cur_level].face_conn[nepf * count_nents + k] = vbuffer[idx];
1146  }
1147  ent_buffer[i] = level_mesh[cur_level].start_face + count_nents;
1148  count_nents += 1;
1149  }
1150 
1151  // Step 3: Update the local AHF maps
1152  error = update_local_ahf( deg, ftype, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1153 
1154  // Step 4: Add the new vertices to the tracking array
1155  int id;
1156 
1157  for( int i = 0; i < nepf; i++ )
1158  {
1159  // Add the vertices to trackvertsF for fid
1160  for( int j = 0; j < nve; j++ )
1161  {
1162  id = refTemplates[findex][d].vert_on_edges[i][j];
1163  trackvertsF[fid * nepf * nve + nve * i + j] = vbuffer[id];
1164  }
1165 
1166  std::vector< EntityHandle > sibfids;
1167  std::vector< int > sibleids;
1168  std::vector< int > siborient;
1169 
1170  // Add the vertices to trackvertsF for siblings of fid, if any.
1171  error = ahf->get_up_adjacencies_2d( face, i, false, sibfids, &sibleids, &siborient );MB_CHK_ERR( error );
1172 
1173  if( !sibfids.size() ) continue;
1174 
1175  for( int s = 0; s < (int)sibfids.size(); s++ )
1176  {
1177  int sibid;
1178  if( cur_level )
1179  sibid = sibfids[s] - level_mesh[cur_level - 1].start_face;
1180  else
1181  sibid = sibfids[s] - *_infaces.begin();
1182 
1183  if( siborient[s] ) // Same half-edge direction as the current half-edge
1184  {
1185  for( int j = 0; j < nve; j++ )
1186  {
1187  id = refTemplates[findex][d].vert_on_edges[i][j];
1188  trackvertsF[sibid * nepf * nve + nve * sibleids[s] + j] = vbuffer[id];
1189  }
1190  }
1191  else
1192  {
1193  for( int j = 0; j < nve; j++ )
1194  {
1195  id = refTemplates[findex][d].vert_on_edges[i][nve - j - 1];
1196  trackvertsF[sibid * nepf * nve + nve * sibleids[s] + j] = vbuffer[id];
1197  }
1198  }
1199  }
1200  }
1201 
1202  // Step 5: Compute the coordinates of the new vertices, avoids computing more than once via
1203  // the flag_verts array.
1204  std::vector< double > corner_coords( nepf * 3 );
1205  error = get_coordinates( &cur_conn[0], nepf, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
1206 
1207  error = compute_coordinates( cur_level, deg, ftype, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
1208  nverts_prev );MB_CHK_ERR( error );
1209  }
1210 
1211  // Step 6: Update the global maps
1212  error = update_global_ahf( ftype, cur_level, deg );MB_CHK_ERR( error );
1213 
1214  // Step 7: If edges exists, refine them.
1215  if( !_inedges.empty() )
1216  {
1217  error = construct_hm_1D( cur_level, deg, ftype, trackvertsF );MB_CHK_ERR( error );
1218  }
1219 
1220  return MB_SUCCESS;
1221 }
1222 
1224  int deg,
1225  EntityType type,
1226  std::vector< EntityHandle >& trackvertsE,
1227  std::vector< EntityHandle >& trackvertsF )
1228 {
1229  ErrorCode error;
1230 
1231  EntityType ftype = MBTRI;
1232  if( type == MBHEX ) ftype = MBQUAD;
1233 
1234  int d = get_index_from_degree( deg );
1235  int findex = ftype - 1;
1236  int cidx = ahf->get_index_in_lmap( *( _incells.begin() ) );
1237 
1238  int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
1239  int nepc = ahf->lConnMap3D[cidx].num_edges_in_cell;
1240  int nfpc = ahf->lConnMap3D[cidx].num_faces_in_cell;
1241 
1242  int tnv = refTemplates[findex][d].total_new_verts;
1243  int nve = refTemplates[findex][d].nv_edge;
1244  int nvf = refTemplates[findex][d].nv_face;
1245  int vtotal = nepf + tnv;
1246  int etotal = refTemplates[findex][d].total_new_ents;
1247 
1248  std::vector< EntityHandle > vbuffer( vtotal );
1249  std::vector< EntityHandle > ent_buffer( etotal );
1250 
1251  std::vector< EntityHandle > adjents, fconn, cconn;
1252  std::vector< int > leids;
1253  int count_nents = 0;
1254 
1255  int nents_prev, ecount;
1256  if( cur_level )
1257  {
1258  nents_prev = level_mesh[cur_level - 1].num_faces;
1259  ecount = level_mesh[cur_level - 1].num_edges * refTemplates[MBEDGE - 1][d].total_new_ents;
1260  ;
1261  }
1262  else
1263  {
1264  nents_prev = _infaces.size();
1265  ecount = _inedges.size() * refTemplates[MBEDGE - 1][d].total_new_ents;
1266  ;
1267  }
1268 
1269  // Step 1: Create the subentities via refinement of the previous mesh
1270  for( int it = 0; it < nents_prev; it++ )
1271  {
1272  fconn.clear();
1273  cconn.clear();
1274  adjents.clear();
1275  leids.clear();
1276  for( int i = 0; i < vtotal; i++ )
1277  vbuffer[i] = 0;
1278  for( int i = 0; i < etotal; i++ )
1279  ent_buffer[i] = 0;
1280 
1281  // EntityHandle of the working face
1283  if( cur_level )
1284  face = level_mesh[cur_level - 1].start_face + it;
1285  else
1286  face = _infaces[it];
1287 
1288  error = get_connectivity( face, cur_level, fconn );MB_CHK_ERR( error );
1289 
1290  // Add the new handles for old connectivity in the buffer
1291  for( int i = 0; i < (int)fconn.size(); i++ )
1292  {
1293  if( cur_level )
1294  vbuffer[i] = level_mesh[cur_level].start_vertex + ( fconn[i] - level_mesh[cur_level - 1].start_vertex );
1295  else
1296  vbuffer[i] = level_mesh[cur_level].start_vertex + ( fconn[i] - *_inverts.begin() );
1297  }
1298 
1299  // Add handles for vertices on edges and faces from the already refined cell
1300  int fid, lid;
1301  error = ahf->get_up_adjacencies_face_3d( face, adjents, &leids );MB_CHK_ERR( error );
1302 
1303  if( cur_level )
1304  fid = adjents[0] - level_mesh[cur_level - 1].start_cell;
1305  else
1306  fid = _incells.index( adjents[0] );
1307 
1308  lid = leids[0];
1309 
1310  error = get_connectivity( adjents[0], cur_level, cconn );MB_CHK_ERR( error );
1311 
1312  // Find the orientation w.r.t the half-face and then add vertices properly.
1313  std::vector< EntityHandle > fac_conn( nepf );
1314  std::vector< EntityHandle > lfac_conn( nepf );
1315  for( int j = 0; j < nepf; j++ )
1316  {
1317  fac_conn[j] = fconn[j];
1318  int id = ahf->lConnMap3D[cidx].hf2v[lid][j];
1319  lfac_conn[j] = cconn[id];
1320  }
1321 
1322  std::vector< int > le_idx, indices;
1323 
1324  error = reorder_indices( deg, &fac_conn[0], &lfac_conn[0], nepf, le_idx, indices );MB_CHK_ERR( error );
1325 
1326  // Add the existing vertices on edges of the already refined cell to the vbuffer
1327  for( int j = 0; j < nepf; j++ )
1328  {
1329  int id = le_idx[j]; // Corresponding local edge
1330  int idx = ahf->lConnMap3D[cidx].f2leid[lid][id]; // Local edge in the cell
1331 
1332  // Get the orientation of the local edge of the face wrt the corresponding local edge in
1333  // the cell
1334  bool eorient = false;
1335  int fnext = ahf->lConnMap2D[ftype - 2].next[j];
1336  int idx1 = ahf->lConnMap3D[cidx].e2v[idx][0];
1337  int idx2 = ahf->lConnMap3D[cidx].e2v[idx][1];
1338  if( ( fconn[j] == cconn[idx1] ) && ( fconn[fnext] == cconn[idx2] ) ) eorient = true;
1339 
1340  if( eorient )
1341  {
1342  for( int k = 0; k < nve; k++ )
1343  {
1344  int ind = refTemplates[findex][d].vert_on_edges[j][k];
1345  vbuffer[ind] = trackvertsE[fid * nepc * nve + nve * idx + k];
1346  }
1347  }
1348  else
1349  {
1350  for( int k = 0; k < nve; k++ )
1351  {
1352  int ind = refTemplates[findex][d].vert_on_edges[j][nve - k - 1];
1353  vbuffer[ind] = trackvertsE[fid * nepc * nve + nve * idx + k];
1354  }
1355  }
1356  }
1357 
1358  // Add the existing vertices on the face of the refine cell to vbuffer
1359  if( nvf )
1360  {
1361  for( int k = 0; k < nvf; k++ )
1362  {
1363  int ind = refTemplates[findex][d].vert_on_faces[0][k];
1364  vbuffer[ind] = trackvertsF[fid * nfpc * nvf + nvf * lid + indices[k] - 1];
1365  }
1366  }
1367 
1368  // Create the subentities using the template and the vbuffer
1369  for( int i = 0; i < etotal; i++ )
1370  {
1371  for( int k = 0; k < nepf; k++ )
1372  {
1373  int idx = refTemplates[findex][d].ents_conn[i][k];
1374  level_mesh[cur_level].face_conn[nepf * count_nents + k] = vbuffer[idx];
1375  }
1376  ent_buffer[i] = level_mesh[cur_level].start_face + count_nents;
1377  count_nents += 1;
1378  }
1379 
1380  error = update_local_ahf( deg, ftype, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1381 
1382  // Create the interior edges
1383  int id1, id2;
1384 
1385  int ne = intFacEdg[ftype - 2][d].nie;
1386  for( int i = 0; i < ne; i++ )
1387  {
1388  id1 = intFacEdg[ftype - 2][d].ieconn[i][0];
1389  id2 = intFacEdg[ftype - 2][d].ieconn[i][1];
1390  level_mesh[cur_level].edge_conn[2 * ( ecount )] = vbuffer[id1];
1391  level_mesh[cur_level].edge_conn[2 * ( ecount ) + 1] = vbuffer[id2];
1392  ecount += 1;
1393  }
1394  }
1395 
1396  // Step 6: Update the global maps
1397  error = update_global_ahf_2D_sub( cur_level, deg );MB_CHK_ERR( error );
1398 
1399  // Step 7: Update the hf-maps for the edges
1400  error = update_ahf_1D( cur_level );MB_CHK_ERR( error );
1401 
1402  return MB_SUCCESS;
1403 }
1404 
1405 ErrorCode NestedRefine::construct_hm_3D( int cur_level, int deg )
1406 {
1407  ErrorCode error;
1408  EntityType type = mbImpl->type_from_handle( *( _incells.begin() ) );
1409  if( type == MBTET )
1410  {
1411  error = subdivide_tets( cur_level, deg );MB_CHK_ERR( error );
1412  }
1413  else
1414  {
1415  error = subdivide_cells( type, cur_level, deg );MB_CHK_ERR( error );
1416  }
1417 
1418  return MB_SUCCESS;
1419 }
1420 
1421 ErrorCode NestedRefine::subdivide_cells( EntityType type, int cur_level, int deg )
1422 {
1423  ErrorCode error;
1424  int nverts_prev, nents_prev;
1425  if( cur_level )
1426  {
1427  nverts_prev = level_mesh[cur_level - 1].num_verts;
1428  nents_prev = level_mesh[cur_level - 1].num_cells;
1429  }
1430  else
1431  {
1432  nverts_prev = _inverts.size();
1433  nents_prev = _incells.size();
1434  }
1435 
1436  // Create some book-keeping arrays over the parent mesh to avoid introducing duplicate vertices
1437  int cindex = type - 1;
1438  int d = get_index_from_degree( deg );
1439  int ne = refTemplates[cindex][d].nv_edge;
1440  int nvf = refTemplates[cindex][d].nv_face;
1441  int nvtotal = refTemplates[cindex][d].total_new_verts;
1442 
1443  int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
1444  int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
1445  int nepc = ahf->lConnMap3D[index].num_edges_in_cell;
1446  int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
1447 
1448  int vtotal = nvpc + nvtotal;
1449  std::vector< EntityHandle > vbuffer( vtotal );
1450 
1451  std::vector< EntityHandle > trackvertsC_edg( nepc * ne * nents_prev, 0 );
1452  std::vector< EntityHandle > trackvertsC_face( nfpc * nvf * nents_prev, 0 );
1453 
1454  int cur_nverts = level_mesh[cur_level].num_verts;
1455  std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
1456 
1457  int count_nverts = nverts_prev;
1458  int count_ents = 0;
1459  std::vector< EntityHandle > conn, cur_conn;
1460 
1461  // Step 1: Create the subentities via refinement of the previous mesh
1462  for( int cid = 0; cid < nents_prev; cid++ )
1463  {
1464  conn.clear();
1465  cur_conn.clear();
1466  for( int i = 0; i < vtotal; i++ )
1467  vbuffer[i] = 0;
1468 
1469  // EntityHandle of the working cell
1470  EntityHandle cell;
1471  if( cur_level )
1472  cell = level_mesh[cur_level - 1].start_cell + cid;
1473  else
1474  cell = _incells[cid];
1475 
1476  error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
1477 
1478  // Step 1: Add vertices from the current level for the working face that will be used for
1479  // subdivision.
1480  // Add the corners to vbuffer
1481  for( int i = 0; i < (int)conn.size(); i++ )
1482  {
1483  if( cur_level )
1484  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
1485  else
1486  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
1487 
1488  cur_conn.push_back( vbuffer[i] );
1489  }
1490 
1491  // Gather vertices already added to tracking array due to refinement of the sibling cells
1492  for( int i = 0; i < nepc; i++ )
1493  {
1494  for( int j = 0; j < ne; j++ )
1495  {
1496  int idx = refTemplates[cindex][d].vert_on_edges[i][j];
1497  vbuffer[idx] = trackvertsC_edg[cid * nepc * ne + ne * i + j];
1498  }
1499  }
1500 
1501  // Add remaining new vertex handles
1502  for( int i = 0; i < nfpc; i++ )
1503  {
1504  for( int j = 0; j < nvf; j++ )
1505  {
1506  int idx = refTemplates[cindex][d].vert_on_faces[i][j];
1507  vbuffer[idx] = trackvertsC_face[cid * nfpc * nvf + nvf * i + j];
1508  }
1509  }
1510 
1511  // Add the remaining vertex handles to vbuffer for the current level for the working cell
1512  for( int i = 0; i < nvtotal; i++ )
1513  {
1514  if( !vbuffer[i + nvpc] )
1515  {
1516  vbuffer[i + nvpc] = level_mesh[cur_level].start_vertex + count_nverts;
1517  count_nverts += 1;
1518  }
1519  }
1520 
1521  // Step 2: Use the template to obtain the subentities. The coordinates and local ahf maps
1522  // are also constructed. Connectivity of the children
1523  int etotal = refTemplates[type - 1][d].total_new_ents;
1524  std::vector< EntityHandle > ent_buffer( etotal );
1525 
1526  for( int i = 0; i < etotal; i++ )
1527  {
1528  for( int k = 0; k < nvpc; k++ )
1529  {
1530  int idx = refTemplates[type - 1][d].ents_conn[i][k];
1531  level_mesh[cur_level].cell_conn[nvpc * count_ents + k] = vbuffer[idx];
1532  }
1533  ent_buffer[i] = level_mesh[cur_level].start_cell + count_ents;
1534  count_ents += 1;
1535  }
1536 
1537  // Step 3: Update local ahf maps
1538  error = update_local_ahf( deg, type, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1539 
1540  // Step 4: Update tracking information
1541  error = update_tracking_verts( cell, cur_level, deg, trackvertsC_edg, trackvertsC_face, &vbuffer[0] );MB_CHK_ERR( error );
1542 
1543  // Step 5: Coordinates of the new vertices
1544  std::vector< double > corner_coords( nvpc * 3 );
1545  error = get_coordinates( &cur_conn[0], nvpc, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
1546 
1547  error = compute_coordinates( cur_level, deg, type, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
1548  nverts_prev );MB_CHK_ERR( error );
1549  }
1550 
1551  // error = ahf->print_tags(3);
1552 
1553  // Step 6: Update the global maps
1554  error = update_global_ahf( type, cur_level, deg );MB_CHK_ERR( error );
1555 
1556  // Step 7: If edges exists, refine them as well.
1557  if( level_mesh[cur_level].num_edges != 0 )
1558  {
1559  error = construct_hm_1D( cur_level, deg, type, trackvertsC_edg );MB_CHK_ERR( error );
1560  }
1561 
1562  // Step 8: If faces exists, refine them as well.
1563  if( !_infaces.empty() )
1564  {
1565  error = construct_hm_2D( cur_level, deg, type, trackvertsC_edg, trackvertsC_face );MB_CHK_ERR( error );
1566  }
1567 
1568  // error = ahf->print_tags(3);
1569 
1570  return MB_SUCCESS;
1571 }
1572 
1573 ErrorCode NestedRefine::subdivide_tets( int cur_level, int deg )
1574 {
1575  ErrorCode error;
1576  int nverts_prev, nents_prev;
1577  if( cur_level )
1578  {
1579  nverts_prev = level_mesh[cur_level - 1].num_verts;
1580  nents_prev = level_mesh[cur_level - 1].num_cells;
1581  }
1582  else
1583  {
1584  nverts_prev = _inverts.size();
1585  nents_prev = _incells.size();
1586  }
1587 
1588  EntityType type = MBTET;
1589  int cindex = type - 1;
1590  int d = get_index_from_degree( deg );
1591  int ne = refTemplates[cindex][d].nv_edge;
1592  int nvf = refTemplates[cindex][d].nv_face;
1593  int nvtotal = refTemplates[cindex][d].total_new_verts;
1594 
1595  int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
1596  int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
1597  int nepc = ahf->lConnMap3D[index].num_edges_in_cell;
1598  int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
1599 
1600  // Create vertex buffer
1601  int vtotal = nvpc + nvtotal;
1602  std::vector< EntityHandle > vbuffer( vtotal );
1603 
1604  // Create book-keeping arrays over the parent mesh to avoid introducing duplicate vertices
1605  std::vector< EntityHandle > trackvertsC_edg( nepc * ne * nents_prev, 0 );
1606  std::vector< EntityHandle > trackvertsC_face( nfpc * nvf * nents_prev, 0 );
1607 
1608  int cur_nverts = level_mesh[cur_level].num_verts;
1609  std::vector< int > flag_verts( cur_nverts - nverts_prev, 0 );
1610  std::vector< int > cell_patterns( nents_prev, 0 );
1611 
1612  int count_nverts = nverts_prev;
1613  int count_ents = 0;
1614  std::vector< EntityHandle > conn, cur_conn;
1615 
1616  // Step 1: Create the subentities via refinement of the previous mesh
1617  for( int cid = 0; cid < nents_prev; cid++ )
1618  {
1619  conn.clear();
1620  cur_conn.clear();
1621  for( int i = 0; i < vtotal; i++ )
1622  vbuffer[i] = 0;
1623 
1624  // EntityHandle of the working cell
1625  EntityHandle cell;
1626  if( cur_level )
1627  cell = level_mesh[cur_level - 1].start_cell + cid;
1628  else
1629  cell = _incells[cid];
1630 
1631  error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
1632 
1633  // Step 1: Add vertices from the current level for the working face that will be used for
1634  // subdivision.
1635  // Add the corners to vbuffer
1636  for( int i = 0; i < (int)conn.size(); i++ )
1637  {
1638  if( cur_level )
1639  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - level_mesh[cur_level - 1].start_vertex );
1640  else
1641  vbuffer[i] = level_mesh[cur_level].start_vertex + ( conn[i] - *_inverts.begin() );
1642 
1643  cur_conn.push_back( vbuffer[i] );
1644  }
1645 
1646  // Gather vertices already added to tracking array due to refinement of the sibling cells
1647  for( int i = 0; i < nepc; i++ )
1648  {
1649  for( int j = 0; j < ne; j++ )
1650  {
1651  int idx = refTemplates[cindex][d].vert_on_edges[i][j];
1652  vbuffer[idx] = trackvertsC_edg[cid * nepc * ne + ne * i + j];
1653  }
1654  }
1655 
1656  // Add remaining new vertex handles
1657  for( int i = 0; i < nfpc; i++ )
1658  {
1659  for( int j = 0; j < nvf; j++ )
1660  {
1661  int idx = refTemplates[cindex][d].vert_on_faces[i][j];
1662  vbuffer[idx] = trackvertsC_face[cid * nfpc * nvf + nvf * i + j];
1663  }
1664  }
1665 
1666  // Add the remaining vertex handles to vbuffer for the current level for the working cell
1667  for( int i = 0; i < nvtotal; i++ )
1668  {
1669  if( !vbuffer[i + nvpc] )
1670  {
1671  vbuffer[i + nvpc] = level_mesh[cur_level].start_vertex + count_nverts;
1672  count_nverts += 1;
1673  }
1674  }
1675 
1676  // Step 2: Coordinates of the new vertices
1677  std::vector< double > corner_coords( nvpc * 3 );
1678  error = get_coordinates( &cur_conn[0], nvpc, cur_level + 1, &corner_coords[0] );MB_CHK_ERR( error );
1679 
1680  error = compute_coordinates( cur_level, deg, type, &vbuffer[0], vtotal, &corner_coords[0], flag_verts,
1681  nverts_prev );MB_CHK_ERR( error );
1682 
1683  // Step 3: Choose the tet refine pattern to be used for this tet
1684  int diag = find_shortest_diagonal_octahedron( cur_level, deg, &vbuffer[0] );
1685  int pat_id = diag + 2;
1686  cell_patterns[cid] = pat_id;
1687 
1688  // Step 4: Use the template to obtain the subentities. The coordinates and local ahf maps
1689  // are also constructed. Connectivity of the children
1690  int etotal = refTemplates[pat_id][d].total_new_ents;
1691  std::vector< EntityHandle > ent_buffer( etotal );
1692 
1693  for( int i = 0; i < etotal; i++ )
1694  {
1695  for( int k = 0; k < nvpc; k++ )
1696  {
1697  int idx = refTemplates[pat_id][d].ents_conn[i][k];
1698  level_mesh[cur_level].cell_conn[nvpc * count_ents + k] = vbuffer[idx];
1699  }
1700  ent_buffer[i] = level_mesh[cur_level].start_cell + count_ents;
1701  count_ents += 1;
1702  }
1703 
1704  // Step 5: Update local ahf maps
1705  error = update_local_ahf( deg, MBTET, pat_id, &vbuffer[0], &ent_buffer[0], etotal );MB_CHK_ERR( error );
1706 
1707  // Step 6: Update tracking information
1708  error = update_tracking_verts( cell, cur_level, deg, trackvertsC_edg, trackvertsC_face, &vbuffer[0] );MB_CHK_ERR( error );
1709  }
1710 
1711  // Step 7: Update the global maps
1712  // error = update_global_ahf(cur_level, deg, cell_patterns); MB_CHK_ERR(error);
1713  error = update_global_ahf( type, cur_level, deg, &cell_patterns );MB_CHK_ERR( error );
1714 
1715  // Step 8: If edges exists, refine them as well.
1716  if( level_mesh[cur_level].num_edges != 0 )
1717  {
1718  error = construct_hm_1D( cur_level, deg, type, trackvertsC_edg );MB_CHK_ERR( error );
1719  }
1720 
1721  // Step 9: If faces exists, refine them as well.
1722  if( !_infaces.empty() )
1723  {
1724  error = construct_hm_2D( cur_level, deg, type, trackvertsC_edg, trackvertsC_face );MB_CHK_ERR( error );
1725  }
1726 
1727  return MB_SUCCESS;
1728 }
1729 
1731  int deg,
1732  EntityType type,
1733  EntityHandle* vbuffer,
1734  int vtotal,
1735  double* corner_coords,
1736  std::vector< int >& vflag,
1737  int nverts_prev )
1738 {
1739  EntityHandle vstart = level_mesh[cur_level].start_vertex;
1740  int d = get_index_from_degree( deg );
1741 
1742  if( type == MBTRI )
1743  {
1744  double xi, eta, N[3];
1745  int findex = mbImpl->type_from_handle( *( _infaces.begin() ) ) - 1;
1746 
1747  for( int i = 3; i < vtotal; i++ )
1748  {
1749  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1750 
1751  xi = refTemplates[findex][d].vert_nat_coord[i - 3][0];
1752  eta = refTemplates[findex][d].vert_nat_coord[i - 3][1];
1753  N[0] = 1 - xi - eta;
1754  N[1] = xi;
1755  N[2] = eta;
1756 
1757  double x = 0, y = 0, z = 0;
1758  for( int j = 0; j < 3; j++ )
1759  {
1760  x += N[j] * corner_coords[3 * j];
1761  y += N[j] * corner_coords[3 * j + 1];
1762  z += N[j] * corner_coords[3 * j + 2];
1763  }
1764 
1765  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1766  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1767  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1768  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1769  }
1770  }
1771  else if( type == MBQUAD )
1772  {
1773  double xi, eta, N[4];
1774  int findex = mbImpl->type_from_handle( *( _infaces.begin() ) ) - 1;
1775 
1776  for( int i = 4; i < vtotal; i++ )
1777  {
1778  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1779 
1780  xi = refTemplates[findex][d].vert_nat_coord[i - 4][0];
1781  eta = refTemplates[findex][d].vert_nat_coord[i - 4][1];
1782  N[0] = ( 1 - xi ) * ( 1 - eta ) / 4;
1783  N[1] = ( 1 + xi ) * ( 1 - eta ) / 4;
1784  N[2] = ( 1 + xi ) * ( 1 + eta ) / 4, N[3] = ( 1 - xi ) * ( 1 + eta ) / 4;
1785 
1786  double x = 0, y = 0, z = 0;
1787  for( int j = 0; j < 4; j++ )
1788  {
1789  x += N[j] * corner_coords[3 * j];
1790  y += N[j] * corner_coords[3 * j + 1];
1791  z += N[j] * corner_coords[3 * j + 2];
1792  }
1793 
1794  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1795  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1796  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1797  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1798  }
1799  }
1800  else if( type == MBTET )
1801  {
1802  double xi, eta, mu, N[4];
1803  int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
1804 
1805  for( int i = 4; i < vtotal; i++ )
1806  {
1807  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1808 
1809  xi = refTemplates[cindex][d].vert_nat_coord[i - 4][0];
1810  eta = refTemplates[cindex][d].vert_nat_coord[i - 4][1];
1811  mu = refTemplates[cindex][d].vert_nat_coord[i - 4][2];
1812 
1813  N[0] = 1 - xi - eta - mu;
1814  N[1] = xi;
1815  N[2] = eta, N[3] = mu;
1816 
1817  double x = 0, y = 0, z = 0;
1818  for( int j = 0; j < 4; j++ )
1819  {
1820  x += N[j] * corner_coords[3 * j];
1821  y += N[j] * corner_coords[3 * j + 1];
1822  z += N[j] * corner_coords[3 * j + 2];
1823  }
1824 
1825  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1826  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1827  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1828  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1829  }
1830  }
1831  else if( type == MBPRISM )
1832  {
1833  double xi, eta, mu, N[6];
1834  int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
1835 
1836  for( int i = 6; i < vtotal; i++ )
1837  {
1838  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1839 
1840  xi = refTemplates[cindex][d].vert_nat_coord[i - 6][0];
1841  eta = refTemplates[cindex][d].vert_nat_coord[i - 6][1];
1842  mu = refTemplates[cindex][d].vert_nat_coord[i - 6][2];
1843 
1844  N[0] = ( 1 - xi - eta ) * ( 1 - mu ), N[1] = xi * ( 1 - mu ), N[2] = eta * ( 1 - mu ),
1845  N[3] = ( 1 - xi - eta ) * ( 1 + mu ), N[4] = xi * ( 1 + mu ), N[5] = eta * ( 1 + mu );
1846 
1847  double x = 0, y = 0, z = 0;
1848  for( int j = 0; j < 6; j++ )
1849  {
1850  x += N[j] * corner_coords[3 * j];
1851  y += N[j] * corner_coords[3 * j + 1];
1852  z += N[j] * corner_coords[3 * j + 2];
1853  }
1854 
1855  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1856  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1857  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1858  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1859  }
1860  }
1861  else if( type == MBHEX )
1862  {
1863  double xi, eta, mu, N[8];
1864  double d1, d2, d3, s1, s2, s3;
1865  int cindex = mbImpl->type_from_handle( *( _incells.begin() ) ) - 1;
1866 
1867  for( int i = 8; i < vtotal; i++ )
1868  {
1869 
1870  if( vflag[vbuffer[i] - vstart - nverts_prev] ) continue;
1871 
1872  xi = refTemplates[cindex][d].vert_nat_coord[i - 8][0];
1873  eta = refTemplates[cindex][d].vert_nat_coord[i - 8][1];
1874  mu = refTemplates[cindex][d].vert_nat_coord[i - 8][2];
1875 
1876  d1 = 1 - xi;
1877  d2 = 1 - eta;
1878  d3 = 1 - mu;
1879  s1 = 1 + xi;
1880  s2 = 1 + eta;
1881  s3 = 1 + mu;
1882  N[0] = ( d1 * d2 * d3 ) / 8;
1883  N[1] = ( s1 * d2 * d3 ) / 8;
1884  N[2] = ( s1 * s2 * d3 ) / 8;
1885  N[3] = ( d1 * s2 * d3 ) / 8;
1886  N[4] = ( d1 * d2 * s3 ) / 8;
1887  N[5] = ( s1 * d2 * s3 ) / 8;
1888  N[6] = ( s1 * s2 * s3 ) / 8;
1889  N[7] = ( d1 * s2 * s3 ) / 8;
1890 
1891  double x = 0, y = 0, z = 0;
1892  for( int j = 0; j < 8; j++ )
1893  {
1894  x += N[j] * corner_coords[3 * j];
1895  y += N[j] * corner_coords[3 * j + 1];
1896  z += N[j] * corner_coords[3 * j + 2];
1897  }
1898 
1899  level_mesh[cur_level].coordinates[0][vbuffer[i] - vstart] = x;
1900  level_mesh[cur_level].coordinates[1][vbuffer[i] - vstart] = y;
1901  level_mesh[cur_level].coordinates[2][vbuffer[i] - vstart] = z;
1902 
1903  vflag[vbuffer[i] - vstart - nverts_prev] = 1;
1904  }
1905  }
1906  return MB_SUCCESS;
1907 }
1908 /**********************************
1909  * Parallel Communication *
1910  * ********************************/
1911 
1912 #ifdef MOAB_HAVE_MPI
1913 ErrorCode NestedRefine::resolve_shared_ents_parmerge( int level, EntityHandle levelset )
1914 {
1915  // TEMP: Add the adjacencies for MOAB-native DS
1916  // NOTE (VSM): This is expensive since it creates a doubly
1917  // redundant copy of the adjacency data in both MOAB-native
1918  // and AHF. Need to fix this with AHF optimized branch.
1919  ErrorCode error;
1920  ReadUtilIface* read_iface;
1921  error = mbImpl->query_interface( read_iface );MB_CHK_ERR( error );
1922  if( level_mesh[level].num_edges != 0 )
1923  {
1924  error = read_iface->update_adjacencies( level_mesh[level].start_edge, level_mesh[level].num_edges, 2,
1925  level_mesh[level].edge_conn );MB_CHK_ERR( error );
1926  }
1927  if( level_mesh[level].num_faces != 0 )
1928  {
1929  EntityType type = mbImpl->type_from_handle( *( _infaces.begin() ) );
1930  int nvpf = ahf->lConnMap2D[type - 2].num_verts_in_face;
1931  error = read_iface->update_adjacencies( level_mesh[level].start_face, level_mesh[level].num_faces, nvpf,
1932  level_mesh[level].face_conn );MB_CHK_ERR( error );
1933  }
1934  if( level_mesh[level].num_cells != 0 )
1935  {
1936  int index = ahf->get_index_in_lmap( *_incells.begin() );
1937  int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
1938  error = read_iface->update_adjacencies( level_mesh[level].start_cell, level_mesh[level].num_cells, nvpc,
1939  level_mesh[level].cell_conn );MB_CHK_ERR( error );
1940  }
1941 
1942  if( pcomm->size() > 1 )
1943  {
1944 
1945  // get all entities on the rootset
1946  moab::Range vtxs, edgs, facs, elms;
1947  error = mbImpl->get_entities_by_dimension( levelset, 0, vtxs, false );MB_CHK_ERR( error );
1948  error = mbImpl->get_entities_by_dimension( levelset, 1, edgs, false );MB_CHK_ERR( error );
1949  error = mbImpl->get_entities_by_dimension( levelset, 2, facs, false );MB_CHK_ERR( error );
1950  error = mbImpl->get_entities_by_dimension( levelset, 3, elms, false );MB_CHK_ERR( error );
1951 
1952  // set the parallel partition tag data
1953  moab::Tag part_tag;
1954  int partid = pcomm->rank(), dum_id = -1;
1955  error = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, part_tag,
1957  error = mbImpl->tag_set_data( part_tag, &levelset, 1, &partid );MB_CHK_ERR( error );
1958  //
1959  // Now that we have the local piece of the mesh refined consistently,
1960  // call parallel merge instead of resolved_shared to stitch together the meshes
1961  // and to handle parallel communication of remote proc/entity-handle pairs for
1962  // shared + non-owned interfaces entities.
1963  //
1964  // TODO: This needs to be replaced by the following scheme in the future as
1965  // an optimization step.
1966  // > Assign global IDs consistently for shared entities so that parallel
1967  // > resolve shared ents can happen out of the box. This is the fastest option.
1968 
1969  ParallelMergeMesh pm( pcomm, 1e-08 );
1970  error = pm.merge( levelset, true );MB_CHK_ERR( error );
1971 
1972  //
1973  // Parallel Communication complete - all entities resolved
1974  //
1975  }
1976  return MB_SUCCESS;
1977 }
1978 
1979 ErrorCode NestedRefine::resolve_shared_ents_opt( EntityHandle* hm_set, int num_levels )
1980 {
1981  assert( pcomm->size() > 1 );
1982 
1983  ErrorCode error;
1984 
1985  // Step 1A: Pre-processing: setting the parallel partition tag data
1986  for( int i = 0; i < num_levels; i++ )
1987  {
1988  Tag part_tag;
1989  int partid = pcomm->rank(), dum_id = -1;
1990  error = mbImpl->tag_get_handle( "PARALLEL_PARTITION", 1, moab::MB_TYPE_INTEGER, part_tag,
1992  error = mbImpl->tag_set_data( part_tag, &hm_set[i], 1, &partid );MB_CHK_ERR( error );
1993  }
1994 
1995  // Step 1B: Pre-processing: gather all shared entities and list entities shared with each
1996  // sharing processor
1997 
1998  // All shared processors
1999  std::set< unsigned int > shprocs;
2000  error = pcomm->get_comm_procs( shprocs );MB_CHK_ERR( error );
2001 
2002  std::vector< int > sharedprocs;
2003  for( std::set< unsigned int >::iterator it = shprocs.begin(); it != shprocs.end(); it++ )
2004  sharedprocs.push_back( *it );
2005  int nprocs = sharedprocs.size();
2006 
2007  // Create buffer variables storing the entities to be sent and received
2008  std::vector< std::vector< int > > nsharedEntsperproc( nprocs );
2009  std::vector< std::vector< EntityHandle > > localBuffs( nprocs );
2010  std::vector< std::vector< EntityHandle > > remlocalBuffs( nprocs );
2011  std::vector< std::vector< EntityHandle > > remoteBuffs( nprocs );
2012 
2013  int i;
2014  Range sharedentities;
2015 
2016  for( i = 0; i < nprocs; i++ )
2017  {
2018  // List of shared entities at the coarsest level
2019  sharedentities.clear();
2020  error = pcomm->get_shared_entities( sharedprocs[i], sharedentities, -1, true );MB_CHK_ERR( error );
2021 
2022  // Get the list shared edges and vertices that are not part of the shared edges
2023  Range allEnts;
2024  error = collect_shared_entities_by_dimension( sharedentities, allEnts );MB_CHK_ERR( error );
2025 
2026  Range V0, E0, F0;
2027  V0 = allEnts.subset_by_dimension( 0 );
2028  E0 = allEnts.subset_by_dimension( 1 );
2029  F0 = allEnts.subset_by_dimension( 2 );
2030 
2031  // Step 2A: Prepare msg to be sent:
2032  //
2033  // FList = <F, FC, FE> where F, FC and FE are vectors containing the faces, their
2034  // connectivities and edges. F = <F_0, F_1, F_2, ..., F_L> where
2035  // F_0 is the list of shared faces at the coarsest level,
2036  // F_i are the list of children faces at subsequent levels.
2037  // FC vector contains the connectivities of F_i's and FE vector contains the bounding edges
2038  // of F_i's.
2039  //
2040  // EList = <E, EC> where E and EC are vectors containing the edges and their connectivities.
2041  // E = <E_0, E_1, ..., E_L> where
2042  // E_0 is the list of shared edges at the coarsest level and
2043  // E_i are the list of children edges at subsequent levels.
2044  // The EC vector contains the connectivities of E_i's.
2045  //
2046  // VList = <V> = <V_0, V_1, ...., V_L> where V_0 are shared vertices at the coarsest level
2047  // and V_i are the duplicates in the subsequent levels.
2048 
2049  std::vector< EntityHandle > locFList, remFList;
2050  std::vector< EntityHandle > locEList, remEList;
2051  std::vector< EntityHandle > locVList, remVList;
2052 
2053  // collect faces
2054  if( !F0.empty() )
2055  {
2056  error = collect_FList( sharedprocs[i], F0, locFList, remFList );MB_CHK_ERR( error );
2057  }
2058 
2059  // collect edges
2060  if( !E0.empty() )
2061  {
2062  error = collect_EList( sharedprocs[i], E0, locEList, remEList );MB_CHK_ERR( error );
2063  }
2064 
2065  // collect vertices
2066  if( !V0.empty() )
2067  {
2068  error = collect_VList( sharedprocs[i], V0, locVList, remVList );MB_CHK_ERR( error );
2069  }
2070 
2071  // Step 2B: Add data to NR local buffer to be sent
2072  std::vector< int > msgsz;
2073  msgsz.push_back( F0.size() );
2074  msgsz.push_back( locFList.size() );
2075  msgsz.push_back( E0.size() );
2076  msgsz.push_back( locEList.size() );
2077  msgsz.push_back( V0.size() );
2078  msgsz.push_back( locVList.size() );
2079  nsharedEntsperproc[i].insert( nsharedEntsperproc[i].end(), msgsz.begin(), msgsz.end() );
2080 
2081  if( !F0.empty() )
2082  {
2083  localBuffs[i].insert( localBuffs[i].end(), locFList.begin(), locFList.end() );
2084  remlocalBuffs[i].insert( remlocalBuffs[i].end(), remFList.begin(), remFList.end() );
2085  }
2086  if( !E0.empty() )
2087  {
2088  localBuffs[i].insert( localBuffs[i].end(), locEList.begin(), locEList.end() );
2089  remlocalBuffs[i].insert( remlocalBuffs[i].end(), remEList.begin(), remEList.end() );
2090  }
2091  if( !V0.empty() )
2092  {
2093  localBuffs[i].insert( localBuffs[i].end(), locVList.begin(), locVList.end() );
2094  remlocalBuffs[i].insert( remlocalBuffs[i].end(), remVList.begin(), remVList.end() );
2095  }
2096  }
2097 
2098  // Step 3: Send and receive the remote collection of child ents
2099  error = pcomm->send_recv_entities( sharedprocs, nsharedEntsperproc, remlocalBuffs, remoteBuffs );MB_CHK_ERR( error );
2100 
2101  // Step 5: Resolve shared child entities and update parallel tags
2102  std::multimap< EntityHandle, int > rprocs;
2103  std::multimap< EntityHandle, EntityHandle > rhandles;
2104 
2105  error = decipher_remote_handles( sharedprocs, nsharedEntsperproc, localBuffs, remoteBuffs, rprocs, rhandles );MB_CHK_ERR( error );
2106 
2107  // Step 6: Update pcomm tags
2108  error = update_parallel_tags( rprocs, rhandles );MB_CHK_ERR( error );
2109 
2110  return MB_SUCCESS;
2111 }
2112 
2113 ErrorCode NestedRefine::collect_shared_entities_by_dimension( Range sharedEnts, Range& allEnts )
2114 {
2115  ErrorCode error;
2116 
2117  Range F0, E0, V0, E0all, V0all;
2118  std::vector< EntityHandle > ents;
2119 
2120  F0 = sharedEnts.subset_by_dimension( 2 );
2121  E0all = sharedEnts.subset_by_dimension( 1 );
2122  V0all = sharedEnts.subset_by_dimension( 0 );
2123 
2124  if( !F0.empty() )
2125  {
2126  Range edges, verts;
2127 
2128  for( Range::iterator it = F0.begin(); it != F0.end(); it++ )
2129  {
2130  ents.clear();
2131  error = ahf->get_adjacencies( *it, 1, ents );MB_CHK_ERR( error );
2132  std::copy( ents.begin(), ents.end(), range_inserter( edges ) );
2133  ents.clear();
2134  error = mbImpl->get_connectivity( &( *it ), 0, ents );MB_CHK_ERR( error );
2135  std::copy( ents.begin(), ents.end(), range_inserter( verts ) );
2136  }
2137 
2138  E0 = subtract( E0all, edges );
2139  if( !E0.empty() )
2140  {
2141  for( Range::iterator it = E0.begin(); it != E0.end(); it++ )
2142  {
2143  ents.clear();
2144  error = mbImpl->get_connectivity( &( *it ), 0, ents );MB_CHK_ERR( error );
2145  std::copy( ents.begin(), ents.end(), range_inserter( verts ) );
2146  }
2147  }
2148  V0 = subtract( V0all, verts );
2149  }
2150  else if( !E0all.empty() )
2151  {
2152  Range verts;
2153  for( Range::iterator it = E0all.begin(); it != E0all.end(); it++ )
2154  {
2155  ents.clear();
2156  error = mbImpl->get_connectivity( &( *it ), 1, ents );MB_CHK_ERR( error );
2157  std::copy( ents.begin(), ents.end(), range_inserter( verts ) );
2158  }
2159  E0 = E0all;
2160  V0 = subtract( V0all, verts );
2161  }
2162  else if( !V0all.empty() )
2163  {
2164  V0 = V0all;
2165  }
2166  else
2167  MB_SET_ERR( MB_FAILURE, "Trying to pack unsupported sub-entities for shared interface entities" );
2168 
2169  if( !F0.empty() ) std::copy( F0.begin(), F0.end(), range_inserter( allEnts ) );
2170  if( !E0.empty() ) std::copy( E0.begin(), E0.end(), range_inserter( allEnts ) );
2171  if( !V0.empty() ) std::copy( V0.begin(), V0.end(), range_inserter( allEnts ) );
2172 
2173  return MB_SUCCESS;
2174 }
2175 
2176 ErrorCode NestedRefine::collect_FList( int to_proc,
2177  Range faces,
2178  std::vector< EntityHandle >& FList,
2179  std::vector< EntityHandle >& RList )
2180 {
2181  ErrorCode error;
2182 
2183  FList.clear();
2184  std::vector< EntityHandle > F, FC, FE, lF, lFC, lFE, rF, rFC, rFE;
2185  std::vector< EntityHandle > childEnts, conn, fedges;
2186 
2187  for( Range::iterator it = faces.begin(); it != faces.end(); it++ )
2188  {
2189  EntityHandle face = *it;
2190  conn.clear();
2191  fedges.clear();
2192  error = mbImpl->get_connectivity( &face, 1, conn );MB_CHK_ERR( error );
2193  error = ahf->get_face_edges( *it, fedges );MB_CHK_ERR( error );
2194 
2195  // add local handles
2196  lF.push_back( *it );
2197  lFC.insert( lFC.end(), conn.begin(), conn.end() );
2198  lFE.insert( lFE.end(), fedges.begin(), fedges.end() );
2199 
2200  // replace local handles with remote handles on to_proc
2201  EntityHandle rval = 0;
2202  error = pcomm->get_remote_handles( &face, &rval, 1, to_proc );MB_CHK_ERR( error );
2203  rF.push_back( rval );
2204 
2205  for( int i = 0; i < (int)conn.size(); i++ )
2206  {
2207  error = pcomm->get_remote_handles( &conn[i], &rval, 1, to_proc );MB_CHK_ERR( error );
2208  rFC.push_back( rval );
2209 
2210  error = pcomm->get_remote_handles( &fedges[i], &rval, 1, to_proc );MB_CHK_ERR( error );
2211  rFE.push_back( rval );
2212  }
2213  }
2214 
2215  for( int l = 0; l < nlevels; l++ )
2216  {
2217  for( Range::iterator it = faces.begin(); it != faces.end(); it++ )
2218  {
2219  childEnts.clear();
2220  error = parent_to_child( *it, 0, l + 1, childEnts );MB_CHK_ERR( error );
2221 
2222  for( int i = 0; i < (int)childEnts.size(); i++ )
2223  {
2224  conn.clear();
2225  fedges.clear();
2226  error = mbImpl->get_connectivity( &childEnts[i], 1, conn );MB_CHK_ERR( error );
2227  error = ahf->get_face_edges( childEnts[i], fedges );MB_CHK_ERR( error );
2228 
2229  F.push_back( childEnts[i] );
2230  FC.insert( FC.end(), conn.begin(), conn.end() );
2231  FE.insert( FE.end(), fedges.begin(), fedges.end() );
2232  }
2233  }
2234  }
2235 
2236  FList.insert( FList.end(), lF.begin(), lF.end() );
2237  FList.insert( FList.end(), F.begin(), F.end() );
2238  FList.insert( FList.end(), lFC.begin(), lFC.end() );
2239  FList.insert( FList.end(), FC.begin(), FC.end() );
2240  FList.insert( FList.end(), lFE.begin(), lFE.end() );
2241  FList.insert( FList.end(), FE.begin(), FE.end() );
2242 
2243  RList.insert( RList.end(), rF.begin(), rF.end() );
2244  RList.insert( RList.end(), F.begin(), F.end() );
2245  RList.insert( RList.end(), rFC.begin(), rFC.end() );
2246  RList.insert( RList.end(), FC.begin(), FC.end() );
2247  RList.insert( RList.end(), rFE.begin(), rFE.end() );
2248  RList.insert( RList.end(), FE.begin(), FE.end() );
2249 
2250  return MB_SUCCESS;
2251 }
2252 
2253 ErrorCode NestedRefine::collect_EList( int to_proc,
2254  Range edges,
2255  std::vector< EntityHandle >& EList,
2256  std::vector< EntityHandle >& RList )
2257 {
2258  ErrorCode error;
2259  EList.clear();
2260  std::vector< EntityHandle > E, EC, lE, lEC, rE, rEC;
2261  std::vector< EntityHandle > childEnts, conn;
2262 
2263  // Add the edges and their connectivities at the coarsest level first.
2264  for( Range::iterator it = edges.begin(); it != edges.end(); it++ )
2265  {
2266  EntityHandle edg = *it;
2267  conn.clear();
2268  error = mbImpl->get_connectivity( &edg, 1, conn );MB_CHK_ERR( error );
2269 
2270  // add local handles
2271  lE.push_back( edg );
2272  lEC.insert( lEC.end(), conn.begin(), conn.end() );
2273 
2274  // replace local handles with remote handle on to_proc
2275  EntityHandle rval = 0;
2276  error = pcomm->get_remote_handles( &edg, &rval, 1, to_proc );MB_CHK_ERR( error );
2277  rE.push_back( rval );
2278  error = pcomm->get_remote_handles( &conn[0], &rval, 1, to_proc );MB_CHK_ERR( error );
2279  rEC.push_back( rval );
2280  error = pcomm->get_remote_handles( &conn[1], &rval, 1, to_proc );MB_CHK_ERR( error );
2281  rEC.push_back( rval );
2282  }
2283 
2284  // Add the edges and their connectivities at subsequent levels.
2285  for( int l = 0; l < nlevels; l++ )
2286  {
2287  for( Range::iterator it = edges.begin(); it != edges.end(); it++ )
2288  {
2289  childEnts.clear();
2290  error = parent_to_child( *it, 0, l + 1, childEnts );MB_CHK_ERR( error );
2291 
2292  for( int i = 0; i < (int)childEnts.size(); i++ )
2293  {
2294  conn.clear();
2295  error = mbImpl->get_connectivity( &childEnts[i], 1, conn );MB_CHK_ERR( error );
2296  E.push_back( childEnts[i] );
2297  EC.insert( EC.end(), conn.begin(), conn.end() );
2298  }
2299  }
2300  }
2301 
2302  EList.insert( EList.end(), lE.begin(), lE.end() );
2303  EList.insert( EList.end(), E.begin(), E.end() );
2304  EList.insert( EList.end(), lEC.begin(), lEC.end() );
2305  EList.insert( EList.end(), EC.begin(), EC.end() );
2306 
2307  RList.insert( RList.end(), rE.begin(), rE.end() );
2308  RList.insert( RList.end(), E.begin(), E.end() );
2309  RList.insert( RList.end(), rEC.begin(), rEC.end() );
2310  RList.insert( RList.end(), EC.begin(), EC.end() );
2311 
2312  return MB_SUCCESS;
2313 }
2314 
2315 ErrorCode NestedRefine::collect_VList( int to_proc,
2316  Range verts,
2317  std::vector< EntityHandle >& VList,
2318  std::vector< EntityHandle >& RList )
2319 {
2320  ErrorCode error;
2321  std::vector< EntityHandle > V, lV, rV;
2322  VList.clear();
2323 
2324  // Add the vertices at the coarsest level first.
2325  for( Range::iterator it = verts.begin(); it != verts.end(); it++ )
2326  {
2327  EntityHandle v = *it;
2328 
2329  lV.push_back( v );
2330  EntityHandle rval = 0;
2331  error = pcomm->get_remote_handles( &v, &rval, 1, to_proc );MB_CHK_ERR( error );
2332 
2333  rV.push_back( rval );
2334  }
2335 
2336  // Add the vertices at the subsequent levels .
2337  for( int l = 0; l < nlevels; l++ )
2338  {
2339  for( Range::iterator it = verts.begin(); it != verts.end(); it++ )
2340  {
2341  EntityHandle dupvert = 0;
2342  error = get_vertex_duplicates( *it, l + 1, dupvert );MB_CHK_ERR( error );
2343  V.push_back( dupvert );
2344  }
2345  }
2346 
2347  // local vertex handles at the coarsest level
2348  VList.insert( VList.end(), lV.begin(), lV.end() );
2349  VList.insert( VList.end(), V.begin(), V.end() );
2350 
2351  // remote vertex handles at the coarsest level
2352  RList.insert( RList.end(), rV.begin(), rV.end() );
2353  RList.insert( RList.end(), V.begin(), V.end() );
2354 
2355  return MB_SUCCESS;
2356 }
2357 
2358 ErrorCode NestedRefine::decipher_remote_handles( std::vector< int >& sharedprocs,
2359  std::vector< std::vector< int > >& auxinfo,
2360  std::vector< std::vector< EntityHandle > >& localbuffers,
2361  std::vector< std::vector< EntityHandle > >& remotebuffers,
2362  std::multimap< EntityHandle, int >& remProcs,
2363  std::multimap< EntityHandle, EntityHandle >& remHandles )
2364 {
2365  ErrorCode error;
2366 
2367  int i;
2368  int nprocs = sharedprocs.size();
2369 
2370  for( i = 0; i < nprocs; i++ )
2371  {
2372  std::vector< int > msgsz;
2373  for( int j = 0; j < (int)auxinfo[i].size(); j++ )
2374  {
2375  msgsz.push_back( auxinfo[i][j] );
2376  }
2377 
2378  if( msgsz[0] != 0 ) // Faces
2379  {
2380  // Get the local and remote face handles from the buffers
2381  std::vector< EntityHandle > LFList, RFList;
2382  LFList.insert( LFList.end(), localbuffers[i].begin(), localbuffers[i].begin() + msgsz[1] );
2383  RFList.insert( RFList.end(), remotebuffers[i].begin(), remotebuffers[i].begin() + msgsz[1] );
2384 
2385  error = decipher_remote_handles_face( sharedprocs[i], msgsz[0], LFList, RFList, remProcs, remHandles );MB_CHK_ERR( error );
2386 
2387  if( msgsz[2] != 0 ) // Edges
2388  {
2389  std::vector< EntityHandle > LEList, REList;
2390  LEList.insert( LEList.end(), localbuffers[i].begin() + msgsz[1] + 1,
2391  localbuffers[i].begin() + msgsz[1] + msgsz[3] );
2392  REList.insert( REList.end(), remotebuffers[i].begin() + msgsz[1] + 1,
2393  remotebuffers[i].begin() + msgsz[1] + msgsz[3] );
2394 
2395  error = decipher_remote_handles_edge( sharedprocs[i], msgsz[2], LEList, REList, remProcs, remHandles );MB_CHK_ERR( error );
2396 
2397  if( msgsz[4] != 0 ) // Vertices
2398  {
2399  // Get the local and remote face handles from the buffers
2400  std::vector< EntityHandle > LVList, RVList;
2401  LVList.insert( LVList.end(), localbuffers[i].begin() + msgsz[1] + msgsz[3] + 1,
2402  localbuffers[i].begin() + msgsz[1] + msgsz[3] + msgsz[5] );
2403  RVList.insert( RVList.end(), remotebuffers[i].begin() + msgsz[1] + msgsz[3] + 1,
2404  remotebuffers[i].begin() + msgsz[1] + msgsz[3] + msgsz[5] );
2405 
2406  error = decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs,
2407  remHandles );MB_CHK_ERR( error );
2408  }
2409  }
2410  else if( msgsz[4] != 0 ) // Vertices
2411  {
2412  // Get the local and remote face handles from the buffers
2413  std::vector< EntityHandle > LVList, RVList;
2414  LVList.insert( LVList.end(), localbuffers[i].begin() + msgsz[1] + 1,
2415  localbuffers[i].begin() + msgsz[1] + msgsz[5] );
2416  RVList.insert( RVList.end(), remotebuffers[i].begin() + msgsz[1] + 1,
2417  remotebuffers[i].begin() + msgsz[1] + msgsz[5] );
2418 
2419  error =
2420  decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs, remHandles );MB_CHK_ERR( error );
2421  }
2422  }
2423 
2424  else if( msgsz[2] != 0 ) // Edges
2425  {
2426  // Get the local and remote face handles from the buffers
2427  std::vector< EntityHandle > LEList, REList;
2428  LEList.insert( LEList.end(), localbuffers[i].begin(), localbuffers[i].begin() + msgsz[3] );
2429  REList.insert( REList.end(), remotebuffers[i].begin(), remotebuffers[i].begin() + msgsz[3] );
2430 
2431  error = decipher_remote_handles_edge( sharedprocs[i], msgsz[2], LEList, REList, remProcs, remHandles );MB_CHK_ERR( error );
2432 
2433  if( msgsz[4] != 0 ) // Vertices
2434  {
2435  // Get the local and remote face handles from the buffers
2436  std::vector< EntityHandle > LVList, RVList;
2437  LVList.insert( LVList.end(), localbuffers[i].begin() + msgsz[3] + 1,
2438  localbuffers[i].begin() + msgsz[3] + msgsz[5] );
2439  RVList.insert( RVList.end(), remotebuffers[i].begin() + msgsz[3] + 1,
2440  remotebuffers[i].begin() + msgsz[3] + msgsz[5] );
2441 
2442  error =
2443  decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs, remHandles );MB_CHK_ERR( error );
2444  }
2445  }
2446 
2447  else if( msgsz[4] != 0 ) // Vertices
2448  {
2449  // Get the local and remote face handles from the buffers
2450  std::vector< EntityHandle > LVList, RVList;
2451  LVList.insert( LVList.end(), localbuffers[i].begin(), localbuffers[i].end() );
2452  RVList.insert( RVList.end(), remotebuffers[i].begin(), remotebuffers[i].end() );
2453 
2454  error = decipher_remote_handles_vertex( sharedprocs[i], msgsz[4], LVList, RVList, remProcs, remHandles );MB_CHK_ERR( error );
2455  }
2456  else
2457  MB_SET_ERR( MB_FAILURE, "Trying to decipher entities other than verts, edges, faces" );
2458  }
2459 
2460  return MB_SUCCESS;
2461 }
2462 
2463 ErrorCode NestedRefine::decipher_remote_handles_face( int shared_proc,
2464  int numfaces,
2465  std::vector< EntityHandle >& localFaceList,
2466  std::vector< EntityHandle >& remFaceList,
2467  std::multimap< EntityHandle, int >& remProcs,
2468  std::multimap< EntityHandle, EntityHandle >& remHandles )
2469 {
2470  ErrorCode error;
2471 
2472  for( int i = 0; i < numfaces; i++ )
2473  {
2474  // find local and remote handles of the coarsest face
2475  EntityHandle Lface = localFaceList[i];
2476  int Rface_idx =
2477  ( std::find( remFaceList.begin(), remFaceList.begin() + numfaces - 1, Lface ) ) - remFaceList.begin();
2478 
2479  // get connectivities of the local and remote coarsest faces
2480  std::vector< EntityHandle > Lface_conn, Rface_conn;
2481  error = get_data_from_buff( 2, 1, 0, i, numfaces, localFaceList, Lface_conn );MB_CHK_ERR( error );
2482  error = get_data_from_buff( 2, 1, 0, Rface_idx, numfaces, remFaceList, Rface_conn );MB_CHK_ERR( error );
2483 
2484  // find the combination difference between local and remote coarsest face
2485  std::vector< int > cmap;
2486  int comb = 0;
2487  int nvF = (int)Lface_conn.size();
2488  error = reorder_indices( &Lface_conn[0], &Rface_conn[0], nvF, &cmap[0], comb );MB_CHK_ERR( error );
2489 
2490  // go into loop over all levels
2491  std::vector< EntityHandle > lchildents, lparents;
2492  std::vector< EntityHandle > lcents, rcents;
2493  int lidx, ridx;
2494 
2495  for( int l = 0; l < nlevels; l++ )
2496  {
2497  lchildents.clear();
2498  lparents.clear();
2499  lcents.clear();
2500  rcents.clear();
2501 
2502  // obtain children at the current level
2503  error = get_data_from_buff( 2, 0, l + 1, i, numfaces, localFaceList, lchildents );MB_CHK_ERR( error );
2504 
2505  // obtain parents at the previous level
2506  if( l == 0 )
2507  {
2508  lparents.push_back( Lface );
2509  }
2510  else
2511  {
2512  error = get_data_from_buff( 2, 0, l, i, numfaces, localFaceList, lparents );MB_CHK_ERR( error );
2513  }
2514 
2515  //#children at the previous level and the current level
2516  EntityType ftype = mbImpl->type_from_handle( Lface );
2517  int d = get_index_from_degree( level_dsequence[l] );
2518  int nch = refTemplates[ftype - 1][d].total_new_ents;
2519  std::vector< int > fmap;
2520  error = reorder_indices( level_dsequence[l], nvF, comb, &fmap[0] );MB_CHK_ERR( error );
2521 
2522  // loop over all the lparents
2523  for( int j = 0; j < (int)lparents.size(); j++ )
2524  {
2525  // list local childrent at the current level
2526  lidx = std::find( localFaceList.begin(), localFaceList.end(), lparents[j] ) - localFaceList.begin();
2527  error = get_data_from_buff( 2, 0, l + 1, lidx, numfaces, localFaceList, lcents );MB_CHK_ERR( error );
2528 
2529  // find the corresponding remote of lparent and its children
2530  EntityHandle rparent = 0;
2531  error = check_for_parallelinfo( lparents[j], shared_proc, remHandles, remProcs, rparent );MB_CHK_ERR( error );
2532  ridx = std::find( remFaceList.begin(), remFaceList.end(), rparent ) - remFaceList.begin();
2533  error = get_data_from_buff( 2, 0, l + 1, ridx, numfaces, remFaceList, rcents );MB_CHK_ERR( error );
2534 
2535  // match up local face with remote handles according to cmap
2536  std::vector< EntityHandle > lconn, rconn, ledg, redg;
2537  for( int k = 0; k < nch; k++ )
2538  {
2539  lconn.clear();
2540  rconn.clear();
2541  ledg.clear();
2542  redg.clear();
2543 
2544  // matching the local children face handles with their remotes
2545  bool found = check_for_parallelinfo( lcents[k], shared_proc, remProcs );
2546  if( !found )
2547  {
2548  remProcs.insert( std::pair< EntityHandle, int >( lcents[k], shared_proc ) );
2549  remHandles.insert( std::pair< EntityHandle, EntityHandle >( lcents[k], rcents[fmap[k]] ) );
2550  }
2551 
2552  // find indices of the matched child face
2553  lidx = std::find( localFaceList.begin(), localFaceList.end(), lcents[k] ) - localFaceList.begin();
2554  ridx = std::find( remFaceList.begin(), remFaceList.end(), rcents[fmap[k]] ) - remFaceList.begin();
2555 
2556  // find bounding edges and connectivity of the matched child face
2557  error = get_data_from_buff( 2, 2, l + 1, lidx, numfaces, localFaceList, ledg );MB_CHK_ERR( error );
2558  error = get_data_from_buff( 2, 1, l + 1, lidx, numfaces, localFaceList, lconn );MB_CHK_ERR( error );
2559  error = get_data_from_buff( 2, 2, l + 1, ridx, numfaces, remFaceList, redg );MB_CHK_ERR( error );
2560  error = get_data_from_buff( 2, 1, l + 1, ridx, numfaces, remFaceList, rconn );MB_CHK_ERR( error );
2561 
2562  // now match the handles of the bounding edges and the vertices using
2563  // combination difference
2564  for( int m = 0; m < (int)ledg.size(); m++ )
2565  {
2566  found = check_for_parallelinfo( ledg[m], shared_proc, remProcs );
2567 
2568  if( !found )
2569  {
2570  remProcs.insert( std::pair< EntityHandle, int >( ledg[m], shared_proc ) );
2571  remHandles.insert( std::pair< EntityHandle, EntityHandle >( ledg[m], redg[cmap[m]] ) );
2572  }
2573 
2574  found = check_for_parallelinfo( lconn[m], shared_proc, remProcs );
2575 
2576  if( !found )
2577  {
2578  remProcs.insert( std::pair< EntityHandle, int >( lconn[m], shared_proc ) );
2579  remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[m], rconn[cmap[m]] ) );
2580  }
2581  }
2582  }
2583  }
2584  }
2585  }
2586 
2587  return MB_SUCCESS;
2588 }
2589 
2590 ErrorCode NestedRefine::decipher_remote_handles_edge( int shared_proc,
2591  int numedges,
2592  std::vector< EntityHandle >& localEdgeList,
2593  std::vector< EntityHandle >& remEdgeList,
2594  std::multimap< EntityHandle, int >& remProcs,
2595  std::multimap< EntityHandle, EntityHandle >& remHandles )
2596 {
2597  ErrorCode error;
2598 
2599  for( int i = 0; i < numedges; i++ )
2600  {
2601  EntityHandle Ledge = localEdgeList[i];
2602  int Redge_idx =
2603  ( std::find( remEdgeList.begin(), remEdgeList.begin() + numedges - 1, Ledge ) ) - remEdgeList.begin();
2604 
2605  std::vector< EntityHandle > Ledge_conn, Redge_conn;
2606  error = get_data_from_buff( 1, 1, 0, i, numedges, localEdgeList, Ledge_conn );MB_CHK_ERR( error );
2607  error = get_data_from_buff( 1, 1, 0, Redge_idx, numedges, remEdgeList, Redge_conn );MB_CHK_ERR( error );
2608 
2609  bool orient = true;
2610  if( ( Ledge_conn[0] == Redge_conn[1] ) && ( Ledge_conn[1] == Redge_conn[0] ) ) orient = false;
2611 
2612  if( orient ) assert( ( Ledge_conn[0] == Redge_conn[0] ) && ( Ledge_conn[1] == Redge_conn[1] ) );
2613 
2614  std::vector< EntityHandle > lchildEdgs, rchildEdgs, lconn, rconn;
2615  for( int l = 0; l < nlevels; l++ )
2616  {
2617  lchildEdgs.clear();
2618  rchildEdgs.clear();
2619  error = get_data_from_buff( 1, 0, l + 1, i, numedges, localEdgeList, lchildEdgs );MB_CHK_ERR( error );
2620  error = get_data_from_buff( 1, 0, l + 1, Redge_idx, numedges, remEdgeList, rchildEdgs );MB_CHK_ERR( error );
2621 
2622  int nchd = lchildEdgs.size();
2623  if( orient )
2624  {
2625  for( int j = 0; j < nchd; j++ )
2626  {
2627  // match entityhandles of child edges
2628  bool found = check_for_parallelinfo( lchildEdgs[j], shared_proc, remProcs );
2629  if( !found )
2630  {
2631  remProcs.insert( std::pair< EntityHandle, int >( lchildEdgs[j], shared_proc ) );
2632  remHandles.insert( std::pair< EntityHandle, EntityHandle >( lchildEdgs[j], rchildEdgs[j] ) );
2633  }
2634 
2635  // match entityhandles of child vertices
2636  lconn.clear();
2637  rconn.clear();
2638 
2639  int lidx =
2640  std::find( localEdgeList.begin(), localEdgeList.end(), lchildEdgs[j] ) - localEdgeList.begin();
2641  int ridx = std::find( remEdgeList.begin(), remEdgeList.end(), rchildEdgs[j] ) - remEdgeList.begin();
2642 
2643  error = get_data_from_buff( 1, 1, l + 1, lidx, numedges, localEdgeList, lconn );MB_CHK_ERR( error );
2644  error = get_data_from_buff( 1, 1, l + 1, ridx, numedges, remEdgeList, rconn );MB_CHK_ERR( error );
2645 
2646  found = check_for_parallelinfo( lconn[0], shared_proc, remProcs );
2647  if( !found )
2648  {
2649  remProcs.insert( std::pair< EntityHandle, int >( lconn[0], shared_proc ) );
2650  remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[0], rconn[0] ) );
2651  }
2652  found = check_for_parallelinfo( lconn[1], shared_proc, remProcs );
2653  if( !found )
2654  {
2655  remProcs.insert( std::pair< EntityHandle, int >( lconn[1], shared_proc ) );
2656  remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[1], rconn[0] ) );
2657  }
2658  }
2659  }
2660 
2661  else
2662  {
2663  for( int j = 0; j < nchd; j++ )
2664  {
2665  // match entityhandles of child edges
2666  bool found = check_for_parallelinfo( lchildEdgs[j], shared_proc, remProcs );
2667 
2668  if( !found )
2669  {
2670  remProcs.insert( std::pair< EntityHandle, int >( lchildEdgs[j], shared_proc ) );
2671  remHandles.insert(
2672  std::pair< EntityHandle, EntityHandle >( lchildEdgs[j], rchildEdgs[nchd - j - 1] ) );
2673  }
2674 
2675  // match entityhandles of child vertices
2676  lconn.clear();
2677  rconn.clear();
2678 
2679  int lidx =
2680  std::find( localEdgeList.begin(), localEdgeList.end(), lchildEdgs[j] ) - localEdgeList.begin();
2681  int ridx = std::find( remEdgeList.begin(), remEdgeList.end(), rchildEdgs[nchd - j - 1] ) -
2682  remEdgeList.begin();
2683 
2684  error = get_data_from_buff( 1, 1, l + 1, lidx, numedges, localEdgeList, lconn );MB_CHK_ERR( error );
2685  error = get_data_from_buff( 1, 1, l + 1, ridx, numedges, remEdgeList, rconn );MB_CHK_ERR( error );
2686  found = check_for_parallelinfo( lconn[0], shared_proc, remProcs );
2687 
2688  if( !found )
2689  {
2690  remProcs.insert( std::pair< EntityHandle, int >( lconn[0], shared_proc ) );
2691  remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[0], rconn[1] ) );
2692  }
2693 
2694  found = check_for_parallelinfo( lconn[1], shared_proc, remProcs );
2695  if( !found )
2696  {
2697  remProcs.insert( std::pair< EntityHandle, int >( lconn[1], shared_proc ) );
2698  remHandles.insert( std::pair< EntityHandle, EntityHandle >( lconn[1], rconn[1] ) );
2699  }
2700  }
2701  }
2702  }
2703  }
2704 
2705  return MB_SUCCESS;
2706 }
2707 
2708 ErrorCode NestedRefine::decipher_remote_handles_vertex( int shared_proc,
2709  int numverts,
2710  std::vector< EntityHandle >& localVertexList,
2711  std::vector< EntityHandle >& remVertexList,
2712  std::multimap< EntityHandle, int >& remProcs,
2713  std::multimap< EntityHandle, EntityHandle >& remHandles )
2714 {
2715  // LVList = <V> where V = <LV0, LV1, ..,LVL>
2716  // RVList = <V> where V = <LV0', RV1, ..,RVL>, LV0' is the local handles of coarsest vertices
2717  // but in the order in which they appear on the remote/shared proc
2718 
2719  ErrorCode error;
2720 
2721  for( int i = 0; i < numverts; i++ )
2722  {
2723  EntityHandle v = localVertexList[i];
2724  int Rvert_idx =
2725  ( std::find( remVertexList.begin(), remVertexList.begin() + numverts - 1, v ) ) - remVertexList.begin();
2726 
2727  std::vector< EntityHandle > lverts, rverts;
2728  for( int l = 0; l < nlevels; l++ )
2729  {
2730  error = get_data_from_buff( 0, 0, l + 1, i, numverts, localVertexList, lverts );MB_CHK_ERR( error );
2731  error = get_data_from_buff( 0, 0, l + 1, Rvert_idx, numverts, remVertexList, rverts );MB_CHK_ERR( error );
2732 
2733  bool found = check_for_parallelinfo( lverts[0], shared_proc, remProcs );
2734  if( !found )
2735  {
2736  remProcs.insert( std::pair< EntityHandle, int >( lverts[0], shared_proc ) );
2737  remHandles.insert( std::pair< EntityHandle, EntityHandle >( lverts[0], rverts[0] ) );
2738  }
2739  }
2740  }
2741 
2742  return MB_SUCCESS;
2743 }
2744 
2745 ErrorCode NestedRefine::update_parallel_tags( std::multimap< EntityHandle, int >& remProcs,
2746  std::multimap< EntityHandle, EntityHandle >& remHandles )
2747 {
2748  ErrorCode error;
2749 
2750  std::vector< int > rprocs;
2751  std::vector< EntityHandle > rhandles;
2752 
2753  std::multimap< EntityHandle, int >::iterator it;
2754  std::pair< std::multimap< EntityHandle, int >::iterator, std::multimap< EntityHandle, int >::iterator > it_procs;
2755  std::pair< std::multimap< EntityHandle, EntityHandle >::iterator,
2756  std::multimap< EntityHandle, EntityHandle >::iterator >
2757  it_handles;
2758 
2759  it = remProcs.begin();
2760  while( it != remProcs.end() )
2761  {
2762  rprocs.clear();
2763  rhandles.clear();
2764 
2765  EntityHandle entity = it->first;
2766  it_procs = remProcs.equal_range( entity );
2767  it_handles = remHandles.equal_range( entity );
2768 
2769  for( std::multimap< EntityHandle, int >::iterator pit = it_procs.first; pit != it_procs.second; pit++ )
2770  rprocs.push_back( pit->second );
2771  for( std::multimap< EntityHandle, EntityHandle >::iterator pit = it_handles.first; pit != it_handles.second;
2772  pit++ )
2773  rhandles.push_back( pit->second );
2774 
2775  error = pcomm->update_remote_data( entity, rprocs, rhandles );MB_CHK_ERR( error );
2776 
2777  it = remProcs.upper_bound( it->first );
2778  }
2779 
2780  return MB_SUCCESS;
2781 }
2782 
2783 ErrorCode NestedRefine::get_data_from_buff( int listtype,
2784  int datatype,
2785  int level,
2786  int entity_index,
2787  int nentities,
2788  std::vector< EntityHandle >& buffer,
2789  std::vector< EntityHandle >& data )
2790 {
2791  /**
2792  * listtype = 2, 1, 0 for FList (faces), EList (edge) and VList (vertices), respectively
2793  * datatype = 0(entityhandles of entities at level 0 and their children upto nlevels, in the
2794  *case of vertices its the duplicate vertices at subsequent levels), = 1(connectivity for the
2795  *above), = 2(subentities, only case is edges from FList) entity_index = integer index of the
2796  *entity in the buffer nentities = number of entities at the coarsest level i.e., |F0| or |E0|
2797  *or |V0|
2798  *
2799  * Two types of queries:
2800  * 1) Given an entity, find its children at a particular level.
2801  * 2) Given any entity at any level , find its connectivity (or bounding edges)
2802  *
2803  **/
2804 
2805  data.clear();
2806 
2807  if( listtype == 2 ) // FList
2808  {
2809  // FList = <F, FC, FE> where F = <F0, F1, ..,FL>, FC = <FC0, FC1, .., FCL> and FE = <FE0,
2810  // FE1, .., FEL>
2811  if( datatype == 0 ) // requesting child entities
2812  {
2813  EntityType ftype = mbImpl->type_from_handle( buffer[0] );
2814 
2815  int start, end, toadd, prev;
2816  start = end = entity_index;
2817  toadd = prev = nentities;
2818  for( int i = 0; i < level; i++ )
2819  {
2820  int d = get_index_from_degree( level_dsequence[i] );
2821  int nch = refTemplates[ftype - 1][d].total_new_ents;
2822  start = start * nch;
2823  end = end * nch + nch - 1;
2824  if( i < level - 1 )
2825  {
2826  prev = prev * nch;
2827  toadd += prev;
2828  }
2829  }
2830 
2831  start += toadd;
2832  end += toadd;
2833 
2834  int num_child = end - start + 1;
2835  data.reserve( num_child );
2836 
2837  for( int i = start; i <= end; i++ )
2838  {
2839  EntityHandle child = buffer[i];
2840  data.push_back( child );
2841  }
2842  }
2843  else if( datatype == 1 ) // requesting connectivity of an entity
2844  {
2845  EntityType ftype = mbImpl->type_from_handle( buffer[0] );
2846  int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
2847 
2848  int toadd = nentities, prev = nentities;
2849  for( int i = 0; i < nlevels; i++ )
2850  {
2851  int d = get_index_from_degree( level_dsequence[i] );
2852  int nch = refTemplates[ftype - 1][d].total_new_ents;
2853  prev = prev * nch;
2854  toadd += prev;
2855  }
2856 
2857  for( int i = 0; i < nepf; i++ )
2858  data.push_back( buffer[toadd + nepf * entity_index + i] );
2859  }
2860  else if( datatype == 2 ) // requesting bounding edges of an entity
2861  {
2862  EntityType ftype = mbImpl->type_from_handle( buffer[0] );
2863  int nepf = ahf->lConnMap2D[ftype - 2].num_verts_in_face;
2864 
2865  int toadd = nentities, prev = nentities;
2866  for( int i = 0; i < nlevels; i++ )
2867  {
2868  int d = get_index_from_degree( level_dsequence[i] );
2869  int nch = refTemplates[ftype - 1][d].total_new_ents;
2870  prev = prev * nch;
2871  toadd += prev;
2872  }
2873  toadd += toadd * nepf;
2874 
2875  for( int i = 0; i < nepf; i++ )
2876  data.push_back( buffer[toadd + nepf * entity_index + i] );
2877  }
2878  else
2879  MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer for faces" );
2880  }
2881  else if( listtype == 1 ) // EList
2882  {
2883  // EList = <E, EC> where E = <E0, E1, ..,EL> and EC = <EC0, EC1, .., ECL>
2884  if( datatype == 0 ) // requesting child entities
2885  {
2886  int start, end, toadd, prev;
2887  start = end = entity_index;
2888  toadd = prev = nentities;
2889  for( int i = 0; i < level; i++ )
2890  {
2891  int d = get_index_from_degree( level_dsequence[i] );
2892  int nch = refTemplates[MBEDGE - 1][d].total_new_ents;
2893  start = start * nch;
2894  end = end * nch + nch - 1;
2895 
2896  if( i < level - 1 )
2897  {
2898  prev = prev * nch;
2899  toadd += prev;
2900  }
2901  }
2902 
2903  start += toadd;
2904  end += toadd;
2905 
2906  int num_child = end - start + 1;
2907  data.reserve( num_child );
2908 
2909  for( int i = start; i <= end; i++ )
2910  {
2911  EntityHandle child = buffer[i];
2912  data.push_back( child );
2913  }
2914  }
2915  else if( datatype == 1 ) // requesting connectivities of child entities
2916  {
2917  int toadd = nentities, prev = nentities;
2918  for( int i = 0; i < nlevels; i++ )
2919  {
2920  int d = get_index_from_degree( level_dsequence[i] );
2921  int nch = refTemplates[MBEDGE - 1][d].total_new_ents;
2922  prev = prev * nch;
2923  toadd += prev;
2924  }
2925 
2926  data.push_back( buffer[toadd + 2 * entity_index] );
2927  data.push_back( buffer[toadd + 2 * entity_index + 1] );
2928  }
2929  else
2930  MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer for edges" );
2931  }
2932  else if( listtype == 0 ) // VList
2933  {
2934  // VList = <V> where V = <V0, V1, ..,VL>
2935  if( datatype == 0 )
2936  {
2937  int idx = level * nentities + entity_index;
2938  data.push_back( buffer[idx] );
2939  }
2940  else
2941  MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer for vertices" );
2942  }
2943  else
2944  MB_SET_ERR( MB_FAILURE, "Requesting invalid info from buffer" );
2945 
2946  return MB_SUCCESS;
2947 }
2948 
2949 bool NestedRefine::check_for_parallelinfo( EntityHandle entity, int proc, std::multimap< EntityHandle, int >& remProcs )
2950 {
2951  bool found = false;
2952 
2953  std::pair< std::multimap< EntityHandle, int >::iterator, std::multimap< EntityHandle, int >::iterator > it_hes;
2954  it_hes = remProcs.equal_range( entity );
2955 
2956  for( std::multimap< EntityHandle, int >::iterator it = it_hes.first; it != it_hes.second; ++it )
2957  {
2958  if( it->second == proc )
2959  {
2960  found = true;
2961  break;
2962  }
2963  }
2964 
2965  return found;
2966 }
2967 
2968 ErrorCode NestedRefine::check_for_parallelinfo( EntityHandle entity,
2969  int proc,
2970  std::multimap< EntityHandle, EntityHandle >& remHandles,
2971  std::multimap< EntityHandle, int >& remProcs,
2972  EntityHandle& rhandle )
2973 {
2974  // shared procs for given entity
2975  std::pair< std::multimap< EntityHandle, int >::iterator, std::multimap< EntityHandle, int >::iterator > it_ps;
2976  it_ps = remProcs.equal_range( entity );
2977 
2978  // shared remote handles for given entity
2979  std::pair< std::multimap< EntityHandle, EntityHandle >::iterator,
2980  std::multimap< EntityHandle, EntityHandle >::iterator >
2981  it_hs;
2982  it_hs = remHandles.equal_range( entity );
2983 
2984  std::multimap< EntityHandle, int >::iterator itp;
2985  std::multimap< EntityHandle, EntityHandle >::iterator ith;
2986 
2987  for( itp = it_ps.first, ith = it_hs.first; itp != it_ps.second; ++itp, ++ith )
2988  {
2989  if( itp->second == proc )
2990  {
2991  rhandle = ith->second;
2992  break;
2993  }
2994  }
2995 
2996  return MB_SUCCESS;
2997 }
2998 
2999 #endif
3000 
3001 /**********************************
3002  * Update AHF maps *
3003  * ********************************/
3004 
3006  EntityType type,
3007  int pat_id,
3008  EntityHandle* vbuffer,
3009  EntityHandle* ent_buffer,
3010  int etotal )
3011 {
3012  ErrorCode error;
3013  int nhf = 0, nv = 0, total_new_verts = 0;
3014  int d = get_index_from_degree( deg );
3015 
3016  // Get the number of half-facets
3017  if( type == MBEDGE )
3018  {
3019  nhf = 2;
3020  nv = 2;
3021  total_new_verts = refTemplates[0][d].total_new_verts;
3022  }
3023  else if( type == MBTRI || type == MBQUAD )
3024  {
3025  nhf = ahf->lConnMap2D[type - 2].num_verts_in_face;
3026  nv = nhf;
3027  total_new_verts = refTemplates[pat_id][d].total_new_verts;
3028  }
3029  else if( type == MBTET || type == MBHEX )
3030  {
3031  int index = ahf->get_index_in_lmap( *_incells.begin() );
3032  nhf = ahf->lConnMap3D[index].num_faces_in_cell;
3033  nv = ahf->lConnMap3D[index].num_verts_in_cell;
3034  total_new_verts = refTemplates[pat_id][d].total_new_verts;
3035  }
3036 
3037  std::vector< EntityHandle > ent;
3038  std::vector< int > lid;
3039 
3040  // Update the vertex to half-facet map
3041  for( int i = 0; i < total_new_verts; i++ )
3042  {
3043  ent.clear();
3044  lid.clear();
3045  EntityHandle vid = vbuffer[i + nv];
3046  error = ahf->get_incident_map( type, vid, ent, lid );MB_CHK_ERR( error );
3047 
3048  if( ent[0] ) continue;
3049 
3050  int id = refTemplates[pat_id][d].v2hf[i + nv][0] - 1;
3051  ent[0] = ent_buffer[id];
3052  lid[0] = refTemplates[pat_id][d].v2hf[i + nv][1];
3053 
3054  error = ahf->set_incident_map( type, vid, ent, lid );MB_CHK_ERR( error );
3055  }
3056 
3057  // Update the sibling half-facet map
3058  for( int i = 0; i < etotal; i++ )
3059  {
3060  std::vector< EntityHandle > sib_entids( nhf );
3061  std::vector< int > sib_lids( nhf );
3062 
3063  error = ahf->get_sibling_map( type, ent_buffer[i], &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3064 
3065  for( int l = 0; l < nhf; l++ )
3066  {
3067  if( sib_entids[l] ) continue;
3068 
3069  // Fill out the sibling values
3070  int id = refTemplates[pat_id][d].ents_opphfs[i][2 * l];
3071  if( id )
3072  {
3073  sib_entids[l] = ent_buffer[id - 1];
3074  sib_lids[l] = refTemplates[pat_id][d].ents_opphfs[i][2 * l + 1];
3075  }
3076  else
3077  {
3078  sib_entids[l] = 0;
3079  sib_lids[l] = 0;
3080  }
3081  }
3082 
3083  error = ahf->set_sibling_map( type, ent_buffer[i], &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3084 
3085  for( int l = 0; l < nhf; l++ )
3086  {
3087  if( sib_entids[l] )
3088  {
3089 
3090  EntityHandle set_entid = ent_buffer[i];
3091  int set_lid = l;
3092 
3093  error = ahf->set_sibling_map( type, sib_entids[l], sib_lids[l], set_entid, set_lid );MB_CHK_ERR( error );
3094  }
3095  }
3096  }
3097  return MB_SUCCESS;
3098 }
3099 
3101  EntityType type,
3102  EntityHandle* vbuffer,
3103  EntityHandle* ent_buffer,
3104  int etotal )
3105 {
3106  ErrorCode error;
3107  assert( type != MBTET );
3108  error = update_local_ahf( deg, type, type - 1, vbuffer, ent_buffer, etotal );MB_CHK_ERR( error );
3109 
3110  return MB_SUCCESS;
3111 }
3112 
3113 ErrorCode NestedRefine::update_global_ahf( EntityType type, int cur_level, int deg, std::vector< int >* pattern_ids )
3114 {
3115 
3116  ErrorCode error;
3117 
3118  // Get the number of half-facets and number of children of each type
3119  if( type == MBEDGE )
3120  {
3121  assert( pattern_ids == NULL );
3122  error = update_global_ahf_1D( cur_level, deg );MB_CHK_ERR( error );
3123  }
3124  else if( type == MBTRI || type == MBQUAD )
3125  {
3126  assert( pattern_ids == NULL );
3127  error = update_global_ahf_2D( cur_level, deg );MB_CHK_ERR( error );
3128  }
3129  else if( type == MBHEX )
3130  {
3131  assert( pattern_ids == NULL );
3132  error = update_global_ahf_3D( cur_level, deg );MB_CHK_ERR( error );
3133  }
3134  else if( type == MBTET )
3135  {
3136  assert( pattern_ids != NULL );
3137  error = update_global_ahf_3D( cur_level, deg, pattern_ids );MB_CHK_ERR( error );
3138  }
3139  else
3140  MB_SET_ERR( MB_NOT_IMPLEMENTED, "Requesting AHF update for an unsupported mesh entity type" );
3141 
3142  return MB_SUCCESS;
3143 }
3144 
3145 /*ErrorCode NestedRefine::update_global_ahf(int cur_level, int deg, std::vector<int> &pattern_ids)
3146 {
3147  ErrorCode error;
3148  error = update_global_ahf_3D(cur_level, deg, pattern_ids); MB_CHK_ERR(error);
3149 
3150  return MB_SUCCESS;
3151 }*/
3152 
3154 {
3155  ErrorCode error;
3156  int d = get_index_from_degree( deg );
3157  int nhf, nchilds, nverts_prev, nents_prev;
3158  nhf = 2;
3159  nchilds = refTemplates[0][d].total_new_ents;
3160  if( cur_level )
3161  {
3162  nverts_prev = level_mesh[cur_level - 1].num_verts;
3163  nents_prev = level_mesh[cur_level - 1].num_edges;
3164  }
3165  else
3166  {
3167  nverts_prev = _inverts.size();
3168  nents_prev = _inedges.size();
3169  }
3170 
3171  std::vector< EntityHandle > inci_ent, child_ents;
3172  std::vector< int > inci_lid, child_lids;
3173 
3174  // Update the vertex to half-facet maps for duplicate vertices
3175  for( int i = 0; i < nverts_prev; i++ )
3176  {
3177  inci_ent.clear();
3178  inci_lid.clear();
3179  child_ents.clear();
3180  child_lids.clear();
3181 
3182  // Vertex id in the previous mesh and the current one
3183  EntityHandle vid;
3184  if( cur_level )
3185  vid = level_mesh[cur_level - 1].start_vertex + i;
3186  else
3187  vid = _inverts[i];
3188  EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
3189 
3190  // Get the incident half-vert in the previous mesh
3191  error = ahf->get_incident_map( MBEDGE, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3192 
3193  // Obtain the corresponding incident child in the current mesh
3194  int lvid = get_local_vid( vid, inci_ent[0], cur_level - 1 );
3195  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3196  int chid = refTemplates[0][d].v2hf[lvid][0] - 1;
3197 
3198  int pid;
3199  if( cur_level )
3200  pid = inci_ent[0] - level_mesh[cur_level - 1].start_edge;
3201  else
3202  pid = inci_ent[0] - *_inedges.begin();
3203 
3204  int ind = nchilds * pid;
3205 
3206  child_ents.push_back( level_mesh[cur_level].start_edge + ind + chid );
3207  child_lids.push_back( refTemplates[0][d].v2hf[lvid][1] );
3208 
3209  error = ahf->set_incident_map( MBEDGE, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3210  }
3211 
3212  // Update the sibling half-facet maps across entities
3213  for( int i = 0; i < nents_prev; i++ )
3214  {
3215  EntityHandle ent;
3216  if( cur_level )
3217  ent = level_mesh[cur_level - 1].start_edge + i;
3218  else
3219  ent = _inedges[i];
3220 
3221  std::vector< EntityHandle > sib_entids( nhf );
3222  std::vector< int > sib_lids( nhf );
3223 
3224  error = ahf->get_sibling_map( MBEDGE, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3225 
3226  int id, idx;
3227 
3228  for( int l = 0; l < nhf; l++ )
3229  {
3230  if( !sib_entids[l] ) continue;
3231 
3232  // Find the child incident on the half-facet
3233  id = refTemplates[0][d].ents_on_pent[l][1] - 1;
3234  idx = nchilds * i;
3235  EntityHandle child_ent = level_mesh[cur_level].start_edge + idx + id;
3236  int ch_lid = l;
3237 
3238  // Find the sibling of the child
3239  std::vector< EntityHandle > sib_childs( nhf );
3240  std::vector< int > sib_chlids( nhf );
3241 
3242  error = ahf->get_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
3243 
3244  // If the sibling already exists, dont do anything
3245  if( sib_childs[ch_lid] ) continue;
3246 
3247  // Get the correponding child of the sibling of the current parent
3248  int psib;
3249  if( cur_level )
3250  psib = sib_entids[l] - level_mesh[cur_level - 1].start_edge;
3251  else
3252  psib = sib_entids[l] - *_inedges.begin();
3253 
3254  int plid = sib_lids[l];
3255 
3256  id = refTemplates[0][d].ents_on_pent[plid][1] - 1;
3257  idx = nchilds * psib;
3258 
3259  EntityHandle psib_child = level_mesh[cur_level].start_edge + idx + id;
3260  int psib_chlid = plid;
3261 
3262  // Set the siblings
3263  sib_childs[ch_lid] = psib_child;
3264  sib_chlids[ch_lid] = psib_chlid;
3265 
3266  error = ahf->set_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
3267  }
3268  }
3269 
3270  return MB_SUCCESS;
3271 }
3272 
3274 {
3275  ErrorCode error;
3276  int d = get_index_from_degree( deg );
3277  int nhf, nchilds, nents_prev;
3278  nhf = 2;
3279  nchilds = refTemplates[0][d].total_new_ents;
3280  if( cur_level )
3281  {
3282  nents_prev = level_mesh[cur_level - 1].num_edges;
3283  }
3284  else
3285  {
3286  nents_prev = _inedges.size();
3287  }
3288 
3289  // Update the sibling half-facet maps across entities
3290 
3291  std::vector< EntityHandle > conn;
3292  for( int i = 0; i < nents_prev; i++ )
3293  {
3294  EntityHandle ent;
3295  if( cur_level )
3296  ent = level_mesh[cur_level - 1].start_edge + i;
3297  else
3298  ent = _inedges[i];
3299 
3300  // Set incident hv maps
3301  conn.clear();
3302  error = get_connectivity( ent, cur_level, conn );MB_CHK_ERR( error );
3303 
3304  std::vector< EntityHandle > inci_ent, child_ents;
3305  std::vector< int > inci_lid, child_lids;
3306  for( int j = 0; j < 2; j++ )
3307  {
3308  inci_ent.clear();
3309  inci_lid.clear();
3310  child_ents.clear();
3311  child_lids.clear();
3312 
3313  // Get the entityhandle of the vertex from previous level in the current level
3314  EntityHandle cur_vid;
3315  if( cur_level )
3316  cur_vid = level_mesh[cur_level].start_vertex + ( conn[j] - level_mesh[cur_level - 1].start_vertex );
3317  else
3318  cur_vid = level_mesh[cur_level].start_vertex + ( conn[j] - *_inverts.begin() );
3319 
3320  // Obtain the incident half-facet. If exists, then no need to assign another
3321  error = ahf->get_incident_map( MBEDGE, cur_vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3322  if( inci_ent[0] != 0 ) continue;
3323 
3324  // Get the incident half-facet on the old vertex
3325  error = ahf->get_incident_map( MBEDGE, conn[j], inci_ent, inci_lid );MB_CHK_ERR( error );
3326 
3327  // Obtain the corresponding incident child in the current mesh
3328  int lvid = get_local_vid( conn[j], inci_ent[0], cur_level - 1 );
3329  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3330  int chid = refTemplates[0][d].v2hf[lvid][0] - 1;
3331 
3332  int pid;
3333  if( cur_level )
3334  pid = inci_ent[0] - level_mesh[cur_level - 1].start_edge;
3335  else
3336  pid = inci_ent[0] - *_inedges.begin();
3337 
3338  int ind = nchilds * pid;
3339 
3340  child_ents.push_back( level_mesh[cur_level].start_edge + ind + chid );
3341  child_lids.push_back( refTemplates[0][d].v2hf[lvid][1] );
3342 
3343  error = ahf->set_incident_map( MBEDGE, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3344  }
3345 
3346  std::vector< EntityHandle > sib_entids( nhf );
3347  std::vector< int > sib_lids( nhf );
3348 
3349  error = ahf->get_sibling_map( MBEDGE, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3350 
3351  int id, idx;
3352 
3353  for( int l = 0; l < nhf; l++ )
3354  {
3355  if( !sib_entids[l] ) continue;
3356 
3357  // Find the child incident on the half-facet
3358  id = refTemplates[0][d].ents_on_pent[l][1] - 1;
3359  idx = nchilds * i;
3360  EntityHandle child_ent = level_mesh[cur_level].start_edge + idx + id;
3361  int ch_lid = l;
3362 
3363  // Find the sibling of the child
3364  std::vector< EntityHandle > sib_childs( nhf );
3365  std::vector< int > sib_chlids( nhf );
3366 
3367  error = ahf->get_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
3368 
3369  // If the sibling already exists, dont do anything
3370  if( sib_childs[ch_lid] ) continue;
3371 
3372  // Get the correponding child of the sibling of the current parent
3373  int psib;
3374  if( cur_level )
3375  psib = sib_entids[l] - level_mesh[cur_level - 1].start_edge;
3376  else
3377  psib = sib_entids[l] - *_inedges.begin();
3378 
3379  int plid = sib_lids[l];
3380 
3381  id = refTemplates[0][d].ents_on_pent[plid][1] - 1;
3382  idx = nchilds * psib;
3383 
3384  EntityHandle psib_child = level_mesh[cur_level].start_edge + idx + id;
3385  int psib_chlid = plid;
3386 
3387  // Set the siblings
3388  sib_childs[ch_lid] = psib_child;
3389  sib_chlids[ch_lid] = psib_chlid;
3390 
3391  error = ahf->set_sibling_map( MBEDGE, child_ent, &sib_childs[0], &sib_chlids[0], nhf );MB_CHK_ERR( error );
3392  }
3393  }
3394 
3395  return MB_SUCCESS;
3396 }
3397 
3399 {
3400  ErrorCode error;
3401  error = ahf->determine_sibling_halfverts( level_mesh[cur_level].verts, level_mesh[cur_level].edges );MB_CHK_ERR( error );
3402 
3404 
3405  return MB_SUCCESS;
3406 }
3407 
3409 {
3410  ErrorCode error;
3411 
3412  EntityType type = mbImpl->type_from_handle( *_infaces.begin() );
3413  int nhf, nchilds, nverts_prev, nents_prev;
3414 
3415  nhf = ahf->lConnMap2D[type - 2].num_verts_in_face;
3416  int d = get_index_from_degree( deg );
3417  nchilds = refTemplates[type - 1][d].total_new_ents;
3418 
3419  if( cur_level )
3420  {
3421  nverts_prev = level_mesh[cur_level - 1].num_verts;
3422  nents_prev = level_mesh[cur_level - 1].num_faces;
3423  }
3424  else
3425  {
3426  nverts_prev = _inverts.size();
3427  nents_prev = _infaces.size();
3428  }
3429 
3430  std::vector< EntityHandle > inci_ent, child_ents;
3431  std::vector< int > inci_lid, child_lids;
3432 
3433  // Update the vertex to half-edge maps for old/duplicate vertices
3434  for( int i = 0; i < nverts_prev; i++ )
3435  {
3436  inci_ent.clear();
3437  inci_lid.clear();
3438  child_ents.clear();
3439  child_lids.clear();
3440 
3441  // Vertex id in the previous mesh
3442  EntityHandle vid;
3443  if( cur_level )
3444  vid = level_mesh[cur_level - 1].start_vertex + i;
3445  else
3446  vid = _inverts[i];
3447  EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
3448 
3449  // Get the incident half-vert in the previous mesh
3450  error = ahf->get_incident_map( type, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3451 
3452  // Obtain the corresponding incident child in the current mesh
3453  for( int j = 0; j < (int)inci_ent.size(); j++ )
3454  {
3455  int lvid = get_local_vid( vid, inci_ent[j], cur_level - 1 );
3456  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3457  int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
3458 
3459  int pid;
3460  if( cur_level )
3461  pid = inci_ent[j] - level_mesh[cur_level - 1].start_face;
3462  else
3463  pid = inci_ent[j] - *_infaces.begin();
3464 
3465  int ind = nchilds * pid;
3466 
3467  child_ents.push_back( level_mesh[cur_level].start_face + ind + chid );
3468  child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
3469  }
3470  error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3471  }
3472 
3473  EntityHandle fedge[2];
3474 
3475  // Update the sibling half-facet maps across entities
3476  for( int i = 0; i < nents_prev; i++ )
3477  {
3478  EntityHandle ent;
3479  if( cur_level )
3480  ent = level_mesh[cur_level - 1].start_face + i;
3481  else
3482  ent = _infaces[i];
3483 
3484  std::vector< EntityHandle > fid_conn;
3485  error = get_connectivity( ent, cur_level, fid_conn );
3486  if( MB_SUCCESS != error ) return error;
3487 
3488  std::vector< EntityHandle > sib_entids( nhf );
3489  std::vector< int > sib_lids( nhf );
3490 
3491  error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3492 
3493  int id, idx;
3494 
3495  for( int l = 0; l < nhf; l++ )
3496  {
3497  if( !sib_entids[l] ) continue;
3498 
3499  int nidx = ahf->lConnMap2D[type - 2].next[l];
3500  fedge[0] = fid_conn[l];
3501  fedge[1] = fid_conn[nidx];
3502 
3503  EntityHandle sfid = sib_entids[l];
3504  int slid = sib_lids[l];
3505 
3506  std::vector< EntityHandle > conn;
3507  error = get_connectivity( sfid, cur_level, conn );
3508  if( MB_SUCCESS != error ) return error;
3509 
3510  bool orient = true;
3511  nidx = ahf->lConnMap2D[type - 2].next[slid];
3512  if( ( fedge[1] == conn[slid] ) && ( fedge[0] == conn[nidx] ) ) orient = false;
3513 
3514  if( orient ) assert( ( fedge[0] == conn[slid] ) && ( fedge[1] == conn[nidx] ) );
3515 
3516  // Find the childrens incident on the half-facet
3517  int nch = refTemplates[type - 1][d].ents_on_pent[l][0];
3518  idx = nchilds * i;
3519 
3520  // Loop over all the incident childrens
3521  for( int k = 0; k < nch; k++ )
3522  {
3523  id = refTemplates[type - 1][d].ents_on_pent[l][k + 1] - 1;
3524  EntityHandle child_ent = level_mesh[cur_level].start_face + idx + id;
3525  int child_lid = l;
3526 
3527  // Find the sibling of the child
3528  EntityHandle child_sibent;
3529  int child_siblid;
3530  error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
3531 
3532  if( child_sibent != 0 ) continue;
3533 
3534  // Get the correponding child of the sibling of the current parent
3535  int psib;
3536  if( cur_level )
3537  psib = sfid - level_mesh[cur_level - 1].start_face;
3538  else
3539  psib = sfid - *_infaces.begin();
3540 
3541  int plid = slid;
3542 
3543  if( orient )
3544  id = refTemplates[type - 1][d].ents_on_pent[plid][k + 1] - 1;
3545  else
3546  id = refTemplates[type - 1][d].ents_on_pent[plid][nch - k] - 1;
3547 
3548  int sidx = nchilds * psib;
3549 
3550  EntityHandle psib_child = level_mesh[cur_level].start_face + sidx + id;
3551  int psib_chlid = plid;
3552 
3553  // Set the siblings
3554  error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
3555  }
3556  }
3557  }
3558 
3559  return MB_SUCCESS;
3560 }
3561 
3563 {
3564  ErrorCode error;
3565  int d = get_index_from_degree( deg );
3566  EntityType type = mbImpl->type_from_handle( *_infaces.begin() );
3567  int nhf, nchilds, nents_prev;
3568  nhf = ahf->lConnMap2D[type - 2].num_verts_in_face;
3569  nchilds = refTemplates[type - 1][d].total_new_ents;
3570 
3571  if( cur_level )
3572  nents_prev = level_mesh[cur_level - 1].num_faces;
3573  else
3574  nents_prev = _infaces.size();
3575 
3576  EntityHandle fedge[2];
3577 
3578  // Update the sibling half-facet maps across entities
3579  for( int i = 0; i < nents_prev; i++ )
3580  {
3581  EntityHandle ent;
3582  if( cur_level )
3583  ent = level_mesh[cur_level - 1].start_face + i;
3584  else
3585  ent = _infaces[i];
3586 
3587  std::vector< EntityHandle > fid_conn;
3588  error = get_connectivity( ent, cur_level, fid_conn );
3589  if( MB_SUCCESS != error ) return error;
3590 
3591  std::vector< EntityHandle > inci_ent, child_ents;
3592  std::vector< int > inci_lid, child_lids;
3593 
3594  // Set incident half-edges
3595  for( int j = 0; j < nhf; j++ )
3596  {
3597  inci_ent.clear();
3598  inci_lid.clear();
3599  child_ents.clear();
3600  child_lids.clear();
3601  EntityHandle cur_vid;
3602  if( cur_level )
3603  cur_vid = level_mesh[cur_level].start_vertex + ( fid_conn[j] - level_mesh[cur_level - 1].start_vertex );
3604  else
3605  cur_vid = level_mesh[cur_level].start_vertex + ( fid_conn[j] - *_inverts.begin() );
3606 
3607  // Obtain the incident half-facet. If exists, then no need to assign another
3608  error = ahf->get_incident_map( type, cur_vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3609  if( inci_ent[0] != 0 ) continue;
3610 
3611  // Get the incident half-facet on the old vertex
3612  error = ahf->get_incident_map( type, fid_conn[j], inci_ent, inci_lid );MB_CHK_ERR( error );
3613 
3614  // Obtain the corresponding incident child in the current mesh
3615  for( int k = 0; k < (int)inci_ent.size(); k++ )
3616  {
3617  int lvid = get_local_vid( fid_conn[j], inci_ent[k], cur_level - 1 );
3618  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3619  int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
3620 
3621  int pid;
3622  if( cur_level )
3623  pid = inci_ent[k] - level_mesh[cur_level - 1].start_face;
3624  else
3625  pid = inci_ent[k] - *_infaces.begin();
3626 
3627  int ind = nchilds * pid;
3628 
3629  child_ents.push_back( level_mesh[cur_level].start_face + ind + chid );
3630  child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
3631  }
3632 
3633  error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3634  }
3635 
3636  // Set sibling half-edges
3637  std::vector< EntityHandle > sib_entids( nhf );
3638  std::vector< int > sib_lids( nhf );
3639 
3640  error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3641 
3642  int id, idx;
3643 
3644  for( int l = 0; l < nhf; l++ )
3645  {
3646  if( !sib_entids[l] ) continue;
3647 
3648  int nidx = ahf->lConnMap2D[type - 2].next[l];
3649  fedge[0] = fid_conn[l];
3650  fedge[1] = fid_conn[nidx];
3651 
3652  EntityHandle sfid = sib_entids[l];
3653  int slid = sib_lids[l];
3654 
3655  std::vector< EntityHandle > conn;
3656  error = get_connectivity( sfid, cur_level, conn );MB_CHK_ERR( error );
3657 
3658  assert( (int)conn.size() > nidx && (int)conn.size() > slid );
3659 
3660  bool orient = true;
3661  nidx = ahf->lConnMap2D[type - 2].next[slid];
3662  if( ( fedge[1] == conn[slid] ) && ( fedge[0] == conn[nidx] ) ) orient = false;
3663 
3664  if( orient ) assert( ( fedge[0] == conn[slid] ) && ( fedge[1] == conn[nidx] ) );
3665 
3666  // Find the childrens incident on the half-facet
3667  int nch = refTemplates[type - 1][d].ents_on_pent[l][0];
3668  idx = nchilds * i;
3669 
3670  // Loop over all the incident childrens
3671  for( int k = 0; k < nch; k++ )
3672  {
3673  id = refTemplates[type - 1][d].ents_on_pent[l][k + 1] - 1;
3674  EntityHandle child_ent = level_mesh[cur_level].start_face + idx + id;
3675  int child_lid = l;
3676 
3677  // Find the sibling of the child
3678  EntityHandle child_sibent;
3679  int child_siblid;
3680  error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
3681 
3682  if( child_sibent != 0 ) continue;
3683 
3684  // Get the correponding child of the sibling of the current parent
3685  int psib;
3686  if( cur_level )
3687  psib = sfid - level_mesh[cur_level - 1].start_face;
3688  else
3689  psib = sfid - *_infaces.begin();
3690 
3691  int plid = slid;
3692 
3693  if( orient )
3694  id = refTemplates[type - 1][d].ents_on_pent[plid][k + 1] - 1;
3695  else
3696  id = refTemplates[type - 1][d].ents_on_pent[plid][nch - k] - 1;
3697 
3698  int sidx = nchilds * psib;
3699 
3700  EntityHandle psib_child = level_mesh[cur_level].start_face + sidx + id;
3701  int psib_chlid = plid;
3702 
3703  // Set the siblings
3704  error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
3705  }
3706  }
3707  }
3708 
3709  return MB_SUCCESS;
3710 }
3711 
3712 ErrorCode NestedRefine::update_global_ahf_3D( int cur_level, int deg, std::vector< int >* pattern_ids )
3713 {
3714  ErrorCode error;
3715  int nvpc, ne, nhf, nchilds, nverts_prev, nents_prev;
3716 
3717  EntityType type = mbImpl->type_from_handle( *_incells.begin() );
3718  int index = ahf->get_index_in_lmap( *_incells.begin() );
3719  int d = get_index_from_degree( deg );
3720 
3721  nhf = ahf->lConnMap3D[index].num_faces_in_cell;
3722  ne = ahf->lConnMap3D[index].num_edges_in_cell;
3723  nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
3724  nchilds = refTemplates[type - 1][d].total_new_ents;
3725 
3726  if( cur_level )
3727  {
3728  nverts_prev = level_mesh[cur_level - 1].num_verts;
3729  nents_prev = level_mesh[cur_level - 1].num_cells;
3730  }
3731  else
3732  {
3733  nverts_prev = _inverts.size();
3734  nents_prev = _incells.size();
3735  }
3736 
3737  std::vector< EntityHandle > inci_ent, child_ents;
3738  std::vector< int > inci_lid, child_lids;
3739 
3740  // Step 1: Update the V2HF maps for old/duplicate vertices
3741  for( int i = 0; i < nverts_prev; i++ )
3742  {
3743  inci_ent.clear();
3744  inci_lid.clear();
3745  child_ents.clear();
3746  child_lids.clear();
3747 
3748  // Vertex id in the previous mesh
3749  EntityHandle vid;
3750  if( cur_level )
3751  vid = level_mesh[cur_level - 1].start_vertex + i;
3752  else
3753  vid = _inverts[i];
3754  EntityHandle cur_vid = level_mesh[cur_level].start_vertex + i;
3755 
3756  // Get the incident half-vert in the previous mesh
3757  error = ahf->get_incident_map( type, vid, inci_ent, inci_lid );MB_CHK_ERR( error );
3758 
3759  // Obtain the corresponding incident child in the current mesh
3760  for( int j = 0; j < (int)inci_ent.size(); j++ )
3761  {
3762  int lvid = get_local_vid( vid, inci_ent[j], cur_level - 1 );
3763  if( lvid < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex ix " );
3764  int chid = refTemplates[type - 1][d].v2hf[lvid][0] - 1;
3765 
3766  int pid;
3767  if( cur_level )
3768  pid = inci_ent[j] - level_mesh[cur_level - 1].start_cell;
3769  else
3770  pid = inci_ent[j] - *_incells.begin();
3771 
3772  int ind = nchilds * pid;
3773 
3774  // EntityHandle child_ent = level_mesh[cur_level].start_cell + ind+chid ;
3775  // int child_lid = refTemplates[type-1][d].v2hf[lvid][1];
3776  child_ents.push_back( level_mesh[cur_level].start_cell + ind + chid );
3777  child_lids.push_back( refTemplates[type - 1][d].v2hf[lvid][1] );
3778  }
3779 
3780  error = ahf->set_incident_map( type, cur_vid, child_ents, child_lids );MB_CHK_ERR( error );
3781  }
3782 
3783  // error = ahf->determine_incident_halffaces( level_mesh[cur_level].cells);MB_CHK_ERR(error);
3784 
3785  // Step 2: Update SIBHFS maps
3786  for( int i = 0; i < nents_prev; i++ )
3787  {
3788  EntityHandle ent;
3789  if( cur_level )
3790  ent = level_mesh[cur_level - 1].start_cell + i;
3791  else
3792  ent = _incells[i];
3793 
3794  std::vector< EntityHandle > sib_entids( nhf );
3795  std::vector< int > sib_lids( nhf );
3796 
3797  error = ahf->get_sibling_map( type, ent, &sib_entids[0], &sib_lids[0], nhf );MB_CHK_ERR( error );
3798 
3799  int id, idx;
3800 
3801  for( int l = 0; l < nhf; l++ )
3802  {
3803 
3804  if( !sib_entids[l] ) continue;
3805 
3806  // Get the number of children incident on this half-face
3807  int pat_id;
3808  if( type == MBTET )
3809  pat_id = ( *pattern_ids )[i];
3810  else
3811  pat_id = type - 1;
3812  int nch = refTemplates[pat_id][d].ents_on_pent[l][0];
3813 
3814  // Get the order of children indices incident on this half-face
3815  std::vector< int > id_sib( nch );
3816  for( int k = 0; k < nch; k++ )
3817  id_sib[k] = 0;
3818 
3819  error = reorder_indices( cur_level, deg, ent, l, sib_entids[l], sib_lids[l], 1, &id_sib[0] );MB_CHK_ERR( error );
3820 
3821  // Get the parent index of the sibling cell
3822  int psib;
3823  if( cur_level )
3824  psib = sib_entids[l] - level_mesh[cur_level - 1].start_cell;
3825  else
3826  psib = sib_entids[l] - *_incells.begin();
3827 
3828  int plid = sib_lids[l];
3829  int sidx = nchilds * psib;
3830  int sibpat_id;
3831  if( type == MBTET )
3832  sibpat_id = ( *pattern_ids )[psib];
3833  else
3834  sibpat_id = type - 1;
3835 
3836  // Loop over all the childs incident on the working half-face
3837  idx = nchilds * i;
3838 
3839  for( int k = 0; k < nch; k++ )
3840  {
3841  id = refTemplates[pat_id][d].ents_on_pent[l][k + 1] - 1;
3842  EntityHandle child_ent = level_mesh[cur_level].start_cell + idx + id;
3843  int child_lid = l;
3844 
3845  // Find the sibling of the working child
3846  EntityHandle child_sibent;
3847  int child_siblid;
3848  error = ahf->get_sibling_map( type, child_ent, child_lid, child_sibent, child_siblid );MB_CHK_ERR( error );
3849 
3850  if( child_sibent != 0 ) continue;
3851 
3852  // Get the correponding child of the sibling of the current parent
3853  // We have already computed the order the children on incident corresponding to the
3854  // working half-face
3855  id = refTemplates[sibpat_id][d].ents_on_pent[plid][id_sib[k]] - 1;
3856 
3857  EntityHandle psib_child = level_mesh[cur_level].start_cell + sidx + id;
3858  int psib_chlid = plid;
3859 
3860  // Set the siblings of children incident on current half-face
3861  error = ahf->set_sibling_map( type, child_ent, child_lid, psib_child, psib_chlid );MB_CHK_ERR( error );
3862 
3863  // Set the sibling of the sibling of the children to the children
3864  error = ahf->set_sibling_map( type, psib_child, psib_chlid, child_ent, child_lid );MB_CHK_ERR( error );
3865  }
3866  }
3867 
3868  // Loop over edges to check if there are any non-manifold edges. If there are then the v2hfs
3869  // map should be updated for the new vertices on it.
3870  const EntityHandle* conn;
3871  error = mbImpl->get_connectivity( ent, conn, nvpc );MB_CHK_ERR( error );
3872 
3873  int nv = refTemplates[type - 1][d].nv_edge; //#verts on each edge
3874  for( int l = 0; l < ne; l++ )
3875  {
3876  id = ahf->lConnMap3D[index].e2v[l][0];
3877  EntityHandle v_start = conn[id];
3878  id = ahf->lConnMap3D[index].e2v[l][1];
3879  EntityHandle v_end = conn[id];
3880 
3881  bool visited = false;
3882 
3883  std::vector< EntityHandle > inci_ent1, inci_ent2;
3884  std::vector< int > inci_lid1, inci_lid2;
3885  error = ahf->get_incident_map( type, v_start, inci_ent1, inci_lid1 );MB_CHK_ERR( error );
3886  error = ahf->get_incident_map( type, v_end, inci_ent2, inci_lid2 );MB_CHK_ERR( error );
3887 
3888  if( inci_ent1.size() > 1 && inci_ent2.size() > 1 )
3889  {
3890  std::vector< EntityHandle > cell_comps;
3891  std::vector< int > leid_comps;
3892 
3893  error = ahf->get_up_adjacencies_edg_3d_comp( ent, l, cell_comps, &leid_comps );MB_CHK_ERR( error );
3894 
3895  int ncomps = cell_comps.size();
3896  std::vector< EntityHandle > edgverts;
3897  std::vector< EntityHandle > compchildents( nv * ncomps );
3898  std::vector< int > compchildlfids( nv * ncomps );
3899 
3900  for( int s = 0; s < nv * ncomps; s++ )
3901  {
3902  compchildents[s] = 0;
3903  compchildlfids[s] = 0;
3904  }
3905 
3906  for( int j = 0; j < (int)cell_comps.size(); j++ )
3907  {
3908  int ind;
3909  if( cur_level )
3910  ind = level_mesh[cur_level - 1].cells.index( cell_comps[j] );
3911  else
3912  ind = _incells.index( cell_comps[j] );
3913 
3914  for( int k = 0; k < nv; k++ )
3915  {
3916  int chid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k] - 1;
3917  int lfid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k + 1];
3918  int lvid = refTemplates[type - 1][d].ents_on_vedge[leid_comps[j]][3 * k + 2];
3919 
3920  EntityHandle childcell = level_mesh[cur_level].start_cell + ind * nchilds + chid;
3921 
3922  const EntityHandle* econn;
3923  error = mbImpl->get_connectivity( childcell, econn, nvpc );MB_CHK_ERR( error );
3924 
3925  EntityHandle vert = econn[lvid];
3926 
3927  if( ahf->check_nonmanifold_vertices( type, vert ) )
3928  {
3929  visited = true;
3930  break;
3931  }
3932 
3933  if( edgverts.empty() )
3934  {
3935  edgverts.push_back( vert );
3936  compchildents[0] = childcell;
3937  compchildlfids[0] = lfid;
3938  }
3939  else
3940  {
3941  std::vector< EntityHandle >::iterator it;
3942  it = find( edgverts.begin(), edgverts.end(), vert );
3943  int indx = it - edgverts.begin();
3944 
3945  if( it == edgverts.end() )
3946  {
3947  edgverts.push_back( vert );
3948  compchildents[k * ncomps] = childcell;
3949  compchildlfids[k * ncomps] = lfid;
3950  }
3951  else
3952  {
3953  compchildents[indx * ncomps + j] = childcell;
3954  compchildlfids[indx * ncomps + j] = lfid;
3955  }
3956  }
3957  }
3958  }
3959 
3960  if( visited )
3961  {
3962  break;
3963  }
3964 
3965  // Set the incident half-facet map
3966  for( int k = 0; k < nv; k++ )
3967  {
3968  std::vector< EntityHandle > set_childents;
3969  std::vector< int > set_childlfids;
3970  for( int j = 0; j < ncomps; j++ )
3971  {
3972  set_childents.push_back( compchildents[k * ncomps + j] );
3973  set_childlfids.push_back( compchildlfids[k * ncomps + j] );
3974  }
3975 
3976  error = ahf->set_incident_map( type, edgverts[k], set_childents, set_childlfids );MB_CHK_ERR( error );
3977  }
3978  }
3979  }
3980  }
3981 
3982  return MB_SUCCESS;
3983 }
3984 
3986  int deg,
3987  int lfid,
3988  int leid,
3989  std::vector< int >& child_ids,
3990  std::vector< int >& child_lvids )
3991 {
3992  int index = ahf->get_index_in_lmap( *_incells.begin() );
3993  int d = get_index_from_degree( deg );
3994 
3995  // int lv0 = ahf->lConnMap3D[index].e2v[leid][0];
3996  // int lv1 = ahf->lConnMap3D[index].e2v[leid][1];
3997  int nvpc = ahf->lConnMap3D[index].num_verts_in_cell;
3998 
3999  int nv = refTemplates[type - 1][d].nv_edge;
4000  int nch = refTemplates[type - 1][d].ents_on_pent[lfid][0];
4001 
4002  for( int i = 0; i < nch; i++ )
4003  {
4004  int id = refTemplates[type - 1][d].ents_on_pent[lfid][i + 1] - 1;
4005  for( int j = 0; j < nvpc; j++ )
4006  {
4007  int lv = refTemplates[type - 1][d].ents_conn[id][j];
4008  for( int k = 0; k < nv; k++ )
4009  {
4010  if( lv == refTemplates[type - 1][d].vert_on_edges[leid][k] )
4011  {
4012  child_ids.push_back( id );
4013  child_lvids.push_back( j );
4014  }
4015  }
4016  }
4017  }
4018 
4019  return MB_SUCCESS;
4020 }
4021 
4022 /* **********************************
4023  * * Boundary Functions *
4024  ************************************/
4025 
4027 {
4028  ErrorCode error;
4029  EntityHandle sibents[27];
4030  int siblids[27];
4031  std::vector< EntityHandle > ent;
4032  std::vector< int > lid;
4033 
4034  int nhf;
4035  if( elementype == MBEDGE )
4036  nhf = 2;
4037  else if( ( elementype == MBTRI ) || ( elementype == MBQUAD ) )
4039  else if( ( elementype == MBTET ) || ( elementype == MBHEX ) )
4040  {
4041  int idx = ahf->get_index_in_lmap( *_incells.begin() );
4042  nhf = ahf->lConnMap3D[idx].num_faces_in_cell;
4043  }
4044  else
4045  MB_SET_ERR( MB_FAILURE, "Requesting vertex boundary information for an unsupported entity type" );
4046 
4048  error = ahf->get_sibling_map( elementype, ent[0], &sibents[0], &siblids[0], nhf );MB_CHK_ERR( error );
4049 
4050  return ( sibents[lid[0]] == 0 );
4051 }
4052 
4054 {
4055  ErrorCode error;
4056  bool is_border = false;
4057  if( meshdim == 1 ) // The edge has a vertex on the boundary in the curve mesh
4058  {
4059  EntityHandle sibents[2];
4060  int siblids[2];
4061  error = ahf->get_sibling_map( MBEDGE, entity, &sibents[0], &siblids[0], 2 );MB_CHK_ERR( error );
4062  for( int i = 0; i < 2; i++ )
4063  {
4064  if( sibents[i] == 0 )
4065  {
4066  is_border = true;
4067  break;
4068  }
4069  }
4070  }
4071  else if( meshdim == 2 ) // The edge is on the boundary of the 2d mesh
4072  {
4073  std::vector< EntityHandle > adjents;
4074  error = ahf->get_up_adjacencies_2d( entity, adjents );MB_CHK_ERR( error );
4075  if( adjents.size() == 1 ) is_border = true;
4076  }
4077  else if( meshdim == 3 ) // The edge is part of a face on the boundary of the 3d mesh
4078  {
4079  std::vector< EntityHandle > adjents;
4080  std::vector< int > leids;
4081  error = ahf->get_up_adjacencies_edg_3d( entity, adjents, &leids );MB_CHK_ERR( error );
4082  assert( !adjents.empty() );
4083 
4084  int index = ahf->get_index_in_lmap( adjents[0] );
4085  int nhf = ahf->lConnMap3D[index].num_faces_in_cell;
4086 
4087  for( int i = 0; i < (int)adjents.size(); i++ )
4088  {
4089  EntityHandle sibents[6];
4090  int siblids[6];
4091  error = ahf->get_sibling_map( elementype, adjents[0], &sibents[0], &siblids[0], nhf );MB_CHK_ERR( error );
4092  for( int k = 0; k < 2; k++ )
4093  {
4094  int hf = ahf->lConnMap3D[index].e2hf[leids[0]][k];
4095  if( sibents[hf] == 0 )
4096  {
4097  is_border = true;
4098  break;
4099  }
4100  }
4101  }
4102  }
4103  return is_border;
4104 }
4105 
4107 {
4108  ErrorCode error;
4109  bool is_border = false;
4110 
4111  if( meshdim == 1 )
4112  MB_SET_ERR( MB_FAILURE, "Requesting boundary information for a face entity type on a curve mesh" );
4113  else if( meshdim == 2 ) // The face has a local edge on the boundary of the 2d mesh
4114  {
4115  EntityHandle sibents[4];
4116  int siblids[4];
4117  int nepf = ahf->lConnMap2D[elementype - 2].num_verts_in_face;
4118  error = ahf->get_sibling_map( elementype, entity, &sibents[0], &siblids[0], nepf );MB_CHK_ERR( error );
4119 
4120  for( int i = 0; i < nepf; i++ )
4121  {
4122  if( sibents[i] == 0 )
4123  {
4124  is_border = true;
4125  break;
4126  }
4127  }
4128  }
4129  else if( meshdim == 3 ) // The face lies on the boundary of the 3d mesh
4130  {
4131  std::vector< EntityHandle > adjents;
4132  error = ahf->get_up_adjacencies_face_3d( entity, adjents );MB_CHK_ERR( error );
4133  if( adjents.size() == 1 ) is_border = true;
4134  }
4135  return is_border;
4136 }
4137 
4139 {
4140  if( meshdim != 3 )
4141  MB_SET_ERR( MB_FAILURE, "Requesting boundary information for a cell entity type on a curve or surface mesh" );
4142 
4143  bool is_border = false;
4144  int index = ahf->get_index_in_lmap( *_incells.begin() );
4145  int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
4146  EntityHandle sibents[6];
4147  int siblids[6];
4148 
4149  ErrorCode error = ahf->get_sibling_map( elementype, entity, &sibents[0], &siblids[0], nfpc );MB_CHK_ERR( error );
4150 
4151  for( int i = 0; i < nfpc; i++ )
4152  {
4153  if( sibents[i] == 0 )
4154  {
4155  is_border = true;
4156  break;
4157  }
4158  }
4159  return is_border;
4160 }
4161 
4162 /* **********************************
4163  * Helper Functions *
4164  ************************************/
4166 {
4167  ErrorCode error;
4168 
4169  if( cur_level )
4170  {
4171  int nverts_prev = level_mesh[cur_level - 1].num_verts;
4172  for( int i = 0; i < nverts_prev; i++ )
4173  {
4174  level_mesh[cur_level].coordinates[0][i] = level_mesh[cur_level - 1].coordinates[0][i];
4175  level_mesh[cur_level].coordinates[1][i] = level_mesh[cur_level - 1].coordinates[1][i];
4176  level_mesh[cur_level].coordinates[2][i] = level_mesh[cur_level - 1].coordinates[2][i];
4177  }
4178  }
4179  else // Copy the vertices from the input mesh
4180  {
4181  int nverts_in = _inverts.size();
4182  std::vector< double > vcoords( 3 * nverts_in );
4183  error = mbImpl->get_coords( _inverts, &vcoords[0] );MB_CHK_ERR( error );
4184 
4185  for( int i = 0; i < nverts_in; i++ )
4186  {
4187  level_mesh[cur_level].coordinates[0][i] = vcoords[3 * i];
4188  level_mesh[cur_level].coordinates[1][i] = vcoords[3 * i + 1];
4189  level_mesh[cur_level].coordinates[2][i] = vcoords[3 * i + 2];
4190  }
4191  }
4192  return MB_SUCCESS;
4193  // To add: Map from old vertices to new duplicates: NOT NEEDED
4194 }
4195 
4197  int cur_level,
4198  int deg,
4199  std::vector< EntityHandle >& trackvertsC_edg,
4200  std::vector< EntityHandle >& trackvertsC_face,
4201  EntityHandle* vbuffer )
4202 {
4203  // The vertices in the vbuffer are added to appropriate edges and faces of cells that are
4204  // incident on the working cell.
4205  ErrorCode error;
4206 
4207  EntityHandle cstart_prev;
4208  if( cur_level )
4209  cstart_prev = level_mesh[cur_level - 1].start_cell;
4210  else
4211  cstart_prev = *_incells.begin();
4212 
4213  EntityType cell_type = mbImpl->type_from_handle( cstart_prev );
4214  int cindex = cell_type - 1;
4215  int d = get_index_from_degree( deg );
4216 
4217  int nve = refTemplates[cindex][d].nv_edge;
4218  int nvf = refTemplates[cindex][d].nv_face;
4219 
4220  int index = ahf->get_index_in_lmap( *( _incells.begin() ) );
4221  int nepc = ahf->lConnMap3D[index].num_edges_in_cell;
4222  int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
4223 
4224  // Step 1: Add the vertices on an edge of the working cell to tracking array of incident cells.
4225  for( int i = 0; i < nepc; i++ )
4226  {
4227  // Add the vertices to edges of the current cell
4228  for( int j = 0; j < nve; j++ )
4229  {
4230  int id = refTemplates[cindex][d].vert_on_edges[i][j];
4231  int idx = cid - cstart_prev;
4232  int aid = idx * nve * nepc + nve * i + j;
4233 
4234  if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
4235  }
4236 
4237  // Obtain all the incident cells
4238  std::vector< EntityHandle > inc_cids;
4239  std::vector< int > inc_leids, inc_orient;
4240 
4241  error = ahf->get_up_adjacencies_edg_3d( cid, i, inc_cids, &inc_leids, &inc_orient );MB_CHK_ERR( error );
4242 
4243  if( inc_cids.size() == 1 ) continue;
4244 
4245  // Add the vertices to the edges of the incident cells
4246  for( int k = 0; k < (int)inc_cids.size(); k++ )
4247  {
4248  if( inc_cids[k] == cid ) continue;
4249 
4250  int idx = inc_cids[k] - cstart_prev;
4251 
4252  if( inc_orient[k] ) // Same edge direction as the current edge
4253  {
4254  for( int j = 0; j < nve; j++ )
4255  {
4256  int id = refTemplates[cindex][d].vert_on_edges[i][j];
4257  int aid = idx * nve * nepc + nve * inc_leids[k] + j;
4258 
4259  if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
4260  }
4261  }
4262  else
4263  {
4264  for( int j = 0; j < nve; j++ )
4265  {
4266  int id = refTemplates[cindex][d].vert_on_edges[i][nve - j - 1];
4267  int aid = idx * nve * nepc + nve * inc_leids[k] + j;
4268 
4269  if( !trackvertsC_edg[aid] ) trackvertsC_edg[aid] = vbuffer[id];
4270  }
4271  }
4272  }
4273  }
4274 
4275  // Step 2: Add the vertices on a face of the working cell to tracking array of incident cells.
4276  if( nvf )
4277  {
4278 
4279  for( int i = 0; i < nfpc; i++ )
4280  {
4281  // Add vertices to the tracking array of vertices on faces for the current cell
4282  std::vector< EntityHandle > face_vbuf( nvf, 0 );
4283  for( int j = 0; j < nvf; j++ )
4284  {
4285  int id = refTemplates[cindex][d].vert_on_faces[i][j];
4286  int idx = cid - cstart_prev;
4287  int aid = idx * nvf * nfpc + nvf * i + j;
4288 
4289  if( !trackvertsC_face[aid] ) trackvertsC_face[aid] = vbuffer[id];
4290 
4291  face_vbuf[j] = vbuffer[id];
4292  }
4293 
4294  // Obtain all the incident cells
4295  std::vector< EntityHandle > sib_cids;
4296  std::vector< int > sib_lfids;
4297  error = ahf->get_up_adjacencies_face_3d( cid, i, sib_cids, &sib_lfids );MB_CHK_ERR( error );
4298 
4299  if( sib_cids.size() == 1 ) continue;
4300 
4301  // Reorder the vertex local ids incident on the half-face
4302  std::vector< int > id_sib( nvf );
4303  for( int k = 0; k < nvf; k++ )
4304  id_sib[k] = 0;
4305 
4306  error = reorder_indices( cur_level, deg, sib_cids[1], sib_lfids[1], cid, i, 0, &id_sib[0] );MB_CHK_ERR( error );
4307 
4308  // Add vertices to the tracking array of vertices on faces for the sibling cell of the
4309  // current cell
4310  for( int j = 0; j < nvf; j++ )
4311  {
4312  int idx = sib_cids[1] - cstart_prev;
4313  int aid = idx * nvf * nfpc + nvf * sib_lfids[1] + j;
4314 
4315  if( !trackvertsC_face[aid] ) trackvertsC_face[aid] = face_vbuf[id_sib[j] - 1];
4316  }
4317  }
4318  }
4319  return MB_SUCCESS;
4320 }
4321 
4323  int deg,
4324  EntityHandle cell,
4325  int lfid,
4326  EntityHandle sib_cell,
4327  int sib_lfid,
4328  int index,
4329  int* id_sib )
4330 {
4331  // Reorders the indices of either vertices or children cell local ids to match with order of the
4332  // given cell and a local face. index = 0 : vertices,
4333  // = 1 : face
4334 
4335  assert( deg == 2 || deg == 3 );
4336 
4337  ErrorCode error;
4338  int idx = ahf->get_index_in_lmap( *_incells.begin() );
4339  int nvF = ahf->lConnMap3D[idx].hf2v_num[lfid];
4340  int nco = permutation[nvF - 3].num_comb;
4341 
4342  if( !index && ( ( nvF == 3 && deg == 3 ) || ( nvF == 4 && deg == 2 ) ) )
4343  {
4344  id_sib[0] = 1;
4345  }
4346  else
4347  {
4348  // Get connectivity of the cell and its sibling cell
4349  std::vector< EntityHandle > conn, sib_conn;
4350  error = get_connectivity( cell, cur_level, conn );MB_CHK_ERR( error );
4351 
4352  error = get_connectivity( sib_cell, cur_level, sib_conn );MB_CHK_ERR( error );
4353 
4354  // Get the connectivity of the local face in the cell and its sibling
4355  std::vector< EntityHandle > lface( nvF );
4356  std::vector< EntityHandle > lface_sib( nvF );
4357  for( int i = 0; i < nvF; i++ )
4358  {
4359  int id = ahf->lConnMap3D[idx].hf2v[lfid][i];
4360  lface[i] = conn[id];
4361 
4362  id = ahf->lConnMap3D[idx].hf2v[sib_lfid][i];
4363  lface_sib[i] = sib_conn[id];
4364  }
4365 
4366  // Find the combination
4367  int c = 0;
4368  for( int i = 0; i < nco; i++ )
4369  {
4370  int count = 0;
4371  for( int j = 0; j < nvF; j++ )
4372  {
4373  int id = permutation[nvF - 3].comb[i][j];
4374  if( lface[j] == lface_sib[id] ) count += 1;
4375  }
4376 
4377  if( count == nvF )
4378  {
4379  c = i;
4380  break;
4381  }
4382  }
4383 
4384  if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
4385 
4386  // Get the ordered indices
4387  if( ( ( !index ) && ( nvF == 4 ) && ( deg == 3 ) ) || ( deg == 2 ) )
4388  {
4389  for( int i = 0; i < 4; i++ )
4390  id_sib[i] = permutation[nvF - 3].porder2[c][i];
4391  }
4392  else
4393  {
4394  for( int i = 0; i < 9; i++ )
4395  id_sib[i] = permutation[nvF - 3].porder3[c][i];
4396  }
4397  }
4398 
4399  return MB_SUCCESS;
4400 }
4401 
4403  EntityHandle* face1_conn,
4404  EntityHandle* face2_conn,
4405  int nvF,
4406  std::vector< int >& lemap,
4407  std::vector< int >& vidx,
4408  int* leorient )
4409 {
4410  // Given the connectivities of two faces, get the permuted indices w.r.t first face.
4411  // Step 1: First find the orientation
4412  int nco = permutation[nvF - 3].num_comb;
4413  int c = 0;
4414  for( int i = 0; i < nco; i++ )
4415  {
4416  int count = 0;
4417  for( int j = 0; j < nvF; j++ )
4418  {
4419  int id = permutation[nvF - 3].comb[i][j];
4420  if( face1_conn[j] == face2_conn[id] ) count += 1;
4421  }
4422 
4423  if( count == nvF )
4424  {
4425  c = i;
4426  break;
4427  }
4428  }
4429 
4430  if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
4431 
4432  // Add the corresponding local edges
4433  lemap.reserve( nvF );
4434  for( int i = 0; i < nvF; i++ )
4435  {
4436  lemap.push_back( permutation[nvF - 3].lemap[c][i] );
4437  }
4438  if( leorient ) leorient[0] = permutation[nvF - 3].orient[c];
4439 
4440  if( nvF == 3 && deg == 2 ) return MB_SUCCESS;
4441 
4442  if( ( nvF == 3 && deg == 3 ) || ( nvF == 4 && deg == 2 ) )
4443  {
4444  vidx.push_back( 1 );
4445  }
4446  else if( nvF == 4 && deg == 3 )
4447  {
4448  for( int i = 0; i < 4; i++ )
4449  vidx.push_back( permutation[nvF - 3].porder2[c][i] );
4450  }
4451 
4452  return MB_SUCCESS;
4453 }
4454 
4455 ErrorCode NestedRefine::reorder_indices( int deg, int nvF, int comb, int* childfid_map )
4456 {
4457  // Given connectivities of two faces and a degree, get the permuted indices of the children
4458  // faces w.r.t first face.
4459 
4460  assert( deg == 2 || deg == 3 );
4461 
4462  // Get the ordered indices
4463  if( deg == 2 )
4464  {
4465  for( int i = 0; i < 4; i++ )
4466  childfid_map[i] = permutation[nvF - 3].porder2[comb][i];
4467  }
4468  else
4469  {
4470  for( int i = 0; i < 9; i++ )
4471  childfid_map[i] = permutation[nvF - 3].porder3[comb][i];
4472  }
4473 
4474  return MB_SUCCESS;
4475 }
4476 
4478  EntityHandle* face2_conn,
4479  int nvF,
4480  int* conn_map,
4481  int& comb,
4482  int* orient )
4483 {
4484  // Given connectivities of two faces and a degree, get the permuted indices of the children
4485  // faces w.r.t first face.
4486 
4487  // Step 1: First find the combination
4488  int nco = permutation[nvF - 3].num_comb;
4489  int c = 0;
4490  for( int i = 0; i < nco; i++ )
4491  {
4492  int count = 0;
4493  for( int j = 0; j < nvF; j++ )
4494  {
4495  int id = permutation[nvF - 3].comb[i][j];
4496  if( face1_conn[j] == face2_conn[id] ) count += 1;
4497  }
4498 
4499  if( count == nvF )
4500  {
4501  c = i;
4502  break;
4503  }
4504  }
4505 
4506  if( c > nco ) MB_SET_ERR( MB_FAILURE, "Getting a combination number more than currently supported" );
4507 
4508  comb = c;
4509 
4510  if( orient ) orient[0] = permutation[nvF - 3].orient[c];
4511 
4512  for( int j = 0; j < nvF; j++ )
4513  {
4514  conn_map[j] = permutation[nvF - 3].comb[c][j];
4515  }
4516 
4517  return MB_SUCCESS;
4518 }
4519 
4520 ErrorCode NestedRefine::count_subentities( EntityHandle set, int cur_level, int* nedges, int* nfaces )
4521 {
4522  ErrorCode error;
4523 
4524  if( cur_level >= 0 )
4525  {
4526  Range edges, faces, cells;
4527 
4528  error = mbImpl->get_entities_by_dimension( set, 1, edges );MB_CHK_ERR( error );
4529 
4530  error = mbImpl->get_entities_by_dimension( set, 2, faces );MB_CHK_ERR( error );
4531 
4532  error = mbImpl->get_entities_by_dimension( set, 3, cells );MB_CHK_ERR( error );
4533 
4534  error = ahf->count_subentities( edges, faces, cells, nedges, nfaces );MB_CHK_ERR( error );
4535  }
4536  else
4537  {
4538  error = ahf->count_subentities( _inedges, _infaces, _incells, nedges, nfaces );MB_CHK_ERR( error );
4539  }
4540 
4541  return MB_SUCCESS;
4542 }
4543 
4544 ErrorCode NestedRefine::get_octahedron_corner_coords( int cur_level, int deg, EntityHandle* vbuffer, double* ocoords )
4545 {
4546  int lid[6] = { 0, 0, 0, 0, 0, 0 };
4547 
4548  if( deg == 2 )
4549  {
4550  lid[0] = 5;
4551  lid[1] = 8;
4552  lid[2] = 9;
4553  lid[3] = 6;
4554  lid[4] = 4;
4555  lid[5] = 7;
4556  }
4557  else if( deg == 3 )
4558  {
4559  lid[0] = 19;
4560  lid[1] = 16;
4561  lid[2] = 18;
4562  lid[3] = 9;
4563  lid[4] = 4;
4564  lid[5] = 10;
4565  }
4566 
4567  EntityHandle vstart = level_mesh[cur_level].start_vertex;
4568 
4569  for( int i = 0; i < 6; i++ )
4570  {
4571  EntityHandle vid = vbuffer[lid[i]];
4572  ocoords[3 * i] = level_mesh[cur_level].coordinates[0][vid - vstart];
4573  ocoords[3 * i + 1] = level_mesh[cur_level].coordinates[1][vid - vstart];
4574  ocoords[3 * i + 2] = level_mesh[cur_level].coordinates[2][vid - vstart];
4575  }
4576 
4577  return MB_SUCCESS;
4578 }
4579 
4580 int NestedRefine::find_shortest_diagonal_octahedron( int cur_level, int deg, EntityHandle* vbuffer )
4581 {
4582  ErrorCode error;
4583  double coords[18];
4584  error = get_octahedron_corner_coords( cur_level, deg, vbuffer, coords );
4585  if( error != MB_SUCCESS ) MB_SET_ERR( MB_FAILURE, "Error in obtaining octahedron corner coordinates" );
4586 
4587  int diag_map[6] = { 1, 3, 2, 4, 5, 0 };
4588  double length = std::numeric_limits< double >::max();
4589 
4590  int diag = 0;
4591  double x, y, z;
4592  x = y = z = 0;
4593 
4594  for( int d = 0; d < 3; d++ )
4595  {
4596  int id1 = diag_map[2 * d];
4597  int id2 = diag_map[2 * d + 1];
4598  x = coords[3 * id1] - coords[3 * id2];
4599  y = coords[3 * id1 + 1] - coords[3 * id2 + 1];
4600  z = coords[3 * id1 + 2] - coords[3 * id2 + 2];
4601  double dlen = sqrt( x * x + y * y + z * z );
4602  if( dlen < length )
4603  {
4604  length = dlen;
4605  diag = d + 1;
4606  }
4607  }
4608 
4609  return diag;
4610 }
4611 
4613 {
4614  ErrorCode error;
4615  // Given a vertex, find its local id in the given entity
4616  std::vector< EntityHandle > conn;
4617 
4618  error = get_connectivity( ent, level + 1, conn );
4619  if( error != MB_SUCCESS ) MB_SET_ERR( MB_FAILURE, "Error in getting connectivity of the requested entity" );
4620 
4621  int lid = -1;
4622  for( int i = 0; i < (int)conn.size(); i++ )
4623  {
4624  if( conn[i] == vid )
4625  {
4626  lid = i;
4627  break;
4628  }
4629  }
4630  if( lid < 0 ) MB_SET_ERR( MB_FAILURE, "Error in getting local vertex id in the given entity" );
4631  return lid;
4632 }
4633 
4635 {
4636  int d = deg_index.find( degree )->second;
4637  return d;
4638 }
4639 
4640 /*
4641 ErrorCode NestedRefine::print_maps_1D(int level)
4642 {
4643  ErrorCode error;
4644  int nv, ne;
4645  nv = level_mesh[level].num_verts;
4646  ne = level_mesh[level].num_edges;
4647 
4648  EntityHandle start_edge = level_mesh[level].start_edge;
4649 
4650  //V2HV
4651  std::cout<<"<V2HV_EID, V2HV_LVID>"<<std::endl;
4652  for (int i=0; i<nv; i++)
4653  {
4654  EntityHandle eid=0;
4655  int lvid=0;
4656  EntityHandle vid = level_mesh[level].start_vertex+i;
4657  error = ahf->get_incident_map(MBEDGE, vid, eid, lvid); MB_CHK_ERR(error);
4658 
4659  std::cout<<"For vertex = "<<vid<<"::Incident halfvertex "<<eid<<" "<<lvid<<std::endl;
4660  }
4661 
4662  //SIBHVS
4663  std::cout<<"start_edge = "<<start_edge<<std::endl;
4664  std::cout<<"<SIBHVS_EID,SIBHVS_LVID>"<<std::endl;
4665  for (int i=0; i<ne; i++)
4666  {
4667  EntityHandle ent = start_edge+i;
4668 
4669  EntityHandle eid[2]; int lvid[2];
4670  error = ahf->get_sibling_map(MBEDGE, ent, &eid[0], &lvid[0], 2); MB_CHK_ERR(error);
4671  std::cout<<"<"<<eid[0]<<","<<lvid[0]<<">"<<" "<<"<"<<eid[1]<<","<<lvid[1]<<">"<<std::endl;
4672  }
4673 
4674  return MB_SUCCESS;
4675 }
4676 
4677 ErrorCode NestedRefine::print_maps_2D(int level, EntityType type)
4678 {
4679  ErrorCode error;
4680  int nv, nf;
4681  nv = level_mesh[level].num_verts;
4682  nf = level_mesh[level].num_faces;
4683 
4684  EntityHandle start_face = level_mesh[level].start_face;
4685 
4686  //V2HV
4687  std::cout<<"<V2HE_FID, V2HE_LEID>"<<std::endl;
4688  for (int i=0; i<nv; i++)
4689  {
4690  EntityHandle fid=0;
4691  int leid=0;
4692  EntityHandle vid = level_mesh[level].start_vertex+i;
4693  error = ahf->get_incident_map(type, vid, fid, leid); MB_CHK_ERR(error);
4694 
4695  std::cout<<"For vertex = "<<vid<<"::Incident halfedge "<<fid<<" "<<leid<<std::endl;
4696  }
4697 
4698  //SIBHES
4699  std::cout<<"start_face = "<<start_face<<std::endl;
4700  std::cout<<"<SIBHES_FID,SIBHES_LEID>"<<std::endl;
4701  EntityType ftype = mbImpl->type_from_handle(*_infaces.begin());
4702  int nepf = ahf->lConnMap2D[ftype-2].num_verts_in_face;
4703 
4704  EntityHandle *fid = new EntityHandle[nepf];
4705  int *leid = new int[nepf];
4706 
4707  for (int i=0; i<nf; i++)
4708  {
4709  for (int j=0; j<nepf; j++)
4710  {
4711  fid[j] = 0;
4712  leid[j] = 0;
4713  }
4714 
4715  EntityHandle ent = start_face+i;
4716  error = ahf->get_sibling_map(type, ent, fid, leid, nepf); MB_CHK_ERR(error);
4717 
4718  for (int j=0; j<nepf; j++){
4719  std::cout<<"<"<<fid[j]<<","<<leid[j]<<">"<<" ";
4720  }
4721  std::cout<<std::endl;
4722  }
4723 
4724  delete [] fid;
4725  delete [] leid;
4726 
4727  return MB_SUCCESS;
4728 }
4729 
4730 ErrorCode NestedRefine::print_maps_3D(int level, EntityType type)
4731 {
4732  ErrorCode error;
4733  int nv, nc;
4734  nv = level_mesh[level].num_verts;
4735  nc = level_mesh[level].num_cells;
4736  EntityHandle start_cell = level_mesh[level].start_cell;
4737 
4738  //V2HF
4739  std::cout<<"<V2HF_CID, V2HF_LFID>"<<std::endl;
4740  for (int i=0; i<nv; i++)
4741  {
4742  EntityHandle cid=0;
4743  int lfid=0;
4744  EntityHandle vid = level_mesh[level].start_vertex+i;
4745  error = ahf->get_incident_map(type, vid, cid, lfid); MB_CHK_ERR(error);
4746 
4747  std::cout<<"For vertex = "<<vid<<"::Incident halfface "<<cid<<" "<<lfid<<std::endl;
4748  }
4749 
4750  //SIBHFS
4751  std::cout<<"start_cell = "<<start_cell<<std::endl;
4752  std::cout<<"<SIBHFS_CID,SIBHFS_LFID>"<<std::endl;
4753  int index = ahf->get_index_in_lmap(start_cell);
4754  int nfpc = ahf->lConnMap3D[index].num_faces_in_cell;
4755 
4756  EntityHandle *cid = new EntityHandle[nfpc];
4757  int *lfid = new int[nfpc];
4758  for (int i=0; i<nc; i++)
4759  {
4760  for (int k=0; k<nfpc; k++)
4761  {
4762  cid[k] = 0;
4763  lfid[k] = 0;
4764  }
4765 
4766  EntityHandle ent = start_cell+i;
4767  error = ahf->get_sibling_map(type, ent, cid, lfid, nfpc); MB_CHK_ERR(error);
4768 
4769  for (int j=0; j<nfpc; j++){
4770  std::cout<<"<"<<cid[j]<<","<<lfid[j]<<">"<<" ";
4771  }
4772  std::cout<<std::endl;
4773  }
4774 
4775  delete [] cid;
4776  delete [] lfid;
4777 
4778  return MB_SUCCESS;
4779 }
4780 
4781 */
4782 } // namespace moab