 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "moab/NestedRefine.hpp"
2 #include "moab/NestedRefineTemplates.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"
11 #include "moab/ParallelMergeMesh.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
24 NestedRefine::NestedRefine( Core* impl, ParallelComm* comm, EntityHandle rset )
25 : mbImpl( impl ), pcomm( comm ), _rset( rset )
26 {
27 ErrorCode error;
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
43 NestedRefine::~NestedRefine()
44 {
45 #ifdef MOAB_HAVE_AHF
46 ahf = NULL;
47 #else
48 delete ahf;
49 #endif
50 delete tm;
51 }
52
53 ErrorCode NestedRefine::initialize()
54 {
55 ErrorCode error;
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
71 error = ahf->initialize();MB_CHK_ERR( error );
72 error = ahf->get_entity_ranges( _inverts, _inedges, _infaces, _incells );MB_CHK_ERR( error );
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;
95 elementype = MBEDGE;
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
114 ErrorCode NestedRefine::generate_mesh_hierarchy( int num_level,
115 int* level_degrees,
116 std::vector< EntityHandle >& level_sets,
117 bool optimize )
118 {
119 assert( num_level > 0 );
120 nlevels = num_level;
121
122 ErrorCode error;
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 {
155 ErrorCode error;
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 {
216 ErrorCode error;
217 error = mbImpl->get_coords( verts, num_verts, coords );MB_CHK_ERR( error );
218 }
219
220 return MB_SUCCESS;
221 }
222
223 ErrorCode NestedRefine::get_adjacencies( const EntityHandle source_entity,
224 const unsigned int target_dimension,
225 std::vector< EntityHandle >& target_entities )
226
227 {
228 ErrorCode error;
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
285 ErrorCode NestedRefine::parent_to_child( EntityHandle parent,
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 {
323 int d = get_index_from_degree( level_dsequence[i] );
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 {
334 EntityHandle child;
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
348 ErrorCode NestedRefine::vertex_to_entities_up( EntityHandle vertex,
349 int vert_level,
350 int parent_level,
351 std::vector< EntityHandle >& incident_entities )
352 {
353 assert( vert_level > parent_level );
354 ErrorCode error;
355
356 // Step 1: Get the incident entities at the current level
357 std::vector< EntityHandle > inents;
358 if( meshdim == 1 )
359 {
360 error = ahf->get_up_adjacencies_1d( vertex, inents );MB_CHK_ERR( error );
361 }
362 else if( meshdim == 2 )
363 {
364 error = ahf->get_up_adjacencies_vert_2d( vertex, inents );MB_CHK_ERR( error );
365 }
366 else if( meshdim == 3 )
367 {
368 error = ahf->get_up_adjacencies_vert_3d( vertex, inents );MB_CHK_ERR( error );
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
388 ErrorCode NestedRefine::vertex_to_entities_down( EntityHandle vertex,
389 int vert_level,
390 int child_level,
391 std::vector< EntityHandle >& incident_entities )
392 {
393 assert( vert_level < child_level );
394 ErrorCode error;
395
396 // Step 1: Get the incident entities at the current level
397 std::vector< EntityHandle > inents;
398 if( meshdim == 1 )
399 {
400 error = ahf->get_up_adjacencies_1d( vertex, inents );MB_CHK_ERR( error );
401 }
402 else if( meshdim == 2 )
403 {
404 error = ahf->get_up_adjacencies_vert_2d( vertex, inents );MB_CHK_ERR( error );
405 }
406 else if( meshdim == 3 )
407 {
408 error = ahf->get_up_adjacencies_vert_3d( vertex, inents );MB_CHK_ERR( error );
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
425 ErrorCode NestedRefine::get_vertex_duplicates( EntityHandle vertex, int level, EntityHandle& dupvertex )
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
435 bool NestedRefine::is_entity_on_boundary( const EntityHandle& entity )
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 {
456 ErrorCode error;
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;
465 error = pcomm->exchange_tags( GLOBAL_ID_TAG_NAME, empty_range );MB_CHK_ERR( error );
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
496 ErrorCode NestedRefine::update_special_tags( int level, EntityHandle& lset )
497 {
498 assert( level > 0 && level < nlevels + 1 );
499
500 ErrorCode error;
501 std::vector< Tag > mtags( 3 );
502
503 error = mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mtags[0] );MB_CHK_ERR( error );
504 error = mbImpl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mtags[1] );MB_CHK_ERR( error );
505 error = mbImpl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mtags[2] );MB_CHK_ERR( error );
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( ¢s[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 {
566 ErrorCode error;
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
624 ErrorCode NestedRefine::create_hm_storage_single_level( EntityHandle* set, int cur_level, int estL[4] )
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;
644 error = mbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, gidtag );MB_CHK_ERR( error );
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
700 error = ahf->update_entity_ranges( *set );MB_CHK_ERR( error );
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 {
717 ErrorCode error;
718
719 Tag gidtag;
720 error = mbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, gidtag );MB_CHK_ERR( error );
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
746 error = copy_vertices_from_prev_level( l );MB_CHK_ERR( error );
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 }
778 timeall.tm_total = timeall.tm_refine + timeall.tm_resolve;
779
780 return MB_SUCCESS;
781 }
782
783 ErrorCode NestedRefine::construct_hm_entities( int cur_level, int deg )
784 {
785 ErrorCode error;
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 {
806 ErrorCode error;
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
833 EntityHandle 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
917 ErrorCode NestedRefine::construct_hm_1D( int cur_level,
918 int deg,
919 EntityType type,
920 std::vector< EntityHandle >& trackverts )
921 {
922 ErrorCode error;
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
967 EntityHandle edge;
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
1096 EntityHandle 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
1223 ErrorCode NestedRefine::construct_hm_2D( int cur_level,
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
1282 EntityHandle 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
1730 ErrorCode NestedRefine::compute_coordinates( int cur_level,
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,
1956 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id );MB_CHK_ERR( error );
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,
1991 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE, &dum_id );MB_CHK_ERR( error );
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
3005 ErrorCode NestedRefine::update_local_ahf( int deg,
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
3100 ErrorCode NestedRefine::update_local_ahf( int deg,
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
3153 ErrorCode NestedRefine::update_global_ahf_1D( int cur_level, int deg )
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
3273 ErrorCode NestedRefine::update_global_ahf_1D_sub( int cur_level, int deg )
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
3398 ErrorCode NestedRefine::update_ahf_1D( int cur_level )
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
3403 error = ahf->determine_incident_halfverts( level_mesh[cur_level].edges );MB_CHK_ERR( error );
3404
3405 return MB_SUCCESS;
3406 }
3407
3408 ErrorCode NestedRefine::update_global_ahf_2D( int cur_level, int deg )
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
3562 ErrorCode NestedRefine::update_global_ahf_2D_sub( int cur_level, int deg )
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
3985 ErrorCode NestedRefine::get_lid_inci_child( EntityType type,
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
4026 bool NestedRefine::is_vertex_on_boundary( const EntityHandle& vertex )
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 ) )
4038 nhf = ahf->lConnMap2D[elementype - 2].num_verts_in_face;
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
4047 error = ahf->get_incident_map( elementype, vertex, ent, lid );MB_CHK_ERR( error );
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
4053 bool NestedRefine::is_edge_on_boundary( const EntityHandle& entity )
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
4106 bool NestedRefine::is_face_on_boundary( const EntityHandle& entity )
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
4138 bool NestedRefine::is_cell_on_boundary( const EntityHandle& entity )
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 ************************************/
4165 ErrorCode NestedRefine::copy_vertices_from_prev_level( int cur_level )
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
4196 ErrorCode NestedRefine::update_tracking_verts( EntityHandle cid,
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
4322 ErrorCode NestedRefine::reorder_indices( int cur_level,
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
4402 ErrorCode NestedRefine::reorder_indices( int deg,
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
4477 ErrorCode NestedRefine::reorder_indices( EntityHandle* face1_conn,
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
4612 int NestedRefine::get_local_vid( EntityHandle vid, EntityHandle ent, int level )
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
4634 int NestedRefine::get_index_from_degree( int degree )
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
1.9.1.