Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
HalfFacetRep.cpp
Go to the documentation of this file.
1 /**
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 #include "moab/HalfFacetRep.hpp"
17 #include "Internals.hpp"
18 #include <iostream>
19 #include <cassert>
20 #include <vector>
21 #include <map>
22 #include "MBTagConventions.hpp"
23 #include "moab/ScdInterface.hpp"
24 #ifdef MOAB_HAVE_MPI
25 #include "moab/ParallelComm.hpp"
26 #endif
27 
28 namespace moab
29 {
30 
31 inline EntityHandle CREATE_HALFFACET( const unsigned lid, const EntityID id )
32 {
33  assert( id <= MB_END_ID && lid < 6 );
34  return ( ( (HFacet)lid ) << MB_ID_WIDTH ) | id;
35 }
36 inline EntityID FID_FROM_HALFFACET( HFacet handle )
37 {
38  return ( handle & MB_ID_MASK );
39 }
40 inline int LID_FROM_HALFFACET( HFacet handle )
41 {
42  return static_cast< int >( handle >> MB_ID_WIDTH );
43 }
44 
45 HalfFacetRep::HalfFacetRep( Core* impl, ParallelComm* comm, moab::EntityHandle rset, bool filter_ghosts )
46  : thismeshtype( CURVE ), mb( impl ), pcomm( comm ), _rset( rset ), _filterghost( filter_ghosts )
47 {
48  assert( NULL != impl );
49  mInitAHFmaps = false;
50  chk_mixed = false;
51  is_mixed = false;
52 }
53 
54 HalfFacetRep::~HalfFacetRep() {}
55 
56 MESHTYPE HalfFacetRep::get_mesh_type( int nverts, int nedges, int nfaces, int ncells )
57 {
58  MESHTYPE mesh_type = CURVE;
59 
60  if( nverts && nedges && ( !nfaces ) && ( !ncells ) )
61  mesh_type = CURVE;
62  else if( nverts && !nedges && nfaces && !ncells )
63  mesh_type = SURFACE;
64  else if( nverts && nedges && nfaces && !ncells )
65  mesh_type = SURFACE_MIXED;
66  else if( nverts && !nedges && !nfaces && ncells )
67  mesh_type = VOLUME;
68  else if( nverts && nedges && !nfaces && ncells )
69  mesh_type = VOLUME_MIXED_1;
70  else if( nverts && !nedges && nfaces && ncells )
71  mesh_type = VOLUME_MIXED_2;
72  else if( nverts && nedges && nfaces && ncells )
73  mesh_type = VOLUME_MIXED;
74 
75  return mesh_type;
76 }
77 
78 const HalfFacetRep::adj_matrix HalfFacetRep::adjMatrix[7] = {
79  // Stores the adjacency matrix for each mesh type.
80  // CURVE
81  { { { 0, 1, 0, 0 }, { 1, 1, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
82 
83  // SURFACE
84  { { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 1, 0 }, { 0, 0, 0, 0 } } },
85 
86  // SURFACE_MIXED
87  { { { 0, 1, 1, 0 }, { 1, 1, 1, 0 }, { 1, 1, 1, 0 }, { 0, 0, 0, 0 } } },
88 
89  // VOLUME
90  { { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 1 } } },
91 
92  // VOLUME_MIXED_1
93  { { { 0, 1, 0, 1 }, { 1, 1, 0, 1 }, { 0, 0, 0, 0 }, { 1, 1, 0, 1 } } },
94 
95  // VOLUME_MIXED_2
96  { { { 0, 0, 1, 1 }, { 0, 0, 0, 0 }, { 1, 0, 1, 1 }, { 1, 0, 1, 1 } } },
97 
98  // VOLUME_MIXED
99  { { { 0, 1, 1, 1 }, { 1, 1, 1, 1 }, { 1, 1, 1, 1 }, { 1, 1, 1, 1 } } } };
100 
101 int HalfFacetRep::get_index_for_meshtype( MESHTYPE mesh_type )
102 {
103  int index = 0;
104  switch( mesh_type )
105  {
106  case CURVE:
107  index = 0;
108  break;
109  case SURFACE:
110  index = 1;
111  break;
112  case SURFACE_MIXED:
113  index = 2;
114  break;
115  case VOLUME:
116  index = 3;
117  break;
118  case VOLUME_MIXED_1:
119  index = 4;
120  break;
121  case VOLUME_MIXED_2:
122  index = 5;
123  break;
124  case VOLUME_MIXED:
125  index = 6;
126  break;
127  }
128  return index;
129 }
130 
131 bool HalfFacetRep::check_mixed_entity_type()
132 {
133  if( !chk_mixed )
134  {
135  chk_mixed = true;
136 
137  ErrorCode error;
138  Range felems, celems;
139 
140  error = mb->get_entities_by_dimension( this->_rset, 2, felems );MB_CHK_ERR( error );
141 
142  if( felems.size() )
143  {
144  Range tris, quad, poly;
145  tris = felems.subset_by_type( MBTRI );
146  quad = felems.subset_by_type( MBQUAD );
147  poly = felems.subset_by_type( MBPOLYGON );
148  if( ( tris.size() && quad.size() ) || ( tris.size() && poly.size() ) || ( quad.size() && poly.size() ) )
149  is_mixed = true;
150  if( poly.size() ) is_mixed = true;
151 
152  if( is_mixed ) return is_mixed;
153  }
154 
155  error = mb->get_entities_by_dimension( this->_rset, 3, celems );MB_CHK_ERR( error );
156  if( celems.size() )
157  {
158  Range tet, pyr, prism, hex, polyhed;
159  tet = celems.subset_by_type( MBTET );
160  pyr = celems.subset_by_type( MBPYRAMID );
161  prism = celems.subset_by_type( MBPRISM );
162  hex = celems.subset_by_type( MBHEX );
163  polyhed = celems.subset_by_type( MBPOLYHEDRON );
164  if( ( tet.size() && pyr.size() ) || ( tet.size() && prism.size() ) || ( tet.size() && hex.size() ) ||
165  ( tet.size() && polyhed.size() ) || ( pyr.size() && prism.size() ) || ( pyr.size() && hex.size() ) ||
166  ( pyr.size() && polyhed.size() ) || ( prism.size() && hex.size() ) ||
167  ( prism.size() && polyhed.size() ) || ( hex.size() && polyhed.size() ) )
168  is_mixed = true;
169 
170  if( polyhed.size() ) is_mixed = true;
171  }
172 
173  ScdInterface* scdi = NULL;
174  error = mb->query_interface( scdi );MB_CHK_ERR( error );
175  if( scdi )
176  {
177  Range boxes;
178  error = scdi->find_boxes( boxes );MB_CHK_ERR( error );
179 
180  if( !boxes.empty() ) is_mixed = true;
181  }
182  }
183  return is_mixed;
184 }
185 
186 /*******************************************************
187  * initialize *
188  ******************************************************/
189 
190 ErrorCode HalfFacetRep::initialize()
191 {
192  ErrorCode error;
193 
194  if( !mInitAHFmaps )
195  {
196  mInitAHFmaps = true;
197 #ifdef MOAB_HAVE_MPI
198  if( pcomm && _filterghost )
199  {
200  moab::Range _averts, _aedgs, _afacs, _acels;
201  error = mb->get_entities_by_dimension( this->_rset, 0, _averts, true );MB_CHK_ERR( error );
202  error = mb->get_entities_by_dimension( this->_rset, 1, _aedgs, true );MB_CHK_ERR( error );
203  error = mb->get_entities_by_dimension( this->_rset, 2, _afacs, true );MB_CHK_ERR( error );
204  error = mb->get_entities_by_dimension( this->_rset, 3, _acels, true );MB_CHK_ERR( error );
205 
206  // filter based on parallel status
207  error = pcomm->filter_pstatus( _averts, PSTATUS_GHOST, PSTATUS_NOT, -1, &_verts );MB_CHK_ERR( error );
208  error = pcomm->filter_pstatus( _aedgs, PSTATUS_GHOST, PSTATUS_NOT, -1, &_edges );MB_CHK_ERR( error );
209  error = pcomm->filter_pstatus( _afacs, PSTATUS_GHOST, PSTATUS_NOT, -1, &_faces );MB_CHK_ERR( error );
210  error = pcomm->filter_pstatus( _acels, PSTATUS_GHOST, PSTATUS_NOT, -1, &_cells );MB_CHK_ERR( error );
211  }
212  else
213  {
214  error = mb->get_entities_by_dimension( this->_rset, 0, _verts, true );MB_CHK_ERR( error );
215  error = mb->get_entities_by_dimension( this->_rset, 1, _edges, true );MB_CHK_ERR( error );
216  error = mb->get_entities_by_dimension( this->_rset, 2, _faces, true );MB_CHK_ERR( error );
217  error = mb->get_entities_by_dimension( this->_rset, 3, _cells, true );MB_CHK_ERR( error );
218  }
219 #else
220  error = mb->get_entities_by_dimension( this->_rset, 0, _verts, true );MB_CHK_ERR( error );
221  error = mb->get_entities_by_dimension( this->_rset, 1, _edges, true );MB_CHK_ERR( error );
222  error = mb->get_entities_by_dimension( this->_rset, 2, _faces, true );MB_CHK_ERR( error );
223  error = mb->get_entities_by_dimension( this->_rset, 3, _cells, true );MB_CHK_ERR( error );
224 
225 #endif
226 
227  int nverts = _verts.size();
228  int nedges = _edges.size();
229  int nfaces = _faces.size();
230  int ncells = _cells.size();
231 
232  MESHTYPE mesh_type = get_mesh_type( nverts, nedges, nfaces, ncells );
233  thismeshtype = mesh_type;
234 
235  // Initialize mesh type specific maps
236  if( thismeshtype == CURVE )
237  {
238  error = init_curve();MB_CHK_ERR( error );
239  }
240  else if( thismeshtype == SURFACE )
241  {
242  error = init_surface();MB_CHK_ERR( error );
243  }
244  else if( thismeshtype == SURFACE_MIXED )
245  {
246  error = init_curve();MB_CHK_ERR( error );
247  error = init_surface();MB_CHK_ERR( error );
248  }
249  else if( thismeshtype == VOLUME )
250  {
251  error = init_volume();MB_CHK_ERR( error );
252  }
253  else if( thismeshtype == VOLUME_MIXED_1 )
254  {
255  error = init_curve();MB_CHK_ERR( error );
256  error = init_volume();MB_CHK_ERR( error );
257  }
258  else if( thismeshtype == VOLUME_MIXED_2 )
259  {
260  error = init_surface();MB_CHK_ERR( error );
261  error = init_volume();MB_CHK_ERR( error );
262  }
263  else if( thismeshtype == VOLUME_MIXED )
264  {
265  error = init_curve();MB_CHK_ERR( error );
266  error = init_surface();MB_CHK_ERR( error );
267  error = init_volume();MB_CHK_ERR( error );
268  }
269  }
270  return MB_SUCCESS;
271 }
272 
273 ErrorCode HalfFacetRep::deinitialize()
274 {
275  return MB_SUCCESS;
276 }
277 
278 ErrorCode HalfFacetRep::init_curve()
279 {
280  ErrorCode error;
281 
282  int nv = ID_FROM_HANDLE( *( _verts.end() - 1 ) );
283  int ne = ID_FROM_HANDLE( *( _edges.end() - 1 ) );
284 
285  v2hv.resize( nv, 0 );
286  sibhvs.resize( ne * 2, 0 );
287 
288  error = determine_sibling_halfverts( _verts, _edges );MB_CHK_ERR( error );
289  error = determine_incident_halfverts( _edges );MB_CHK_ERR( error );
290 
291  return MB_SUCCESS;
292 }
293 
294 ErrorCode HalfFacetRep::init_surface()
295 {
296  ErrorCode error;
297  EntityType ftype = mb->type_from_handle( *_faces.begin() );
298  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
299 
300  int nv = ID_FROM_HANDLE( *( _verts.end() - 1 ) );
301  int nf = ID_FROM_HANDLE( *( _faces.end() - 1 ) );
302 
303  v2he.resize( nv, 0 );
304  sibhes.resize( nf * nepf, 0 );
305 
306  // Construct ahf maps
307  error = determine_sibling_halfedges( _faces );MB_CHK_ERR( error );
308  error = determine_incident_halfedges( _faces );MB_CHK_ERR( error );
309 
310  // Initialize queues for storing face and local id's during local search
311  for( int i = 0; i < MAXSIZE; i++ )
312  {
313  queue_fid[i] = 0;
314  queue_lid[i] = 0;
315  trackfaces[i] = 0;
316  }
317 
318  return MB_SUCCESS;
319 }
320 
321 ErrorCode HalfFacetRep::init_volume()
322 {
323  ErrorCode error;
324 
325  // Initialize std::map between cell-types and their index in lConnMap3D
326  cell_index[MBTET] = 0;
327  cell_index[MBPYRAMID] = 1;
328  cell_index[MBPRISM] = 2;
329  cell_index[MBHEX] = 3;
330 
331  int index = get_index_in_lmap( *_cells.begin() );
332  int nfpc = lConnMap3D[index].num_faces_in_cell;
333  int nv = ID_FROM_HANDLE( *( _verts.end() - 1 ) );
334  int nc = ID_FROM_HANDLE( *( _cells.end() - 1 ) );
335  ;
336 
337  v2hf.resize( nv, 0 );
338  sibhfs.resize( nc * nfpc, 0 );
339 
340  // Construct the maps
341  error = determine_sibling_halffaces( _cells );MB_CHK_ERR( error );
342  error = determine_incident_halffaces( _cells );MB_CHK_ERR( error );
343 
344  // Initialize queues for storing face and local id's during local search
345  for( int i = 0; i < MAXSIZE; i++ )
346  {
347  Stkcells[i] = 0;
348  cellq[i] = 0;
349  trackcells[i] = 0;
350  }
351 
352  return MB_SUCCESS;
353 }
354 
355 //////////////////////////////////////////////////
356 ErrorCode HalfFacetRep::print_tags( int dim )
357 {
358  if( dim == 1 )
359  {
360  EntityHandle start_edge = *_edges.begin();
361  std::cout << "start_edge = " << start_edge << std::endl;
362  std::cout << "<SIBHVS_EID,SIBHVS_LVID>" << std::endl;
363 
364  for( Range::iterator i = _edges.begin(); i != _edges.end(); ++i )
365  {
366  EntityHandle eid[2];
367  int lvid[2];
368  int eidx = ID_FROM_HANDLE( *i ) - 1;
369  HFacet hf1 = sibhvs[2 * eidx];
370  HFacet hf2 = sibhvs[2 * eidx + 1];
371  eid[0] = fid_from_halfacet( hf1, MBEDGE );
372  eid[1] = fid_from_halfacet( hf2, MBEDGE );
373  lvid[0] = lid_from_halffacet( hf1 );
374  lvid[1] = lid_from_halffacet( hf2 );
375  std::cout << "Entity = " << *i << " :: <" << eid[0] << "," << lvid[0] << ">"
376  << " "
377  << "<" << eid[1] << "," << lvid[1] << ">" << std::endl;
378  }
379 
380  std::cout << "<V2HV_EID, V2HV_LVID>" << std::endl;
381 
382  for( Range::iterator i = _verts.begin(); i != _verts.end(); ++i )
383  {
384  int vidx = ID_FROM_HANDLE( *i ) - 1;
385  HFacet hf = v2hv[vidx];
386  EntityHandle eid = fid_from_halfacet( hf, MBEDGE );
387  int lvid = lid_from_halffacet( hf );
388  std::cout << "Vertex = " << *i << " :: <" << eid << "," << lvid << ">" << std::endl;
389  }
390  }
391  else if( dim == 2 )
392  {
393  EntityType ftype = mb->type_from_handle( *_faces.begin() );
394  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
395  EntityHandle start_face = *_faces.begin();
396  std::cout << "start_face = " << start_face << std::endl;
397  std::cout << "<SIBHES_FID,SIBHES_LEID>" << std::endl;
398 
399  for( Range::iterator i = _faces.begin(); i != _faces.end(); ++i )
400  {
401  int fidx = ID_FROM_HANDLE( *i ) - 1;
402  std::cout << "Entity = " << *i;
403  for( int j = 0; j < nepf; j++ )
404  {
405  HFacet hf = sibhes[nepf * fidx + j];
406  EntityHandle sib = fid_from_halfacet( hf, ftype );
407  int lid = lid_from_halffacet( hf );
408  std::cout << " :: <" << sib << "," << lid << ">"
409  << " ";
410  }
411  std::cout << std::endl;
412  }
413 
414  std::cout << "<V2HE_FID, V2HE_LEID>" << std::endl;
415 
416  for( Range::iterator i = _verts.begin(); i != _verts.end(); ++i )
417  {
418  int vidx = ID_FROM_HANDLE( *i ) - 1;
419  HFacet hf = v2he[vidx];
420  EntityHandle fid = fid_from_halfacet( hf, ftype );
421  int lid = lid_from_halffacet( hf );
422  std::cout << "Vertex = " << *i << " :: <" << fid << "," << lid << ">" << std::endl;
423  }
424  }
425  else if( dim == 3 )
426  {
427  EntityType ctype = mb->type_from_handle( *_cells.begin() );
428 
429  int index = get_index_in_lmap( *_cells.begin() );
430  int nfpc = lConnMap3D[index].num_faces_in_cell;
431  EntityHandle start_cell = *_cells.begin();
432  std::cout << "start_cell = " << start_cell << std::endl;
433  std::cout << "<SIBHES_CID,SIBHES_LFID>" << std::endl;
434 
435  for( Range::iterator i = _cells.begin(); i != _cells.end(); ++i )
436  {
437  int cidx = ID_FROM_HANDLE( *i ) - 1;
438  std::cout << "Entity = " << *i;
439  for( int j = 0; j < nfpc; j++ )
440  {
441  HFacet hf = sibhfs[nfpc * cidx + j];
442  EntityHandle sib = fid_from_halfacet( hf, ctype );
443  int lid = lid_from_halffacet( hf );
444  std::cout << " :: <" << sib << "," << lid << ">"
445  << " ";
446  }
447  std::cout << std::endl;
448  }
449 
450  std::cout << "<V2HF_CID, V2HF_LFID>" << std::endl;
451  EntityHandle cid;
452  int lid;
453 
454  for( Range::iterator i = _verts.begin(); i != _verts.end(); ++i )
455  {
456  int vidx = ID_FROM_HANDLE( *i ) - 1;
457  HFacet hf = v2hf[vidx];
458 
459  if( hf == 0 && ( v2hfs.find( *i ) != v2hfs.end() ) )
460  {
461  std::pair< std::multimap< EntityHandle, HFacet >::iterator,
462  std::multimap< EntityHandle, HFacet >::iterator >
463  it_hfs;
464  it_hfs = v2hfs.equal_range( *i );
465 
466  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hfs.first; it != it_hfs.second; ++it )
467  {
468  cid = fid_from_halfacet( it->second, ctype );
469  lid = lid_from_halffacet( hf );
470 
471  std::cout << "Vertex = " << *i << " :: <" << cid << "," << lid << ">" << std::endl;
472  }
473  }
474  else
475  {
476  cid = fid_from_halfacet( hf, ctype );
477  lid = lid_from_halffacet( hf );
478  std::cout << "Vertex = " << *i << " :: <" << cid << "," << lid << ">" << std::endl;
479  }
480  }
481  }
482  return MB_SUCCESS;
483 }
484 
485 /**********************************************************
486  * User interface for adjacency functions *
487  ********************************************************/
488 
489 ErrorCode HalfFacetRep::get_adjacencies( const EntityHandle source_entity,
490  const unsigned int target_dimension,
491  std::vector< EntityHandle >& target_entities )
492 {
493 
494  ErrorCode error;
495 
496  unsigned int source_dimension = mb->dimension_from_handle( source_entity );
497  assert( ( source_dimension <= target_dimension ) || ( source_dimension > target_dimension ) );
498 
499  if( mInitAHFmaps == false )
500  {
501  error = initialize();MB_CHK_ERR( error );
502  }
503 
504  int mindex = get_index_for_meshtype( thismeshtype );
505  int adj_possible = adjMatrix[mindex].val[source_dimension][target_dimension];
506 
507  if( adj_possible )
508  {
509  if( source_dimension < target_dimension )
510  {
511  error = get_up_adjacencies( source_entity, target_dimension, target_entities );MB_CHK_ERR( error );
512  }
513  else if( source_dimension == target_dimension )
514  {
515  error = get_neighbor_adjacencies( source_entity, target_entities );MB_CHK_ERR( error );
516  }
517  else
518  {
519  error = get_down_adjacencies( source_entity, target_dimension, target_entities );MB_CHK_ERR( error );
520  }
521  }
522  else
523  return MB_SUCCESS;
524 
525  return MB_SUCCESS;
526 }
527 
528 ErrorCode HalfFacetRep::get_up_adjacencies( EntityHandle ent,
529  int out_dim,
530  std::vector< EntityHandle >& adjents,
531  std::vector< int >* lids )
532 {
533  ErrorCode error;
534  int in_dim = mb->dimension_from_handle( ent );
535  assert( ( in_dim >= 0 && in_dim <= 2 ) && ( out_dim > in_dim ) );
536 
537  if( in_dim == 0 )
538  {
539  if( out_dim == 1 )
540  {
541  error = get_up_adjacencies_1d( ent, adjents, lids );MB_CHK_ERR( error );
542  }
543  else if( out_dim == 2 )
544  {
545  error = get_up_adjacencies_vert_2d( ent, adjents );MB_CHK_ERR( error );
546  }
547  else if( out_dim == 3 )
548  {
549  error = get_up_adjacencies_vert_3d( ent, adjents );MB_CHK_ERR( error );
550  }
551  }
552 
553  else if( ( in_dim == 1 ) && ( out_dim == 2 ) )
554  {
555  error = get_up_adjacencies_2d( ent, adjents, lids );MB_CHK_ERR( error );
556  }
557  else if( ( in_dim == 1 ) && ( out_dim == 3 ) )
558  {
559  error = get_up_adjacencies_edg_3d( ent, adjents, lids );MB_CHK_ERR( error );
560  }
561  else if( ( in_dim == 2 ) && ( out_dim == 3 ) )
562  {
563  error = get_up_adjacencies_face_3d( ent, adjents, lids );MB_CHK_ERR( error );
564  }
565  return MB_SUCCESS;
566 }
567 
568 ErrorCode HalfFacetRep::get_neighbor_adjacencies( EntityHandle ent, std::vector< EntityHandle >& adjents )
569 {
570  ErrorCode error;
571  int in_dim = mb->dimension_from_handle( ent );
572  assert( in_dim >= 1 && in_dim <= 3 );
573 
574  if( in_dim == 1 )
575  {
576  error = get_neighbor_adjacencies_1d( ent, adjents );MB_CHK_ERR( error );
577  }
578 
579  else if( in_dim == 2 )
580  {
581  error = get_neighbor_adjacencies_2d( ent, adjents );MB_CHK_ERR( error );
582  }
583  else if( in_dim == 3 )
584  {
585  error = get_neighbor_adjacencies_3d( ent, adjents );MB_CHK_ERR( error );
586  }
587  return MB_SUCCESS;
588 }
589 
590 ErrorCode HalfFacetRep::get_down_adjacencies( EntityHandle ent, int out_dim, std::vector< EntityHandle >& adjents )
591 {
592  ErrorCode error;
593  int in_dim = mb->dimension_from_handle( ent );
594  assert( ( in_dim >= 2 && in_dim <= 3 ) && ( out_dim < in_dim ) );
595 
596  if( ( in_dim == 2 ) && ( out_dim == 1 ) )
597  {
598  error = get_down_adjacencies_2d( ent, adjents );MB_CHK_ERR( error );
599  }
600  else if( ( in_dim == 3 ) && ( out_dim == 1 ) )
601  {
602  error = get_down_adjacencies_edg_3d( ent, adjents );MB_CHK_ERR( error );
603  }
604  else if( ( in_dim == 3 ) && ( out_dim == 2 ) )
605  {
606  error = get_down_adjacencies_face_3d( ent, adjents );MB_CHK_ERR( error );
607  }
608  return MB_SUCCESS;
609 }
610 
611 ErrorCode HalfFacetRep::count_subentities( Range& edges, Range& faces, Range& cells, int* nedges, int* nfaces )
612 {
613  ErrorCode error;
614  if( edges.size() && !faces.size() && !cells.size() )
615  {
616  nedges[0] = edges.size();
617  nfaces[0] = 0;
618  }
619  else if( faces.size() && !cells.size() )
620  {
621  nedges[0] = find_total_edges_2d( faces );
622  nfaces[0] = 0;
623  }
624  else if( cells.size() )
625  {
626  error = find_total_edges_faces_3d( cells, nedges, nfaces );MB_CHK_ERR( error );
627  }
628  return MB_SUCCESS;
629 }
630 
631 /********************************************************
632  * 1D: sibhvs, v2hv, incident and neighborhood queries *
633  *********************************************************/
634 ErrorCode HalfFacetRep::determine_sibling_halfverts( Range& verts, Range& edges )
635 {
636  ErrorCode error;
637 
638  // Step 1: Create an index list storing the starting position for each vertex
639  int nv = verts.size();
640  std::vector< int > is_index( nv + 1 );
641  for( int i = 0; i < nv + 1; i++ )
642  is_index[i] = 0;
643 
644  for( Range::iterator eid = edges.begin(); eid != edges.end(); ++eid )
645  {
646  const EntityHandle* conn;
647  int num_conn = 0;
648  error = mb->get_connectivity( *eid, conn, num_conn, true );MB_CHK_ERR( error );
649 
650  int index = verts.index( conn[0] );
651  is_index[index + 1] += 1;
652  index = verts.index( conn[1] );
653  is_index[index + 1] += 1;
654  }
655  is_index[0] = 0;
656 
657  for( int i = 0; i < nv; i++ )
658  is_index[i + 1] = is_index[i] + is_index[i + 1];
659 
660  // Step 2: Define two arrays v2hv_eid, v2hv_lvid storing every half-facet on a vertex
661  std::vector< EntityHandle > v2hv_map_eid( 2 * edges.size() );
662  std::vector< int > v2hv_map_lvid( 2 * edges.size() );
663 
664  for( Range::iterator eid = edges.begin(); eid != edges.end(); ++eid )
665  {
666  const EntityHandle* conn;
667  int num_conn = 0;
668  error = mb->get_connectivity( *eid, conn, num_conn, true );MB_CHK_ERR( error );
669 
670  for( int j = 0; j < 2; j++ )
671  {
672  int v = verts.index( conn[j] );
673  v2hv_map_eid[is_index[v]] = *eid;
674  v2hv_map_lvid[is_index[v]] = j;
675  is_index[v] += 1;
676  }
677  }
678 
679  for( int i = nv - 2; i >= 0; i-- )
680  is_index[i + 1] = is_index[i];
681  is_index[0] = 0;
682 
683  // Step 3: Fill up sibling half-verts map
684  for( Range::iterator vid = verts.begin(); vid != verts.end(); ++vid )
685  {
686  int v = verts.index( *vid );
687  int last = is_index[v + 1] - 1;
688  if( last > is_index[v] )
689  {
690  EntityHandle prev_eid = v2hv_map_eid[last];
691  int prev_lvid = v2hv_map_lvid[last];
692 
693  for( int i = is_index[v]; i <= last; i++ )
694  {
695  EntityHandle cur_eid = v2hv_map_eid[i];
696  int cur_lvid = v2hv_map_lvid[i];
697 
698  int pidx = ID_FROM_HANDLE( prev_eid ) - 1;
699  sibhvs[2 * pidx + prev_lvid] = create_halffacet( cur_eid, cur_lvid );
700 
701  prev_eid = cur_eid;
702  prev_lvid = cur_lvid;
703  }
704  }
705  }
706 
707  return MB_SUCCESS;
708 }
709 
710 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
711 ErrorCode HalfFacetRep::determine_incident_halfverts( Range& edges )
712 {
713  ErrorCode error;
714 
715  for( Range::iterator e_it = edges.begin(); e_it != edges.end(); ++e_it )
716  {
717  EntityHandle cur_eid = *e_it;
718  const EntityHandle* conn;
719  int num_conn = 0;
720  error = mb->get_connectivity( *e_it, conn, num_conn, true );MB_CHK_ERR( error );
721 
722  for( int i = 0; i < 2; ++i )
723  {
724  EntityHandle v = conn[i];
725  int vidx = ID_FROM_HANDLE( v ) - 1;
726 
727  HFacet hf = v2hv[vidx];
728  EntityHandle eid = fid_from_halfacet( hf, MBEDGE );
729  if( eid == 0 )
730  {
731  v2hv[vidx] = create_halffacet( cur_eid, i );
732  }
733  }
734  }
735 
736  return MB_SUCCESS;
737 }
738 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
739 ErrorCode HalfFacetRep::get_up_adjacencies_1d( EntityHandle vid,
740  std::vector< EntityHandle >& adjents,
741  std::vector< int >* lvids )
742 {
743  adjents.clear();
744  adjents.reserve( 20 );
745 
746  if( lvids != NULL ) lvids->reserve( 20 );
747 
748  int vidx = ID_FROM_HANDLE( vid ) - 1;
749  HFacet hf = v2hv[vidx];
750 
751  EntityHandle start_eid = fid_from_halfacet( hf, MBEDGE );
752  int start_lid = lid_from_halffacet( hf );
753 
754  EntityHandle eid = start_eid;
755  int lid = start_lid;
756 
757  if( eid != 0 )
758  {
759  adjents.push_back( eid );
760  if( lvids != NULL ) lvids->push_back( lid );
761 
762  while( eid != 0 )
763  {
764  int eidx = ID_FROM_HANDLE( eid ) - 1;
765  HFacet shf = sibhvs[2 * eidx + lid];
766  eid = fid_from_halfacet( shf, MBEDGE );
767  lid = lid_from_halffacet( shf );
768 
769  if( ( !eid ) || ( eid == start_eid ) ) break;
770 
771  adjents.push_back( eid );
772  if( lvids != NULL ) lvids->push_back( lid );
773  }
774  }
775 
776  return MB_SUCCESS;
777 }
778 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
779 ErrorCode HalfFacetRep::get_neighbor_adjacencies_1d( EntityHandle eid, std::vector< EntityHandle >& adjents )
780 {
781  adjents.clear();
782  adjents.reserve( 20 );
783 
784  EntityHandle sibhv_eid;
785  int sibhv_lid;
786  int eidx = ID_FROM_HANDLE( eid ) - 1;
787 
788  for( int lid = 0; lid < 2; ++lid )
789  {
790  HFacet shf = sibhvs[2 * eidx + lid];
791  sibhv_eid = fid_from_halfacet( shf, MBEDGE );
792  sibhv_lid = lid_from_halffacet( shf );
793 
794  if( sibhv_eid != 0 )
795  {
796  adjents.push_back( sibhv_eid );
797 
798  eidx = ID_FROM_HANDLE( sibhv_eid ) - 1;
799  HFacet nhf = sibhvs[2 * eidx + sibhv_lid];
800  EntityHandle hv_eid = fid_from_halfacet( nhf, MBEDGE );
801  int hv_lid = lid_from_halffacet( nhf );
802 
803  while( hv_eid != 0 )
804  {
805  if( hv_eid != eid ) adjents.push_back( hv_eid );
806 
807  eidx = ID_FROM_HANDLE( hv_eid ) - 1;
808  HFacet hf = sibhvs[2 * eidx + hv_lid];
809  EntityHandle edge = fid_from_halfacet( hf, MBEDGE );
810  if( edge == sibhv_eid ) break;
811 
812  hv_eid = edge;
813  hv_lid = lid_from_halffacet( hf );
814  }
815  }
816  }
817 
818  return MB_SUCCESS;
819 }
820 
821 /*******************************************************
822  * 2D: sibhes, v2he, incident and neighborhood queries *
823  ********************************************************/
824 const HalfFacetRep::LocalMaps2D HalfFacetRep::lConnMap2D[2] = {
825  // Triangle
826  { 3, { 1, 2, 0 }, { 2, 0, 1 } },
827  // Quad
828  { 4, { 1, 2, 3, 0 }, { 3, 0, 1, 2 } } };
829 
830 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
831 ErrorCode HalfFacetRep::determine_sibling_halfedges( Range& faces )
832 {
833  ErrorCode error;
834  EntityHandle start_face = *faces.begin();
835  EntityType ftype = mb->type_from_handle( start_face );
836  int nfaces = faces.size();
837  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
838 
839  // Step 1: Create an index list storing the starting position for each vertex
840  int nv = _verts.size();
841  std::vector< int > is_index( nv + 1 );
842  for( int i = 0; i < nv + 1; i++ )
843  is_index[i] = 0;
844 
845  int index;
846 
847  for( Range::iterator fid = faces.begin(); fid != faces.end(); ++fid )
848  {
849  const EntityHandle* conn;
850  error = mb->get_connectivity( *fid, conn, nepf, true );MB_CHK_ERR( error );
851 
852  for( int i = 0; i < nepf; i++ )
853  {
854  index = _verts.index( conn[i] );
855  is_index[index + 1] += 1;
856  }
857  }
858  is_index[0] = 0;
859 
860  for( int i = 0; i < nv; i++ )
861  is_index[i + 1] = is_index[i] + is_index[i + 1];
862 
863  // Step 2: Define two arrays v2hv_eid, v2hv_lvid storing every half-facet on a vertex
864  std::vector< EntityHandle > v2nv( nepf * nfaces );
865  std::vector< EntityHandle > v2he_map_fid( nepf * nfaces );
866  std::vector< int > v2he_map_leid( nepf * nfaces );
867 
868  for( Range::iterator fid = faces.begin(); fid != faces.end(); ++fid )
869  {
870  const EntityHandle* conn;
871  error = mb->get_connectivity( *fid, conn, nepf, true );MB_CHK_ERR( error );
872 
873  for( int j = 0; j < nepf; j++ )
874  {
875  int v = _verts.index( conn[j] );
876  int nidx = lConnMap2D[ftype - 2].next[j];
877  v2nv[is_index[v]] = conn[nidx];
878  v2he_map_fid[is_index[v]] = *fid;
879  v2he_map_leid[is_index[v]] = j;
880  is_index[v] += 1;
881  }
882  }
883 
884  for( int i = nv - 2; i >= 0; i-- )
885  is_index[i + 1] = is_index[i];
886  is_index[0] = 0;
887 
888  // Step 3: Fill up sibling half-verts map
889  for( Range::iterator fid = faces.begin(); fid != faces.end(); ++fid )
890  {
891  const EntityHandle* conn;
892  error = mb->get_connectivity( *fid, conn, nepf, true );MB_CHK_ERR( error );
893 
894  int fidx = ID_FROM_HANDLE( *fid ) - 1;
895  for( int k = 0; k < nepf; k++ )
896  {
897  HFacet hf = sibhes[nepf * fidx + k];
898  EntityHandle sibfid = fid_from_halfacet( hf, ftype );
899 
900  if( sibfid != 0 ) continue;
901 
902  int nidx = lConnMap2D[ftype - 2].next[k];
903  int v = _verts.index( conn[k] );
904  int vn = _verts.index( conn[nidx] );
905 
906  EntityHandle first_fid = *fid;
907  int first_leid = k;
908 
909  EntityHandle prev_fid = *fid;
910  int prev_leid = k;
911 
912  for( index = is_index[vn]; index <= is_index[vn + 1] - 1; index++ )
913  {
914  if( v2nv[index] == conn[k] )
915  {
916  EntityHandle cur_fid = v2he_map_fid[index];
917  int cur_leid = v2he_map_leid[index];
918 
919  int pidx = ID_FROM_HANDLE( prev_fid ) - 1;
920  sibhes[nepf * pidx + prev_leid] = create_halffacet( cur_fid, cur_leid );
921 
922  prev_fid = cur_fid;
923  prev_leid = cur_leid;
924  }
925  }
926 
927  for( index = is_index[v]; index <= is_index[v + 1] - 1; index++ )
928  {
929  if( ( v2nv[index] == conn[nidx] ) && ( v2he_map_fid[index] != *fid ) )
930  {
931 
932  EntityHandle cur_fid = v2he_map_fid[index];
933  int cur_leid = v2he_map_leid[index];
934 
935  int pidx = ID_FROM_HANDLE( prev_fid ) - 1;
936  sibhes[nepf * pidx + prev_leid] = create_halffacet( cur_fid, cur_leid );
937 
938  prev_fid = cur_fid;
939  prev_leid = cur_leid;
940  }
941  }
942 
943  if( prev_fid != first_fid )
944  {
945 
946  int pidx = ID_FROM_HANDLE( prev_fid ) - 1;
947  sibhes[nepf * pidx + prev_leid] = create_halffacet( first_fid, first_leid );
948  }
949  }
950  }
951 
952  return MB_SUCCESS;
953 }
954 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
955 ErrorCode HalfFacetRep::determine_incident_halfedges( Range& faces )
956 {
957  ErrorCode error;
958  EntityType ftype = mb->type_from_handle( *faces.begin() );
959  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
960 
961  std::vector< char > markEdges( nepf * faces.size(), 0 );
962 
963  for( Range::iterator it = faces.begin(); it != faces.end(); ++it )
964  {
965  EntityHandle fid = *it;
966  const EntityHandle* conn;
967  error = mb->get_connectivity( fid, conn, nepf, true );MB_CHK_ERR( error );
968 
969  for( int i = 0; i < nepf; ++i )
970  {
971  EntityHandle v = conn[i];
972  int vidx = ID_FROM_HANDLE( v ) - 1;
973  HFacet hf = v2he[vidx];
974 
975  if( hf == 0 && ( v2hes.empty() || ( v2hes.find( v ) == v2hes.end() ) ) )
976  {
977  // This is the first time a half-facet is assigned to a vertex.
978  HFacet nwhf = 0;
979  error = mark_halfedges( v, fid, i, faces, markEdges, nwhf );MB_CHK_ERR( error );
980 
981  if( nwhf == 0 ) nwhf = create_halffacet( fid, i );
982 
983  v2he[vidx] = nwhf;
984  }
985  else if( hf != 0 && !markEdges[nepf * faces.index( fid ) + i] )
986  {
987  // This is the first time a non-manifold vertex is encountered. Copy the existing he
988  // in v2he[v] to the multimap.
989  v2hes.insert( std::pair< EntityHandle, HFacet >( v, hf ) );
990  HFacet nwhf = 0;
991  error = mark_halfedges( v, fid, i, faces, markEdges, nwhf );MB_CHK_ERR( error );
992 
993  if( nwhf == 0 ) nwhf = create_halffacet( fid, i );
994 
995  v2hes.insert( std::pair< EntityHandle, HFacet >( v, nwhf ) );
996  v2he[vidx] = 0;
997  }
998  else if( hf == 0 && ( !v2hes.empty() ) && ( v2hes.find( v ) != v2hes.end() ) &&
999  !markEdges[nepf * faces.index( fid ) + i] )
1000  {
1001  // This is check if reached if the vertex is non-manifold and has encountered a
1002  // half-facet to a new component.
1003  HFacet nwhf = 0;
1004  error = mark_halfedges( v, fid, i, faces, markEdges, nwhf );MB_CHK_ERR( error );
1005 
1006  if( nwhf == 0 ) nwhf = create_halffacet( fid, i );
1007 
1008  v2hes.insert( std::pair< EntityHandle, HFacet >( v, nwhf ) );
1009  }
1010  }
1011  }
1012 
1013  // error = print_tags(2);
1014 
1015  return MB_SUCCESS;
1016 }
1017 
1018 ///////////////////////////////////////////////////////////////////
1019 ErrorCode HalfFacetRep::mark_halfedges( EntityHandle vid,
1020  EntityHandle he_fid,
1021  int he_lid,
1022  Range& faces,
1023  std::vector< char >& markHEdgs,
1024  HFacet& bnd_hf )
1025 {
1026  ErrorCode error;
1027  EntityType ftype = mb->type_from_handle( he_fid );
1028  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
1029 
1030  int qsize = 0, count = -1;
1031  int num_qvals = 0;
1032 
1033  error = gather_halfedges( vid, he_fid, he_lid, &qsize, &count );MB_CHK_ERR( error );
1034 
1035  while( num_qvals < qsize )
1036  {
1037  EntityHandle curfid = queue_fid[num_qvals];
1038  int curlid = queue_lid[num_qvals];
1039  num_qvals += 1;
1040 
1041  int fidx = ID_FROM_HANDLE( curfid ) - 1;
1042 
1043  const EntityHandle* conn;
1044  error = mb->get_connectivity( curfid, conn, nepf, true );MB_CHK_ERR( error );
1045 
1046  if( !markHEdgs[nepf * faces.index( curfid ) + curlid] && ( conn[curlid] == vid ) )
1047  {
1048  markHEdgs[nepf * faces.index( curfid ) + curlid] = 1;
1049  HFacet hf = sibhes[nepf * fidx + curlid];
1050  EntityHandle sibfid = fid_from_halfacet( hf, ftype );
1051  if( sibfid == 0 ) bnd_hf = create_halffacet( curfid, curlid );
1052  }
1053 
1054  EntityHandle he2_fid = 0;
1055  int he2_lid = 0;
1056  error = another_halfedge( vid, curfid, curlid, &he2_fid, &he2_lid );MB_CHK_ERR( error );
1057 
1058  if( !markHEdgs[nepf * faces.index( curfid ) + he2_lid] && ( conn[he2_lid] == vid ) )
1059  {
1060  markHEdgs[nepf * faces.index( curfid ) + he2_lid] = 1;
1061  HFacet hf = sibhes[nepf * fidx + he2_lid];
1062  EntityHandle sibfid = fid_from_halfacet( hf, ftype );
1063  if( sibfid == 0 ) bnd_hf = create_halffacet( he2_fid, he2_lid );
1064  }
1065 
1066  bool val = find_match_in_array( he2_fid, trackfaces, count );
1067 
1068  if( val ) continue;
1069 
1070  count += 1;
1071  trackfaces[count] = he2_fid;
1072 
1073  error = get_up_adjacencies_2d( he2_fid, he2_lid, &qsize, &count );MB_CHK_ERR( error );
1074  }
1075 
1076  // Change the visited faces to false, also empty the queue
1077  for( int i = 0; i <= qsize; i++ )
1078  {
1079  queue_fid[i] = 0;
1080  queue_lid[i] = 0;
1081  }
1082 
1083  for( int i = 0; i <= count; i++ )
1084  trackfaces[i] = 0;
1085  return MB_SUCCESS;
1086 }
1087 
1088 ///////////////////////////////////////////////////////////////////
1089 ErrorCode HalfFacetRep::get_up_adjacencies_vert_2d( EntityHandle vid, std::vector< EntityHandle >& adjents )
1090 {
1091  ErrorCode error;
1092  EntityType ftype = mb->type_from_handle( *_faces.begin() );
1093 
1094  int vidx = ID_FROM_HANDLE( vid ) - 1;
1095  HFacet hf = v2he[vidx];
1096 
1097  std::vector< EntityHandle > start_fids;
1098  std::vector< int > start_lids;
1099 
1100  if( hf == 0 && ( v2hes.find( vid ) != v2hes.end() ) )
1101  {
1102  std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
1103  it_hes;
1104  it_hes = v2hes.equal_range( vid );
1105 
1106  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
1107  {
1108  start_fids.push_back( fid_from_halfacet( it->second, ftype ) );
1109  start_lids.push_back( lid_from_halffacet( it->second ) );
1110  }
1111  }
1112  else if( hf != 0 )
1113  {
1114  start_fids.push_back( fid_from_halfacet( hf, ftype ) );
1115  start_lids.push_back( lid_from_halffacet( hf ) );
1116  }
1117 
1118  if( start_fids.empty() ) return MB_SUCCESS;
1119 
1120  int qsize = 0, count = -1;
1121  int num_qvals = 0;
1122 
1123  adjents.reserve( (int)start_fids.size() );
1124 
1125  for( int i = 0; i < (int)start_fids.size(); i++ )
1126  {
1127  adjents.push_back( start_fids[i] );
1128  error = gather_halfedges( vid, start_fids[i], start_lids[i], &qsize, &count );MB_CHK_ERR( error );
1129  }
1130 
1131  while( num_qvals < qsize )
1132  {
1133  EntityHandle curfid = queue_fid[num_qvals];
1134  int curlid = queue_lid[num_qvals];
1135  num_qvals += 1;
1136 
1137  EntityHandle he2_fid = 0;
1138  int he2_lid = 0;
1139  error = another_halfedge( vid, curfid, curlid, &he2_fid, &he2_lid );MB_CHK_ERR( error );
1140 
1141  bool val = find_match_in_array( he2_fid, trackfaces, count );
1142 
1143  if( val ) continue;
1144 
1145  count += 1;
1146  trackfaces[count] = he2_fid;
1147 
1148  error = get_up_adjacencies_2d( he2_fid, he2_lid, &qsize, &count );MB_CHK_ERR( error );
1149 
1150  adjents.push_back( he2_fid );
1151  }
1152 
1153  // Change the visited faces to false, also empty the queue
1154  for( int i = 0; i <= qsize; i++ )
1155  {
1156  queue_fid[i] = 0;
1157  queue_lid[i] = 0;
1158  }
1159 
1160  for( int i = 0; i <= count; i++ )
1161  trackfaces[i] = 0;
1162 
1163  return MB_SUCCESS;
1164 }
1165 
1166 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1167 ErrorCode HalfFacetRep::get_up_adjacencies_2d( EntityHandle eid,
1168  std::vector< EntityHandle >& adjents,
1169  std::vector< int >* leids )
1170 {
1171 
1172  // Given an explicit edge eid, find the incident faces.
1173  ErrorCode error;
1174  EntityHandle he_fid = 0;
1175  int he_lid = 0;
1176 
1177  // Step 1: Given an explicit edge, find a corresponding half-edge from the surface mesh.
1178  bool found = find_matching_halfedge( eid, &he_fid, &he_lid );
1179 
1180  // Step 2: If there is a corresponding half-edge, collect all sibling half-edges and store the
1181  // incident faces.
1182  if( found )
1183  {
1184  error = get_up_adjacencies_2d( he_fid, he_lid, true, adjents, leids );MB_CHK_ERR( error );
1185  }
1186 
1187  return MB_SUCCESS;
1188 }
1189 
1190 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1191 ErrorCode HalfFacetRep::get_up_adjacencies_2d( EntityHandle fid,
1192  int leid,
1193  bool add_inent,
1194  std::vector< EntityHandle >& adj_ents,
1195  std::vector< int >* adj_leids,
1196  std::vector< int >* adj_orients )
1197 {
1198  // Given an implicit half-edge <fid, leid>, find the incident half-edges.
1199  ErrorCode error;
1200 
1201  EntityType ftype = mb->type_from_handle( fid );
1202  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
1203 
1204  if( !fid ) return MB_FAILURE;
1205  adj_ents.reserve( 20 );
1206 
1207  bool local_id = false;
1208  bool orient = false;
1209  if( adj_leids != NULL )
1210  {
1211  local_id = true;
1212  adj_leids->reserve( 20 );
1213  }
1214  if( adj_orients != NULL )
1215  {
1216  orient = true;
1217  adj_orients->reserve( 20 );
1218  }
1219 
1220  if( add_inent )
1221  {
1222  adj_ents.push_back( fid );
1223  if( local_id ) adj_leids->push_back( leid );
1224  }
1225 
1226  EntityHandle fedge[2] = { 0, 0 };
1227 
1228  if( orient )
1229  {
1230  // get connectivity and match their directions
1231  const EntityHandle* fid_conn;
1232  error = mb->get_connectivity( fid, fid_conn, nepf, true );MB_CHK_ERR( error );
1233 
1234  int nidx = lConnMap2D[ftype - 2].next[leid];
1235  fedge[0] = fid_conn[leid];
1236  fedge[1] = fid_conn[nidx];
1237  }
1238 
1239  int fidx = ID_FROM_HANDLE( fid ) - 1;
1240  HFacet hf = sibhes[nepf * fidx + leid];
1241  EntityHandle curfid = fid_from_halfacet( hf, ftype );
1242  int curlid = lid_from_halffacet( hf );
1243 
1244  while( ( curfid != fid ) && ( curfid != 0 ) )
1245  { // Should not go into the loop when no sibling exists
1246  adj_ents.push_back( curfid );
1247 
1248  if( local_id ) adj_leids->push_back( curlid );
1249 
1250  if( orient )
1251  {
1252  // get connectivity and match their directions
1253  const EntityHandle* conn;
1254  error = mb->get_connectivity( curfid, conn, nepf, true );MB_CHK_ERR( error );
1255 
1256  int nidx = lConnMap2D[ftype - 2].next[curlid];
1257 
1258  if( ( fedge[0] == conn[curlid] ) && ( fedge[1] == conn[nidx] ) )
1259  adj_orients->push_back( 1 );
1260  else if( ( fedge[1] == conn[curlid] ) && ( fedge[0] == conn[nidx] ) )
1261  adj_orients->push_back( 0 );
1262  }
1263 
1264  int cidx = ID_FROM_HANDLE( curfid ) - 1;
1265  hf = sibhes[nepf * cidx + curlid];
1266  curfid = fid_from_halfacet( hf, ftype );
1267  curlid = lid_from_halffacet( hf );
1268  }
1269 
1270  return MB_SUCCESS;
1271 }
1272 
1273 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1274 ErrorCode HalfFacetRep::get_up_adjacencies_2d( EntityHandle fid, int lid, int* qsize, int* count )
1275 {
1276  EntityType ftype = mb->type_from_handle( fid );
1277  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
1278 
1279  int fidx = ID_FROM_HANDLE( fid ) - 1;
1280  HFacet hf = sibhes[nepf * fidx + lid];
1281  EntityHandle curfid = fid_from_halfacet( hf, ftype );
1282  int curlid = lid_from_halffacet( hf );
1283 
1284  if( curfid == 0 )
1285  {
1286  int index = 0;
1287  bool found_ent = find_match_in_array( fid, queue_fid, qsize[0] - 1, true, &index );
1288 
1289  if( ( !found_ent ) || ( ( found_ent ) && ( queue_lid[index] != lid ) ) )
1290  {
1291  queue_fid[qsize[0]] = fid;
1292  queue_lid[qsize[0]] = lid;
1293  qsize[0] += 1;
1294  }
1295  }
1296 
1297  while( ( curfid != fid ) && ( curfid != 0 ) )
1298  {
1299  bool val = find_match_in_array( curfid, trackfaces, count[0] );
1300 
1301  if( !val )
1302  {
1303  queue_fid[qsize[0]] = curfid;
1304  queue_lid[qsize[0]] = curlid;
1305  qsize[0] += 1;
1306  }
1307 
1308  int cidx = ID_FROM_HANDLE( curfid ) - 1;
1309  hf = sibhes[nepf * cidx + curlid];
1310  curfid = fid_from_halfacet( hf, ftype );
1311  curlid = lid_from_halffacet( hf );
1312  }
1313 
1314  return MB_SUCCESS;
1315 }
1316 
1317 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1318 bool HalfFacetRep::find_matching_halfedge( EntityHandle eid, EntityHandle* hefid, int* helid )
1319 {
1320  ErrorCode error;
1321  EntityType ftype = mb->type_from_handle( *_faces.begin() );
1322 
1323  const EntityHandle* conn;
1324  int num_conn = 0;
1325  error = mb->get_connectivity( eid, conn, num_conn, true );MB_CHK_ERR( error );
1326 
1327  EntityHandle vid = conn[0];
1328  int vidx = ID_FROM_HANDLE( conn[0] ) - 1;
1329  HFacet hf = v2he[vidx];
1330 
1331  if( hf == 0 )
1332  {
1333  vidx = ID_FROM_HANDLE( conn[1] ) - 1;
1334  hf = v2he[vidx];
1335 
1336  if( hf == 0 ) // The edge is either a dangling edge or attached to two non-manifold
1337  // vertices
1338  return MB_SUCCESS;
1339 
1340  vid = conn[1];
1341  }
1342 
1343  EntityHandle fid = fid_from_halfacet( hf, ftype );
1344  int lid = lid_from_halffacet( hf );
1345 
1346  bool found = false;
1347  int qsize = 0, count = -1;
1348 
1349  error = gather_halfedges( vid, fid, lid, &qsize, &count );MB_CHK_ERR( error );
1350 
1351  found = collect_and_compare( vid, conn, &qsize, &count, hefid, helid );MB_CHK_ERR( error );
1352 
1353  // Change the visited faces to false
1354  for( int i = 0; i < qsize; i++ )
1355  {
1356  queue_fid[i] = 0;
1357  queue_lid[i] = 0;
1358  }
1359 
1360  for( int i = 0; i <= count; i++ )
1361  trackfaces[i] = 0;
1362 
1363  return found;
1364 }
1365 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1366 ErrorCode HalfFacetRep::gather_halfedges( EntityHandle vid, EntityHandle he_fid, int he_lid, int* qsize, int* count )
1367 {
1368  ErrorCode error;
1369  EntityHandle he2_fid = 0;
1370  int he2_lid = 0;
1371 
1372  error = another_halfedge( vid, he_fid, he_lid, &he2_fid, &he2_lid );MB_CHK_ERR( error );
1373 
1374  queue_fid[*qsize] = he_fid;
1375  queue_lid[*qsize] = he_lid;
1376  *qsize += 1;
1377 
1378  queue_fid[*qsize] = he2_fid;
1379  queue_lid[*qsize] = he2_lid;
1380  *qsize += 1;
1381 
1382  *count += 1;
1383  trackfaces[*count] = he_fid;
1384 
1385  error = get_up_adjacencies_2d( he_fid, he_lid, qsize, count );MB_CHK_ERR( error );
1386  error = get_up_adjacencies_2d( he2_fid, he2_lid, qsize, count );MB_CHK_ERR( error );
1387 
1388  return MB_SUCCESS;
1389 }
1390 
1391 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1392 ErrorCode HalfFacetRep::another_halfedge( EntityHandle vid,
1393  EntityHandle he_fid,
1394  int he_lid,
1395  EntityHandle* he2_fid,
1396  int* he2_lid )
1397 {
1398  ErrorCode error;
1399  EntityType ftype = mb->type_from_handle( he_fid );
1400  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
1401 
1402  const EntityHandle* conn;
1403  error = mb->get_connectivity( he_fid, conn, nepf, true );MB_CHK_ERR( error );
1404 
1405  *he2_fid = he_fid;
1406  if( conn[he_lid] == vid )
1407  *he2_lid = lConnMap2D[ftype - 2].prev[he_lid];
1408  else
1409  *he2_lid = lConnMap2D[ftype - 2].next[he_lid];
1410 
1411  return MB_SUCCESS;
1412 }
1413 
1414 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1415 bool HalfFacetRep::collect_and_compare( const EntityHandle vid,
1416  const EntityHandle* edg_vert,
1417  int* qsize,
1418  int* count,
1419  EntityHandle* he_fid,
1420  int* he_lid )
1421 {
1422  ErrorCode error;
1423  EntityType ftype = mb->type_from_handle( *_faces.begin() );
1424  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
1425 
1426  bool found = false;
1427  int num_qvals = 0, counter = 0;
1428 
1429  while( num_qvals < *qsize && counter < MAXSIZE )
1430  {
1431  EntityHandle curfid = queue_fid[num_qvals];
1432  int curlid = queue_lid[num_qvals];
1433  num_qvals += 1;
1434 
1435  const EntityHandle* conn;
1436  error = mb->get_connectivity( curfid, conn, nepf, true );MB_CHK_ERR( error );
1437 
1438  int id = lConnMap2D[ftype - 2].next[curlid];
1439  if( ( ( conn[curlid] == edg_vert[0] ) && ( conn[id] == edg_vert[1] ) ) ||
1440  ( ( conn[curlid] == edg_vert[1] ) && ( conn[id] == edg_vert[0] ) ) )
1441  {
1442  *he_fid = curfid;
1443  *he_lid = curlid;
1444  found = true;
1445  break;
1446  }
1447 
1448  bool val = find_match_in_array( curfid, trackfaces, count[0] );
1449 
1450  if( val ) continue;
1451 
1452  count[0] += 1;
1453  trackfaces[*count] = curfid;
1454 
1455  EntityHandle he2_fid;
1456  int he2_lid;
1457  error = another_halfedge( vid, curfid, curlid, &he2_fid, &he2_lid );MB_CHK_ERR( error );
1458  error = get_up_adjacencies_2d( he2_fid, he2_lid, qsize, count );MB_CHK_ERR( error );
1459 
1460  counter += 1;
1461  }
1462  return found;
1463 }
1464 
1465 ///////////////////////////////////////////////////////////////////////////////////////////////////
1466 ErrorCode HalfFacetRep::get_neighbor_adjacencies_2d( EntityHandle fid, std::vector< EntityHandle >& adjents )
1467 {
1468  ErrorCode error;
1469 
1470  if( fid != 0 )
1471  {
1472  EntityType ftype = mb->type_from_handle( fid );
1473  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
1474 
1475  for( int lid = 0; lid < nepf; ++lid )
1476  {
1477  error = get_up_adjacencies_2d( fid, lid, false, adjents );MB_CHK_ERR( error );
1478  }
1479  }
1480 
1481  return MB_SUCCESS;
1482 }
1483 
1484 /////////////////////////////////////////////////////////////////////////////////////////////////
1485 ErrorCode HalfFacetRep::get_down_adjacencies_2d( EntityHandle fid, std::vector< EntityHandle >& adjents )
1486 {
1487  // Returns explicit edges, if any, of the face
1488  ErrorCode error;
1489  adjents.reserve( 10 );
1490  EntityType ftype = mb->type_from_handle( fid );
1491  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
1492 
1493  const EntityHandle* conn;
1494  error = mb->get_connectivity( fid, conn, nepf, true );MB_CHK_ERR( error );
1495 
1496  std::vector< EntityHandle > temp;
1497 
1498  // Loop over 2 vertices
1499  for( int i = 0; i < 2; i++ )
1500  {
1501  // Choose the two adjacent vertices for a triangle, and two opposite vertices for a quad
1502  int l;
1503  if( ftype == MBTRI )
1504  l = i;
1505  else
1506  l = 2 * i;
1507 
1508  // Get the current, next and prev vertices
1509  int nidx = lConnMap2D[ftype - 2].next[l];
1510  int pidx = lConnMap2D[ftype - 2].prev[l];
1511  EntityHandle v = conn[l];
1512  EntityHandle vnext = conn[nidx];
1513  EntityHandle vprev = conn[pidx];
1514 
1515  // Get incident edges on v
1516  error = get_up_adjacencies_1d( v, temp );MB_CHK_ERR( error );
1517 
1518  // Loop over the incident edges and check if its end vertices match those in the face
1519  for( int k = 0; k < (int)temp.size(); k++ )
1520  {
1521  const EntityHandle* econn;
1522  int num_conn = 0;
1523  error = mb->get_connectivity( temp[k], econn, num_conn, true );MB_CHK_ERR( error );
1524 
1525  if( ( econn[0] == v && econn[1] == vnext ) || ( econn[0] == v && econn[1] == vprev ) ||
1526  ( econn[0] == vnext && econn[1] == v ) || ( econn[0] == vprev && econn[1] == v ) )
1527  {
1528  bool found = false;
1529  for( int j = 0; j < (int)adjents.size(); j++ )
1530  {
1531  if( adjents[j] == temp[k] )
1532  {
1533  found = true;
1534  break;
1535  }
1536  }
1537  if( !found ) adjents.push_back( temp[k] );
1538  }
1539  }
1540  }
1541 
1542  return MB_SUCCESS;
1543 }
1544 
1545 ////////////////////////////////////////////////////////////////////////////////////////////////
1546 int HalfFacetRep::find_total_edges_2d( Range& faces )
1547 {
1548  ErrorCode error;
1549  EntityType ftype = mb->type_from_handle( *faces.begin() );
1550  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
1551  int nfaces = faces.size();
1552 
1553  int total_edges = nepf * nfaces;
1554 
1555  std::vector< int > trackF( total_edges, 0 );
1556  std::vector< EntityHandle > adj_fids;
1557  std::vector< int > adj_lids;
1558 
1559  for( Range::iterator f = faces.begin(); f != faces.end(); ++f )
1560  {
1561  for( int l = 0; l < nepf; l++ )
1562  {
1563 
1564  adj_fids.clear();
1565  adj_lids.clear();
1566 
1567  int id = nepf * ( faces.index( *f ) ) + l;
1568  if( !trackF[id] )
1569  {
1570  error = get_up_adjacencies_2d( *f, l, false, adj_fids, &adj_lids );MB_CHK_ERR( error );
1571 
1572  total_edges -= adj_fids.size();
1573 
1574  for( int i = 0; i < (int)adj_fids.size(); i++ )
1575  trackF[nepf * ( faces.index( adj_fids[i] ) ) + adj_lids[i]] = 1;
1576  };
1577  };
1578  };
1579 
1580  return total_edges;
1581 }
1582 
1583 ErrorCode HalfFacetRep::get_face_edges( EntityHandle fid, std::vector< EntityHandle >& edges )
1584 {
1585  ErrorCode error;
1586  edges.clear();
1587 
1588  EntityType ftype = mb->type_from_handle( fid );
1589  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
1590 
1591  std::vector< EntityHandle > conn;
1592  error = mb->get_connectivity( &fid, 1, conn );MB_CHK_ERR( error );
1593 
1594  for( int i = 0; i < nepf; i++ )
1595  {
1596  EntityHandle v0 = conn[i];
1597  EntityHandle v1 = conn[lConnMap2D[ftype - 2].next[i]];
1598 
1599  std::vector< EntityHandle > e0, e1, ecom;
1600  error = get_up_adjacencies_1d( v0, e0 );MB_CHK_ERR( error );
1601  error = get_up_adjacencies_1d( v1, e1 );MB_CHK_ERR( error );
1602 
1603  std::sort( e0.begin(), e0.end() );
1604  std::sort( e1.begin(), e1.end() );
1605  std::set_intersection( e0.begin(), e0.end(), e1.begin(), e1.end(), std::back_inserter( ecom ) );
1606  assert( ecom.size() == 1 || ecom.size() == 0 );
1607  if( ecom.size() == 0 )
1608  edges.push_back( 0 );
1609  else
1610  edges.push_back( ecom[0] );
1611  }
1612 
1613  return MB_SUCCESS;
1614 }
1615 
1616 /*******************************************************
1617  * 3D: sibhfs, v2hf, incident and neighborhood queries *
1618  ********************************************************/
1619 
1620 int HalfFacetRep::get_index_in_lmap( EntityHandle cid )
1621 {
1622  EntityType type = mb->type_from_handle( cid );
1623  int index = cell_index.find( type )->second;
1624  return index;
1625 }
1626 
1627 const HalfFacetRep::LocalMaps3D HalfFacetRep::lConnMap3D[4] = {
1628  // Tet
1629  { 4,
1630  6,
1631  4,
1632  { 3, 3, 3, 3 },
1633  { { 0, 1, 3 }, { 1, 2, 3 }, { 2, 0, 3 }, { 0, 2, 1 } },
1634  { 3, 3, 3, 3 },
1635  { { 0, 2, 3 }, { 0, 1, 3 }, { 1, 2, 3 }, { 0, 1, 2 } },
1636  { { 0, 1 }, { 1, 2 }, { 2, 0 }, { 0, 3 }, { 1, 3 }, { 2, 3 } },
1637  { { 3, 0 }, { 3, 1 }, { 3, 2 }, { 0, 2 }, { 0, 1 }, { 1, 2 } },
1638  { { 0, 4, 3 }, { 1, 5, 4 }, { 2, 3, 5 }, { 2, 1, 0 } },
1639  { { -1, 0, 2, 3 }, { 0, -1, 1, 4 }, { 2, 1, -1, 5 }, { 3, 4, 5, -1 } },
1640  { 3, 0, 1, 2 },
1641  { 0, 1 },
1642  { { 3, 1, 2, 3 }, { 2, 2, 3 }, { 1, 3 } } },
1643 
1644  // Pyramid: Note: In MOAB pyramid follows the CGNS convention. Look up src/MBCNArrays.hpp
1645  { 5,
1646  8,
1647  5,
1648  { 4, 3, 3, 3, 3 },
1649  { { 0, 3, 2, 1 }, { 0, 1, 4 }, { 1, 2, 4 }, { 2, 3, 4 }, { 3, 0, 4 } },
1650  { 3, 3, 3, 3, 4 },
1651  { { 0, 1, 4 }, { 0, 1, 2 }, { 0, 2, 3 }, { 0, 3, 4 }, { 1, 2, 3, 4 } },
1652  { { 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 0 }, { 0, 4 }, { 1, 4 }, { 2, 4 }, { 3, 4 } },
1653  { { 0, 1 }, { 0, 2 }, { 0, 3 }, { 0, 4 }, { 1, 4 }, { 1, 2 }, { 2, 3 }, { 3, 4 } },
1654  { { 3, 2, 1, 0 }, { 0, 5, 4 }, { 1, 6, 5 }, { 2, 7, 6 }, { 3, 4, 7 } },
1655  { { -1, 0, -1, 3, 4 }, { 0, -1, 1, -1, 5 }, { -1, 1, -1, 2, 6 }, { 3, -1, 2, -1, 7 }, { 4, 5, 6, 7, -1 } },
1656  { 3, 4, 2, 0 },
1657  { 0, 4 },
1658  { { 4, 0, 1, 2, 3 }, { 2, 1, 3 }, { 2, 1, 3 } } },
1659 
1660  // Prism
1661  { 6,
1662  9,
1663  5,
1664  { 4, 4, 4, 3, 3 },
1665  { { 0, 1, 4, 3 }, { 1, 2, 5, 4 }, { 0, 3, 5, 2 }, { 0, 2, 1 }, { 3, 4, 5 } },
1666  { 3, 3, 3, 3, 3, 3 },
1667  { { 0, 2, 3 }, { 0, 1, 3 }, { 1, 2, 3 }, { 0, 2, 4 }, { 0, 1, 4 }, { 1, 4, 2 } },
1668  { { 0, 1 }, { 1, 2 }, { 2, 0 }, { 0, 3 }, { 1, 4 }, { 2, 5 }, { 3, 4 }, { 4, 5 }, { 5, 3 } },
1669  { { 0, 3 }, { 1, 3 }, { 2, 3 }, { 0, 2 }, { 0, 1 }, { 1, 2 }, { 0, 4 }, { 1, 4 }, { 2, 4 } },
1670  { { 0, 4, 6, 3 }, { 1, 5, 7, 4 }, { 2, 3, 8, 5 }, { 2, 1, 0 }, { 6, 7, 8 } },
1671  { { -1, 0, 2, 3, -1, -1 },
1672  { 0, -1, 1, -1, 4, -1 },
1673  { 2, 1, -1, -1, -1, 5 },
1674  { 3, -1, -1, -1, 6, 8 },
1675  { -1, 4, -1, 6, -1, 7 },
1676  { -1, -1, 5, 8, 7, -1 } },
1677  { 4, 0, 5, 4, 1 },
1678  { 0, 5 },
1679  { { 3, 1, 2, 3 }, { 3, 2, 4, 3 }, { 2, 3, 1 }, { 1, 2 } } },
1680 
1681  // Hex
1682  { 8,
1683  12,
1684  6,
1685  { 4, 4, 4, 4, 4, 4 },
1686  { { 0, 1, 5, 4 }, { 1, 2, 6, 5 }, { 2, 3, 7, 6 }, { 3, 0, 4, 7 }, { 0, 3, 2, 1 }, { 4, 5, 6, 7 } },
1687  { 3, 3, 3, 3, 3, 3, 3, 3 },
1688  { { 0, 3, 4 }, { 0, 1, 4 }, { 1, 2, 4 }, { 2, 3, 4 }, { 0, 3, 5 }, { 0, 1, 5 }, { 1, 2, 5 }, { 2, 3, 5 } },
1689  { { 0, 1 },
1690  { 1, 2 },
1691  { 2, 3 },
1692  { 3, 0 },
1693  { 0, 4 },
1694  { 1, 5 },
1695  { 2, 6 },
1696  { 3, 7 },
1697  { 4, 5 },
1698  { 5, 6 },
1699  { 6, 7 },
1700  { 7, 4 } },
1701  { { 0, 4 },
1702  { 1, 4 },
1703  { 2, 4 },
1704  { 3, 4 },
1705  { 0, 3 },
1706  { 0, 1 },
1707  { 1, 2 },
1708  { 2, 3 },
1709  { 0, 5 },
1710  { 1, 5 },
1711  { 2, 5 },
1712  { 3, 5 } },
1713  { { 0, 5, 8, 4 }, { 1, 6, 9, 5 }, { 2, 7, 10, 6 }, { 3, 4, 11, 7 }, { 3, 2, 1, 0 }, { 8, 9, 10, 11 } },
1714  { { -1, 0, -1, 3, 4, -1, -1, -1 },
1715  { 0, -1, 1, -1, -1, 5, -1, -1 },
1716  { -1, 1, -1, 2, -1, -1, 6, -1 },
1717  { 3, -1, 2, -1, -1, -1, -1, 7 },
1718  { 4, -1, -1, -1, -1, 8, -1, 11 },
1719  { -1, 5, -1, -1, 8, -1, 9, -1 },
1720  { -1, -1, 6, -1, -1, 9, -1, 10 },
1721  { -1, -1, -1, 7, 11, -1, 10, -1 } },
1722  { 4, 0, 2, 5, 7 },
1723  { 0, 6 },
1724  { { 3, 1, 3, 4 }, { 3, 1, 3, 6 }, { 3, 1, 4, 6 }, { 3, 3, 6, 4 } } } };
1725 
1726 //////////////////////////////////////////////////////////////////////////////////
1727 ErrorCode HalfFacetRep::determine_sibling_halffaces( Range& cells )
1728 {
1729  ErrorCode error;
1730  EntityHandle start_cell = *cells.begin();
1731  EntityType ctype = mb->type_from_handle( start_cell );
1732  int index = get_index_in_lmap( start_cell );
1733  int nvpc = lConnMap3D[index].num_verts_in_cell;
1734  int nfpc = lConnMap3D[index].num_faces_in_cell;
1735 
1736  // Step 1: Create an index list storing the starting position for each vertex
1737  int nv = _verts.size();
1738  std::vector< int > is_index( nv + 1 );
1739  for( int i = 0; i < nv + 1; i++ )
1740  is_index[i] = 0;
1741 
1742  int vindex;
1743 
1744  for( Range::iterator cid = cells.begin(); cid != cells.end(); ++cid )
1745  {
1746  const EntityHandle* conn;
1747  error = mb->get_connectivity( *cid, conn, nvpc, true );MB_CHK_ERR( error );
1748 
1749  for( int i = 0; i < nfpc; ++i )
1750  {
1751  int nvF = lConnMap3D[index].hf2v_num[i];
1752  EntityHandle v = 0;
1753  for( int k = 0; k < nvF; k++ )
1754  {
1755  int id = lConnMap3D[index].hf2v[i][k];
1756  if( v <= conn[id] ) v = conn[id];
1757  }
1758  vindex = _verts.index( v );
1759  is_index[vindex + 1] += 1;
1760  }
1761  }
1762  is_index[0] = 0;
1763 
1764  for( int i = 0; i < nv; i++ )
1765  is_index[i + 1] = is_index[i] + is_index[i + 1];
1766 
1767  // Step 2: Define four arrays v2hv_eid, v2hv_lvid storing every half-facet on a vertex
1768  std::vector< EntityHandle > v2oe_v1( is_index[nv] );
1769  std::vector< EntityHandle > v2oe_v2( is_index[nv] );
1770  std::vector< EntityHandle > v2hf_map_cid( is_index[nv] );
1771  std::vector< int > v2hf_map_lfid( is_index[nv] );
1772 
1773  for( Range::iterator cid = cells.begin(); cid != cells.end(); ++cid )
1774  {
1775  const EntityHandle* conn;
1776  error = mb->get_connectivity( *cid, conn, nvpc, true );MB_CHK_ERR( error );
1777 
1778  for( int i = 0; i < nfpc; i++ )
1779  {
1780  int nvF = lConnMap3D[index].hf2v_num[i];
1781  std::vector< EntityHandle > vs( nvF );
1782  EntityHandle vmax = 0;
1783  int lv = -1;
1784  for( int k = 0; k < nvF; k++ )
1785  {
1786  int id = lConnMap3D[index].hf2v[i][k];
1787  vs[k] = conn[id];
1788  if( vmax <= conn[id] )
1789  {
1790  vmax = conn[id];
1791  lv = k;
1792  }
1793  }
1794 
1795  int nidx = lConnMap2D[nvF - 3].next[lv];
1796  int pidx = lConnMap2D[nvF - 3].prev[lv];
1797 
1798  int v = _verts.index( vmax );
1799  v2oe_v1[is_index[v]] = vs[nidx];
1800  v2oe_v2[is_index[v]] = vs[pidx];
1801  v2hf_map_cid[is_index[v]] = *cid;
1802  v2hf_map_lfid[is_index[v]] = i;
1803  is_index[v] += 1;
1804  }
1805  }
1806 
1807  for( int i = nv - 2; i >= 0; i-- )
1808  is_index[i + 1] = is_index[i];
1809  is_index[0] = 0;
1810 
1811  // Step 3: Fill up sibling half-verts map
1812  for( Range::iterator cid = cells.begin(); cid != cells.end(); ++cid )
1813  {
1814  const EntityHandle* conn;
1815  error = mb->get_connectivity( *cid, conn, nvpc, true );MB_CHK_ERR( error );
1816 
1817  int cidx = ID_FROM_HANDLE( *cid ) - 1;
1818  for( int i = 0; i < nfpc; i++ )
1819  {
1820  HFacet hf = sibhfs[nfpc * cidx + i];
1821  EntityHandle sibcid = fid_from_halfacet( hf, ctype );
1822 
1823  if( sibcid != 0 ) continue;
1824 
1825  int nvF = lConnMap3D[index].hf2v_num[i];
1826  std::vector< EntityHandle > vs( nvF );
1827  EntityHandle vmax = 0;
1828  int lv = -1;
1829  for( int k = 0; k < nvF; k++ )
1830  {
1831  int id = lConnMap3D[index].hf2v[i][k];
1832  vs[k] = conn[id];
1833  if( vmax <= conn[id] )
1834  {
1835  vmax = conn[id];
1836  lv = k;
1837  }
1838  }
1839  if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
1840  int nidx = lConnMap2D[nvF - 3].next[lv];
1841  int pidx = lConnMap2D[nvF - 3].prev[lv];
1842 
1843  int v = _verts.index( vmax );
1844  EntityHandle v1 = vs[pidx];
1845  EntityHandle v2 = vs[nidx];
1846 
1847  for( int ind = is_index[v]; ind <= is_index[v + 1] - 1; ind++ )
1848  {
1849  if( ( v2oe_v1[ind] == v1 ) && ( v2oe_v2[ind] == v2 ) )
1850  {
1851  // Map to opposite hf
1852  EntityHandle cur_cid = v2hf_map_cid[ind];
1853  int cur_lfid = v2hf_map_lfid[ind];
1854 
1855  sibhfs[nfpc * cidx + i] = create_halffacet( cur_cid, cur_lfid );
1856 
1857  int scidx = ID_FROM_HANDLE( cur_cid ) - 1;
1858  sibhfs[nfpc * scidx + cur_lfid] = create_halffacet( *cid, i );
1859  }
1860  }
1861  }
1862  }
1863 
1864  return MB_SUCCESS;
1865 }
1866 
1867 ErrorCode HalfFacetRep::determine_incident_halffaces( Range& cells )
1868 {
1869  ErrorCode error;
1870  int index = get_index_in_lmap( *cells.begin() );
1871  int nvpc = lConnMap3D[index].num_verts_in_cell;
1872 
1873  std::multimap< EntityHandle, EntityHandle > comps;
1874  HFacet nwhf;
1875 
1876  for( Range::iterator cid = cells.begin(); cid != cells.end(); ++cid )
1877  {
1878  EntityHandle cell = *cid;
1879  const EntityHandle* conn;
1880  error = mb->get_connectivity( *cid, conn, nvpc, true );MB_CHK_ERR( error );
1881 
1882  for( int i = 0; i < nvpc; ++i )
1883  {
1884  EntityHandle v = conn[i];
1885  int vidx = ID_FROM_HANDLE( v ) - 1;
1886  HFacet hf = v2hf[vidx];
1887 
1888  bool found = find_cell_in_component( v, cell, comps );
1889 
1890  if( hf == 0 && !found && ( v2hfs.empty() || ( v2hfs.find( v ) == v2hfs.end() ) ) )
1891  {
1892  nwhf = 0;
1893  error = add_cells_of_single_component( v, cell, lConnMap3D[index].v2hf[i][0], comps, nwhf );MB_CHK_ERR( error );
1894 
1895  v2hf[vidx] = nwhf;
1896  }
1897  else if( hf != 0 && !found )
1898  {
1899  nwhf = 0;
1900  error = add_cells_of_single_component( v, cell, lConnMap3D[index].v2hf[i][0], comps, nwhf );MB_CHK_ERR( error );
1901 
1902  v2hfs.insert( std::pair< EntityHandle, HFacet >( v, hf ) );
1903  v2hfs.insert( std::pair< EntityHandle, HFacet >( v, nwhf ) );
1904  v2hf[vidx] = 0;
1905  }
1906  else if( hf == 0 && !found && ( !v2hfs.empty() ) && ( v2hfs.find( v ) != v2hfs.end() ) )
1907  {
1908  nwhf = 0;
1909  error = add_cells_of_single_component( v, cell, lConnMap3D[index].v2hf[i][0], comps, nwhf );MB_CHK_ERR( error );
1910  v2hfs.insert( std::pair< EntityHandle, HFacet >( v, nwhf ) );
1911  }
1912  }
1913  }
1914 
1915  // error = print_tags(3);
1916 
1917  return MB_SUCCESS;
1918 }
1919 
1920 ErrorCode HalfFacetRep::determine_border_vertices( Range& cells, Tag isborder )
1921 {
1922  ErrorCode error;
1923  EntityHandle start_cell = *cells.begin();
1924  EntityType ctype = mb->type_from_handle( *cells.begin() );
1925  int index = get_index_in_lmap( start_cell );
1926  int nvpc = lConnMap3D[index].num_verts_in_cell;
1927  int nfpc = lConnMap3D[index].num_faces_in_cell;
1928 
1929  int val = 1;
1930 
1931  for( Range::iterator t = cells.begin(); t != cells.end(); ++t )
1932  {
1933 
1934  const EntityHandle* conn;
1935  error = mb->get_connectivity( *t, conn, nvpc, true );MB_CHK_ERR( error );
1936 
1937  int cidx = ID_FROM_HANDLE( *t ) - 1;
1938  for( int i = 0; i < nfpc; ++i )
1939  {
1940  HFacet hf = sibhfs[nfpc * cidx + i];
1941  EntityHandle sib_cid = fid_from_halfacet( hf, ctype );
1942 
1943  if( sib_cid == 0 )
1944  {
1945  int nvF = lConnMap3D[index].hf2v_num[i];
1946 
1947  for( int j = 0; j < nvF; ++j )
1948  {
1949  int ind = lConnMap3D[index].hf2v[i][j];
1950  error = mb->tag_set_data( isborder, &conn[ind], 1, &val );MB_CHK_ERR( error );
1951  }
1952  }
1953  }
1954  }
1955 
1956  return MB_SUCCESS;
1957 }
1958 
1959 /////////////////////////////////////////////////////////////////////////
1960 ErrorCode HalfFacetRep::add_cells_of_single_component( EntityHandle vid,
1961  EntityHandle curcid,
1962  int curlid,
1963  std::multimap< EntityHandle, EntityHandle >& comps,
1964  HFacet& hf )
1965 {
1966  ErrorCode error;
1967  EntityType ctype = mb->type_from_handle( curcid );
1968  int index = get_index_in_lmap( curcid );
1969  int nvpc = lConnMap3D[index].num_verts_in_cell;
1970  int nfpc = lConnMap3D[index].num_faces_in_cell;
1971 
1972  int Stksize = 0, count = -1;
1973  Stkcells[0] = curcid;
1974 
1975  hf = create_halffacet( curcid, curlid );
1976 
1977  EntityHandle cur_cid;
1978  while( Stksize >= 0 )
1979  {
1980  cur_cid = Stkcells[Stksize];
1981  Stksize -= 1;
1982 
1983  bool found = find_match_in_array( cur_cid, trackcells, count );
1984  if( !found )
1985  {
1986  count += 1;
1987  trackcells[count] = cur_cid;
1988 
1989  // Add the current cell
1990  comps.insert( std::pair< EntityHandle, EntityHandle >( vid, cur_cid ) );
1991  }
1992 
1993  // Connectivity of the cell
1994  const EntityHandle* conn;
1995  error = mb->get_connectivity( cur_cid, conn, nvpc, true );MB_CHK_ERR( error );
1996 
1997  // Local id of vid in the cell and the half-faces incident on it
1998  int lv = -1;
1999  for( int i = 0; i < nvpc; ++i )
2000  {
2001  if( conn[i] == vid )
2002  {
2003  lv = i;
2004  break;
2005  }
2006  }
2007  if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
2008  int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
2009  int cidx = ID_FROM_HANDLE( cur_cid ) - 1;
2010 
2011  // Add new cells into the stack
2012  EntityHandle ngb;
2013  HFacet hf_ngb;
2014  for( int i = 0; i < nhf_thisv; ++i )
2015  {
2016  int ind = lConnMap3D[index].v2hf[lv][i];
2017  hf_ngb = sibhfs[nfpc * cidx + ind];
2018  ngb = fid_from_halfacet( hf_ngb, ctype );
2019 
2020  if( ngb )
2021  {
2022  bool found_ent = find_match_in_array( ngb, trackcells, count );
2023 
2024  if( !found_ent )
2025  {
2026  Stksize += 1;
2027  Stkcells[Stksize] = ngb;
2028  }
2029  }
2030  else
2031  hf = create_halffacet( cur_cid, ind );
2032  }
2033  }
2034 
2035  // Change the visited faces to false
2036  for( int i = 0; i < Stksize; i++ )
2037  Stkcells[i] = 0;
2038 
2039  for( int i = 0; i <= count; i++ )
2040  trackcells[i] = 0;
2041 
2042  return MB_SUCCESS;
2043 }
2044 
2045 bool HalfFacetRep::find_cell_in_component( EntityHandle vid,
2046  EntityHandle cell,
2047  std::multimap< EntityHandle, EntityHandle >& comps )
2048 {
2049  bool found = false;
2050 
2051  if( comps.empty() ) return found;
2052 
2053  std::pair< std::multimap< EntityHandle, EntityHandle >::iterator,
2054  std::multimap< EntityHandle, EntityHandle >::iterator >
2055  rit;
2056 
2057  rit = comps.equal_range( vid );
2058 
2059  for( std::multimap< EntityHandle, EntityHandle >::iterator it = rit.first; it != rit.second; ++it )
2060  {
2061  if( it->second == cell )
2062  {
2063  found = true;
2064  break;
2065  }
2066  }
2067 
2068  return found;
2069 }
2070 
2071 ////////////////////////////////////////////////////////////////////////
2072 ErrorCode HalfFacetRep::get_up_adjacencies_vert_3d( EntityHandle vid, std::vector< EntityHandle >& adjents )
2073 {
2074  ErrorCode error;
2075  adjents.reserve( 20 );
2076  EntityType ctype = mb->type_from_handle( *_cells.begin() );
2077 
2078  // Obtain a half-face/s incident on v
2079  int vidx = ID_FROM_HANDLE( vid ) - 1;
2080  HFacet hf = v2hf[vidx];
2081 
2082  std::vector< EntityHandle > start_cells;
2083  if( hf == 0 && ( v2hfs.find( vid ) != v2hfs.end() ) ) // Vertex is non-manifold
2084  {
2085  std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
2086  it_hes;
2087  it_hes = v2hfs.equal_range( vid );
2088 
2089  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
2090  {
2091  start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
2092  }
2093  }
2094  else if( hf != 0 )
2095  start_cells.push_back( fid_from_halfacet( hf, ctype ) );
2096 
2097  if( start_cells.empty() ) return MB_SUCCESS;
2098 
2099  int index = get_index_in_lmap( start_cells[0] );
2100  int nvpc = lConnMap3D[index].num_verts_in_cell;
2101  int nfpc = lConnMap3D[index].num_faces_in_cell;
2102 
2103  for( int i = 0; i < (int)start_cells.size(); i++ )
2104  cellq[i] = start_cells[i];
2105 
2106  int qsize = start_cells.size();
2107  EntityHandle cur_cid;
2108  int num_qvals = 0;
2109 
2110  while( num_qvals < qsize )
2111  {
2112  cur_cid = cellq[num_qvals];
2113  num_qvals += 1;
2114 
2115  // Add the current cell to output adj vector
2116  adjents.push_back( cur_cid );
2117 
2118  // Connectivity of the cell
2119  const EntityHandle* conn;
2120  error = mb->get_connectivity( cur_cid, conn, nvpc, true );MB_CHK_ERR( error );
2121 
2122  // Local id of vid in the cell and the half-faces incident on it
2123  int lv = -1;
2124  for( int i = 0; i < nvpc; ++i )
2125  {
2126  if( conn[i] == vid )
2127  {
2128  lv = i;
2129  break;
2130  }
2131  };
2132  if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
2133  // Number of local half-faces incident on the current vertex
2134  int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
2135  int cidx = ID_FROM_HANDLE( cur_cid ) - 1;
2136 
2137  // Add new cells into the stack
2138  EntityHandle ngb;
2139  for( int i = 0; i < nhf_thisv; ++i )
2140  {
2141  int ind = lConnMap3D[index].v2hf[lv][i];
2142  hf = sibhfs[nfpc * cidx + ind];
2143  ngb = fid_from_halfacet( hf, ctype );
2144 
2145  if( ngb )
2146  {
2147  bool found_ent = find_match_in_array( ngb, cellq, qsize - 1 );
2148 
2149  if( !found_ent )
2150  {
2151  cellq[qsize] = ngb;
2152  qsize += 1;
2153  }
2154  }
2155  }
2156  }
2157 
2158  // Change the visited faces to false
2159  for( int i = 0; i < qsize; i++ )
2160  cellq[i] = 0;
2161 
2162  return MB_SUCCESS;
2163 }
2164 
2165 ErrorCode HalfFacetRep::get_up_adjacencies_edg_3d( EntityHandle eid,
2166  std::vector< EntityHandle >& adjents,
2167  std::vector< int >* leids )
2168 {
2169  ErrorCode error;
2170  EntityType ctype = mb->type_from_handle( *_cells.begin() );
2171  EntityHandle start_cell = *_cells.begin();
2172  int index = get_index_in_lmap( start_cell );
2173  int nvpc = lConnMap3D[index].num_verts_in_cell;
2174  int nfpc = lConnMap3D[index].num_faces_in_cell;
2175 
2176  adjents.reserve( 20 );
2177  if( leids != NULL ) leids->reserve( 20 );
2178 
2179  // Find the edge vertices
2180  const EntityHandle* econn;
2181  int num_conn = 0;
2182  error = mb->get_connectivity( eid, econn, num_conn, true );MB_CHK_ERR( error );
2183 
2184  EntityHandle v_start = econn[0], v_end = econn[1];
2185  int v1idx = ID_FROM_HANDLE( v_start ) - 1;
2186  int v2idx = ID_FROM_HANDLE( v_end ) - 1;
2187 
2188  // Find an half-facets incident to each end vertex of the edge
2189  std::vector< EntityHandle > start_cells;
2190  HFacet hf1 = v2hf[v1idx];
2191  HFacet hf2 = v2hf[v2idx];
2192 
2193  if( hf1 == 0 && !v2hfs.empty() )
2194  {
2195  std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
2196  it_hes;
2197  it_hes = v2hfs.equal_range( v_start );
2198 
2199  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
2200  {
2201  start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
2202  }
2203  }
2204  else if( hf1 != 0 )
2205  {
2206  start_cells.push_back( fid_from_halfacet( hf1, ctype ) );
2207  }
2208 
2209  if( hf2 == 0 && !v2hfs.empty() )
2210  {
2211  std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
2212  it_hes;
2213  it_hes = v2hfs.equal_range( v_end );
2214 
2215  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
2216  {
2217  start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
2218  }
2219  }
2220  else if( hf2 != 0 )
2221  {
2222  start_cells.push_back( fid_from_halfacet( hf2, ctype ) );
2223  }
2224 
2225  if( start_cells.empty() ) return MB_SUCCESS;
2226 
2227  std::sort( start_cells.begin(), start_cells.end() );
2228  std::vector< EntityHandle >::iterator last = std::unique( start_cells.begin(), start_cells.end() );
2229  start_cells.erase( last, start_cells.end() );
2230 
2231  for( int i = 0; i < (int)start_cells.size(); i++ )
2232  cellq[i] = start_cells[i];
2233 
2234  int qsize = start_cells.size();
2235  int num_qvals = 0;
2236 
2237  while( num_qvals < qsize )
2238  {
2239  EntityHandle cell_id = cellq[num_qvals];
2240  num_qvals += 1;
2241 
2242  const EntityHandle* conn;
2243  error = mb->get_connectivity( cell_id, conn, nvpc, true );MB_CHK_ERR( error );
2244 
2245  int lv0 = -1, lv1 = -1, lv = -1;
2246 
2247  // locate v_origin in poped out tet, check if v_end is in
2248  for( int i = 0; i < nvpc; i++ )
2249  {
2250  if( v_start == conn[i] )
2251  {
2252  lv0 = i;
2253  lv = lv0;
2254  }
2255  else if( v_end == conn[i] )
2256  {
2257  lv1 = i;
2258  lv = lv1;
2259  }
2260  }
2261 
2262  if( ( lv0 >= 0 ) && ( lv1 >= 0 ) )
2263  {
2264  adjents.push_back( cell_id );
2265  if( leids != NULL ) leids->push_back( lConnMap3D[index].lookup_leids[lv0][lv1] );
2266  }
2267 
2268  // push back new found unchecked incident tets of v_start
2269  int cidx = ID_FROM_HANDLE( cell_id ) - 1;
2270  if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
2271  int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
2272 
2273  for( int i = 0; i < nhf_thisv; i++ )
2274  {
2275  int ind = lConnMap3D[index].v2hf[lv][i];
2276  HFacet hf = sibhfs[nfpc * cidx + ind];
2277  EntityHandle ngb = fid_from_halfacet( hf, ctype );
2278 
2279  if( ngb )
2280  {
2281  bool found_ent = find_match_in_array( ngb, &cellq[0], qsize - 1 );
2282 
2283  if( !found_ent )
2284  {
2285  cellq[qsize] = ngb;
2286  qsize += 1;
2287  }
2288  }
2289  }
2290  }
2291 
2292  for( int i = 0; i < qsize; i++ )
2293  cellq[i] = 0;
2294 
2295  return MB_SUCCESS;
2296 }
2297 
2298 ErrorCode HalfFacetRep::get_up_adjacencies_edg_3d( EntityHandle cid,
2299  int leid,
2300  std::vector< EntityHandle >& adjents,
2301  std::vector< int >* leids,
2302  std::vector< int >* adj_orients )
2303 {
2304  ErrorCode error;
2305  EntityType ctype = mb->type_from_handle( cid );
2306  int index = get_index_in_lmap( cid );
2307  int nvpc = lConnMap3D[index].num_verts_in_cell;
2308  int nfpc = lConnMap3D[index].num_faces_in_cell;
2309  adjents.clear();
2310  adjents.reserve( 20 );
2311 
2312  if( leids != NULL )
2313  {
2314  leids->clear();
2315  leids->reserve( 20 );
2316  }
2317  if( adj_orients != NULL )
2318  {
2319  adj_orients->clear();
2320  adj_orients->reserve( 20 );
2321  }
2322 
2323  const EntityHandle* econn;
2324  error = mb->get_connectivity( cid, econn, nvpc, true );MB_CHK_ERR( error );
2325 
2326  // Get the end vertices of the edge <cid,leid>
2327  int id = lConnMap3D[index].e2v[leid][0];
2328  EntityHandle v_start = econn[id];
2329  id = lConnMap3D[index].e2v[leid][1];
2330  EntityHandle v_end = econn[id];
2331 
2332  int v1idx = ID_FROM_HANDLE( v_start ) - 1;
2333  int v2idx = ID_FROM_HANDLE( v_end ) - 1;
2334 
2335  // Find an half-facets incident to each end vertex of the edge
2336  std::vector< EntityHandle > start_cells;
2337  HFacet hf1 = v2hf[v1idx];
2338  HFacet hf2 = v2hf[v2idx];
2339 
2340  if( hf1 == 0 && !v2hfs.empty() )
2341  {
2342  std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
2343  it_hes;
2344  it_hes = v2hfs.equal_range( v_start );
2345 
2346  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
2347  {
2348  start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
2349  }
2350  }
2351  else if( hf1 != 0 )
2352  {
2353  start_cells.push_back( fid_from_halfacet( hf1, ctype ) );
2354  }
2355 
2356  if( hf2 == 0 && !v2hfs.empty() )
2357  {
2358  std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
2359  it_hes;
2360  it_hes = v2hfs.equal_range( v_end );
2361 
2362  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
2363  {
2364  start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
2365  }
2366  }
2367  else if( hf2 != 0 )
2368  {
2369  start_cells.push_back( fid_from_halfacet( hf2, ctype ) );
2370  }
2371 
2372  if( start_cells.empty() ) return MB_SUCCESS;
2373 
2374  std::sort( start_cells.begin(), start_cells.end() );
2375  std::vector< EntityHandle >::iterator last = std::unique( start_cells.begin(), start_cells.end() );
2376  start_cells.erase( last, start_cells.end() );
2377 
2378  for( int i = 0; i < (int)start_cells.size(); i++ )
2379  cellq[i] = start_cells[i];
2380 
2381  int qsize = start_cells.size();
2382  int num_qvals = 0;
2383 
2384  while( num_qvals < qsize )
2385  {
2386  EntityHandle cell_id = cellq[num_qvals];
2387  num_qvals += 1;
2388 
2389  const EntityHandle* conn;
2390  error = mb->get_connectivity( cell_id, conn, nvpc, true );MB_CHK_ERR( error );
2391 
2392  int lv0 = -1, lv1 = -1, lv = -1;
2393 
2394  // locate v_origin in poped out tet, check if v_end is in
2395  for( int i = 0; i < nvpc; i++ )
2396  {
2397  if( v_start == conn[i] )
2398  {
2399  lv0 = i;
2400  lv = lv0;
2401  }
2402  else if( v_end == conn[i] )
2403  {
2404  lv1 = i;
2405  lv = lv1;
2406  }
2407  }
2408 
2409  if( ( lv0 >= 0 ) && ( lv1 >= 0 ) )
2410  {
2411  adjents.push_back( cell_id );
2412  if( leids != NULL ) leids->push_back( lConnMap3D[index].lookup_leids[lv0][lv1] );
2413 
2414  if( adj_orients != NULL )
2415  {
2416  int cur_leid = lConnMap3D[index].lookup_leids[lv0][lv1];
2417  int id1 = lConnMap3D[index].e2v[cur_leid][0];
2418  int id2 = lConnMap3D[index].e2v[cur_leid][1];
2419  if( ( v_start == conn[id1] ) && ( v_end == conn[id2] ) )
2420  adj_orients->push_back( 1 );
2421  else if( ( v_start == conn[id2] ) && ( v_end == conn[id1] ) )
2422  adj_orients->push_back( 0 );
2423  }
2424  }
2425  if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
2426  // push back new found unchecked incident tets of v_start
2427  int cidx = ID_FROM_HANDLE( cell_id ) - 1;
2428  int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
2429 
2430  for( int i = 0; i < nhf_thisv; i++ )
2431  {
2432  int ind = lConnMap3D[index].v2hf[lv][i];
2433  HFacet hf = sibhfs[nfpc * cidx + ind];
2434  EntityHandle ngb = fid_from_halfacet( hf, ctype );
2435 
2436  if( ngb )
2437  {
2438  bool found_ent = find_match_in_array( ngb, &cellq[0], qsize - 1 );
2439 
2440  if( !found_ent )
2441  {
2442  cellq[qsize] = ngb;
2443  qsize += 1;
2444  }
2445  }
2446  }
2447  }
2448 
2449  for( int i = 0; i < qsize; i++ )
2450  cellq[i] = 0;
2451 
2452  return MB_SUCCESS;
2453 }
2454 
2455 ErrorCode HalfFacetRep::get_up_adjacencies_edg_3d_comp( EntityHandle cid,
2456  int leid,
2457  std::vector< EntityHandle >& adjents,
2458  std::vector< int >* leids,
2459  std::vector< int >* adj_orients )
2460 {
2461  ErrorCode error;
2462  EntityType ctype = mb->type_from_handle( cid );
2463  int index = get_index_in_lmap( cid );
2464  int nvpc = lConnMap3D[index].num_verts_in_cell;
2465  int nfpc = lConnMap3D[index].num_faces_in_cell;
2466  adjents.clear();
2467  adjents.reserve( 20 );
2468 
2469  if( leids != NULL )
2470  {
2471  leids->clear();
2472  leids->reserve( 20 );
2473  }
2474  if( adj_orients != NULL )
2475  {
2476  adj_orients->clear();
2477  adj_orients->reserve( 20 );
2478  }
2479 
2480  const EntityHandle* econn;
2481  error = mb->get_connectivity( cid, econn, nvpc );MB_CHK_ERR( error );
2482 
2483  // Get the end vertices of the edge <cid,leid>
2484  int id = lConnMap3D[index].e2v[leid][0];
2485  EntityHandle v_start = econn[id];
2486  id = lConnMap3D[index].e2v[leid][1];
2487  EntityHandle v_end = econn[id];
2488 
2489  int v1idx = ID_FROM_HANDLE( v_start ) - 1;
2490  int v2idx = ID_FROM_HANDLE( v_end ) - 1;
2491 
2492  // Find an half-facets incident to each end vertex of the edge
2493  std::vector< EntityHandle > start_cells;
2494  HFacet hf1 = v2hf[v1idx];
2495  HFacet hf2 = v2hf[v2idx];
2496 
2497  if( ( hf1 == 0 ) && ( v2hfs.find( v_start ) != v2hfs.end() ) && ( hf2 == 0 ) &&
2498  ( v2hfs.find( v_end ) != v2hfs.end() ) )
2499  {
2500  std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
2501  it_hes;
2502  it_hes = v2hfs.equal_range( v_start );
2503 
2504  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
2505  {
2506  start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
2507  }
2508  }
2509  else
2510  return MB_SUCCESS;
2511 
2512  if( start_cells.empty() ) return MB_SUCCESS;
2513 
2514  // std::sort(start_cells.begin(), start_cells.end());
2515  // std::vector<EntityHandle>::iterator last = std::unique(start_cells.begin(),
2516  // start_cells.end()); start_cells.erase(last, start_cells.end());
2517 
2518  for( int c = 0; c < (int)start_cells.size(); c++ )
2519  {
2520  cellq[0] = start_cells[c];
2521 
2522  int qsize = 1;
2523  int num_qvals = 0;
2524 
2525  while( num_qvals < qsize )
2526  {
2527  EntityHandle cell_id = cellq[num_qvals];
2528  num_qvals += 1;
2529 
2530  const EntityHandle* conn;
2531  error = mb->get_connectivity( cell_id, conn, nvpc );MB_CHK_ERR( error );
2532 
2533  int lv0 = -1, lv1 = -1, lv = -1;
2534 
2535  // locate v_origin in poped out tet, check if v_end is in
2536  for( int i = 0; i < nvpc; i++ )
2537  {
2538  if( v_start == conn[i] )
2539  {
2540  lv0 = i;
2541  lv = lv0;
2542  }
2543  else if( v_end == conn[i] )
2544  {
2545  lv1 = i;
2546  lv = lv1;
2547  }
2548  }
2549 
2550  if( ( lv0 >= 0 ) && ( lv1 >= 0 ) )
2551  {
2552  adjents.push_back( cell_id );
2553  if( leids != NULL ) leids->push_back( lConnMap3D[index].lookup_leids[lv0][lv1] );
2554 
2555  if( adj_orients != NULL )
2556  {
2557  int cur_leid = lConnMap3D[index].lookup_leids[lv0][lv1];
2558  int id1 = lConnMap3D[index].e2v[cur_leid][0];
2559  int id2 = lConnMap3D[index].e2v[cur_leid][1];
2560  if( ( v_start == conn[id1] ) && ( v_end == conn[id2] ) )
2561  adj_orients->push_back( 1 );
2562  else if( ( v_start == conn[id2] ) && ( v_end == conn[id1] ) )
2563  adj_orients->push_back( 0 );
2564  }
2565 
2566  for( int i = 0; i < qsize; i++ )
2567  cellq[i] = 0;
2568 
2569  break;
2570  }
2571 
2572  // push back new found unchecked incident tets of v_start
2573  int cidx = ID_FROM_HANDLE( cell_id ) - 1;
2574  int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
2575 
2576  for( int i = 0; i < nhf_thisv; i++ )
2577  {
2578  int ind = lConnMap3D[index].v2hf[lv][i];
2579  HFacet hf = sibhfs[nfpc * cidx + ind];
2580  EntityHandle ngb = fid_from_halfacet( hf, ctype );
2581 
2582  if( ngb )
2583  {
2584  bool found_ent = find_match_in_array( ngb, &cellq[0], qsize - 1 );
2585 
2586  if( !found_ent )
2587  {
2588  cellq[qsize] = ngb;
2589  qsize += 1;
2590  }
2591  }
2592  }
2593  }
2594 
2595  for( int i = 0; i < qsize; i++ )
2596  cellq[i] = 0;
2597  }
2598 
2599  return MB_SUCCESS;
2600 }
2601 
2602 ErrorCode HalfFacetRep::get_up_adjacencies_face_3d( EntityHandle fid,
2603  std::vector< EntityHandle >& adjents,
2604  std::vector< int >* lfids )
2605 {
2606  ErrorCode error;
2607 
2608  EntityHandle cid = 0;
2609  int lid = 0;
2610  bool found = find_matching_halfface( fid, &cid, &lid );
2611 
2612  if( found )
2613  {
2614  error = get_up_adjacencies_face_3d( cid, lid, adjents, lfids );MB_CHK_ERR( error );
2615  }
2616 
2617  return MB_SUCCESS;
2618 }
2619 
2620 ErrorCode HalfFacetRep::get_up_adjacencies_face_3d( EntityHandle cid,
2621  int lfid,
2622  std::vector< EntityHandle >& adjents,
2623  std::vector< int >* lfids )
2624 {
2625 
2626  EntityHandle start_cell = *_cells.begin();
2627  EntityType ctype = mb->type_from_handle( start_cell );
2628  int index = get_index_in_lmap( start_cell );
2629  int nfpc = lConnMap3D[index].num_faces_in_cell;
2630 
2631  adjents.reserve( 4 );
2632  adjents.push_back( cid );
2633 
2634  if( lfids != NULL )
2635  {
2636  lfids->reserve( 4 );
2637  lfids->push_back( lfid );
2638  }
2639 
2640  int cidx = ID_FROM_HANDLE( cid ) - 1;
2641  HFacet hf = sibhfs[nfpc * cidx + lfid];
2642  EntityHandle sibcid = fid_from_halfacet( hf, ctype );
2643  int siblid = lid_from_halffacet( hf );
2644 
2645  if( sibcid != 0 )
2646  {
2647  adjents.push_back( sibcid );
2648  if( lfids != NULL ) lfids->push_back( siblid );
2649  }
2650 
2651  return MB_SUCCESS;
2652 }
2653 
2654 bool HalfFacetRep::find_matching_implicit_edge_in_cell( EntityHandle eid,
2655  std::vector< EntityHandle >& cid,
2656  std::vector< int >& leid )
2657 {
2658  ErrorCode error;
2659  EntityType ctype = mb->type_from_handle( *_cells.begin() );
2660  EntityHandle start_cell = *_cells.begin();
2661  int index = get_index_in_lmap( start_cell );
2662  int nvpc = lConnMap3D[index].num_verts_in_cell;
2663  int nfpc = lConnMap3D[index].num_faces_in_cell;
2664 
2665  // Find the edge vertices
2666  const EntityHandle* econn;
2667  int num_conn = 0;
2668  error = mb->get_connectivity( eid, econn, num_conn, true );MB_CHK_ERR( error );
2669 
2670  EntityHandle v_start = econn[0], v_end = econn[1];
2671  int v1idx = ID_FROM_HANDLE( v_start ) - 1;
2672  int v2idx = ID_FROM_HANDLE( v_end ) - 1;
2673 
2674  // Find an incident cell to v_start
2675  std::vector< EntityHandle > start_cells;
2676  HFacet hf1 = v2hf[v1idx];
2677  HFacet hf2 = v2hf[v2idx];
2678 
2679  int ncomps1 = 0, ncomps2 = 0, ncomp;
2680  if( hf1 == 0 && !v2hfs.empty() )
2681  {
2682  std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
2683  it_hes;
2684  it_hes = v2hfs.equal_range( v_start );
2685 
2686  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
2687  {
2688  start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
2689  ncomps1 += 1;
2690  }
2691  }
2692  else if( hf1 != 0 )
2693  {
2694  start_cells.push_back( fid_from_halfacet( hf1, ctype ) );
2695  ncomps1 += 1;
2696  }
2697 
2698  if( hf2 == 0 && !v2hfs.empty() )
2699  {
2700  std::pair< std::multimap< EntityHandle, HFacet >::iterator, std::multimap< EntityHandle, HFacet >::iterator >
2701  it_hes;
2702  it_hes = v2hfs.equal_range( v_end );
2703 
2704  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
2705  {
2706  start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
2707  ncomps2 += 1;
2708  }
2709  }
2710  else if( hf2 != 0 )
2711  {
2712  start_cells.push_back( fid_from_halfacet( hf2, ctype ) );
2713  ncomps2 += 1;
2714  }
2715 
2716  ncomp = std::min( ncomps1, ncomps2 );
2717 
2718  bool found = false;
2719  if( start_cells.empty() ) return found;
2720 
2721  for( int i = 0; i < (int)start_cells.size(); i++ )
2722  cellq[i] = start_cells[i];
2723 
2724  int qsize = start_cells.size();
2725  int num_qvals = 0;
2726 
2727  while( num_qvals < qsize )
2728  {
2729  EntityHandle cell_id = cellq[num_qvals];
2730  num_qvals += 1;
2731 
2732  const EntityHandle* conn;
2733  error = mb->get_connectivity( cell_id, conn, nvpc, true );MB_CHK_ERR( error );
2734 
2735  int lv0 = -1, lv1 = -1, lv = -1;
2736 
2737  // locate v_origin in poped out tet, check if v_end is in
2738  for( int i = 0; i < nvpc; i++ )
2739  {
2740  if( v_start == conn[i] )
2741  {
2742  lv0 = i;
2743  lv = lv0;
2744  }
2745  else if( v_end == conn[i] )
2746  {
2747  lv1 = i;
2748  lv = lv1;
2749  }
2750  }
2751 
2752  if( ( lv0 >= 0 ) && ( lv1 >= 0 ) )
2753  {
2754  found = true;
2755  cid.push_back( cell_id );
2756  leid.push_back( lConnMap3D[index].lookup_leids[lv0][lv1] );
2757 
2758  if( (int)cid.size() == ncomp ) break;
2759  }
2760 
2761  // push back new found unchecked incident tets of v_start
2762  int cidx = ID_FROM_HANDLE( cell_id ) - 1;
2763  int nhf_thisv = lConnMap3D[index].v2hf_num[lv];
2764  if( lv < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
2765  for( int i = 0; i < nhf_thisv; i++ )
2766  {
2767  int ind = lConnMap3D[index].v2hf[lv][i];
2768  HFacet hf = sibhfs[nfpc * cidx + ind];
2769  EntityHandle ngb = fid_from_halfacet( hf, ctype );
2770 
2771  if( ngb )
2772  {
2773  bool found_ent = find_match_in_array( ngb, &cellq[0], qsize - 1 );
2774 
2775  if( !found_ent )
2776  {
2777  cellq[qsize] = ngb;
2778  qsize += 1;
2779  }
2780  }
2781  }
2782  }
2783 
2784  for( int i = 0; i < qsize; i++ )
2785  cellq[i] = 0;
2786 
2787  return found;
2788 }
2789 
2790 bool HalfFacetRep::find_matching_halfface( EntityHandle fid, EntityHandle* cid, int* lid )
2791 {
2792  ErrorCode error;
2793  EntityHandle start_cell = *_cells.begin();
2794  EntityType ctype = mb->type_from_handle( start_cell );
2795  int index = get_index_in_lmap( start_cell );
2796  int nvpc = lConnMap3D[index].num_verts_in_cell;
2797  int nfpc = lConnMap3D[index].num_faces_in_cell;
2798  EntityType ftype = mb->type_from_handle( fid );
2799  int nvF = lConnMap2D[ftype - 2].num_verts_in_face;
2800 
2801  const EntityHandle* fid_verts;
2802  error = mb->get_connectivity( fid, fid_verts, nvF, true );MB_CHK_ERR( error );
2803 
2804  std::vector< EntityHandle > start_cells;
2805  int vidx, locfv0 = -1;
2806  HFacet hf = 0;
2807 
2808  for( int i = 0; i < nvF; i++ )
2809  {
2810  vidx = ID_FROM_HANDLE( fid_verts[i] ) - 1;
2811  hf = v2hf[vidx];
2812  if( hf != 0 )
2813  {
2814  start_cells.push_back( fid_from_halfacet( hf, ctype ) );
2815  locfv0 = i;
2816  break;
2817  }
2818  else if( hf == 0 && v2hfs.find( fid_verts[i] ) != v2hfs.end() )
2819  {
2820  std::pair< std::multimap< EntityHandle, HFacet >::iterator,
2821  std::multimap< EntityHandle, HFacet >::iterator >
2822  it_hfs;
2823  it_hfs = v2hfs.equal_range( fid_verts[i] );
2824 
2825  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hfs.first; it != it_hfs.second; ++it )
2826  {
2827  start_cells.push_back( fid_from_halfacet( it->second, ctype ) );
2828  }
2829  locfv0 = i;
2830  break;
2831  }
2832  }
2833 
2834  if( start_cells.empty() ) return false;
2835 
2836  EntityHandle cur_cid;
2837  bool found = false;
2838 
2839  int Stksize = 0, count = -1;
2840 
2841  for( int i = 0; i < (int)start_cells.size(); i++ )
2842  Stkcells[i] = start_cells[i];
2843 
2844  Stksize = start_cells.size() - 1;
2845 
2846  while( Stksize >= 0 )
2847  {
2848  cur_cid = Stkcells[Stksize];
2849  Stksize -= 1;
2850  count += 1;
2851  trackcells[count] = cur_cid;
2852 
2853  const EntityHandle* conn;
2854  error = mb->get_connectivity( cur_cid, conn, nvpc, true );MB_CHK_ERR( error );
2855 
2856  int lv[4] = { -1, -1, -1, -1 };
2857  int cnt = 0;
2858  for( int i = 0; i < nvpc; i++ )
2859  {
2860  for( int j = 0; j < nvF; j++ )
2861  if( conn[i] == fid_verts[j] )
2862  {
2863  lv[j] = i;
2864  cnt += 1;
2865  }
2866  }
2867  if( cnt == nvF ) // All face verts are part of the cell
2868  {
2869  found = true;
2870  int nhf_thisv = lConnMap3D[index].v2hf_num[lv[locfv0]];
2871  int lfid = -1;
2872  for( int i = 0; i < nhf_thisv; ++i )
2873  {
2874  lfid = lConnMap3D[index].v2hf[lv[locfv0]][i];
2875  int lcnt = 0;
2876  for( int j = 0; j < nvF; ++j )
2877  {
2878  for( int k = 0; k < nvF; ++k )
2879  {
2880  if( lv[k] == lConnMap3D[index].hf2v[lfid][j] ) lcnt += 1;
2881  }
2882  }
2883  if( lcnt == nvF ) break;
2884  }
2885  cid[0] = cur_cid;
2886  lid[0] = lfid;
2887 
2888  break;
2889  }
2890  else
2891  {
2892  // Add other cells that are incident on fid_verts[0]
2893  if( locfv0 < 0 || lv[locfv0] < 0 ) MB_SET_ERR( MB_FAILURE, "did not find local vertex " );
2894  int nhf_thisv = lConnMap3D[index].v2hf_num[lv[locfv0]];
2895  int cidx = ID_FROM_HANDLE( cur_cid ) - 1;
2896 
2897  // Add new cells into the stack
2898  EntityHandle ngb;
2899  for( int i = 0; i < nhf_thisv; ++i )
2900  {
2901  int ind = lConnMap3D[index].v2hf[lv[locfv0]][i];
2902  hf = sibhfs[nfpc * cidx + ind];
2903  ngb = fid_from_halfacet( hf, ctype );
2904 
2905  if( ngb )
2906  {
2907 
2908  bool found_ent = find_match_in_array( ngb, trackcells, count );
2909 
2910  if( !found_ent )
2911  {
2912  Stksize += 1;
2913  Stkcells[Stksize] = ngb;
2914  }
2915  }
2916  }
2917  }
2918  }
2919 
2920  // Change the visited faces to false
2921  for( int i = 0; i < Stksize; i++ )
2922  Stkcells[i] = 0;
2923 
2924  for( int i = 0; i <= count; i++ )
2925  trackcells[i] = 0;
2926 
2927  return found;
2928 }
2929 
2930 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2931 ErrorCode HalfFacetRep::get_neighbor_adjacencies_3d( EntityHandle cid, std::vector< EntityHandle >& adjents )
2932 {
2933  adjents.reserve( 20 );
2934  EntityType ctype = mb->type_from_handle( cid );
2935  int index = get_index_in_lmap( cid );
2936  int nfpc = lConnMap3D[index].num_faces_in_cell;
2937  int cidx = ID_FROM_HANDLE( cid ) - 1;
2938 
2939  if( cid != 0 )
2940  {
2941  for( int lfid = 0; lfid < nfpc; ++lfid )
2942  {
2943  HFacet hf = sibhfs[nfpc * cidx + lfid];
2944  EntityHandle sibcid = fid_from_halfacet( hf, ctype );
2945  if( sibcid != 0 ) adjents.push_back( sibcid );
2946  }
2947  }
2948 
2949  return MB_SUCCESS;
2950 }
2951 
2952 /////////////////////////////////////////////////////////////////////////////////////////////////
2953 ErrorCode HalfFacetRep::get_down_adjacencies_edg_3d( EntityHandle cid, std::vector< EntityHandle >& adjents )
2954 {
2955  // TODO: Try intersection without using std templates
2956  // Returns explicit edges, if any, of the face
2957  ErrorCode error;
2958  adjents.reserve( 20 );
2959  int index = get_index_in_lmap( cid );
2960  int nvpc = lConnMap3D[index].num_verts_in_cell;
2961 
2962  const EntityHandle* conn;
2963  error = mb->get_connectivity( cid, conn, nvpc, true );MB_CHK_ERR( error );
2964 
2965  // Gather all the incident edges on each vertex of the face
2966  int ns = lConnMap3D[index].search_everts[0];
2967  std::vector< EntityHandle > temp;
2968  for( int i = 0; i < ns; i++ )
2969  {
2970  temp.clear();
2971  int lv0 = lConnMap3D[index].search_everts[i + 1];
2972  error = get_up_adjacencies_1d( conn[lv0], temp );MB_CHK_ERR( error );
2973 
2974  int nle = lConnMap3D[index].v2le[i][0];
2975  int count = 0;
2976  for( int j = 0; j < (int)temp.size(); j++ )
2977  {
2978  const EntityHandle* econn;
2979  int nvpe = 0;
2980  error = mb->get_connectivity( temp[j], econn, nvpe, true );MB_CHK_ERR( error );
2981 
2982  for( int k = 0; k < nle; k++ )
2983  {
2984  int lv1 = lConnMap3D[index].v2le[i][k + 1];
2985  if( ( ( econn[0] == conn[lv0] ) && ( econn[1] == conn[lv1] ) ) ||
2986  ( ( econn[1] == conn[lv0] ) && ( econn[0] == conn[lv1] ) ) )
2987  {
2988  adjents.push_back( temp[j] );
2989  count += 1;
2990  }
2991  }
2992  if( count == nle ) break;
2993  }
2994  }
2995  return MB_SUCCESS;
2996 }
2997 
2998 ErrorCode HalfFacetRep::get_down_adjacencies_face_3d( EntityHandle cid, std::vector< EntityHandle >& adjents )
2999 {
3000  // Returns explicit face, if any of the cell
3001  ErrorCode error;
3002  adjents.reserve( 10 );
3003  int index = get_index_in_lmap( cid );
3004  int nvpc = lConnMap3D[index].num_verts_in_cell;
3005  int nfpc = lConnMap3D[index].num_faces_in_cell;
3006 
3007  // Get the connectivity of the input cell
3008  const EntityHandle* conn;
3009  error = mb->get_connectivity( cid, conn, nvpc, true );MB_CHK_ERR( error );
3010 
3011  // Collect all the half-faces of the cell
3012  EntityHandle half_faces[6][4];
3013  for( int i = 0; i < nfpc; i++ )
3014  {
3015  int nvf = lConnMap3D[index].hf2v_num[i];
3016  for( int j = 0; j < nvf; j++ )
3017  {
3018  int ind = lConnMap3D[index].hf2v[i][j];
3019  half_faces[i][j] = conn[ind];
3020  }
3021  }
3022 
3023  // Add two vertices for which the upward adjacencies are computed
3024  int search_verts[2];
3025  search_verts[0] = lConnMap3D[index].search_fverts[0];
3026  search_verts[1] = lConnMap3D[index].search_fverts[1];
3027 
3028  std::vector< EntityHandle > temp;
3029  temp.reserve( 20 );
3030  for( int i = 0; i < 2; i++ )
3031  {
3032  // Get the incident faces on the local vertex
3033  int lv = search_verts[i];
3034  temp.clear();
3035  error = get_up_adjacencies_vert_2d( conn[lv], temp );MB_CHK_ERR( error );
3036 
3037  if( temp.size() == 0 ) continue;
3038 
3039  // Get the half-faces incident on the local vertex and match it with the obtained faces
3040  int nhfthisv = lConnMap3D[index].v2hf_num[lv];
3041  for( int k = 0; k < (int)temp.size(); k++ )
3042  {
3043  const EntityHandle* fid_verts;
3044  int fsize = 0;
3045  error = mb->get_connectivity( temp[k], fid_verts, fsize, true );MB_CHK_ERR( error );
3046 
3047  for( int j = 0; j < nhfthisv; j++ )
3048  {
3049  // Collect all the vertices of this half-face
3050  int idx = lConnMap3D[index].v2hf[lv][j];
3051  int nvF = lConnMap3D[index].hf2v_num[idx];
3052 
3053  if( fsize != nvF ) continue;
3054 
3055  int direct, offset;
3056  bool they_match = CN::ConnectivityMatch( &half_faces[idx][0], &fid_verts[0], nvF, direct, offset );
3057 
3058  if( they_match )
3059  {
3060  bool found = false;
3061  for( int p = 0; p < (int)adjents.size(); p++ )
3062  {
3063  if( adjents[p] == temp[k] )
3064  {
3065  found = true;
3066  break;
3067  }
3068  }
3069  if( !found ) adjents.push_back( temp[k] );
3070  }
3071  }
3072  }
3073  }
3074 
3075  return MB_SUCCESS;
3076 }
3077 ////////////////////////////////////////////////////////////////////////////////
3078 ErrorCode HalfFacetRep::find_total_edges_faces_3d( const Range& cells, int* nedges, int* nfaces )
3079 {
3080  ErrorCode error;
3081  int index = get_index_in_lmap( *cells.begin() );
3082  int nepc = lConnMap3D[index].num_edges_in_cell;
3083  int nfpc = lConnMap3D[index].num_faces_in_cell;
3084  int ncells = cells.size();
3085  int total_edges = nepc * ncells;
3086  int total_faces = nfpc * ncells;
3087 
3088  std::vector< int > trackE( total_edges, 0 );
3089  std::vector< int > trackF( total_faces, 0 );
3090 
3091  std::vector< EntityHandle > inc_cids, sib_cids;
3092  std::vector< int > inc_leids, sib_lfids;
3093 
3094  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
3095  {
3096  // Count edges
3097  for( int i = 0; i < nepc; i++ )
3098  {
3099  inc_cids.clear();
3100  inc_leids.clear();
3101 
3102  int id = nepc * ( cells.index( *it ) ) + i;
3103  if( !trackE[id] )
3104  {
3105  error = get_up_adjacencies_edg_3d( *it, i, inc_cids, &inc_leids );MB_CHK_ERR( error );
3106 
3107  total_edges -= inc_cids.size() - 1;
3108  for( int j = 0; j < (int)inc_cids.size(); j++ )
3109  trackE[nepc * ( cells.index( inc_cids[j] ) ) + inc_leids[j]] = 1;
3110  }
3111  }
3112 
3113  // Count faces
3114  for( int i = 0; i < nfpc; i++ )
3115  {
3116  sib_cids.clear();
3117  sib_lfids.clear();
3118 
3119  int id = nfpc * ( cells.index( *it ) ) + i;
3120  if( !trackF[id] )
3121  {
3122  error = get_up_adjacencies_face_3d( *it, i, sib_cids, &sib_lfids );MB_CHK_ERR( error );
3123 
3124  if( sib_cids.size() == 1 ) continue;
3125 
3126  total_faces -= sib_cids.size() - 1;
3127  trackF[nfpc * ( cells.index( sib_cids[1] ) ) + sib_lfids[1]] = 1;
3128  }
3129  }
3130  }
3131 
3132  nedges[0] = total_edges;
3133  nfaces[0] = total_faces;
3134 
3135  return MB_SUCCESS;
3136 }
3137 
3138 bool HalfFacetRep::find_match_in_array( EntityHandle ent,
3139  EntityHandle* ent_list,
3140  int count,
3141  bool get_index,
3142  int* index )
3143 {
3144  bool found = false;
3145  for( int i = 0; i <= count; i++ )
3146  {
3147  if( ent == ent_list[i] )
3148  {
3149  found = true;
3150  if( get_index ) *index = i;
3151  break;
3152  }
3153  }
3154 
3155  return found;
3156 }
3157 
3158 ErrorCode HalfFacetRep::get_half_facet_in_comp( EntityHandle cid,
3159  int leid,
3160  std::vector< EntityHandle >& ents,
3161  std::vector< int >& lids,
3162  std::vector< int >& lfids )
3163 {
3164  ErrorCode error;
3165  ents.clear();
3166  lids.clear();
3167  EntityType ctype = mb->type_from_handle( cid );
3168  int index = get_index_in_lmap( *_cells.begin() );
3169  int nfpc = lConnMap3D[index].num_faces_in_cell;
3170  int nvpc = lConnMap3D[index].num_verts_in_cell;
3171 
3172  // Get all incident cells
3173  std::vector< EntityHandle > adjents;
3174  std::vector< int > adjlids;
3175  error = get_up_adjacencies_edg_3d( cid, leid, adjents, &adjlids );MB_CHK_ERR( error );
3176 
3177  // Get the end vertices of the edge <cid,leid>
3178  const EntityHandle* econn;
3179  error = mb->get_connectivity( cid, econn, nvpc, true );MB_CHK_ERR( error );
3180  int id = lConnMap3D[index].e2v[leid][0];
3181  EntityHandle vstart = econn[id];
3182  id = lConnMap3D[index].e2v[leid][1];
3183  EntityHandle vend = econn[id];
3184 
3185  std::vector< int > mark( adjents.size(), 0 );
3186  int count = 0;
3187 
3188  for( int k = 0; k < (int)adjents.size(); k++ )
3189  {
3190 
3191  if( mark[k] != 0 ) continue;
3192 
3193  count += 1;
3194  mark[k] = count;
3195 
3196  // Loop over each half-face incident on this local edge
3197  for( int i = 0; i < 2; i++ )
3198  {
3199  EntityHandle cur_cell = adjents[k];
3200  int cur_leid = adjlids[k];
3201 
3202  int lface = i;
3203 
3204  while( true )
3205  {
3206  int cidx = ID_FROM_HANDLE( cur_cell ) - 1;
3207  int lfid = lConnMap3D[index].e2hf[cur_leid][lface];
3208 
3209  HFacet hf = sibhfs[nfpc * cidx + lfid];
3210  cur_cell = fid_from_halfacet( hf, ctype );
3211  lfid = lid_from_halffacet( hf );
3212 
3213  // Check if loop reached starting cell or a boundary
3214  if( ( cur_cell == adjents[k] ) || ( cur_cell == 0 ) ) break;
3215 
3216  const EntityHandle* sib_conn;
3217  error = mb->get_connectivity( cur_cell, sib_conn, nvpc, true );MB_CHK_ERR( error );
3218 
3219  // Find the local edge id wrt to sibhf
3220  int nv_curF = lConnMap3D[index].hf2v_num[lfid];
3221  int lv0 = -1, lv1 = -1, idx = -1;
3222  for( int j = 0; j < nv_curF; j++ )
3223  {
3224  idx = lConnMap3D[index].hf2v[lfid][j];
3225  if( vstart == sib_conn[idx] ) lv0 = idx;
3226  if( vend == sib_conn[idx] ) lv1 = idx;
3227  }
3228 
3229  assert( ( lv0 >= 0 ) && ( lv1 >= 0 ) );
3230  cur_leid = lConnMap3D[index].lookup_leids[lv0][lv1];
3231 
3232  int chk_lfid = lConnMap3D[index].e2hf[cur_leid][0];
3233 
3234  if( lfid == chk_lfid )
3235  lface = 1;
3236  else
3237  lface = 0;
3238 
3239  int ind = std::find( adjents.begin(), adjents.end(), cur_cell ) - adjents.begin();
3240  mark[ind] = count;
3241  }
3242 
3243  // Loop back
3244  if( cur_cell != 0 ) break;
3245  }
3246  }
3247 
3248  // Loop over again to find cells on the boundary
3249  for( int c = 0; c < count; c++ )
3250  {
3251  for( int i = 0; i < (int)adjents.size(); i++ )
3252  {
3253  if( mark[i] == c + 1 )
3254  {
3255  int cidx = ID_FROM_HANDLE( adjents[i] ) - 1;
3256  for( int j = 0; j < nfpc; j++ )
3257  {
3258  HFacet hf = sibhfs[nfpc * cidx + j];
3259  EntityHandle cell = fid_from_halfacet( hf, ctype );
3260  if( cell == 0 )
3261  {
3262  ents.push_back( adjents[i] );
3263  lids.push_back( adjlids[i] );
3264  lfids.push_back( j );
3265  break;
3266  }
3267  }
3268  }
3269  }
3270  }
3271 
3272  return MB_SUCCESS;
3273 }
3274 
3275 ///////////////////////////////////////////////////////////////////////////////////////////
3276 ErrorCode HalfFacetRep::resize_hf_maps( EntityHandle start_vert,
3277  int nverts,
3278  EntityHandle start_edge,
3279  int nedges,
3280  EntityHandle start_face,
3281  int nfaces,
3282  EntityHandle start_cell,
3283  int ncells )
3284 {
3285  int nwsz = 0, insz = 0;
3286  if( nedges )
3287  {
3288  if( ID_FROM_HANDLE( ( *( _edges.end() - 1 ) + 1 ) ) != ID_FROM_HANDLE( start_edge ) )
3289  nwsz = ( ID_FROM_HANDLE( start_edge ) - ID_FROM_HANDLE( *_edges.end() ) + nedges ) * 2;
3290  else
3291  nwsz = nedges * 2;
3292  insz = sibhvs.size();
3293  sibhvs.resize( insz + nwsz, 0 );
3294 
3295  if( v2hv.empty() )
3296  {
3297  if( ( !v2he.empty() ) )
3298  insz = v2he.size();
3299  else if( ( !v2hf.empty() ) )
3300  insz = v2hf.size();
3301  else
3302  MB_SET_ERR( MB_FAILURE, "Trying to resize ahf maps for a mesh with no edges, faces and cells" );
3303  }
3304  else
3305  insz = v2hv.size();
3306 
3307  if( ID_FROM_HANDLE( *( _verts.end() - 1 ) + 1 ) != ID_FROM_HANDLE( start_vert ) )
3308  nwsz = ID_FROM_HANDLE( start_vert ) - ID_FROM_HANDLE( *_verts.end() ) + nverts;
3309  else
3310  nwsz = nverts;
3311  v2hv.resize( insz + nwsz, 0 );
3312  }
3313 
3314  if( nfaces )
3315  {
3316  EntityType ftype = mb->type_from_handle( *_faces.begin() );
3317  int nepf = lConnMap2D[ftype - 2].num_verts_in_face;
3318 
3319  if( ID_FROM_HANDLE( ( *( _faces.end() - 1 ) + 1 ) ) != ID_FROM_HANDLE( start_face ) )
3320  nwsz = ( ID_FROM_HANDLE( start_face ) - ID_FROM_HANDLE( *_faces.end() ) + nfaces ) * nepf;
3321  else
3322  nwsz = nfaces * nepf;
3323  insz = sibhes.size();
3324  sibhes.resize( insz + nwsz, 0 );
3325 
3326  if( ID_FROM_HANDLE( *( _verts.end() - 1 ) + 1 ) != ID_FROM_HANDLE( start_vert ) )
3327  nwsz = ID_FROM_HANDLE( start_vert ) - ID_FROM_HANDLE( *_verts.end() ) + nverts;
3328  else
3329  nwsz = nverts;
3330  insz = v2he.size();
3331  v2he.resize( insz + nwsz, 0 );
3332  }
3333 
3334  if( ncells )
3335  {
3336  int index = get_index_in_lmap( *_cells.begin() );
3337  int nfpc = lConnMap3D[index].num_faces_in_cell;
3338 
3339  if( ID_FROM_HANDLE( ( *( _cells.end() - 1 ) + 1 ) ) != ID_FROM_HANDLE( start_cell ) )
3340  nwsz = ( ID_FROM_HANDLE( start_cell ) - ID_FROM_HANDLE( *_cells.end() ) + ncells ) * nfpc;
3341  else
3342  nwsz = ncells * nfpc;
3343  insz = sibhfs.size();
3344  sibhfs.resize( insz + nwsz, 0 );
3345 
3346  if( ID_FROM_HANDLE( *( _verts.end() - 1 ) + 1 ) != ID_FROM_HANDLE( start_vert ) )
3347  nwsz = ID_FROM_HANDLE( start_vert ) - ID_FROM_HANDLE( *_verts.end() ) + nverts;
3348  else
3349  nwsz = nverts;
3350  insz = v2hf.size();
3351  v2hf.resize( insz + nwsz, 0 );
3352  }
3353 
3354  return MB_SUCCESS;
3355 }
3356 
3357 bool HalfFacetRep::check_nonmanifold_vertices( EntityType type, EntityHandle vid )
3358 {
3359  bool status = false;
3360  if( type == MBTRI || type == MBQUAD )
3361  {
3362  HFacet hf = v2he[ID_FROM_HANDLE( vid ) - 1];
3363  if( hf == 0 && ( v2hes.find( vid ) != v2hes.end() ) ) status = true;
3364  }
3365  else if( type == MBTET || type == MBHEX )
3366  {
3367  HFacet hf = v2hf[ID_FROM_HANDLE( vid ) - 1];
3368  if( hf == 0 && ( v2hfs.find( vid ) != v2hfs.end() ) ) status = true;
3369  }
3370  else
3371  MB_SET_ERR( MB_FAILURE, "Requesting non-manifold vertex checks for either (1) 1D mesh or "
3372  "(2) not-implemented entity types" );
3373 
3374  return status;
3375 }
3376 
3377 ErrorCode HalfFacetRep::get_sibling_map( EntityType type,
3378  EntityHandle ent,
3379  EntityHandle* sib_entids,
3380  int* sib_lids,
3381  int num_halffacets )
3382 {
3383 
3384  if( type == MBEDGE )
3385  {
3386  if( num_halffacets != 2 ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halfvertices." );
3387 
3388  int eidx = ID_FROM_HANDLE( ent ) - 1;
3389  for( int i = 0; i < 2; i++ )
3390  {
3391  HFacet hf = sibhvs[2 * eidx + i];
3392  sib_entids[i] = fid_from_halfacet( hf, MBEDGE );
3393  sib_lids[i] = lid_from_halffacet( hf );
3394  }
3395  }
3396  else if( type == MBTRI || type == MBQUAD )
3397  {
3398  int nepf = lConnMap2D[type - 2].num_verts_in_face;
3399 
3400  if( num_halffacets != nepf ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halfedges." );
3401 
3402  int fidx = ID_FROM_HANDLE( ent ) - 1;
3403  for( int i = 0; i < nepf; i++ )
3404  {
3405  HFacet hf = sibhes[nepf * fidx + i];
3406  sib_entids[i] = fid_from_halfacet( hf, type );
3407  sib_lids[i] = lid_from_halffacet( hf );
3408  }
3409  }
3410  else
3411  {
3412  int idx = get_index_in_lmap( *_cells.begin() );
3413  int nfpc = lConnMap3D[idx].num_faces_in_cell;
3414 
3415  if( num_halffacets != nfpc ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halffaces." );
3416 
3417  int cidx = ID_FROM_HANDLE( ent ) - 1;
3418  for( int i = 0; i < nfpc; i++ )
3419  {
3420  HFacet hf = sibhfs[nfpc * cidx + i];
3421  sib_entids[i] = fid_from_halfacet( hf, type );
3422  sib_lids[i] = lid_from_halffacet( hf );
3423  }
3424  }
3425  return MB_SUCCESS;
3426 }
3427 
3428 ErrorCode HalfFacetRep::get_sibling_map( EntityType type,
3429  EntityHandle ent,
3430  int lid,
3431  EntityHandle& sib_entid,
3432  int& sib_lid )
3433 {
3434 
3435  if( type == MBEDGE )
3436  {
3437  int eidx = ID_FROM_HANDLE( ent ) - 1;
3438  HFacet hf = sibhvs[2 * eidx + lid];
3439  sib_entid = fid_from_halfacet( hf, MBEDGE );
3440  sib_lid = lid_from_halffacet( hf );
3441  }
3442  else if( type == MBTRI || type == MBQUAD )
3443  {
3444  int nepf = lConnMap2D[type - 2].num_verts_in_face;
3445  int fidx = ID_FROM_HANDLE( ent ) - 1;
3446  HFacet hf = sibhes[nepf * fidx + lid];
3447  sib_entid = fid_from_halfacet( hf, type );
3448  sib_lid = lid_from_halffacet( hf );
3449  }
3450  else
3451  {
3452  int idx = get_index_in_lmap( *_cells.begin() );
3453  int nfpc = lConnMap3D[idx].num_faces_in_cell;
3454  int cidx = ID_FROM_HANDLE( ent ) - 1;
3455  HFacet hf = sibhfs[nfpc * cidx + lid];
3456  sib_entid = fid_from_halfacet( hf, type );
3457  sib_lid = lid_from_halffacet( hf );
3458  }
3459  return MB_SUCCESS;
3460 }
3461 
3462 ErrorCode HalfFacetRep::set_sibling_map( EntityType type,
3463  EntityHandle ent,
3464  EntityHandle* set_entids,
3465  int* set_lids,
3466  int num_halffacets )
3467 {
3468 
3469  if( type == MBEDGE )
3470  {
3471  if( num_halffacets != 2 ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halfvertices" );
3472 
3473  int eidx = ID_FROM_HANDLE( ent ) - 1;
3474  for( int i = 0; i < 2; i++ )
3475  {
3476  sibhvs[2 * eidx + i] = create_halffacet( set_entids[i], set_lids[i] );
3477  }
3478  }
3479  else if( type == MBTRI || type == MBQUAD )
3480  {
3481  int nepf = lConnMap2D[type - 2].num_verts_in_face;
3482  if( num_halffacets != nepf ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halfedges." );
3483 
3484  int fidx = ID_FROM_HANDLE( ent ) - 1;
3485  for( int i = 0; i < nepf; i++ )
3486  {
3487  sibhes[nepf * fidx + i] = create_halffacet( set_entids[i], set_lids[i] );
3488  }
3489  }
3490  else
3491  {
3492  int idx = get_index_in_lmap( *_cells.begin() );
3493  int nfpc = lConnMap3D[idx].num_faces_in_cell;
3494  if( num_halffacets != nfpc ) MB_SET_ERR( MB_FAILURE, "Incorrect number of halffaces." );
3495 
3496  int cidx = ID_FROM_HANDLE( ent ) - 1;
3497  for( int i = 0; i < nfpc; i++ )
3498  {
3499  sibhfs[nfpc * cidx + i] = create_halffacet( set_entids[i], set_lids[i] );
3500  }
3501  }
3502 
3503  return MB_SUCCESS;
3504 }
3505 
3506 ErrorCode HalfFacetRep::set_sibling_map( EntityType type,
3507  EntityHandle ent,
3508  int lid,
3509  EntityHandle& set_entid,
3510  int& set_lid )
3511 {
3512 
3513  if( type == MBEDGE )
3514  {
3515  int eidx = ID_FROM_HANDLE( ent ) - 1;
3516  sibhvs[2 * eidx + lid] = create_halffacet( set_entid, set_lid );
3517  }
3518  else if( type == MBTRI || type == MBQUAD )
3519  {
3520  int nepf = lConnMap2D[type - 2].num_verts_in_face;
3521  int fidx = ID_FROM_HANDLE( ent ) - 1;
3522  sibhes[nepf * fidx + lid] = create_halffacet( set_entid, set_lid );
3523  }
3524  else
3525  {
3526  int idx = get_index_in_lmap( *_cells.begin() );
3527  int nfpc = lConnMap3D[idx].num_faces_in_cell;
3528  int cidx = ID_FROM_HANDLE( ent ) - 1;
3529  sibhfs[nfpc * cidx + lid] = create_halffacet( set_entid, set_lid );
3530  }
3531 
3532  return MB_SUCCESS;
3533 }
3534 
3535 ErrorCode HalfFacetRep::get_incident_map( EntityType type,
3536  EntityHandle vid,
3537  std::vector< EntityHandle >& inci_entid,
3538  std::vector< int >& inci_lid )
3539 {
3540  inci_entid.clear();
3541  inci_lid.clear();
3542 
3543  if( type == MBEDGE )
3544  {
3545  HFacet hf = v2hv[ID_FROM_HANDLE( vid ) - 1];
3546  inci_entid.push_back( fid_from_halfacet( hf, type ) );
3547  inci_lid.push_back( lid_from_halffacet( hf ) );
3548  }
3549  else if( type == MBTRI || type == MBQUAD )
3550  {
3551  HFacet hf = v2he[ID_FROM_HANDLE( vid ) - 1];
3552  if( hf == 0 && ( v2hes.find( vid ) != v2hes.end() ) )
3553  {
3554  std::pair< std::multimap< EntityHandle, HFacet >::iterator,
3555  std::multimap< EntityHandle, HFacet >::iterator >
3556  it_hes;
3557  it_hes = v2hes.equal_range( vid );
3558 
3559  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
3560  {
3561  inci_entid.push_back( fid_from_halfacet( it->second, type ) );
3562  inci_lid.push_back( lid_from_halffacet( it->second ) );
3563  }
3564  }
3565  else if( hf != 0 )
3566  {
3567  inci_entid.push_back( fid_from_halfacet( hf, type ) );
3568  inci_lid.push_back( lid_from_halffacet( hf ) );
3569  }
3570  else if( hf == 0 && ( v2hes.find( vid ) == v2hes.end() ) )
3571  {
3572  inci_entid.push_back( fid_from_halfacet( hf, type ) );
3573  inci_lid.push_back( lid_from_halffacet( hf ) );
3574  }
3575  }
3576  else
3577  {
3578  HFacet hf = v2hf[ID_FROM_HANDLE( vid ) - 1];
3579  if( hf == 0 && ( v2hfs.find( vid ) != v2hfs.end() ) )
3580  {
3581  std::pair< std::multimap< EntityHandle, HFacet >::iterator,
3582  std::multimap< EntityHandle, HFacet >::iterator >
3583  it_hes;
3584  it_hes = v2hfs.equal_range( vid );
3585 
3586  for( std::multimap< EntityHandle, HFacet >::iterator it = it_hes.first; it != it_hes.second; ++it )
3587  {
3588  inci_entid.push_back( fid_from_halfacet( it->second, type ) );
3589  inci_lid.push_back( lid_from_halffacet( it->second ) );
3590  }
3591  }
3592  else if( hf != 0 )
3593  {
3594  inci_entid.push_back( fid_from_halfacet( hf, type ) );
3595  inci_lid.push_back( lid_from_halffacet( hf ) );
3596  }
3597  else if( hf == 0 && ( v2hfs.find( vid ) == v2hfs.end() ) )
3598  {
3599  inci_entid.push_back( fid_from_halfacet( hf, type ) );
3600  inci_lid.push_back( lid_from_halffacet( hf ) );
3601  }
3602  }
3603 
3604  return MB_SUCCESS;
3605 }
3606 
3607 ErrorCode HalfFacetRep::set_incident_map( EntityType type,
3608  EntityHandle vid,
3609  std::vector< EntityHandle >& set_entid,
3610  std::vector< int >& set_lid )
3611 {
3612  if( type == MBEDGE )
3613  {
3614  v2hv[ID_FROM_HANDLE( vid ) - 1] = create_halffacet( set_entid[0], set_lid[0] );
3615  }
3616  else if( type == MBTRI || type == MBQUAD )
3617  {
3618  if( set_entid.size() == 1 )
3619  v2he[ID_FROM_HANDLE( vid ) - 1] = create_halffacet( set_entid[0], set_lid[0] );
3620  else
3621  {
3622  HFacet hf = 0;
3623  for( int i = 0; i < (int)set_entid.size(); i++ )
3624  {
3625  hf = create_halffacet( set_entid[i], set_lid[i] );
3626  v2hes.insert( std::pair< EntityHandle, HFacet >( vid, hf ) );
3627  }
3628  }
3629  }
3630  else
3631  {
3632  if( set_entid.size() == 1 )
3633  v2hf[ID_FROM_HANDLE( vid ) - 1] = create_halffacet( set_entid[0], set_lid[0] );
3634  else
3635  {
3636  HFacet hf = v2hf[ID_FROM_HANDLE( vid ) - 1];
3637  if( hf != 0 )
3638  {
3639  v2hf[ID_FROM_HANDLE( vid ) - 1] = 0;
3640  }
3641  for( int i = 0; i < (int)set_entid.size(); i++ )
3642  {
3643  hf = create_halffacet( set_entid[i], set_lid[i] );
3644  v2hfs.insert( std::pair< EntityHandle, HFacet >( vid, hf ) );
3645  }
3646  }
3647  }
3648 
3649  return MB_SUCCESS;
3650 }
3651 
3652 ErrorCode HalfFacetRep::get_entity_ranges( Range& verts, Range& edges, Range& faces, Range& cells )
3653 {
3654  verts = _verts;
3655  edges = _edges;
3656  faces = _faces;
3657  cells = _cells;
3658  return MB_SUCCESS;
3659 }
3660 
3661 ErrorCode HalfFacetRep::update_entity_ranges( EntityHandle fileset )
3662 {
3663  ErrorCode error;
3664 
3665  error = mb->get_entities_by_dimension( fileset, 0, _verts, true );MB_CHK_ERR( error );
3666 
3667  error = mb->get_entities_by_dimension( fileset, 1, _edges, true );MB_CHK_ERR( error );
3668 
3669  error = mb->get_entities_by_dimension( fileset, 2, _faces, true );MB_CHK_ERR( error );
3670 
3671  error = mb->get_entities_by_dimension( fileset, 3, _cells, true );MB_CHK_ERR( error );
3672 
3673  return MB_SUCCESS;
3674 }
3675 
3676 void HalfFacetRep::get_memory_use( unsigned long long& entity_total, unsigned long long& memory_total )
3677 {
3678  entity_total = memory_total = 0;
3679  // 1D
3680  if( !v2hv.empty() ) entity_total += v2hv.capacity() * sizeof( HFacet ) + sizeof( v2hv );
3681  if( !sibhvs.empty() ) entity_total += sibhvs.capacity() * sizeof( HFacet ) + sizeof( sibhvs );
3682 
3683  // 2D
3684  if( !v2he.empty() ) entity_total += v2he.capacity() * sizeof( HFacet ) + sizeof( v2he );
3685  if( !sibhes.empty() ) entity_total += sibhes.capacity() * sizeof( HFacet ) + sizeof( sibhes );
3686 
3687  // 3D
3688  if( !v2hf.empty() ) entity_total += v2hf.capacity() * sizeof( HFacet ) + sizeof( v2hf );
3689  if( !sibhfs.empty() ) entity_total += sibhfs.capacity() * sizeof( HFacet ) + sizeof( sibhfs );
3690 
3691  memory_total = entity_total;
3692 }
3693 
3694 HFacet HalfFacetRep::create_halffacet( EntityHandle handle, int lid )
3695 {
3696  EntityID fid = ID_FROM_HANDLE( handle );
3697  return CREATE_HALFFACET( lid, fid );
3698 }
3699 
3700 EntityHandle HalfFacetRep::fid_from_halfacet( const HFacet facet, EntityType type )
3701 {
3702  EntityID id = FID_FROM_HALFFACET( facet );
3703  EntityHandle handle = 0;
3704  if( id == 0 ) return handle;
3705 
3706  ErrorCode error = mb->handle_from_id( type, id, handle );MB_CHK_ERR( error );
3707  return handle;
3708 }
3709 
3710 int HalfFacetRep::lid_from_halffacet( const HFacet facet )
3711 {
3712  if( facet == 0 )
3713  return 0;
3714  else
3715  return LID_FROM_HALFFACET( facet );
3716 }
3717 
3718 } // namespace moab