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