Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
DualTool.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/DualTool.hpp"
17 #include "moab/Range.hpp"
18 // using Core for call to check_adjacencies
19 #include "moab/Core.hpp"
20 #include "Internals.hpp"
21 #include "MBTagConventions.hpp"
22 #include "moab/Skinner.hpp"
23 #include "moab/Core.hpp"
24 #include "moab/MeshTopoUtil.hpp"
25 #include "AEntityFactory.hpp"
26 #include "moab/CN.hpp"
27 #include <string>
28 #include <algorithm>
29 #include <iostream>
30 #include <sstream>
31 #include <cassert>
32 
33 #define RR \
34  if( MB_SUCCESS != result ) return result
35 #define SWAP( a, b ) \
36  { \
37  EntityHandle tmp_ent = a; \
38  ( a ) = b; \
39  ( b ) = tmp_ent; \
40  }
41 
42 namespace moab
43 {
44 
45 bool debug = false;
46 bool debug_ap = false;
47 
48 //! tag name for dual surfaces
49 const char* DualTool::DUAL_SURFACE_TAG_NAME = "DUAL_SURFACE";
50 
51 //! tag name for dual curves
52 const char* DualTool::DUAL_CURVE_TAG_NAME = "DUAL_CURVE";
53 
54 //! tag name for dual cells
55 const char* DualTool::IS_DUAL_CELL_TAG_NAME = "__IS_DUAL_CELL";
56 
57 //! tag name for dual entities
58 const char* DualTool::DUAL_ENTITY_TAG_NAME = "__DUAL_ENTITY";
59 
60 //! tag name for extra dual entities
61 const char* DualTool::EXTRA_DUAL_ENTITY_TAG_NAME = "__EXTRA_DUAL_ENTITY";
62 
63 //! tag name for graphics point
64 const char* DualTool::DUAL_GRAPHICS_POINT_TAG_NAME = "__DUAL_GRAPHICS_POINT";
65 
66 // const int DualTool::GP_SIZE = 20;
67 
68 DualTool::DualTool( Interface* impl ) : mbImpl( impl )
69 {
70  EntityHandle dum_handle = 0;
71  ErrorCode result;
72 
73  result = mbImpl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dualSurfaceTag,
74  MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle );
75  assert( MB_SUCCESS == result );
76 
77  result = mbImpl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dualCurveTag, MB_TAG_SPARSE | MB_TAG_CREAT,
78  &dum_handle );
79  assert( MB_SUCCESS == result );
80 
81  unsigned int dummy = 0;
82  result = mbImpl->tag_get_handle( IS_DUAL_CELL_TAG_NAME, 1, MB_TYPE_INTEGER, isDualCellTag,
83  MB_TAG_SPARSE | MB_TAG_CREAT, &dummy );
84  assert( MB_SUCCESS == result );
85 
86  result = mbImpl->tag_get_handle( DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, dualEntityTag,
87  MB_TAG_DENSE | MB_TAG_CREAT, &dum_handle );
88  assert( MB_SUCCESS == result );
89 
90  result = mbImpl->tag_get_handle( EXTRA_DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, extraDualEntityTag,
91  MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle );
92  assert( MB_SUCCESS == result );
93 
94  static const char dum_name[CATEGORY_TAG_SIZE] = { 0 };
95  result = mbImpl->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, categoryTag,
96  MB_TAG_SPARSE | MB_TAG_CREAT, dum_name );
97  assert( MB_SUCCESS == result );
98 
99  DualTool::GraphicsPoint dum_pt( 0.0, 0.0, 0.0, -1 );
100  result = mbImpl->tag_get_handle( DUAL_GRAPHICS_POINT_TAG_NAME, sizeof( DualTool::GraphicsPoint ), MB_TYPE_DOUBLE,
101  dualGraphicsPointTag, MB_TAG_DENSE | MB_TAG_CREAT | MB_TAG_BYTES, &dum_pt );
102  assert( MB_SUCCESS == result );
103 
104  globalIdTag = mbImpl->globalId_tag();
105 
106  if( MB_SUCCESS == result )
107  {
108  } // empty statement to get rid of warning.
109 
110  maxHexId = -1;
111 }
112 
113 DualTool::~DualTool() {}
114 
115 //! construct the dual entities for the entire mesh
116 ErrorCode DualTool::construct_dual( EntityHandle* entities, const int num_entities )
117 {
118  // allocate a dual entity for each primal entity in the mesh, starting
119  // with highest dimension and working downward; do each dimension in a separate code
120  // block, since they're all handled slightly differently
121 
122  Range regions, faces, edges, vertices;
123  ErrorCode result;
124 
125  if( NULL == entities || 0 == num_entities )
126  {
127 
128  // first, construct all the aentities, since they're currently needed to
129  // compute the dual
130  result = mbImpl->get_entities_by_dimension( 0, 0, vertices );
131  if( MB_SUCCESS != result ) return result;
132 
133  result = MeshTopoUtil( mbImpl ).construct_aentities( vertices );
134  if( MB_SUCCESS != result ) return result;
135 
136  // get all edges, faces and regions now, so we don't need to filter out dual
137  // entities later
138 
139  result = mbImpl->get_entities_by_dimension( 0, 1, edges );
140  if( MB_SUCCESS != result ) return result;
141  result = mbImpl->get_entities_by_dimension( 0, 2, faces );
142  if( MB_SUCCESS != result ) return result;
143  result = mbImpl->get_entities_by_dimension( 0, 3, regions );
144  if( MB_SUCCESS != result ) return result;
145 
146  // get the max global id for hexes, we'll need for modification ops
147  std::vector< int > gid_vec( regions.size() );
148  result = mbImpl->tag_get_data( globalId_tag(), regions, &gid_vec[0] );
149  if( MB_SUCCESS != result ) return result;
150  maxHexId = -1;
151  Range::iterator rit;
152  unsigned int i;
153  for( rit = regions.begin(), i = 0; rit != regions.end(); ++rit, i++ )
154  {
155  if( gid_vec[i] > maxHexId && mbImpl->type_from_handle( *rit ) == MBHEX ) maxHexId = gid_vec[i];
156  }
157  }
158  else
159  {
160  // get entities of various dimensions adjacent to these
161  result = mbImpl->get_adjacencies( entities, num_entities, 0, true, vertices, Interface::UNION );
162  if( MB_SUCCESS != result ) return result;
163  result = mbImpl->get_adjacencies( entities, num_entities, 1, true, edges, Interface::UNION );
164  if( MB_SUCCESS != result ) return result;
165  result = mbImpl->get_adjacencies( entities, num_entities, 2, true, faces, Interface::UNION );
166  if( MB_SUCCESS != result ) return result;
167  result = mbImpl->get_adjacencies( entities, num_entities, 3, true, regions, Interface::UNION );
168  if( MB_SUCCESS != result ) return result;
169  }
170 
171  Range dual_verts;
172  result = construct_dual_vertices( regions, dual_verts );
173  if( MB_SUCCESS != result || dual_verts.size() != regions.size() ) return result;
174  if( debug ) std::cout << "Constructed " << dual_verts.size() << " dual vertices." << std::endl;
175 
176  // don't really need dual edges, but construct 'em anyway
177  Range dual_edges;
178  result = construct_dual_edges( faces, dual_edges );
179  if( MB_SUCCESS != result || dual_edges.size() != faces.size() ) return result;
180  if( debug ) std::cout << "Constructed " << dual_edges.size() << " dual edges." << std::endl;
181 
182  // construct dual faces
183  Range dual_faces;
184  result = construct_dual_faces( edges, dual_faces );
185  if( MB_SUCCESS != result || dual_faces.size() != edges.size() ) return result;
186  if( debug ) std::cout << "Constructed " << dual_faces.size() << " dual faces." << std::endl;
187 
188  // construct dual cells
189  Range dual_cells;
190  result = construct_dual_cells( vertices, dual_cells );
191  if( MB_SUCCESS != result || dual_cells.size() != vertices.size() ) return result;
192  if( debug ) std::cout << "Constructed " << dual_cells.size() << " dual cells." << std::endl;
193 
194  return MB_SUCCESS;
195 }
196 
197 ErrorCode DualTool::construct_dual_vertices( const Range& all_regions, Range& dual_ents )
198 {
199  if( all_regions.empty() ) return MB_SUCCESS;
200 
201  // make sure they're all regions
202  assert( 3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.begin() ) ) &&
203  3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.rbegin() ) ) );
204 
205  Range::const_iterator rit;
206  EntityHandle dual_ent;
207  ErrorCode tmp_result = MB_SUCCESS;
208  ErrorCode result = MB_SUCCESS;
209 
210  for( rit = all_regions.begin(); rit != all_regions.end(); ++rit )
211  {
212  if( tmp_result != MB_SUCCESS ) result = tmp_result;
213 
214  tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
215  if( MB_SUCCESS == tmp_result && 0 != dual_ent )
216  {
217  dual_ents.insert( dual_ent );
218  continue;
219  }
220  else if( MB_SUCCESS != tmp_result )
221  continue;
222 
223  tmp_result = construct_dual_vertex( *rit, dual_ent, false, true );
224  if( MB_SUCCESS != tmp_result ) continue;
225 
226  // save it in the list of new dual ents
227  dual_ents.insert( dual_ent );
228  }
229 
230  return result;
231 }
232 
233 ErrorCode DualTool::construct_dual_vertex( EntityHandle entity,
234  EntityHandle& dual_ent,
235  const bool extra,
236  const bool add_graphics_pt )
237 {
238  // no dual entity; construct one; first need the avg coordinates
239  unsigned int is_dual = 0x1;
240  double avg_pos[3];
241  ErrorCode result = MeshTopoUtil( mbImpl ).get_average_position( entity, avg_pos );
242  if( MB_SUCCESS != result ) return result;
243 
244  // now construct the new dual entity
245  result = mbImpl->create_vertex( avg_pos, dual_ent );
246  if( MB_SUCCESS != result ) return result;
247 
248  // tag it indicating it's a dual entity
249  result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
250  if( MB_SUCCESS != result ) return result;
251 
252  // tag the primal entity with its dual entity and vica versa
253  if( extra )
254  result = mbImpl->tag_set_data( extraDualEntity_tag(), &( entity ), 1, &dual_ent );
255  else
256  result = mbImpl->tag_set_data( dualEntity_tag(), &( entity ), 1, &dual_ent );
257  if( MB_SUCCESS != result ) return result;
258 
259  result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( entity ) );
260  if( MB_SUCCESS != result ) return result;
261 
262  if( add_graphics_pt )
263  // put a graphics point on that vertex too
264  result = add_graphics_point( dual_ent, avg_pos );
265 
266  return result;
267 }
268 
269 ErrorCode DualTool::add_graphics_point( EntityHandle entity, double* avg_pos )
270 {
271  // add a graphics pt, placed at the same position as the vertex
272  double my_pos[3];
273  ErrorCode result;
274 
275  if( NULL == avg_pos )
276  {
277  result = MeshTopoUtil( mbImpl ).get_average_position( entity, my_pos );
278  if( MB_SUCCESS != result ) return result;
279  }
280  else
281  for( int i = 0; i < 3; i++ )
282  my_pos[i] = avg_pos[i];
283 
284  DualTool::GraphicsPoint dum_pt( my_pos, -1 );
285  result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), &entity, 1, &dum_pt );
286  return result;
287 }
288 
289 ErrorCode DualTool::construct_dual_edges( const Range& all_faces, Range& dual_ents )
290 {
291  if( all_faces.empty() ) return MB_SUCCESS;
292 
293  // make sure they're all faces
294  assert( 2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.begin() ) ) &&
295  2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.rbegin() ) ) );
296 
297  Range::const_iterator rit;
298  EntityHandle dual_ent;
299  unsigned int is_dual = 0x1;
300  ErrorCode tmp_result = MB_SUCCESS;
301  ErrorCode result = MB_SUCCESS;
302 
303  for( rit = all_faces.begin(); rit != all_faces.end(); ++rit )
304  {
305  if( tmp_result != MB_SUCCESS ) result = tmp_result;
306 
307  tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
308  if( MB_SUCCESS == tmp_result && 0 != dual_ent )
309  {
310  dual_ents.insert( dual_ent );
311  continue;
312  }
313 
314  // no dual entity; construct one; get the bounding regions
315  std::vector< EntityHandle > out_ents;
316  tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, out_ents );
317  if( MB_SUCCESS != tmp_result || out_ents.empty() ) continue;
318 
319  // get the dual vertices
320  std::vector< EntityHandle > dual_verts( out_ents.size() );
321  tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &out_ents[0], out_ents.size(), &dual_verts[0] );
322  if( MB_SUCCESS != tmp_result ) continue;
323  assert( dual_verts.size() <= 2 );
324 
325  double avg_pos[3];
326  bool bdy_face = ( dual_verts.size() == 1 ? true : false );
327  if( bdy_face )
328  {
329  // boundary face - make a dual vertex at the face center and put in list
330  tmp_result = construct_dual_vertex( *rit, dual_ent, true, true );
331 
332  // put it on vertex list
333  dual_verts.push_back( dual_ent );
334  }
335 
336  assert( dual_verts.size() == 2 );
337 
338  // now create the dual edge
339  tmp_result = mbImpl->create_element( MBEDGE, &dual_verts[0], 2, dual_ent );
340  if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue;
341 
342  // save it in the list of new dual ents
343  dual_ents.insert( dual_ent );
344 
345  // tag the primal entity with its dual entity and vica versa
346  tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
347  if( MB_SUCCESS != tmp_result ) continue;
348 
349  tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
350  if( MB_SUCCESS != tmp_result ) continue;
351 
352  // tag the edge indicating it's a dual entity
353  tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
354  if( MB_SUCCESS != tmp_result ) continue;
355 
356  // add a graphics point to the edge; position depends on whether it's a
357  // bdy face (mid-pt of dual edge) or not (mid-pt of primal face)
358  if( bdy_face )
359  tmp_result = add_graphics_point( dual_ent );
360  else
361  {
362  // get the face's position
363  tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos );
364  if( MB_SUCCESS != tmp_result ) continue;
365  tmp_result = add_graphics_point( dual_ent, avg_pos );
366  }
367  if( MB_SUCCESS != tmp_result ) continue;
368  }
369 
370  return result;
371 }
372 
373 ErrorCode DualTool::construct_dual_faces( const Range& all_edges, Range& dual_ents )
374 {
375  if( all_edges.empty() ) return MB_SUCCESS;
376 
377  // make sure they're all edges
378  assert( 1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.begin() ) ) &&
379  1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.rbegin() ) ) );
380 
381  Range::const_iterator rit;
382  EntityHandle dual_ent;
383  unsigned int is_dual = 0x1;
384  ErrorCode tmp_result = MB_SUCCESS;
385  ErrorCode result = MB_SUCCESS;
386  Range equiv_edges;
387 #define TRC \
388  if( MB_SUCCESS != tmp_result ) \
389  { \
390  result = tmp_result; \
391  continue; \
392  }
393  for( rit = all_edges.begin(); rit != all_edges.end(); ++rit )
394  {
395 
396  tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
397  if( MB_SUCCESS == tmp_result && 0 != dual_ent )
398  {
399  dual_ents.insert( dual_ent );
400  continue;
401  }
402 
403  // no dual entity; construct one; get the dual vertices bounding the edge in radial order,
404  // then construct the dual face
405  std::vector< EntityHandle > rad_dverts;
406  bool bdy_edge;
407  tmp_result = get_radial_dverts( *rit, rad_dverts, bdy_edge );
408  TRC if( rad_dverts.empty() ) continue;
409 
410  tmp_result = mbImpl->create_element( MBPOLYGON, &rad_dverts[0], rad_dverts.size(), dual_ent );
411  TRC
412 
413  // tag it indicating it's a dual entity, and tag primal/dual with dual/primal
414  tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
415  TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
416  TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
417  TRC
418 
419  // save it in the list of new dual ents
420  dual_ents.insert( dual_ent );
421 
422  // add a graphics point to the cell; position depends on whether it's a
423  // bdy cell (mid-pt of cell's vertices) or not (mid-pt of primal edge)
424  double avg_pos[3];
425  tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos );
426  TRC if( bdy_edge )
427  {
428 
429  // add a new dual edge betw last 2 verts
430  EntityHandle new_edge;
431  tmp_result = mbImpl->create_element( MBEDGE, &rad_dverts[rad_dverts.size() - 2], 2, new_edge );
432  TRC tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &new_edge, 1, &is_dual );
433  TRC
434 
435  // tag the new dual edge with the primal edge as it's dual entity; primal
436  // edge IS NOT likewise tagged, since it's already tagged with the 2cell
437  tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &new_edge, 1, &( *rit ) );
438  TRC
439 
440  // add a graphics pt, position is center of primal edge
441  tmp_result = add_graphics_point( dual_ent );
442  TRC tmp_result = add_graphics_point( new_edge, avg_pos );
443  TRC
444  }
445 
446  else
447  {
448  // if inside, point goes on the 2cell, at primal edge mid-pt
449  tmp_result = add_graphics_point( dual_ent, avg_pos );
450  TRC
451  }
452 
453  // check to see whether we have equiv entities; if we find any, save for later fixup
454  Range dum_edges, dum_poly( dual_ent, dual_ent );
455  tmp_result = mbImpl->get_adjacencies( dum_poly, 1, false, dum_edges );
456  if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result )
457  {
458  // we do - need to add adjacencies to disambiguate; use the primal
459  equiv_edges.merge( dum_edges );
460  }
461  }
462 
463  if( !equiv_edges.empty() ) result = check_dual_equiv_edges( equiv_edges );
464 
465  return result;
466 }
467 
468 ErrorCode DualTool::check_dual_equiv_edges( Range& dual_edges )
469 {
470  // fix equivalent dual edges (i.e. edges whose vertices define multiple edges)
471  // by explicitly adding adjacencies to containing polygons; adjacent polygons
472  // found by going through primal
473  ErrorCode tmp_result, result = MB_SUCCESS;
474 
475  Range all_dedges( dual_edges );
476  // first, go through all dual edges and find equivalent edges (by looking for
477  // up-adjacent edges on the vertices of each edge)
478  for( Range::iterator rit = dual_edges.begin(); rit != dual_edges.end(); ++rit )
479  {
480  Range connect, dum_range( *rit, *rit );
481  tmp_result = mbImpl->get_adjacencies( dum_range, 0, false, connect );
482  if( MB_SUCCESS != tmp_result ) continue;
483  tmp_result = mbImpl->get_adjacencies( connect, 1, false, all_dedges, Interface::UNION );
484  if( MB_SUCCESS != tmp_result ) continue;
485  }
486 
487  // save a copy for checking later
488  Range save_all_2cells;
489 
490  // go through each edge
491  while( !all_dedges.empty() )
492  {
493  EntityHandle this_edge = *all_dedges.begin();
494  all_dedges.erase( all_dedges.begin() );
495 
496  const EntityHandle* connect;
497  int num_connect;
498  result = mbImpl->get_connectivity( this_edge, connect, num_connect );
499  if( MB_SUCCESS != result ) continue;
500 
501  Range dum_edges, verts;
502  verts.insert( connect[0] );
503  verts.insert( connect[1] );
504  tmp_result = mbImpl->get_adjacencies( verts, 1, false, dum_edges );
505  if( MB_SUCCESS != tmp_result )
506  {
507  result = tmp_result;
508  continue;
509  }
510  if( dum_edges.size() == 1 )
511  {
512  // not an equiv edge - already removed from list, so just continue
513  continue;
514  }
515 
516  // ok, have an equiv entity - fix by looking through primal
517  // pre-get the primal of these
518  EntityHandle dedge_quad;
519  tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &this_edge, 1, &dedge_quad );
520  if( MB_SUCCESS != tmp_result )
521  {
522  result = tmp_result;
523  continue;
524  }
525 
526  if( MBQUAD == mbImpl->type_from_handle( dedge_quad ) )
527  {
528 
529  // get the primal edges adj to quad
530  Range dum_quad_range( dedge_quad, dedge_quad ), adj_pedges;
531  tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_pedges );
532  if( MB_SUCCESS != tmp_result )
533  {
534  result = tmp_result;
535  continue;
536  }
537  // get the dual 2cells corresponding to those pedges
538  std::vector< EntityHandle > dcells;
539  dcells.resize( adj_pedges.size() );
540  tmp_result = mbImpl->tag_get_data( dualEntity_tag(), adj_pedges, &dcells[0] );
541  if( MB_SUCCESS != tmp_result )
542  {
543  result = tmp_result;
544  continue;
545  }
546  // now add explicit adjacencies from the dedge to those dcells
547  std::vector< EntityHandle >::iterator vit;
548  for( vit = dcells.begin(); vit != dcells.end(); ++vit )
549  {
550  save_all_2cells.insert( *vit );
551 
552  assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) );
553  tmp_result = mbImpl->add_adjacencies( this_edge, &( *vit ), 1, false );
554  if( MB_SUCCESS != tmp_result )
555  {
556  result = tmp_result;
557  continue;
558  }
559  // check that there are really adjacencies and *vit is in them
560  const EntityHandle* adjs;
561  int num_adjs;
562  tmp_result = reinterpret_cast< Core* >( mbImpl )->a_entity_factory()->get_adjacencies( this_edge, adjs,
563  num_adjs );
564  if( NULL == adjs || std::find( adjs, adjs + num_adjs, *vit ) == adjs + num_adjs )
565  std::cout << "Add_adjacencies failed in construct_dual_faces." << std::endl;
566  }
567  }
568  else
569  {
570  // else, have a dual edge representing a bdy edge - tie directly to
571  // dual entity if its dual entity
572  EntityHandle bdy_dcell;
573  tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &dedge_quad, 1, &bdy_dcell );
574  TRC assert( MBPOLYGON == mbImpl->type_from_handle( bdy_dcell ) );
575 
576  tmp_result = mbImpl->add_adjacencies( this_edge, &bdy_dcell, 1, false );
577  if( MB_SUCCESS != tmp_result )
578  {
579  result = tmp_result;
580  continue;
581  }
582  }
583  }
584 
585  // sanity check - look for adj edges again, and check for equiv entities
586  for( Range::iterator vit = save_all_2cells.begin(); vit != save_all_2cells.end(); ++vit )
587  {
588  Range adj_edges, dum_quad_range;
589  dum_quad_range.insert( *vit );
590  assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) );
591  tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_edges );
592  if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result )
593  {
594  std::cout << "Multiple entities returned for polygon " << mbImpl->id_from_handle( *vit ) << "."
595  << std::endl;
596  continue;
597  }
598  }
599  // success!
600  return result;
601 }
602 
603 ErrorCode DualTool::construct_dual_cells( const Range& all_verts, Range& dual_ents )
604 {
605  if( all_verts.empty() ) return MB_SUCCESS;
606 
607  // make sure they're all edges
608  assert( 0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.begin() ) ) &&
609  0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.rbegin() ) ) );
610 
611  Range::const_iterator rit;
612  EntityHandle dual_ent;
613  unsigned int is_dual = 0x1;
614  ErrorCode tmp_result = MB_SUCCESS;
615  ErrorCode result = MB_SUCCESS;
616  std::vector< EntityHandle > edges, dfaces;
617 
618  for( rit = all_verts.begin(); rit != all_verts.end(); ++rit )
619  {
620  if( tmp_result != MB_SUCCESS ) result = tmp_result;
621 
622  tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
623  if( MB_SUCCESS == tmp_result && 0 != dual_ent )
624  {
625  dual_ents.insert( dual_ent );
626  continue;
627  }
628 
629  // no dual entity; construct one; get the edges bounding the vertex
630  edges.clear();
631  dfaces.clear();
632  tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 1, false, edges );
633  if( MB_SUCCESS != tmp_result ) continue;
634 
635  // get the dual faces corresponding to the edges
636  dfaces.resize( edges.size() );
637  tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &edges[0], edges.size(), &dfaces[0] );
638  if( MB_SUCCESS != tmp_result ) continue;
639 
640  // create the dual cell from those faces
641  tmp_result = mbImpl->create_element( MBPOLYHEDRON, &dfaces[0], dfaces.size(), dual_ent );
642  if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue;
643 
644  // save it in the list of new dual ents
645  dual_ents.insert( dual_ent );
646 
647  // tag it indicating it's a dual entity
648  tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
649  if( MB_SUCCESS != tmp_result ) continue;
650 
651  // tag the primal entity with its dual entity and vica versa
652  tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
653  if( MB_SUCCESS != tmp_result ) continue;
654  tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
655  if( MB_SUCCESS != tmp_result ) continue;
656  }
657 
658  return result;
659 }
660 
661 //! given an edge handle, return a list of dual vertices in radial order
662 //! around the edge
663 ErrorCode DualTool::get_radial_dverts( const EntityHandle edge,
664  std::vector< EntityHandle >& rad_dverts,
665  bool& bdy_edge )
666 {
667  rad_dverts.clear();
668 
669  std::vector< EntityHandle > rad_faces, rad_ents;
670  ErrorCode result = MeshTopoUtil( mbImpl ).star_entities( edge, rad_faces, bdy_edge, 0, &rad_ents );
671  if( MB_SUCCESS != result ) return result;
672 
673  if( bdy_edge )
674  {
675  // if we're a bdy edge, change the order back to what DualTool expects
676  rad_ents.push_back( *rad_faces.rbegin() );
677  rad_ents.push_back( *rad_faces.begin() );
678  }
679 
680  rad_dverts.resize( rad_ents.size() );
681  for( unsigned int i = 0; i < rad_ents.size(); i++ )
682  {
683  EntityHandle dual_ent;
684  result = mbImpl->tag_get_data( dualEntity_tag(), &rad_ents[i], 1, &dual_ent );
685  if( !bdy_edge || i < rad_ents.size() - 2 )
686  rad_dverts[i] = dual_ent;
687  else
688  {
689  // fix up this entry
690  assert( mbImpl->type_from_handle( dual_ent ) == MBEDGE );
691 
692  // get connectivity of that edge
693  const EntityHandle* connect;
694  int num_connect;
695  result = mbImpl->get_connectivity( dual_ent, connect, num_connect );
696  if( MB_SUCCESS != result ) return result;
697 
698  // we want the one that's not already on the list; reuse last_face
699  int last_hex = ( i == rad_ents.size() - 1 ? 0 : i - 1 );
700  EntityHandle last_face = ( connect[0] == rad_dverts[last_hex] ? connect[1] : connect[0] );
701  rad_dverts[i] = last_face;
702  }
703  }
704 
705  return result;
706 }
707 
708 //! construct the dual entities for a hex mesh, including dual surfaces & curves
709 ErrorCode DualTool::construct_hex_dual( Range& entities )
710 {
711  std::vector< EntityHandle > evec;
712  std::copy( entities.begin(), entities.end(), std::back_inserter( evec ) );
713  return construct_hex_dual( &evec[0], evec.size() );
714 }
715 
716 //! construct the dual entities for a hex mesh, including dual surfaces & curves
717 ErrorCode DualTool::construct_hex_dual( EntityHandle* entities, const int num_entities )
718 {
719  // really quite simple:
720 
721  // construct the dual...
722  ErrorCode result = construct_dual( entities, num_entities );
723  if( MB_SUCCESS != result )
724  {
725  std::cerr << "Error constructing dual entities for primal entities." << std::endl;
726  return result;
727  }
728 
729  // now traverse to build 1d and 2d hyperplanes
730  result = construct_dual_hyperplanes( 1, entities, num_entities );
731  if( MB_SUCCESS != result )
732  {
733  std::cerr << "Problem traversing 1d hyperplanes." << std::endl;
734  return result;
735  }
736 
737  result = construct_dual_hyperplanes( 2, entities, num_entities );
738  if( MB_SUCCESS != result )
739  {
740  std::cerr << "Problem traversing 2d hyperplanes." << std::endl;
741  return result;
742  }
743 
744  result = construct_hp_parent_child();
745  if( MB_SUCCESS != result )
746  {
747  std::cerr << "Problem constructing parent/child relations between hyperplanes." << std::endl;
748  return result;
749  }
750 
751  // see? simple, just like I said
752  return MB_SUCCESS;
753 }
754 
755 //! get the cells of the dual
756 ErrorCode DualTool::get_dual_entities( const int dim, EntityHandle* entities, const int num_entities, Range& dual_ents )
757 {
758  if( 0 == isDualCell_tag() ) return MB_SUCCESS;
759  if( 0 > dim || 3 < dim ) return MB_INDEX_OUT_OF_RANGE;
760 
761  unsigned int dum = 0x1;
762  const void* dum_ptr = &dum;
763  static EntityType dual_type[] = { MBVERTEX, MBEDGE, MBPOLYGON, MBPOLYHEDRON };
764 
765  Range dim_ents;
766 
767  ErrorCode result;
768 
769  if( 0 == entities || 0 == num_entities )
770  {
771  // just get all the dual entities of this dimension
772  result = mbImpl->get_entities_by_type_and_tag( 0, dual_type[dim], &isDualCellTag, &dum_ptr, 1, dual_ents );
773  }
774  else
775  {
776  // else look for specific dual entities
777  result = mbImpl->get_adjacencies( entities, num_entities, 3 - dim, false, dim_ents, Interface::UNION );
778  if( MB_SUCCESS != result ) return result;
779  std::vector< EntityHandle > dual_ents_vec( dim_ents.size() );
780  result = mbImpl->tag_get_data( dualEntity_tag(), dim_ents, &dual_ents_vec[0] );
781  if( MB_SUCCESS != result ) return result;
782  std::copy( dual_ents_vec.begin(), dual_ents_vec.end(), range_inserter( dual_ents ) );
783  }
784 
785  return result;
786 }
787 
788 //! get the faces of the dual
789 ErrorCode DualTool::get_dual_entities( const int dim,
790  EntityHandle* entities,
791  const int num_entities,
792  std::vector< EntityHandle >& dual_ents )
793 {
794  Range tmp_range;
795  ErrorCode result = get_dual_entities( dim, entities, num_entities, tmp_range );
796  if( MB_SUCCESS != result ) return result;
797 
798  // dual_ents.insert(dual_ents.end(), tmp_range.begin(), tmp_range.end());
799  dual_ents.reserve( dual_ents.size() + tmp_range.size() );
800  for( Range::const_iterator it = tmp_range.begin(); it != tmp_range.end(); ++it )
801  {
802  dual_ents.push_back( *it );
803  }
804  return MB_SUCCESS;
805 }
806 
807 ErrorCode DualTool::get_dual_hyperplanes( const Interface* impl, const int dim, Range& dual_ents )
808 {
809  if( dim != 1 && dim != 2 ) return MB_INDEX_OUT_OF_RANGE;
810 
811  Tag dual_tag;
812  ErrorCode result;
813 
814  if( dim == 1 )
815  result = impl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag );
816  else
817  result = impl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag );
818 
819  if( MB_SUCCESS == result )
820  result = impl->get_entities_by_type_and_tag( 0, MBENTITYSET, &dual_tag, NULL, 1, dual_ents, Interface::UNION );
821 
822  return result;
823 }
824 
825 ErrorCode DualTool::construct_dual_hyperplanes( const int dim, EntityHandle* entities, const int num_entities )
826 {
827  // this function traverses dual faces of input dimension, constructing
828  // dual hyperplanes of them in sets as it goes
829 
830  // check various inputs
831  int num_quads, num_hexes;
832  if(
833  // this function only makes sense for dim == 1 or dim == 2
834  ( dim != 1 && dim != 2 ) ||
835  // should either be quads or hexes around
836  mbImpl->get_number_entities_by_type( 0, MBQUAD, num_quads ) != MB_SUCCESS ||
837  mbImpl->get_number_entities_by_type( 0, MBHEX, num_hexes ) != MB_SUCCESS ||
838  // if we're asking for 1d dual ents, should be quads around
839  ( num_quads == 0 && dim == 1 ) ||
840  // if we're asking for 2d dual ents, should be hexes around
841  ( num_hexes == 0 && dim == 2 ) )
842  return MB_FAILURE;
843 
844  Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() );
845 
846  // two stacks: one completely untreated entities, and the other untreated
847  // entities on the current dual hyperplane
848  std::vector< EntityHandle > tot_untreated;
849 
850  // put dual entities of this dimension on the untreated list
851  ErrorCode result = get_dual_entities( dim, entities, num_entities, tot_untreated );
852  if( MB_SUCCESS != result ) return result;
853 
854  // main part of traversal loop
855  EntityHandle this_ent;
856  EntityHandle this_hp;
857 
858  while( !tot_untreated.empty() )
859  {
860  if( debug && dim == 2 /*(tot_untreated.size()%report == 0)*/ )
861  std::cout << "Untreated list size " << tot_untreated.size() << "." << std::endl;
862 
863  this_ent = tot_untreated.back();
864  tot_untreated.pop_back();
865  result = mbImpl->tag_get_data( hp_tag, &this_ent, 1, &this_hp );
866  if( MB_SUCCESS != result && MB_TAG_NOT_FOUND != result ) return result;
867 
868  // d for this entity having a hyperplane assignment already
869  else if( this_hp != 0 )
870  continue;
871 
872  if( 1 == dim && check_1d_loop_edge( this_ent ) ) continue;
873 
874  // inner loop: traverse the hyperplane 'till we don't have any more
875  result = traverse_hyperplane( hp_tag, this_hp, this_ent );
876  if( MB_SUCCESS != result )
877  {
878  std::cout << "Failed to traverse hyperplane ";
879  if( this_hp )
880  std::cout << mbImpl->id_from_handle( this_hp ) << "." << std::endl;
881  else
882  std::cout << "0." << std::endl;
883  return result;
884  }
885 
886  // ok, now order the edges if it's a chord
887  if( 1 == dim ) order_chord( this_hp );
888  }
889 
890  return MB_SUCCESS;
891 }
892 
893 ErrorCode DualTool::traverse_hyperplane( const Tag hp_tag, EntityHandle& this_hp, EntityHandle this_ent )
894 {
895  Range tmp_star, star, tmp_range, new_hyperplane_ents;
896  std::vector< EntityHandle > hp_untreated;
897  int dim = mbImpl->dimension_from_handle( this_ent );
898  MeshTopoUtil mtu( mbImpl );
899  this_hp = 0;
900  ErrorCode result;
901 
902  unsigned short mark_val = 0x0;
903  Tag mark_tag;
904  result = mbImpl->tag_get_handle( "__hyperplane_mark", 1, MB_TYPE_BIT, mark_tag, MB_TAG_CREAT | MB_TAG_BIT );
905  if( MB_SUCCESS != result ) return result;
906  mark_val = 0x1;
907 
908  while( 0 != this_ent )
909  {
910  EntityHandle tmp_hp = get_dual_hyperplane( this_ent );
911  if( 0 == this_hp && 0 != tmp_hp ) this_hp = tmp_hp;
912 
913  if( 0 == tmp_hp ) new_hyperplane_ents.insert( this_ent );
914 
915  if( debug && hp_untreated.size() % 10 == 0 )
916  std::cout << "Dual surface " << this_hp << ", hp_untreated list size = " << hp_untreated.size() << "."
917  << std::endl;
918 
919  // get the 2nd order adjacencies through lower dimension
920  tmp_range.clear();
921  tmp_star.clear();
922  star.clear();
923  result = mtu.get_bridge_adjacencies( this_ent, dim - 1, dim, star );RR;
924 
925  // get the bridge adjacencies through higher dimension
926  result = mtu.get_bridge_adjacencies( this_ent, dim + 1, dim, tmp_star );RR;
927  tmp_range = subtract( star, tmp_star );
928 
929  for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit )
930  {
931  if( new_hyperplane_ents.find( *rit ) != new_hyperplane_ents.end() ) continue;
932 
933  // check for tag first, 'cuz it's probably faster than checking adjacencies
934  // assign to avoid valgrind warning
935  unsigned short tmp_mark = 0x0;
936  result = mbImpl->tag_get_data( mark_tag, &( *rit ), 1, &tmp_mark );
937  if( MB_SUCCESS == result && mark_val == tmp_mark ) continue;
938 
939  // if it's on the loop, it's not eligible
940  if( 1 == dim && check_1d_loop_edge( *rit ) ) continue;
941 
942  // have one on this hp; just put it on the hp_untreated list for now,
943  // will get tagged and put in the hp set later
944  hp_untreated.push_back( *rit );
945  result = mbImpl->tag_set_data( mark_tag, &( *rit ), 1, &mark_val );
946  if( MB_SUCCESS != result ) return result;
947  }
948 
949  // end of inner loop; get the next this_ent, or set to zero
950  if( hp_untreated.empty() )
951  this_ent = 0;
952  else
953  {
954  this_ent = hp_untreated.back();
955  hp_untreated.pop_back();
956  }
957  }
958 
959  if( debug_ap )
960  {
961  std::string hp_name;
962  if( 2 == dim )
963  hp_name = "sheet";
964  else
965  hp_name = "chord";
966 
967  if( 0 == this_hp )
968  std::cout << "Constructed new " << hp_name << " with ";
969  else
970  {
971  int this_id;
972  result = mbImpl->tag_get_data( globalId_tag(), &this_hp, 1, &this_id );RR;
973  std::cout << "Added to " << hp_name << " " << this_id << " ";
974  }
975  if( dim == 2 )
976  std::cout << "edges:" << std::endl;
977  else
978  std::cout << "quads:" << std::endl;
979  std::vector< EntityHandle > pents( new_hyperplane_ents.size() );
980  result = mbImpl->tag_get_data( dualEntity_tag(), new_hyperplane_ents, &pents[0] );RR;
981  for( std::vector< EntityHandle >::iterator vit = pents.begin(); vit != pents.end(); ++vit )
982  {
983  if( vit != pents.begin() ) std::cout << ", ";
984  std::cout << mbImpl->id_from_handle( *vit );
985  }
986  std::cout << std::endl;
987  }
988 
989  if( 0 == this_hp )
990  {
991  // ok, doesn't have one; make a new hyperplane
992  int new_id = -1;
993  result = construct_new_hyperplane( dim, this_hp, new_id );
994  if( MB_SUCCESS != result ) return result;
995 
996  if( debug_ap )
997  {
998  std::cout << "New ";
999  if( 2 == dim )
1000  std::cout << " sheet ";
1001  else
1002  std::cout << " chord ";
1003  std::cout << new_id << " constructed." << std::endl;
1004  }
1005  }
1006 
1007  // set the hp_val for entities which didn't have one before
1008  std::vector< EntityHandle > hp_tags( new_hyperplane_ents.size() );
1009  std::fill( hp_tags.begin(), hp_tags.end(), this_hp );
1010  result = mbImpl->tag_set_data( hp_tag, new_hyperplane_ents, &hp_tags[0] );
1011  if( MB_SUCCESS != result ) return result;
1012  result = mbImpl->add_entities( this_hp, new_hyperplane_ents );
1013  if( MB_SUCCESS != result ) return result;
1014 
1015  // unmark the entities by removing the tag
1016  result = mbImpl->tag_delete( mark_tag );
1017  if( MB_SUCCESS != result ) return result;
1018 
1019  return MB_SUCCESS;
1020 }
1021 
1022 ErrorCode DualTool::order_chord( EntityHandle chord_set )
1023 {
1024  // re-order the 1cells in the set so they are in order along the chord
1025  // start by finding the vertex dual to a quad
1026  Range verts, one_cells;
1027  ErrorCode result = mbImpl->get_entities_by_dimension( chord_set, 1, one_cells );
1028  if( MB_SUCCESS != result || one_cells.empty() ) return MB_FAILURE;
1029 
1030  result = mbImpl->get_adjacencies( one_cells, 0, false, verts, Interface::UNION );
1031  if( MB_SUCCESS != result || verts.empty() ) return MB_FAILURE;
1032 
1033  EntityHandle last_vert = 0;
1034  for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit )
1035  {
1036  if( TYPE_FROM_HANDLE( get_dual_entity( *rit ) ) == MBQUAD )
1037  {
1038  last_vert = *rit;
1039  break;
1040  }
1041  }
1042  // if there's no vertex owned by a quad, just start with 1st one
1043  if( 0 == last_vert ) last_vert = *verts.begin();
1044 
1045  // now, skip from vertex to vertex, building a list of 1cells
1046  std::vector< EntityHandle > ordered_1cells;
1047  EntityHandle last_1cell = 0;
1048  Range dum1, dum2;
1049  const EntityHandle* connect;
1050  int num_connect;
1051  ErrorCode tmp_result = MB_SUCCESS;
1052  while( ordered_1cells.size() != one_cells.size() )
1053  {
1054  dum1 = one_cells;
1055  result = mbImpl->get_adjacencies( &last_vert, 1, 1, false, dum1 );
1056  if( 0 != last_1cell ) dum1.erase( last_1cell );
1057  // assert(1 == dum1.size());
1058  if( 0 != last_1cell && 1 != dum1.size() )
1059  {
1060  std::cerr << "unexpected size traversing chord." << std::endl;
1061  tmp_result = MB_FAILURE;
1062  }
1063 
1064  last_1cell = *dum1.begin();
1065  ordered_1cells.push_back( last_1cell );
1066  result = mbImpl->get_connectivity( last_1cell, connect, num_connect );RR;
1067  if( last_vert == connect[0] )
1068  last_vert = connect[1];
1069  else
1070  last_vert = connect[0];
1071  }
1072 
1073  // now have the 1cells in order, replace them in the set
1074  if( MB_SUCCESS == tmp_result )
1075  {
1076  result = mbImpl->remove_entities( chord_set, one_cells );RR;
1077  result = mbImpl->add_entities( chord_set, &ordered_1cells[0], ordered_1cells.size() );RR;
1078  }
1079 
1080  return MB_SUCCESS;
1081 }
1082 
1083 ErrorCode DualTool::construct_new_hyperplane( const int dim, EntityHandle& new_hyperplane, int& id )
1084 {
1085  ErrorCode result;
1086  if( 1 == dim )
1087  result = mbImpl->create_meshset( ( MESHSET_ORDERED | MESHSET_TRACK_OWNER ), new_hyperplane );
1088  else
1089  result = mbImpl->create_meshset( ( MESHSET_SET | MESHSET_TRACK_OWNER ), new_hyperplane );
1090  if( MB_SUCCESS != result ) return result;
1091 
1092  if( -1 == id )
1093  {
1094  Range all_hyperplanes;
1095  result = get_dual_hyperplanes( mbImpl, dim, all_hyperplanes );RR;
1096  std::vector< int > gids( all_hyperplanes.size() );
1097  result = mbImpl->tag_get_data( globalIdTag, all_hyperplanes, ( gids.empty() ) ? NULL : &gids[0] );RR;
1098  for( unsigned int i = 0; i < gids.size(); i++ )
1099  if( gids[i] > id ) id = gids[i];
1100  id++;
1101  if( 0 == id ) id++;
1102  }
1103 
1104  result = mbImpl->tag_set_data( globalId_tag(), &new_hyperplane, 1, &id );RR;
1105  Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() );
1106  result = mbImpl->tag_set_data( hp_tag, &new_hyperplane, 1, &new_hyperplane );RR;
1107 
1108  // assign a category name to these sets
1109  static const char dual_category_names[2][CATEGORY_TAG_SIZE] = { "Chord\0", "Sheet\0" };
1110 
1111  result = mbImpl->tag_set_data( categoryTag, &new_hyperplane, 1, dual_category_names[dim - 1] );
1112 
1113  return result;
1114 }
1115 
1116 bool DualTool::check_1d_loop_edge( EntityHandle this_ent )
1117 {
1118  // make sure it's an edge
1119  if( MBEDGE != mbImpl->type_from_handle( this_ent ) ) return false;
1120 
1121  // also has to be a dual entity
1122  unsigned int dum;
1123  ErrorCode result = mbImpl->tag_get_data( isDualCell_tag(), &this_ent, 1, &dum );
1124  if( MB_SUCCESS != result || dum != 0x1 ) return false;
1125 
1126  const EntityHandle* verts;
1127  EntityHandle vert_tags[2];
1128  int num_verts;
1129  result = mbImpl->get_connectivity( this_ent, verts, num_verts );
1130  if( MB_SUCCESS != result ) return false;
1131 
1132  result = mbImpl->tag_get_data( dualEntity_tag(), verts, 2, vert_tags );
1133  if( MB_SUCCESS != result || mbImpl->type_from_handle( vert_tags[0] ) != MBQUAD ||
1134  mbImpl->type_from_handle( vert_tags[1] ) != MBQUAD )
1135  return false;
1136 
1137  else
1138  return true;
1139 }
1140 
1141 ErrorCode DualTool::construct_hp_parent_child()
1142 {
1143  Range dual_surfs, dual_cells, dual_edges;
1144  ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs );
1145  if( MB_SUCCESS != result || dual_surfs.empty() ) return result;
1146  std::vector< EntityHandle > dual_curve_sets;
1147 
1148  for( Range::iterator surf_it = dual_surfs.begin(); surf_it != dual_surfs.end(); ++surf_it )
1149  {
1150  // get all the cells, edges in those cells, and chords for those edges
1151  dual_cells.clear();
1152  result = mbImpl->get_entities_by_handle( *surf_it, dual_cells );
1153  if( MB_SUCCESS != result ) return result;
1154  dual_edges.clear();
1155  result = mbImpl->get_adjacencies( dual_cells, 1, false, dual_edges, Interface::UNION );
1156  if( MB_SUCCESS != result ) return result;
1157  dual_curve_sets.resize( dual_edges.size() );
1158  result = mbImpl->tag_get_data( dualCurve_tag(), dual_edges, &dual_curve_sets[0] );
1159  if( MB_SUCCESS != result ) return result;
1160 
1161  // reuse dual_cells to get unique list of chord sets
1162  dual_cells.clear();
1163  for( unsigned int i = 0; i < dual_edges.size(); i++ )
1164  if( dual_curve_sets[i] != 0 ) dual_cells.insert( dual_curve_sets[i] );
1165 
1166  // now connect up this dual surf with all the 1d ones
1167  for( Range::iterator rit = dual_cells.begin(); rit != dual_cells.end(); ++rit )
1168  {
1169  result = mbImpl->add_parent_child( *surf_it, *rit );
1170  if( MB_SUCCESS != result ) return result;
1171  }
1172  }
1173 
1174  return MB_SUCCESS;
1175 }
1176 
1177 ErrorCode DualTool::get_graphics_points( EntityHandle dual_ent,
1178  std::vector< int >& npts,
1179  std::vector< GraphicsPoint >& points )
1180 {
1181  // shouldn't be a set
1182  assert( MBENTITYSET != mbImpl->type_from_handle( dual_ent ) );
1183 
1184  // get the graphics points comprising the given entity
1185  GraphicsPoint gp_array[DualTool::GP_SIZE];
1186 
1187  ErrorCode result = MB_SUCCESS;
1188 
1189  // switch based on topological dimension
1190  switch( mbImpl->dimension_from_handle( dual_ent ) )
1191  {
1192  case 0:
1193  // just return the vertex point
1194  result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array );
1195  if( MB_SUCCESS == result ) points.push_back( gp_array[0] );
1196 
1197  break;
1198 
1199  case 1:
1200  // get my graphics point then those of my vertices
1201  const EntityHandle* connect;
1202  int num_connect;
1203  result = mbImpl->get_connectivity( dual_ent, connect, num_connect );
1204  if( MB_SUCCESS != result ) break;
1205 
1206  result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, gp_array );
1207  if( MB_SUCCESS == result )
1208  {
1209  points.push_back( gp_array[0] );
1210  points.push_back( gp_array[0] );
1211  points.push_back( gp_array[1] );
1212  result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array );
1213  if( MB_SUCCESS == result ) points[1] = gp_array[0];
1214  }
1215 
1216  npts.push_back( 3 );
1217 
1218  break;
1219 
1220  case 2:
1221  result = get_cell_points( dual_ent, npts, points );
1222  break;
1223  }
1224 
1225  return result;
1226 }
1227 
1228 ErrorCode DualTool::get_cell_points( EntityHandle dual_ent,
1229  std::vector< int >& npts,
1230  std::vector< GraphicsPoint >& points )
1231 {
1232  assert( MBPOLYGON == mbImpl->type_from_handle( dual_ent ) );
1233 
1234  // get the 1cells in this 2cell
1235  Range one_cells;
1236 
1237  Range tc_range;
1238  tc_range.insert( dual_ent );
1239  ErrorCode result = mbImpl->get_adjacencies( tc_range, 1, false, one_cells, Interface::UNION );RR;
1240 
1241  int num_edges = one_cells.size();
1242  std::vector< GraphicsPoint > dum_gps( num_edges + 1 );
1243 
1244  // get graphics points for 0cells and for this cell
1245  result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), one_cells, &dum_gps[0] );RR;
1246  result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, &( dum_gps[num_edges] ) );RR;
1247 
1248  Range::iterator eit;
1249  const EntityHandle* connect;
1250  int num_connect;
1251  GraphicsPoint vert_gps[2];
1252  int i;
1253  for( i = 0, eit = one_cells.begin(); i < num_edges; i++, ++eit )
1254  {
1255  // get the vertices and the graphics points for them
1256  result = mbImpl->get_connectivity( *eit, connect, num_connect );RR;
1257  result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, vert_gps );RR;
1258 
1259  // make the 2 tris corresponding to this edge; don't worry about order
1260  // for now
1261  npts.push_back( 3 );
1262  points.push_back( dum_gps[num_edges] );
1263  points.push_back( vert_gps[0] );
1264  points.push_back( dum_gps[i] );
1265 
1266  npts.push_back( 3 );
1267  points.push_back( dum_gps[num_edges] );
1268  points.push_back( dum_gps[i] );
1269  points.push_back( vert_gps[1] );
1270  }
1271 
1272  return result;
1273 }
1274 
1275 ErrorCode DualTool::get_graphics_points( const Range& in_range,
1276  std::vector< GraphicsPoint >& points,
1277  const bool assign_ids,
1278  const int start_id )
1279 {
1280  // return graphics points on dual entities in in_range or in entities
1281  // in sets in in_range
1282  ErrorCode result;
1283 
1284  // for each dual hyperplane set:
1285  Range::const_iterator rit;
1286 
1287  Range two_cells, all_cells;
1288  for( rit = in_range.begin(); rit != in_range.end(); ++rit )
1289  {
1290  // for each entity:
1291  two_cells.clear();
1292  EntityType this_type = mbImpl->type_from_handle( *rit );
1293  if( MBENTITYSET == this_type )
1294  {
1295  result = mbImpl->get_entities_by_handle( *rit, two_cells );RR;
1296 
1297  std::copy( two_cells.begin(), two_cells.end(), range_inserter( all_cells ) );
1298  }
1299 
1300  else
1301  {
1302  two_cells.insert( *rit );
1303  assert( this_type == MBVERTEX || this_type == MBEDGE || this_type == MBPOLYGON ||
1304  this_type == MBPOLYHEDRON );
1305  }
1306 
1307  result = mbImpl->get_adjacencies( two_cells, 0, false, all_cells, Interface::UNION );RR;
1308  result = mbImpl->get_adjacencies( two_cells, 1, false, all_cells, Interface::UNION );RR;
1309  }
1310 
1311  // get graphics points
1312  points.resize( all_cells.size() );
1313 
1314  result = mbImpl->tag_get_data( dualGraphicsPointTag, all_cells, &points[0] );RR;
1315 
1316  if( assign_ids )
1317  {
1318  int i = start_id;
1319 
1320  for( std::vector< GraphicsPoint >::iterator vit = points.begin(); vit != points.end(); ++vit )
1321  vit->id = i++;
1322 
1323  result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), all_cells, &points[0] );RR;
1324  }
1325 
1326  return result;
1327 }
1328 
1329 EntityHandle DualTool::next_loop_vertex( const EntityHandle last_v,
1330  const EntityHandle this_v,
1331  const EntityHandle dual_surf )
1332 {
1333  // given two vertices, find the next one on the loop; if one is a dual
1334  // surface, then just choose either one for that surface
1335  assert( ( 0 == last_v || mbImpl->type_from_handle( last_v ) == MBVERTEX ) &&
1336  mbImpl->type_from_handle( this_v ) == MBVERTEX && mbImpl->type_from_handle( dual_surf ) == MBENTITYSET );
1337 
1338  // get the connected vertices
1339  MeshTopoUtil tpu( mbImpl );
1340  Range other_verts;
1341  ErrorCode result = tpu.get_bridge_adjacencies( this_v, 1, 0, other_verts );
1342  if( MB_SUCCESS != result || other_verts.empty() ) return 0;
1343 
1344  // if (mbImpl->type_from_handle(last_v) == MBENTITYSET) {
1345  // dual surface, choose either; first get a 2cell on this surface
1346  Range tcells, tcells2, verts;
1347  result = mbImpl->get_entities_by_type( dual_surf, MBPOLYGON, tcells );
1348  if( MB_SUCCESS != result || tcells.empty() ) return 0;
1349 
1350  // ok, pay attention here: first get 2cells common to dual surface and this_v
1351  verts.insert( this_v );
1352  result = mbImpl->get_adjacencies( verts, 2, false, tcells );
1353  if( MB_SUCCESS != result || tcells.empty() ) return 0;
1354 
1355  // next get vertices common to both 2cells and subtract from other_verts; also
1356  // remove last_v if it's non-zero
1357  verts.clear();
1358  result = mbImpl->get_adjacencies( tcells, 0, false, verts );
1359  if( MB_SUCCESS != result || verts.empty() ) return 0;
1360 
1361  Range tmp_verts = subtract( other_verts, verts );
1362  other_verts.swap( tmp_verts );
1363  if( 0 != last_v ) other_verts.erase( last_v );
1364 
1365  // now get intersection of remaining vertices and 2 2cells vertices
1366  // look at each one successively; should find one, maybe not on both
1367  tmp_verts = other_verts;
1368  Range tmp_faces( *tcells.begin(), *tcells.begin() );
1369  result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts );
1370  if( MB_SUCCESS == result && !tmp_verts.empty() ) return *tmp_verts.begin();
1371  tmp_faces.clear();
1372  tmp_faces.insert( *tcells.rbegin() );
1373  result = mbImpl->get_adjacencies( tmp_faces, 0, false, other_verts );
1374  if( MB_SUCCESS == result && !other_verts.empty() ) return *other_verts.begin();
1375 
1376  // if we got here, there isn't any
1377  return 0;
1378 }
1379 
1380 EntityHandle DualTool::get_dual_hyperplane( const EntityHandle ncell )
1381 {
1382  // get the sheet or chord it's in
1383  std::vector< EntityHandle > adj_sets;
1384  ErrorCode result = mbImpl->get_adjacencies( &ncell, 1, 4, false, adj_sets );
1385  if( MB_SUCCESS != result ) return 0;
1386 
1387  EntityHandle dum_set;
1388  for( std::vector< EntityHandle >::iterator vit = adj_sets.begin(); vit != adj_sets.end(); ++vit )
1389  {
1390  if( mbImpl->tag_get_data( dualCurve_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND ||
1391  mbImpl->tag_get_data( dualSurface_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND )
1392  return *vit;
1393  }
1394 
1395  return 0;
1396 }
1397 
1398 //! set the dual surface or curve for an entity
1399 ErrorCode DualTool::set_dual_surface_or_curve( EntityHandle entity,
1400  const EntityHandle dual_hyperplane,
1401  const int dual_entity_dimension )
1402 {
1403  if( 1 == dual_entity_dimension )
1404  mbImpl->tag_set_data( dualCurve_tag(), &entity, 1, &dual_hyperplane );
1405  else if( 2 == dual_entity_dimension )
1406  mbImpl->tag_set_data( dualSurface_tag(), &entity, 1, &dual_hyperplane );
1407  else
1408  return MB_INDEX_OUT_OF_RANGE;
1409 
1410  return MB_SUCCESS;
1411 }
1412 
1413 //! return the corresponding dual entity
1414 EntityHandle DualTool::get_dual_entity( const EntityHandle this_ent ) const
1415 {
1416  EntityHandle dual_ent;
1417  ErrorCode result = mbImpl->tag_get_data( dualEntity_tag(), &this_ent, 1, &dual_ent );
1418  if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result )
1419  return 0;
1420  else
1421  return dual_ent;
1422 }
1423 
1424 //! return the corresponding dual entity
1425 EntityHandle DualTool::get_extra_dual_entity( const EntityHandle this_ent )
1426 {
1427  EntityHandle dual_ent;
1428  ErrorCode result = mbImpl->tag_get_data( extraDualEntity_tag(), &this_ent, 1, &dual_ent );
1429  if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result )
1430  return 0;
1431  else
1432  return dual_ent;
1433 }
1434 
1435 ErrorCode DualTool::atomic_pillow( EntityHandle odedge, EntityHandle& quad1, EntityHandle& quad2 )
1436 {
1437  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1438 
1439  if( debug_ap )
1440  {
1441  Range sets;
1442  Tag ms_tag;
1443 
1444  ErrorCode result = mbImpl->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, ms_tag );
1445  if( MB_SUCCESS == result )
1446  {
1447  result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &ms_tag, NULL, 1, sets );
1448  if( MB_SUCCESS == result ) result = mbImpl->delete_entities( sets );
1449  }
1450  }
1451 
1452  std::cout << "-AP(";
1453  print_cell( odedge );
1454  std::cout << ")" << std::endl;
1455 
1456  // perform an atomic pillow operation around dedge
1457 
1458  // grab the quad before deleting the odedge
1459  quad1 = get_dual_entity( odedge );
1460  assert( 0 != quad1 );
1461 
1462  // 0. get star 2cells around odedge (before odedge changes) and 3cells around
1463  // those 2cells (going to delete the 2cells, therefore need to delete the 3cells
1464  // that depend on those too)
1465  MeshTopoUtil mtu( mbImpl );
1466  Range star_cells, tmp_cells;
1467  ErrorCode result = mbImpl->get_adjacencies( &odedge, 1, 2, false, star_cells );RR;
1468  result = mbImpl->get_adjacencies( star_cells, 3, false, tmp_cells, Interface::UNION );RR;
1469  star_cells.merge( tmp_cells );
1470  star_cells.insert( odedge );
1471 
1472  // tear down the dual entities which will be modified by the ap first
1473  result = delete_dual_entities( star_cells );RR;
1474 
1475  // now change the quad to an ap
1476  std::vector< EntityHandle > verts;
1477  result = mbImpl->get_connectivity( &quad1, 1, verts );RR;
1478 
1479  // get average position of vertices
1480  double coords[12], avg[3] = { 0.0, 0.0, 0.0 };
1481  result = mbImpl->get_coords( &verts[0], verts.size(), coords );RR;
1482  for( int i = 0; i < 4; i++ )
1483  {
1484  avg[0] += coords[3 * i];
1485  avg[1] += coords[3 * i + 1];
1486  avg[2] += coords[3 * i + 2];
1487  }
1488  for( int i = 0; i < 3; i++ )
1489  avg[i] *= 0.25;
1490 
1491  // for each position, get a corresponding position 1/2 way to avg
1492  double new_coords[12];
1493  for( int i = 0; i < 4; i++ )
1494  {
1495  new_coords[3 * i] = avg[0] + .5 * ( coords[3 * i] - avg[0] );
1496  new_coords[3 * i + 1] = avg[1] + .5 * ( coords[3 * i + 1] - avg[1] );
1497  new_coords[3 * i + 2] = avg[2] + .5 * ( coords[3 * i + 2] - avg[2] );
1498  }
1499 
1500  // make the 4 new vertices; store in vector long enough for hex connectivity
1501  for( int i = 0; i < 4; i++ )
1502  {
1503  verts.push_back( 0 );
1504  result = mbImpl->create_vertex( &new_coords[3 * i], verts[4 + i] );RR;
1505  }
1506 
1507  // get the hexes connected to the quad
1508  Range hexes;
1509  result = mbImpl->get_adjacencies( &quad1, 1, 3, false, hexes );RR;
1510  assert( hexes.size() <= 2 );
1511 
1512  // remove any explicit adjacency from the first hex, since that'll get connected
1513  // to the new quad; add adjacency between quad and second hex, if there is a 2nd
1514  result = mbImpl->remove_adjacencies( quad1, &( *hexes.begin() ), 1 );RR;
1515  if( hexes.size() == 2 )
1516  {
1517  result = mbImpl->add_adjacencies( quad1, &( *hexes.rbegin() ), 1, false );RR;
1518  }
1519 
1520  // create the new, inner quad, and make it explicitly adjacent to 1st hex;
1521  // make the connectivity of this quad same as the original one
1522  std::vector< EntityHandle > tmp_verts;
1523  std::copy( verts.begin(), verts.end(), std::back_inserter( tmp_verts ) );
1524 
1525  result = mbImpl->create_element( MBQUAD, &tmp_verts[0], 4, quad2 );RR;
1526  result = mbImpl->add_adjacencies( quad2, &( *hexes.begin() ), 1, false );RR;
1527 
1528  // reverse the connectivity of the 1st hex
1529  std::reverse( verts.begin(), verts.begin() + 4 );
1530  std::reverse( verts.begin() + 4, verts.end() );
1531 
1532  // now make two inner hexes; note connectivity array is flipped for the two hexes
1533  EntityHandle new_hexes[2];
1534  result = mbImpl->create_element( MBHEX, &verts[0], 8, new_hexes[0] );RR;
1535  result = mbImpl->create_element( MBHEX, &tmp_verts[0], 8, new_hexes[1] );RR;
1536 
1537  // set the global id tag on the new hexes
1538  int new_hex_ids[2] = { maxHexId + 1, maxHexId + 2 };
1539  maxHexId += 2;
1540  result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 2, new_hex_ids );
1541  if( MB_SUCCESS != result ) return result;
1542 
1543  // by definition, quad1 is adj to new_hexes[0]
1544  result = mbImpl->add_adjacencies( quad1, &new_hexes[0], 1, false );RR;
1545  result = mbImpl->add_adjacencies( quad2, &new_hexes[1], 1, false );RR;
1546 
1547  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1548 
1549  // now update the dual
1550  result = construct_hex_dual( &new_hexes[0], 2 );RR;
1551 
1552  // get the new dual surface, by getting one of the edges between the center
1553  // and outer vertex rings
1554  Range new_edge;
1555  verts[1] = verts[4];
1556  result = mbImpl->get_adjacencies( &verts[0], 2, 1, false, new_edge );
1557  if( MB_SUCCESS != result || new_edge.size() != 1 ) return result;
1558 
1559  return MB_SUCCESS;
1560 }
1561 
1562 //! effect reverse atomic pillow operation
1563 ErrorCode DualTool::rev_atomic_pillow( EntityHandle pillow, Range& chords )
1564 {
1565  // get the dual entities associated with elements in the pillow; go through
1566  // the elements instead of the pillow sheet so you get all of them, not just
1567  // the ones on the sheet
1568  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1569 
1570  std::cout << "-AP(";
1571  print_cell( pillow );
1572  std::cout << ")" << std::endl;
1573 
1574  Range dverts;
1575  ErrorCode result = get_dual_entities( pillow, NULL, NULL, &dverts, NULL, NULL );
1576  if( MB_SUCCESS != result ) return result;
1577  assert( 2 == dverts.size() );
1578 
1579  EntityHandle hexes[2];
1580  result = mbImpl->tag_get_data( dualEntity_tag(), dverts, hexes );RR;
1581  assert( hexes[0] != 0 && hexes[1] != 0 );
1582 
1583  std::vector< EntityHandle > dcells[4];
1584  Range pcells[4];
1585  std::copy( hexes, hexes + 2, range_inserter( pcells[3] ) );
1586  std::copy( dverts.begin(), dverts.end(), std::back_inserter( dcells[0] ) );
1587  for( int dim = 0; dim <= 2; dim++ )
1588  {
1589  result = mbImpl->get_adjacencies( hexes, 2, dim, false, pcells[dim], Interface::UNION );RR;
1590  dcells[3 - dim].resize( pcells[dim].size() );
1591  result = mbImpl->tag_get_data( dualEntity_tag(), pcells[dim], &dcells[3 - dim][0] );RR;
1592  }
1593 
1594  // delete the dual entities which are part of the original pillow
1595  result = mbImpl->delete_entities( &pillow, 1 );
1596  if( MB_SUCCESS != result ) return result;
1597 
1598  result = mbImpl->delete_entities( chords );
1599  if( MB_SUCCESS != result ) return result;
1600 
1601  for( int i = 3; i >= 0; i-- )
1602  {
1603  result = delete_dual_entities( &dcells[i][0], dcells[i].size() );RR;
1604  }
1605 
1606  // delete the primal entities inserted by the ap; be careful to get the right
1607  // faces, edges and vertices
1608  Range del_faces, del_edges, del_verts, tmp_faces, tmp_verts;
1609  // faces are the shared 5 and the 1 other one with greater handle (which
1610  // we know will be later in the range)
1611  result = mbImpl->get_adjacencies( hexes, 2, 2, false, del_faces );RR;
1612  assert( 5 == del_faces.size() );
1613  std::copy( pcells[2].begin(), pcells[2].end(), range_inserter( tmp_faces ) );
1614  tmp_faces = subtract( tmp_faces, del_faces );
1615  del_faces.insert( *tmp_faces.rbegin() );
1616  result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts );RR;
1617  std::copy( pcells[0].begin(), pcells[0].end(), range_inserter( del_verts ) );
1618  del_verts = subtract( del_verts, tmp_verts );
1619  assert( 4 == del_verts.size() );
1620  result = mbImpl->get_adjacencies( del_verts, 1, false, del_edges, Interface::UNION );RR;
1621  assert( 8 == del_edges.size() );
1622 
1623  result = mbImpl->delete_entities( hexes, 2 );RR;
1624  result = mbImpl->delete_entities( del_faces );RR;
1625  result = mbImpl->delete_entities( del_edges );RR;
1626  result = mbImpl->delete_entities( del_verts );RR;
1627 
1628  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1629 
1630  // recompute the dual for the hexes on either side of the quad affected
1631  // by the ap removal
1632  Range tmp_hexes;
1633  result = mbImpl->get_adjacencies( tmp_verts, 3, false, tmp_hexes, Interface::UNION );RR;
1634  result = construct_hex_dual( tmp_hexes );RR;
1635 
1636  return MB_SUCCESS;
1637 }
1638 
1639 ErrorCode DualTool::delete_dual_entities( EntityHandle* entities, int num_entities )
1640 {
1641  Range tmp_ents;
1642  std::copy( entities, entities + num_entities, range_inserter( tmp_ents ) );
1643  return delete_dual_entities( tmp_ents );
1644 }
1645 
1646 ErrorCode DualTool::delete_dual_entities( Range& entities )
1647 {
1648  if( entities.empty() ) return delete_whole_dual();
1649 
1650  EntityHandle null_entity = 0;
1651  ErrorCode result;
1652  Range ents_to_delete;
1653 
1654  while( !entities.empty() )
1655  {
1656  EntityHandle this_entity = entities.pop_back();
1657 
1658  // reset the primal's dual entity
1659  EntityHandle primal = get_dual_entity( this_entity );
1660  if( get_dual_entity( primal ) == this_entity )
1661  {
1662  result = mbImpl->tag_set_data( dualEntity_tag(), &primal, 1, &null_entity );RR;
1663  }
1664  EntityHandle extra = get_extra_dual_entity( primal );
1665  if( 0 != extra )
1666  {
1667  result = mbImpl->tag_set_data( extraDualEntity_tag(), &primal, 1, &null_entity );RR;
1668  }
1669 
1670  ents_to_delete.insert( this_entity );
1671 
1672  // check for extra dual entities
1673  if( mbImpl->type_from_handle( this_entity ) == MBPOLYGON )
1674  {
1675  // for 2cell, might be a loop edge
1676  Range loop_edges;
1677  result = mbImpl->get_adjacencies( &this_entity, 1, 1, false, loop_edges );
1678  for( Range::iterator rit = loop_edges.begin(); rit != loop_edges.end(); ++rit )
1679  if( check_1d_loop_edge( *rit ) ) entities.insert( *rit );
1680  }
1681  else if( extra && extra != this_entity )
1682  // just put it on the list; primal for which we're extra has already been
1683  // reset to not point to extra entity
1684  ents_to_delete.insert( extra );
1685  }
1686 
1687  // now delete the entities (sheets and chords will be updated automatically)
1688  return mbImpl->delete_entities( ents_to_delete );
1689 }
1690 
1691 void DualTool::print_cell( EntityHandle cell )
1692 {
1693  const EntityHandle* connect;
1694  int num_connect;
1695  ErrorCode result = mbImpl->get_connectivity( cell, connect, num_connect );
1696  if( MB_SUCCESS != result ) return;
1697  bool first = true;
1698  EntityHandle primals[20];
1699  std::vector< int > ids;
1700 
1701  assert( num_connect < 20 );
1702  result = mbImpl->tag_get_data( dualEntityTag, connect, num_connect, primals );
1703  if( MB_SUCCESS != result ) return;
1704  ids.resize( num_connect );
1705  result = mbImpl->tag_get_data( globalIdTag, primals, num_connect, &ids[0] );
1706  if( MB_SUCCESS != result ) return;
1707  for( int i = 0; i < num_connect; i++ )
1708  {
1709  if( !first ) std::cout << "-";
1710  EntityType this_type = mbImpl->type_from_handle( primals[i] );
1711  if( this_type == MBHEX )
1712  std::cout << "h";
1713  else if( this_type == MBQUAD )
1714  std::cout << "f";
1715  else
1716  std::cout << "u";
1717 
1718  if( ids[i] != 0 )
1719  std::cout << ids[i];
1720  else
1721  std::cout << mbImpl->id_from_handle( primals[i] );
1722 
1723  first = false;
1724  }
1725 }
1726 
1727 ErrorCode DualTool::face_open_collapse( EntityHandle ocl, EntityHandle ocr )
1728 {
1729  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1730 
1731  MeshTopoUtil mtu( mbImpl );
1732 
1733  std::cout << "OC(";
1734  print_cell( ocl );
1735  std::cout << ")-(";
1736  print_cell( ocr );
1737  std::cout << ")" << std::endl;
1738 
1739  // get the primal entities we're dealing with
1740  EntityHandle split_quads[2] = { 0 }, split_edges[3] = { 0 }, split_nodes[2] = { 0 }, other_edges[6] = { 0 },
1741  other_nodes[6] = { 0 };
1742  Range hexes;
1743  ErrorCode result = foc_get_ents( ocl, ocr, split_quads, split_edges, split_nodes, hexes, other_edges, other_nodes );RR;
1744 
1745  // get star entities around edges, separated into halves
1746  std::vector< EntityHandle > star_dp1[2], star_dp2[2];
1747  result = foc_get_stars( split_quads, split_edges, star_dp1, star_dp2 );RR;
1748 
1749  if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) )
1750  return MB_TYPE_OUT_OF_RANGE;
1751 
1752  result = foc_delete_dual( split_quads, split_edges, hexes );
1753  if( MB_SUCCESS != result ) return result;
1754 
1755  EntityHandle new_quads[2], new_edges[3], new_nodes[2];
1756  result = split_pair_nonmanifold( split_quads, split_edges, split_nodes, star_dp1, star_dp2, other_edges,
1757  other_nodes, new_quads, new_edges, new_nodes );
1758  if( MB_SUCCESS != result ) return result;
1759 
1760  // now merge entities, the C of foc
1761  EntityHandle keepit, deleteit;
1762 #define MIN( a, b ) ( ( a ) < ( b ) ? ( a ) : ( b ) )
1763 #define MAX( a, b ) ( ( a ) > ( b ) ? ( a ) : ( b ) )
1764 #define KEEP_DELETE( a, b, c, d ) \
1765  { \
1766  ( c ) = MIN( a, b ); \
1767  ( d ) = MAX( a, b ); \
1768  }
1769 
1770  // find how many shared edges there were
1771  int num_shared_edges = ( split_edges[2] ? 3 : ( split_edges[1] ? 2 : 1 ) );
1772 
1773  // first the node(s)
1774  for( int i = 0; i < 3 - num_shared_edges; i++ )
1775  {
1776  KEEP_DELETE( other_nodes[2 + 2 * i], other_nodes[3 + 2 * i], keepit, deleteit );
1777  result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
1778  }
1779 
1780  // now the edges
1781  for( int i = 0; i < 4 - num_shared_edges; i++ )
1782  {
1783  KEEP_DELETE( other_edges[2 * i], other_edges[2 * i + 1], keepit, deleteit );
1784  result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
1785  }
1786 
1787  // now the faces
1788  KEEP_DELETE( split_quads[0], split_quads[1], keepit, deleteit );
1789  result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
1790 
1791  result = mbImpl->merge_entities( new_quads[0], new_quads[1], false, true );RR;
1792 
1793  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1794 
1795  // reconstruct dual
1796  result = construct_hex_dual( hexes );
1797  if( MB_SUCCESS != result ) return result;
1798 
1799  return check_dual_adjs();
1800 }
1801 
1802 ErrorCode DualTool::foc_get_ents( EntityHandle ocl,
1803  EntityHandle ocr,
1804  EntityHandle* split_quads,
1805  EntityHandle* split_edges,
1806  EntityHandle* split_nodes,
1807  Range& hexes,
1808  EntityHandle* other_edges,
1809  EntityHandle* other_nodes )
1810 {
1811  // get the entities used for foc; ocl and ocr are dual 1-cells
1812  // representing quads to be split; returned from this function:
1813  // quads[2] - 2 quads to be split
1814  // split_edges[2] - edge(s) to be split (2nd is 0 if only one)
1815  // split_node - node to be split, if any (otherwise 0)
1816  // hexes - connected hexes to split_edges
1817  // other_edges[0], [1] - edges in quads[0] and [1] sharing node with
1818  // one end of split_edges[0]
1819  // other_edges[2], [3] - other end of split_edges[0] (or [1] if 2
1820  // split_edges)
1821  // other_edges[4], [5] - edges in quads[0], [1] opposite to split_edges[0]
1822  // other_nodes[0], [1] - nodes on other_edges[0], [1] not shared with
1823  // split_edges[0]
1824  // other_nodes[2], [3] - nodes on other_edges[2], [3] not shared with
1825  // split_edges[0] (if 2 split_edges, there's only 1 opposite node
1826  // in each split quad)
1827  // (for diagram, see Tim's notes from 11/12/07)
1828 
1829  split_quads[0] = get_dual_entity( ocl );
1830  split_quads[1] = get_dual_entity( ocr );
1831  if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) )
1832  return MB_TYPE_OUT_OF_RANGE;
1833 
1834  Range common_edges;
1835  ErrorCode result = mbImpl->get_adjacencies( split_quads, 2, 1, false, common_edges );
1836  if( MB_SUCCESS != result ) return result;
1837 
1838  if( common_edges.empty() ) return MB_FAILURE;
1839  for( unsigned int i = 0; i < common_edges.size(); i++ )
1840  split_edges[i] = common_edges[i];
1841 
1842  MeshTopoUtil mtu( mbImpl );
1843 
1844  if( common_edges.size() == 3 )
1845  {
1846  // find other (non-shared) edges
1847  for( int i = 0; i < 2; i++ )
1848  {
1849  Range tmp_edges;
1850  result = mbImpl->get_adjacencies( &split_quads[i], 1, 1, false, tmp_edges );
1851  if( MB_SUCCESS != result ) return result;
1852  tmp_edges = subtract( tmp_edges, common_edges );
1853  assert( tmp_edges.size() == 1 );
1854  other_edges[i] = *tmp_edges.begin();
1855  }
1856  assert( other_edges[0] && other_edges[1] && other_edges[0] != other_edges[1] );
1857 
1858  // arrange common edges so middle is in middle
1859  result = mtu.opposite_entity( split_quads[0], other_edges[0], split_edges[1] );RR;
1860  common_edges.erase( split_edges[1] );
1861  split_edges[0] = *common_edges.begin();
1862  split_edges[2] = *common_edges.rbegin();
1863  common_edges.insert( split_edges[1] );
1864 
1865  // get split nodes and other nodes
1866  split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 );
1867  split_nodes[1] = mtu.common_entity( split_edges[2], split_edges[1], 0 );
1868  other_nodes[0] = mtu.common_entity( split_edges[0], other_edges[0], 0 );
1869  other_nodes[1] = mtu.common_entity( split_edges[2], other_edges[1], 0 );
1870 
1871  assert( other_nodes[0] && other_nodes[1] && split_nodes[0] && split_nodes[1] );
1872  assert( split_edges[0] && split_edges[1] && split_edges[2] && split_edges[0] != split_edges[1] &&
1873  split_edges[1] != split_edges[2] && split_edges[0] != split_edges[2] );
1874  }
1875  else if( common_edges.size() == 2 )
1876  {
1877  // split node is shared by split edges
1878  split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 );
1879  if( 0 == split_nodes[0] ) return MB_FAILURE;
1880  // first two other nodes are on split edges opposite split node
1881  result = mtu.opposite_entity( split_edges[0], split_nodes[0], other_nodes[0] );RR;
1882  result = mtu.opposite_entity( split_edges[1], split_nodes[0], other_nodes[1] );RR;
1883  // over split quads:
1884  for( int i = 0; i < 2; i++ )
1885  {
1886  // 1st other edge is opposite second split edge
1887  result = mtu.opposite_entity( split_quads[i], split_edges[1], other_edges[i] );RR;
1888  // 2nd other edge is opposite first split edge
1889  result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR;
1890  // last other node is opposite split node on split quad
1891  result = mtu.opposite_entity( split_quads[i], split_nodes[0], other_nodes[2 + i] );RR;
1892  }
1893  }
1894  else
1895  {
1896  const EntityHandle* connect;
1897  int num_connect;
1898  result = mbImpl->get_connectivity( split_edges[0], connect, num_connect );
1899  if( MB_SUCCESS != result ) return result;
1900  // other_nodes[0], [1] are on split edge
1901  other_nodes[0] = connect[0];
1902  other_nodes[1] = connect[1];
1903 
1904  // for each of the split quads
1905  for( int i = 0; i < 2; i++ )
1906  {
1907  // get the other edge on the split quad adj to node 0 on the split edge, by getting
1908  // edges adj to split quad and node and removing split edge; that's other_edge[i]
1909  Range tmp_range1, tmp_range2;
1910  tmp_range1.insert( connect[0] );
1911  tmp_range1.insert( split_quads[i] );
1912  result = mbImpl->get_adjacencies( tmp_range1, 1, false, tmp_range2 );
1913  if( MB_SUCCESS != result ) return result;
1914  tmp_range2.erase( split_edges[0] );
1915  assert( tmp_range2.size() == 1 );
1916  other_edges[i] = *tmp_range2.begin();
1917  // get edge connected to other node on split edge & split quad; that's
1918  // opposite prev other_edges on the split quad; that's other_edges[4+i]
1919  result = mtu.opposite_entity( split_quads[i], other_edges[i], other_edges[4 + i] );RR;
1920  // get the edge on the split quad opposite the split edge; that's other_edges[2+i]
1921  result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR;
1922  // get nodes on other side of split quad from split edge, by getting common
1923  // node between top/bottom edge and opposite edge
1924  other_nodes[2 + i] = mtu.common_entity( other_edges[i], other_edges[2 + i], 0 );
1925  other_nodes[4 + i] = mtu.common_entity( other_edges[4 + i], other_edges[2 + i], 0 );
1926  if( 0 == other_nodes[2 + i] || 0 == other_nodes[4 + i] ) return MB_FAILURE;
1927  }
1928  }
1929 
1930  result = mbImpl->get_adjacencies( split_edges, common_edges.size(), 3, false, hexes, Interface::UNION );
1931  if( MB_SUCCESS != result ) return result;
1932 
1933  assert( "split node not in other_nodes" && other_nodes[0] != split_nodes[0] && other_nodes[0] != split_nodes[1] &&
1934  other_nodes[1] != split_nodes[0] && other_nodes[1] != split_nodes[1] );
1935  assert( "each split node on an end of a split edge" && mtu.common_entity( other_nodes[0], split_edges[0], 0 ) &&
1936  ( ( ( split_edges[2] && mtu.common_entity( other_nodes[1], split_edges[2], 0 ) ) ||
1937  ( split_edges[1] && mtu.common_entity( other_nodes[1], split_edges[1], 0 ) ) ||
1938  mtu.common_entity( other_nodes[1], split_edges[0], 0 ) ) ) );
1939  assert( "opposite other edges meet at an other node" &&
1940  ( mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[0] ||
1941  ( split_edges[2] && mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[1] ) ) &&
1942  ( split_edges[2] ||
1943  ( split_edges[1] && mtu.common_entity( other_edges[2], other_edges[3], 0 ) == other_nodes[1] ) ||
1944  mtu.common_entity( other_edges[4], other_edges[5], 0 ) == other_nodes[1] ) );
1945 
1946  return MB_SUCCESS;
1947 }
1948 
1949 ErrorCode DualTool::split_pair_nonmanifold( EntityHandle* split_quads,
1950  EntityHandle* split_edges,
1951  EntityHandle* split_nodes,
1952  std::vector< EntityHandle >* star_dp1,
1953  std::vector< EntityHandle >* star_dp2,
1954  EntityHandle* /*other_edges*/,
1955  EntityHandle* /*other_nodes*/,
1956  EntityHandle* new_quads,
1957  EntityHandle* new_edges,
1958  EntityHandle* new_nodes )
1959 {
1960 
1961  // if there's a bdy in the star around the shared edge(s), get the quads on that
1962  // bdy so we know which edges to merge after the split-nonmanifold
1963  MeshTopoUtil mtu( mbImpl );
1964  ErrorCode result;
1965 
1966  // get which star the split faces are in, and choose the other one
1967  int new_side = -1;
1968  if( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() )
1969  new_side = 1;
1970  else if( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() )
1971  new_side = 0;
1972  assert( -1 != new_side );
1973  if( -1 == new_side ) return MB_FAILURE;
1974 
1975  //=============== split faces
1976 
1977  for( int i = 0; i < 2; i++ )
1978  {
1979  // get a hex in star_dp2[new_side] that's adj to this split quad, to tell
1980  // mtu which one the new quad should go with; there should be at least one,
1981  // if we have any hexes connected to the split quad
1982  EntityHandle gowith_hex = 0;
1983  for( std::vector< EntityHandle >::iterator vit = star_dp2[new_side].begin(); vit != star_dp2[new_side].end();
1984  ++vit )
1985  {
1986  if( mtu.common_entity( *vit, split_quads[i], 2 ) )
1987  {
1988  gowith_hex = *vit;
1989  break;
1990  }
1991  }
1992  assert( 0 != gowith_hex );
1993 
1994  // split manifold each of the split_quads, and put the results on the merge list
1995  result =
1996  mtu.split_entities_manifold( split_quads + i, 1, new_quads + i, NULL, ( gowith_hex ? &gowith_hex : NULL ) );RR;
1997  }
1998 
1999  // make ranges of faces which need to be explicitly adj to old, new
2000  // edge; faces come from stars and new_quads (which weren't in the stars);
2001  // new_quads go with side j, which does not have split quads
2002  Range tmp_addl_faces[2], addl_faces[2];
2003  for( int i = 0; i < 2; i++ )
2004  {
2005  std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_faces[i] ) );
2006  tmp_addl_faces[new_side].insert( new_quads[i] );
2007  }
2008 #ifndef NDEBUG
2009  bool cond1 = ( "split_quads on 1, new_quads on 0" &&
2010  tmp_addl_faces[0].find( split_quads[0] ) == tmp_addl_faces[0].end() &&
2011  tmp_addl_faces[0].find( split_quads[1] ) == tmp_addl_faces[0].end() &&
2012  tmp_addl_faces[1].find( split_quads[0] ) != tmp_addl_faces[1].end() &&
2013  tmp_addl_faces[1].find( split_quads[1] ) != tmp_addl_faces[1].end() &&
2014  tmp_addl_faces[0].find( new_quads[0] ) != tmp_addl_faces[0].end() &&
2015  tmp_addl_faces[0].find( new_quads[1] ) != tmp_addl_faces[0].end() &&
2016  tmp_addl_faces[1].find( new_quads[0] ) == tmp_addl_faces[1].end() &&
2017  tmp_addl_faces[1].find( new_quads[1] ) == tmp_addl_faces[1].end() ),
2018  cond2 = ( "split_quads on 0, new_quads on 1" &&
2019  tmp_addl_faces[0].find( split_quads[0] ) != tmp_addl_faces[0].end() &&
2020  tmp_addl_faces[0].find( split_quads[1] ) != tmp_addl_faces[0].end() &&
2021  tmp_addl_faces[1].find( split_quads[0] ) == tmp_addl_faces[1].end() &&
2022  tmp_addl_faces[1].find( split_quads[1] ) == tmp_addl_faces[1].end() &&
2023  tmp_addl_faces[0].find( new_quads[0] ) == tmp_addl_faces[0].end() &&
2024  tmp_addl_faces[0].find( new_quads[1] ) == tmp_addl_faces[0].end() &&
2025  tmp_addl_faces[1].find( new_quads[0] ) != tmp_addl_faces[1].end() &&
2026  tmp_addl_faces[1].find( new_quads[1] ) != tmp_addl_faces[1].end() );
2027 
2028  assert( cond1 || cond2 );
2029 #endif
2030 
2031  //=============== split edge(s)
2032  for( int j = 0; j < 3; j++ )
2033  {
2034  if( !split_edges[j] ) break;
2035 
2036  // filter add'l faces to only those adj to split_edges[j]
2037  addl_faces[0] = tmp_addl_faces[0];
2038  addl_faces[1] = tmp_addl_faces[1];
2039  for( int i = 0; i < 2; i++ )
2040  {
2041  result = mbImpl->get_adjacencies( &split_edges[j], 1, 2, false, addl_faces[i] );RR;
2042  }
2043 
2044  // split edge
2045  result = mtu.split_entity_nonmanifold( split_edges[j], addl_faces[1 - new_side], addl_faces[new_side],
2046  new_edges[j] );RR;
2047  }
2048 
2049  //=============== split node(s)
2050 
2051  for( int j = 0; j < 2; j++ )
2052  {
2053  if( !split_nodes[j] ) break;
2054 
2055  // if we're splitting multiple edges, there might be other edges that have the split
2056  // node; also need to know which side they're on
2057  Range tmp_addl_edges[2];
2058  result = foc_get_addl_ents( star_dp1, star_dp2, split_edges, split_nodes[j], tmp_addl_edges );RR;
2059 
2060  // also, we need to know which of the split/new edges go
2061  // with the split/new node; new edges go with side 0, split with 1
2062  for( int i = 0; i < 3; i++ )
2063  {
2064  if( !split_edges[i] ) break;
2065  tmp_addl_edges[new_side].insert( new_edges[i] );
2066  tmp_addl_edges[1 - new_side].insert( split_edges[i] );
2067  }
2068 
2069  // same for star faces and hexes
2070  for( int i = 0; i < 2; i++ )
2071  {
2072  std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_edges[i] ) );
2073  std::copy( star_dp2[i].begin(), star_dp2[i].end(), range_inserter( tmp_addl_edges[i] ) );
2074  }
2075 
2076  // finally, new quads
2077  for( int i = 0; i < 2; i++ )
2078  tmp_addl_edges[new_side].insert( new_quads[i] );
2079 
2080  // filter the entities, keeping only the ones adjacent to this node
2081  Range addl_edges[2];
2082  for( int i = 0; i < 2; i++ )
2083  {
2084  for( Range::reverse_iterator rit = tmp_addl_edges[i].rbegin(); rit != tmp_addl_edges[i].rend(); ++rit )
2085  {
2086  if( mtu.common_entity( *rit, split_nodes[j], 0 ) ) addl_edges[i].insert( *rit );
2087  }
2088  }
2089 
2090  // now split the node too
2091  result = mtu.split_entity_nonmanifold( split_nodes[j], addl_edges[1 - new_side], addl_edges[new_side],
2092  new_nodes[j] );RR;
2093  }
2094 
2095  return MB_SUCCESS;
2096 }
2097 
2098 ErrorCode DualTool::foc_get_addl_ents( std::vector< EntityHandle >* star_dp1,
2099  std::vector< EntityHandle >* /*star_dp2*/,
2100  EntityHandle* split_edges,
2101  EntityHandle split_node,
2102  Range* addl_ents )
2103 {
2104  // if we're splitting 2 edges, there might be other edges that have the split
2105  // node; also need to know which side they're on
2106 
2107  // algorithm: for a given star_dp1 (faces) on a side:
2108  // - get all edges adj to all faces -> R1
2109  // - get all edges adj to split_node -> R2
2110  // - R3 = R1 & R2 (edges connected to split_node & adj to a star face)
2111  // - R3 -= split_edges (take split edges off addl_ents)
2112 
2113  Range R2;
2114  MeshTopoUtil mtu( mbImpl );
2115  ErrorCode result = mbImpl->get_adjacencies( &split_node, 1, 1, false, R2 );RR;
2116  Range::iterator rit;
2117 
2118  for( int i = 0; i < 2; i++ )
2119  {
2120  Range R1, R3;
2121  result = mbImpl->get_adjacencies( &star_dp1[i][0], star_dp1[i].size(), 1, false, R1, Interface::UNION );RR;
2122  R3 = intersect( R1, R2 );
2123  for( int j = 0; j < 3; j++ )
2124  if( split_edges[j] ) R3.erase( split_edges[j] );
2125  addl_ents[i].merge( R3 );
2126  }
2127 
2128  return MB_SUCCESS;
2129 }
2130 
2131 ErrorCode DualTool::foc_get_stars( EntityHandle* split_quads,
2132  EntityHandle* split_edges,
2133  std::vector< EntityHandle >* star_dp1,
2134  std::vector< EntityHandle >* star_dp2 )
2135 {
2136  bool on_bdy = false, on_bdy_tmp;
2137  ErrorCode result;
2138  MeshTopoUtil mtu( mbImpl );
2139 
2140  // get the star around the split_edge
2141  std::vector< EntityHandle > qstar, hstar;
2142  unsigned int qpos = 0;
2143 
2144  for( int i = 0; i < 3; i++ )
2145  {
2146  if( !split_edges[i] ) break;
2147 
2148  // get the star around this split edge
2149  unsigned int qpos_tmp = 0;
2150  std::vector< EntityHandle > qstar_tmp, hstar_tmp;
2151  result = mtu.star_entities( split_edges[i], qstar_tmp, on_bdy_tmp, 0, &hstar_tmp );RR;
2152  // if we're on the bdy, add a null to the hex star too
2153  if( on_bdy_tmp )
2154  {
2155  assert( hstar_tmp.size() == qstar_tmp.size() - 1 );
2156  hstar_tmp.push_back( 0 );
2157  on_bdy = true;
2158  }
2159 
2160  // get the position of first split quad in star
2161  while( qpos_tmp < qstar_tmp.size() && qstar_tmp[qpos_tmp] != split_quads[0] )
2162  qpos_tmp++;
2163  if( qpos_tmp == qstar_tmp.size() ) return MB_FAILURE;
2164 
2165  bool forward;
2166  // 1st iteration is forward by definition
2167  if( 0 == i ) forward = true;
2168 
2169  // need to be careful about direction on later iters
2170  else if( hstar[qpos] == hstar_tmp[qpos_tmp] )
2171  forward = true;
2172  else if( hstar[qpos] == hstar_tmp[( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size()] &&
2173  hstar_tmp[qpos_tmp] == hstar[( qpos + qstar.size() - 1 ) % qstar.size()] )
2174  forward = false;
2175  else
2176  return MB_FAILURE;
2177 
2178  if( forward )
2179  {
2180  // 1st half of star
2181  // save hex right after split_quad[0] first
2182  star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
2183  qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
2184  while( qstar_tmp[qpos_tmp] != split_quads[1] )
2185  {
2186  star_dp1[0].push_back( qstar_tmp[qpos_tmp] );
2187  star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
2188  qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
2189  }
2190  // 2nd half of star
2191  // save hex right after split_quad[1] first
2192  star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
2193  qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
2194  while( qstar_tmp[qpos_tmp] != split_quads[0] )
2195  {
2196  star_dp1[1].push_back( qstar_tmp[qpos_tmp] );
2197  star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
2198  qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
2199  }
2200  }
2201  else
2202  {
2203  // go in reverse - take prev hex instead of current
2204  // one, and step in reverse
2205 
2206  // save hex right after split_quad[0] first
2207  qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
2208  star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
2209  while( qstar_tmp[qpos_tmp] != split_quads[1] )
2210  {
2211  star_dp1[0].push_back( qstar_tmp[qpos_tmp] );
2212  qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
2213  star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
2214  }
2215  // 2nd half of star
2216  // save hex right after split_quad[1] first
2217  qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
2218  star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
2219  while( qstar_tmp[qpos_tmp] != split_quads[0] )
2220  {
2221  star_dp1[1].push_back( qstar_tmp[qpos_tmp] );
2222  qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
2223  star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
2224  }
2225  }
2226 
2227  if( 0 == i )
2228  {
2229  // if we're on the first iteration, save results and continue, other iters
2230  // get compared to this one
2231  qstar.swap( qstar_tmp );
2232  hstar.swap( hstar_tmp );
2233  on_bdy = on_bdy_tmp;
2234  qpos = qpos_tmp;
2235  }
2236  }
2237 
2238  // split quads go on list with NULLs, if any, otherwise on 2nd
2239  if( on_bdy )
2240  {
2241  if( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) != star_dp2[0].end() )
2242  {
2243  // remove *all* the zeros
2244  star_dp2[0].erase( std::remove( star_dp2[0].begin(), star_dp2[0].end(), 0 ), star_dp2[0].end() );
2245  // put the split quads on this half
2246  star_dp1[0].push_back( split_quads[0] );
2247  star_dp1[0].push_back( split_quads[1] );
2248  }
2249  else
2250  {
2251  star_dp2[1].erase( std::remove( star_dp2[1].begin(), star_dp2[1].end(), 0 ), star_dp2[1].end() );
2252  // put the split quads on this half
2253  star_dp1[1].push_back( split_quads[0] );
2254  star_dp1[1].push_back( split_quads[1] );
2255  }
2256  }
2257  else
2258  {
2259  star_dp1[1].push_back( split_quads[0] );
2260  star_dp1[1].push_back( split_quads[1] );
2261  }
2262 
2263  // some error checking
2264  if( !( ( ( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) == star_dp1[0].end() &&
2265  std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) == star_dp1[0].end() &&
2266  std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() &&
2267  std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) != star_dp1[1].end() ) ||
2268 
2269  ( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) == star_dp1[1].end() &&
2270  std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) == star_dp1[1].end() &&
2271  std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() &&
2272  std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) != star_dp1[0].end() ) ) ) )
2273  {
2274  std::cerr << "foc_get_stars: both split quads should be on the same star list half and not "
2275  << "on the other, failed" << std::endl;
2276  return MB_FAILURE;
2277  }
2278 
2279  if( !( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) == star_dp2[0].end() &&
2280  std::find( star_dp2[1].begin(), star_dp2[1].end(), 0 ) == star_dp2[1].end() ) )
2281  {
2282  std::cerr << "foc_get_stars: no NULLs on the hstar lists, failed";
2283  return MB_FAILURE;
2284  }
2285 
2286  return MB_SUCCESS;
2287 }
2288 
2289 ErrorCode DualTool::foc_delete_dual( EntityHandle* split_quads, EntityHandle* split_edges, Range& hexes )
2290 {
2291  // special delete dual procedure, because in some cases we need to delete
2292  // a sheet too since it'll get merged into another
2293 
2294  // figure out whether we'll need to delete a sheet
2295  EntityHandle sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[0] ) );
2296  if( split_edges[1] ) sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[1] ) );
2297  EntityHandle chordl = get_dual_hyperplane( get_dual_entity( split_quads[0] ) );
2298  EntityHandle chordr = get_dual_hyperplane( get_dual_entity( split_quads[1] ) );
2299  assert( 0 != sheet1 && 0 != chordl && 0 != chordr );
2300  Range parentsl, parentsr;
2301  ErrorCode result = mbImpl->get_parent_meshsets( chordl, parentsl );
2302  if( MB_SUCCESS != result ) return result;
2303  result = mbImpl->get_parent_meshsets( chordr, parentsr );
2304  if( MB_SUCCESS != result ) return result;
2305  parentsl.erase( sheet1 );
2306  parentsr.erase( sheet1 );
2307 
2308  // before deciding which one to delete, collect the other cells which must
2309  // be deleted, and all the chords/sheets they're on
2310  Range adj_ents, dual_ents, cells1or2;
2311  for( int i = 0; i < 3; i++ )
2312  {
2313  result = mbImpl->get_adjacencies( hexes, i, false, adj_ents, Interface::UNION );
2314  if( MB_SUCCESS != result ) return result;
2315  }
2316 
2317  // cache any adjacent hexes, for rebuilding the dual later
2318  result = mbImpl->get_adjacencies( adj_ents, 3, false, hexes, Interface::UNION );
2319  if( MB_SUCCESS != result ) return result;
2320 
2321  for( Range::iterator rit = adj_ents.begin(); rit != adj_ents.end(); ++rit )
2322  {
2323  EntityHandle this_ent = get_dual_entity( *rit );
2324  dual_ents.insert( this_ent );
2325  int dim = mbImpl->dimension_from_handle( this_ent );
2326  if( 1 == dim || 2 == dim ) cells1or2.insert( this_ent );
2327  }
2328 
2329  Range dual_hps;
2330  for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit )
2331  dual_hps.insert( get_dual_hyperplane( *rit ) );
2332 
2333  result = delete_dual_entities( dual_ents );
2334  if( MB_SUCCESS != result ) return result;
2335 
2336  // now decide which sheet to delete (to be merged into the other)
2337  EntityHandle sheet_delete = 0;
2338  if( is_blind( *parentsl.begin() ) )
2339  sheet_delete = *parentsl.begin();
2340  else if( is_blind( *parentsr.begin() ) )
2341  sheet_delete = *parentsr.begin();
2342  else
2343  {
2344  // neither is blind, take the one with fewer cells
2345  Range tmp_ents;
2346  int numl, numr;
2347  result = mbImpl->get_number_entities_by_handle( *parentsl.begin(), numl );
2348  if( MB_SUCCESS != result ) return result;
2349  result = mbImpl->get_number_entities_by_handle( *parentsr.begin(), numr );
2350  if( MB_SUCCESS != result ) return result;
2351  sheet_delete = ( numl > numr ? *parentsr.begin() : *parentsl.begin() );
2352  }
2353  assert( 0 != sheet_delete );
2354 
2355  // after deleting cells, check for empty chords & sheets, and delete those too
2356  for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
2357  {
2358  Range tmp_ents;
2359  result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
2360  if( MB_SUCCESS != result ) return result;
2361  if( tmp_ents.empty() )
2362  {
2363  result = mbImpl->delete_entities( &( *rit ), 1 );
2364  if( MB_SUCCESS != result ) return result;
2365  }
2366  else if( *rit == sheet_delete )
2367  {
2368  // delete the sheet
2369  result = mbImpl->delete_entities( &( *rit ), 1 );
2370  if( MB_SUCCESS != result ) return result;
2371  }
2372  }
2373 
2374  // now just to be safe, add the hexes bridge-adjacent across vertices
2375  // to the hexes we already have
2376  Range tmp_hexes;
2377  MeshTopoUtil mtu( mbImpl );
2378  for( Range::iterator rit = hexes.begin(); rit != hexes.end(); ++rit )
2379  {
2380  result = mtu.get_bridge_adjacencies( *rit, 0, 3, tmp_hexes );
2381  if( MB_SUCCESS != result ) return result;
2382  }
2383  hexes.merge( tmp_hexes );
2384 
2385  return MB_SUCCESS;
2386 }
2387 
2388 //! returns true if all vertices are dual to hexes (not faces)
2389 bool DualTool::is_blind( const EntityHandle chord_or_sheet )
2390 {
2391  // must be an entity set
2392  if( TYPE_FROM_HANDLE( chord_or_sheet ) != MBENTITYSET ) return false;
2393 
2394  // get the vertices
2395  Range verts, ents;
2396  ErrorCode result = mbImpl->get_entities_by_handle( chord_or_sheet, ents );
2397  if( MB_SUCCESS != result || ents.empty() ) return false;
2398 
2399  result = mbImpl->get_adjacencies( ents, 0, false, verts, Interface::UNION );
2400  if( MB_SUCCESS != result || verts.empty() ) return false;
2401 
2402  for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit )
2403  {
2404  // get dual entity for this vertex
2405  EntityHandle dual_ent = get_dual_entity( *rit );
2406  if( 0 == dual_ent ) continue;
2407  if( TYPE_FROM_HANDLE( dual_ent ) == MBQUAD ) return false;
2408  }
2409 
2410  // if none of the vertices' duals were quads, chord_or_sheet must be blind
2411  return true;
2412 }
2413 
2414 //! given a 1-cell and a chord, return the neighboring vertices on the
2415 //! chord, in the same order as the 1-cell's vertices
2416 ErrorCode DualTool::get_opposite_verts( const EntityHandle middle_edge, const EntityHandle chord, EntityHandle* verts )
2417 {
2418  // get the edges on the chord, in order, and move to middle_edge
2419  std::vector< EntityHandle > chord_edges;
2420  const EntityHandle* connect;
2421  int num_connect;
2422 
2423  ErrorCode result = mbImpl->get_entities_by_handle( chord, chord_edges );RR;
2424  std::vector< EntityHandle >::iterator vit = std::find( chord_edges.begin(), chord_edges.end(), middle_edge );
2425  result = mbImpl->get_connectivity( middle_edge, connect, num_connect );RR;
2426 
2427  if(
2428  // middle_edge isn't on this chord
2429  vit == chord_edges.end() ||
2430  // chord only has 1 edge
2431  chord_edges.size() == 1 ||
2432  // middle_edge is at beginning or end and chord isn't blind
2433  ( ( vit == chord_edges.begin() || vit == chord_edges.end() - 1 ) && !is_blind( chord ) ) )
2434  return MB_FAILURE;
2435 
2436  else if( chord_edges.size() == 2 )
2437  {
2438  // easier if it's a 2-edge blind chord, just get vertices in opposite order
2439  verts[0] = connect[1];
2440  verts[1] = connect[0];
2441  return MB_SUCCESS;
2442  }
2443 
2444  // get vertices with the prev edge & subtract vertices of 1-cell
2445  if( vit == chord_edges.begin() )
2446  vit = chord_edges.end() - 1;
2447  else
2448  --vit;
2449  Range dum_connect, middle_connect;
2450  result = mbImpl->get_connectivity( &middle_edge, 1, middle_connect );RR;
2451  result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR;
2452  dum_connect = subtract( dum_connect, middle_connect );
2453  if( dum_connect.size() != 1 )
2454  {
2455  std::cerr << "Trouble traversing chord." << std::endl;
2456  return MB_FAILURE;
2457  }
2458 
2459  // put in verts[0]
2460  verts[0] = *dum_connect.begin();
2461 
2462  // same with prev edge
2463  ++vit;
2464  if( vit == chord_edges.end() ) vit = chord_edges.begin();
2465  ++vit;
2466  dum_connect.clear();
2467  result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR;
2468  dum_connect = subtract( dum_connect, middle_connect );
2469  if( dum_connect.size() != 1 )
2470  {
2471  std::cerr << "Trouble traversing chord." << std::endl;
2472  return MB_FAILURE;
2473  }
2474 
2475  // put in verts[1]
2476  verts[1] = *dum_connect.begin();
2477 
2478  // if verts[0] and 1st vertex of 1cell don't have common edge, switch verts
2479  MeshTopoUtil mtu( mbImpl );
2480  if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) )
2481  {
2482  EntityHandle dum_h = verts[0];
2483  verts[0] = verts[1];
2484  verts[1] = dum_h;
2485  }
2486 
2487  if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) )
2488  {
2489  std::cerr << "Trouble traversing chord." << std::endl;
2490  return MB_FAILURE;
2491  }
2492 
2493  return MB_SUCCESS;
2494 }
2495 
2496 ErrorCode DualTool::get_dual_entities( const EntityHandle dual_ent,
2497  Range* dcells,
2498  Range* dedges,
2499  Range* dverts,
2500  Range* dverts_loop,
2501  Range* dedges_loop )
2502 {
2503  ErrorCode result = MB_SUCCESS;
2504 
2505  if( NULL != dcells )
2506  {
2507  result = mbImpl->get_entities_by_type( dual_ent, MBPOLYGON, *dcells );
2508  if( MB_SUCCESS != result ) return result;
2509  }
2510 
2511  if( NULL != dedges )
2512  {
2513  if( NULL != dcells )
2514  result = mbImpl->get_adjacencies( *dcells, 1, false, *dedges, Interface::UNION );
2515  else
2516  result = mbImpl->get_entities_by_type( dual_ent, MBEDGE, *dedges );
2517 
2518  if( MB_SUCCESS != result ) return result;
2519  }
2520 
2521  if( NULL != dverts )
2522  {
2523  if( NULL != dcells )
2524  result = mbImpl->get_adjacencies( *dcells, 0, false, *dverts, Interface::UNION );
2525  else if( NULL != dedges )
2526  result = mbImpl->get_adjacencies( *dedges, 0, false, *dverts, Interface::UNION );
2527  else
2528  {
2529  Range all_ents;
2530  result = mbImpl->get_entities_by_handle( dual_ent, all_ents );RR;
2531  result = mbImpl->get_adjacencies( all_ents, 0, false, *dverts, Interface::UNION );
2532  }
2533 
2534  if( MB_SUCCESS != result ) return result;
2535  }
2536 
2537  if( NULL != dverts_loop && NULL != dverts )
2538  {
2539  static std::vector< EntityHandle > dual_ents;
2540  dual_ents.resize( dverts->size() );
2541  result = mbImpl->tag_get_data( dualEntity_tag(), *dverts, &dual_ents[0] );
2542  if( MB_SUCCESS != result ) return result;
2543  Range::iterator rit;
2544  unsigned int i;
2545  for( rit = dverts->begin(), i = 0; rit != dverts->end(); ++rit, i++ )
2546  if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBQUAD ) dverts_loop->insert( *rit );
2547  }
2548 
2549  if( NULL != dedges_loop && NULL != dedges )
2550  {
2551  static std::vector< EntityHandle > dual_ents;
2552  dual_ents.resize( dedges->size() );
2553  result = mbImpl->tag_get_data( dualEntity_tag(), *dedges, &dual_ents[0] );
2554  if( MB_SUCCESS != result ) return result;
2555  Range::iterator rit;
2556  unsigned int i;
2557  for( rit = dedges->begin(), i = 0; rit != dedges->end(); ++rit, i++ )
2558  if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBEDGE ) dedges_loop->insert( *rit );
2559  }
2560 
2561  return result;
2562 }
2563 
2564 ErrorCode DualTool::list_entities( const EntityHandle* entities, const int num_entities ) const
2565 {
2566  Range temp_range;
2567  ErrorCode result;
2568  if( NULL == entities && 0 == num_entities )
2569  return mbImpl->list_entities( entities, num_entities );
2570 
2571  else if( NULL == entities && 0 < num_entities )
2572  {
2573 
2574  // list all entities of all types
2575  std::cout << std::endl;
2576  for( EntityType this_type = MBVERTEX; this_type < MBMAXTYPE; this_type++ )
2577  {
2578  result = mbImpl->get_entities_by_type( 0, this_type, temp_range );
2579  if( MB_SUCCESS != result ) return result;
2580  }
2581  }
2582 
2583  else
2584  {
2585  std::copy( entities, entities + num_entities, range_inserter( temp_range ) );
2586  }
2587 
2588  return list_entities( temp_range );
2589 }
2590 
2591 ErrorCode DualTool::list_entities( const Range& entities ) const
2592 {
2593  // now print each entity, listing the dual information first then calling Interface to do
2594  // the rest
2595  ErrorCode result = MB_SUCCESS, tmp_result;
2596  for( Range::const_iterator iter = entities.begin(); iter != entities.end(); ++iter )
2597  {
2598  EntityType this_type = TYPE_FROM_HANDLE( *iter );
2599  std::cout << CN::EntityTypeName( this_type ) << " " << ID_FROM_HANDLE( *iter ) << ":" << std::endl;
2600 
2601  EntityHandle dual_ent = get_dual_entity( *iter );
2602  if( 0 != dual_ent )
2603  {
2604  std::cout << "Dual to " << CN::EntityTypeName( mbImpl->type_from_handle( dual_ent ) ) << " "
2605  << mbImpl->id_from_handle( dual_ent ) << std::endl;
2606  }
2607 
2608  if( TYPE_FROM_HANDLE( *iter ) == MBENTITYSET )
2609  {
2610  EntityHandle chord = 0, sheet = 0;
2611  int id;
2612  result = mbImpl->tag_get_data( dualCurve_tag(), &( *iter ), 1, &chord );
2613  if( MB_SUCCESS != result ) return result;
2614  result = mbImpl->tag_get_data( dualSurface_tag(), &( *iter ), 1, &sheet );
2615  if( MB_SUCCESS != result ) return result;
2616  result = mbImpl->tag_get_data( globalId_tag(), &( *iter ), 1, &id );
2617  if( MB_SUCCESS != result ) return result;
2618 
2619  if( 0 != chord ) std::cout << "(Dual chord " << id << ")" << std::endl;
2620  if( 0 != sheet ) std::cout << "(Dual sheet " << id << ")" << std::endl;
2621  }
2622 
2623  tmp_result = mbImpl->list_entity( *iter );
2624  if( MB_SUCCESS != tmp_result ) result = tmp_result;
2625  }
2626 
2627  return result;
2628 }
2629 
2630 ErrorCode DualTool::face_shrink( EntityHandle odedge )
2631 {
2632  // some preliminary checking
2633  if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE;
2634 
2635  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
2636 
2637  std::cout << "FS(";
2638  print_cell( odedge );
2639  std::cout << ")" << std::endl;
2640 
2641  EntityHandle quads[4], hexes[2];
2642  std::vector< EntityHandle > connects[4], side_quads[2];
2643 
2644  // get the quads along the chord through the 2 hexes, and the vertices
2645  // for those quads
2646  ErrorCode result = fs_get_quads( odedge, quads, hexes, connects );
2647  if( MB_SUCCESS != result ) return result;
2648 
2649  // flip/rotate connect arrays so they align & are same sense
2650  result = fs_check_quad_sense( hexes[0], quads[0], connects );
2651  if( MB_SUCCESS != result ) return result;
2652 
2653  // get the quad loops along the "side" surfaces
2654  result = fs_get_quad_loops( hexes, connects, side_quads );
2655  if( MB_SUCCESS != result ) return result;
2656 
2657  // ok, done with setup; now delete dual entities affected by this operation,
2658  // which is all the entities adjacent to vertices of dual edge
2659  Range adj_verts, adj_edges, dual_ents, cells1or2;
2660  MeshTopoUtil mtu( mbImpl );
2661  result = mtu.get_bridge_adjacencies( odedge, 0, 1, adj_edges );
2662  if( MB_SUCCESS != result ) return result;
2663  result = mbImpl->get_adjacencies( adj_edges, 0, false, adj_verts, Interface::UNION );
2664  if( MB_SUCCESS != result ) return result;
2665  for( int i = 1; i <= 3; i++ )
2666  {
2667  result = mbImpl->get_adjacencies( adj_verts, i, false, dual_ents, Interface::UNION );
2668  if( MB_SUCCESS != result ) return result;
2669  }
2670 
2671  // before deleting dual, grab the 1- and 2-cells
2672  for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit )
2673  {
2674  int dim = mbImpl->dimension_from_handle( *rit );
2675  if( 1 == dim || 2 == dim ) cells1or2.insert( *rit );
2676  }
2677  Range dual_hps;
2678  for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit )
2679  dual_hps.insert( get_dual_hyperplane( *rit ) );
2680 
2681  dual_ents.insert( odedge );
2682  result = delete_dual_entities( dual_ents );
2683  if( MB_SUCCESS != result ) return result;
2684 
2685  // after deleting cells, check for empty chords & sheets, and delete those too
2686  for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
2687  {
2688  Range tmp_ents;
2689  result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
2690  if( MB_SUCCESS != result ) return result;
2691  if( tmp_ents.empty() )
2692  {
2693  result = mbImpl->delete_entities( &( *rit ), 1 );
2694  if( MB_SUCCESS != result ) return result;
2695  }
2696  }
2697 
2698  // remove any explicit adjacencies between side quads and hexes; don't check
2699  // for error, since there might not be adjacencies
2700  for( int i = 0; i < 4; i++ )
2701  {
2702  for( int j = 0; j < 2; j++ )
2703  {
2704  result = mbImpl->remove_adjacencies( side_quads[j][i], &hexes[j], 1 );
2705  }
2706  }
2707 
2708  // make inner ring of vertices
2709  // get centroid of quad2
2710  double q2coords[12], avg[3] = { 0.0, 0.0, 0.0 };
2711  result = mbImpl->get_coords( &connects[1][0], 4, q2coords );
2712  if( MB_SUCCESS != result ) return result;
2713  for( int i = 0; i < 4; i++ )
2714  {
2715  avg[0] += q2coords[3 * i];
2716  avg[1] += q2coords[3 * i + 1];
2717  avg[2] += q2coords[3 * i + 2];
2718  }
2719  avg[0] *= .25;
2720  avg[1] *= .25;
2721  avg[2] *= .25;
2722  // position new vertices
2723  connects[3].resize( 4 );
2724  for( int i = 0; i < 4; i++ )
2725  {
2726  q2coords[3 * i] = avg[0] + .25 * ( q2coords[3 * i] - avg[0] );
2727  q2coords[3 * i + 1] = avg[1] + .25 * ( q2coords[3 * i + 1] - avg[1] );
2728  q2coords[3 * i + 2] = avg[2] + .25 * ( q2coords[3 * i + 2] - avg[2] );
2729  result = mbImpl->create_vertex( &q2coords[3 * i], connects[3][i] );
2730  if( MB_SUCCESS != result ) return result;
2731  }
2732 
2733  // ok, now have the 4 connectivity arrays for 4 quads; construct hexes
2734  EntityHandle hconnect[8], new_hexes[4];
2735  int new_hex_ids[4];
2736 
2737  for( int i = 0; i < 4; i++ )
2738  {
2739  int i1 = i, i2 = ( i + 1 ) % 4;
2740  hconnect[0] = connects[0][i1];
2741  hconnect[1] = connects[0][i2];
2742  hconnect[2] = connects[3][i2];
2743  hconnect[3] = connects[3][i1];
2744 
2745  hconnect[4] = connects[1][i1];
2746  hconnect[5] = connects[1][i2];
2747  hconnect[6] = connects[2][i2];
2748  hconnect[7] = connects[2][i1];
2749 
2750  result = mbImpl->create_element( MBHEX, hconnect, 8, new_hexes[i] );
2751  if( MB_SUCCESS != result ) return result;
2752 
2753  // test for equiv entities from the side quads, and make explicit adjacencies
2754  // if there are any
2755  for( int j = 0; j < 2; j++ )
2756  {
2757  if( mtu.equivalent_entities( side_quads[j][i] ) )
2758  {
2759  result = mbImpl->add_adjacencies( side_quads[j][i], &new_hexes[i], 1, false );
2760  if( MB_SUCCESS != result ) return result;
2761  }
2762  }
2763 
2764  new_hex_ids[i] = ++maxHexId;
2765  }
2766 
2767  // set the global id tag on the new hexes
2768  result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 4, new_hex_ids );
2769  if( MB_SUCCESS != result ) return result;
2770 
2771  // now fixup other two hexes; start by getting hex through quads 0, 1
2772  // make this first hex switch to the other side, to make the dual look like
2773  // a hex push
2774  int tmp_ids[2];
2775  result = mbImpl->tag_get_data( globalId_tag(), hexes, 2, tmp_ids );
2776  if( MB_SUCCESS != result ) return result;
2777 
2778  result = mbImpl->delete_entities( hexes, 2 );
2779  if( MB_SUCCESS != result ) return result;
2780  result = mbImpl->delete_entities( &quads[1], 1 );
2781  if( MB_SUCCESS != result ) return result;
2782  for( int i = 0; i < 4; i++ )
2783  {
2784  hconnect[i] = connects[3][i];
2785  hconnect[4 + i] = connects[2][i];
2786  }
2787  result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[0] );
2788  if( MB_SUCCESS != result ) return result;
2789 
2790  for( int i = 0; i < 4; i++ )
2791  {
2792  hconnect[i] = connects[0][i];
2793  hconnect[4 + i] = connects[3][i];
2794  }
2795  result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[1] );
2796  if( MB_SUCCESS != result ) return result;
2797 
2798  // check for and fix any explicit adjacencies on either end quad
2799  if( mtu.equivalent_entities( quads[0] ) ) mbImpl->add_adjacencies( quads[0], &hexes[1], 1, false );
2800  if( mtu.equivalent_entities( quads[2] ) ) mbImpl->add_adjacencies( quads[2], &hexes[0], 1, false );
2801 
2802  // re-set the global ids for the hexes to what they were
2803  result = mbImpl->tag_set_data( globalId_tag(), hexes, 2, tmp_ids );
2804  if( MB_SUCCESS != result ) return result;
2805 
2806  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
2807 
2808  // now update the dual
2809  Range tmph;
2810  result = mtu.get_bridge_adjacencies( hexes[0], 0, 3, tmph );
2811  if( MB_SUCCESS != result ) return result;
2812  result = mtu.get_bridge_adjacencies( hexes[1], 0, 3, tmph );
2813  if( MB_SUCCESS != result ) return result;
2814  tmph.insert( hexes[1] );
2815  result = construct_hex_dual( tmph );
2816  if( MB_SUCCESS != result ) return result;
2817 
2818  return result;
2819 }
2820 
2821 ErrorCode DualTool::fs_get_quad_loops( EntityHandle* hexes,
2822  std::vector< EntityHandle >* connects,
2823  std::vector< EntityHandle >* side_quads )
2824 {
2825  for( int i = 0; i < 4; i++ )
2826  {
2827  for( int j = 0; j < 2; j++ )
2828  {
2829  Range adj_ents, dum_quads;
2830  adj_ents.insert( hexes[j] );
2831  adj_ents.insert( connects[j][i] );
2832  adj_ents.insert( connects[j][( i + 1 ) % 4] );
2833  adj_ents.insert( connects[j + 1][i] );
2834  adj_ents.insert( connects[j + 1][( i + 1 ) % 4] );
2835 
2836  ErrorCode result = mbImpl->get_adjacencies( adj_ents, 2, false, dum_quads );
2837  if( MB_SUCCESS != result ) return result;
2838  assert( 1 == dum_quads.size() );
2839  side_quads[j].push_back( *dum_quads.begin() );
2840  }
2841  }
2842 
2843  return MB_SUCCESS;
2844 }
2845 
2846 ErrorCode DualTool::fs_check_quad_sense( EntityHandle hex0, EntityHandle quad0, std::vector< EntityHandle >* connects )
2847 {
2848  // check sense of 0th quad wrt hex; since sense is out of element,
2849  // switch if quad is NOT reversed wrt hex
2850  int dum1, dum2, sense = 0;
2851  ErrorCode result = mbImpl->side_number( hex0, quad0, dum1, sense, dum2 );
2852  if( MB_SUCCESS != result ) return result;
2853  assert( 0 != sense );
2854  if( 1 == sense )
2855  {
2856  // just switch sense of this one; others will get switched next
2857  EntityHandle dum = connects[0][0];
2858  connects[0][0] = connects[0][2];
2859  connects[0][2] = dum;
2860  }
2861 
2862  // check sense of 1st, 2nd quads, rotate if necessary to align connect arrays
2863  int index0 = -1, index2 = -1, sense0 = 0, sense2 = 0;
2864  MeshTopoUtil mtu( mbImpl );
2865  for( int i = 0; i < 4; i++ )
2866  {
2867  if( 0 != mtu.common_entity( connects[0][0], connects[1][i], 1 ) )
2868  {
2869  index0 = i;
2870  if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 1 ) % 4], 1 ) )
2871  sense0 = 1;
2872  else if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 4 - 1 ) % 4], 1 ) )
2873  sense0 = -1;
2874  break;
2875  }
2876  }
2877 
2878  assert( index0 != -1 && sense0 != 0 );
2879 
2880  if( sense0 == -1 )
2881  {
2882  EntityHandle dumh = connects[1][0];
2883  connects[1][0] = connects[1][2];
2884  connects[1][2] = dumh;
2885  if( index0 % 2 == 0 ) index0 = ( index0 + 2 ) % 4;
2886  }
2887 
2888  if( index0 != 0 )
2889  {
2890  std::vector< EntityHandle > tmpc;
2891  for( int i = 0; i < 4; i++ )
2892  tmpc.push_back( connects[1][( index0 + i ) % 4] );
2893  connects[1].swap( tmpc );
2894  }
2895 
2896  for( int i = 0; i < 4; i++ )
2897  {
2898  if( 0 != mtu.common_entity( connects[1][0], connects[2][i], 1 ) )
2899  {
2900  index2 = i;
2901  if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 1 ) % 4], 1 ) )
2902  sense2 = 1;
2903  else if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 4 - 1 ) % 4], 1 ) )
2904  sense2 = -1;
2905  break;
2906  }
2907  }
2908 
2909  assert( index2 != -1 && sense2 != 0 );
2910 
2911  if( sense2 == -1 )
2912  {
2913  EntityHandle dumh = connects[2][0];
2914  connects[2][0] = connects[2][2];
2915  connects[2][2] = dumh;
2916  if( index2 % 2 == 0 ) index2 = ( index2 + 2 ) % 4;
2917  }
2918 
2919  if( index2 != 0 )
2920  {
2921  std::vector< EntityHandle > tmpc;
2922  for( int i = 0; i < 4; i++ )
2923  tmpc.push_back( connects[2][( index2 + i ) % 4] );
2924  connects[2].swap( tmpc );
2925  }
2926 
2927  return MB_SUCCESS;
2928 }
2929 
2930 //! effect reverse face shrink operation
2931 ErrorCode DualTool::rev_face_shrink( EntityHandle odedge )
2932 {
2933  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
2934 
2935  // some preliminary checking
2936  if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE;
2937 
2938  std::cout << "-FS(";
2939  print_cell( odedge );
2940  std::cout << ")" << std::endl;
2941 
2942  EntityHandle quads[4], hexes[2];
2943  std::vector< EntityHandle > connects[4], side_quads[2];
2944 
2945  // get three quads (shared quad & 2 end quads), hexes, and quad
2946  // connects
2947  ErrorCode result = fs_get_quads( odedge, quads, hexes, connects );
2948  if( MB_SUCCESS != result ) return result;
2949 
2950  // adjust sense & rotation so they're aligned, together & wrt first
2951  // hex
2952  result = fs_check_quad_sense( hexes[0], quads[0], connects );
2953  if( MB_SUCCESS != result ) return result;
2954 
2955  result = fsr_get_fourth_quad( connects, side_quads );
2956  if( MB_SUCCESS != result )
2957  {
2958  std::cout << "Can't do -FS here, two hexes must be adjacent to ring of 4 hexes." << std::endl;
2959  return result;
2960  }
2961 
2962  Range adj_ents, outer_hexes, all_adjs;
2963 
2964  // first get the entities connected to interior 4 verts
2965  for( int i = 1; i <= 3; i++ )
2966  {
2967  result = mbImpl->get_adjacencies( &connects[1][0], 4, i, false, adj_ents, Interface::UNION );
2968  if( MB_SUCCESS != result ) return result;
2969  }
2970 
2971  // next get all entities adjacent to those; these will have their dual
2972  // entities deleted
2973  for( int i = 0; i < 3; i++ )
2974  {
2975  result = mbImpl->get_adjacencies( adj_ents, i, false, all_adjs, Interface::UNION );
2976  if( MB_SUCCESS != result ) return result;
2977  }
2978 
2979  // get the dual entities and delete them
2980  Range dual_ents, dual_hps;
2981  for( Range::iterator rit = all_adjs.begin(); rit != all_adjs.end(); ++rit )
2982  {
2983  EntityHandle this_ent = get_dual_entity( *rit );
2984  dual_ents.insert( this_ent );
2985  }
2986 
2987  // before deleting dual, grab the 1- and 2-cells
2988  for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit )
2989  {
2990  int dim = mbImpl->dimension_from_handle( *rit );
2991  if( 1 == dim || 2 == dim ) dual_hps.insert( get_dual_hyperplane( *rit ) );
2992  }
2993 
2994  result = delete_dual_entities( dual_ents );
2995  if( MB_SUCCESS != result ) return result;
2996 
2997  // after deleting cells, check for empty chords & sheets, and delete those too
2998  for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
2999  {
3000  Range tmp_ents;
3001  result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
3002  if( MB_SUCCESS != result ) return result;
3003  if( tmp_ents.empty() )
3004  {
3005  result = mbImpl->delete_entities( &( *rit ), 1 );
3006  if( MB_SUCCESS != result ) return result;
3007  }
3008  }
3009 
3010  // before re-connecting two hexes, check for existing quad on 4th quad vertices;
3011  // if there is a quad there, need to add explicit adjs to any adj hexes, since
3012  // by definition there'll be another quad on those vertices
3013  bool need_explicit = false;
3014  Range adj_quads;
3015  result = mbImpl->get_adjacencies( &connects[3][0], 4, 2, false, adj_quads );
3016  if( MB_MULTIPLE_ENTITIES_FOUND == result || !adj_quads.empty() )
3017  {
3018  // there's already a quad for these 4 vertices; by definition,
3019  // we'll be creating equivalent entities, so that original quad
3020  // needs explicit adj's to its bounding elements
3021  need_explicit = true;
3022  for( Range::iterator rit = adj_quads.begin(); rit != adj_quads.end(); ++rit )
3023  {
3024  Range adj_hexes;
3025  result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, adj_hexes );RR;
3026  result = mbImpl->add_adjacencies( *rit, adj_hexes, false );RR;
3027  }
3028  }
3029 
3030  // re-connect the two hexes
3031  std::vector< EntityHandle > new_connect;
3032  std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) );
3033  std::copy( connects[2].begin(), connects[2].end(), std::back_inserter( new_connect ) );
3034  result = mbImpl->set_connectivity( hexes[0], &new_connect[0], 8 );
3035  if( MB_SUCCESS != result ) return result;
3036 
3037  new_connect.clear();
3038  std::copy( connects[0].begin(), connects[0].end(), std::back_inserter( new_connect ) );
3039  std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) );
3040  result = mbImpl->set_connectivity( hexes[1], &new_connect[0], 8 );
3041  if( MB_SUCCESS != result ) return result;
3042 
3043  // test for equiv entities from the side quads, and make explicit
3044  // adjacencies if there are any
3045  MeshTopoUtil mtu( mbImpl );
3046  for( int j = 0; j < 2; j++ )
3047  {
3048  for( int i = 0; i < 4; i++ )
3049  {
3050  if( mtu.equivalent_entities( side_quads[j][i] ) )
3051  {
3052  result = mbImpl->add_adjacencies( side_quads[j][i], &hexes[j], 1, false );
3053  if( MB_SUCCESS != result ) return result;
3054  }
3055  }
3056  }
3057 
3058  // remove hexes we want to keep
3059  adj_ents.erase( hexes[0] );
3060  adj_ents.erase( hexes[1] );
3061 
3062  // delete the other interior entities
3063  result = mbImpl->delete_entities( adj_ents );
3064  if( MB_SUCCESS != result ) return result;
3065 
3066  EntityHandle new_quad;
3067  result = mbImpl->create_element( MBQUAD, &connects[3][0], 4, new_quad );RR;
3068  if( need_explicit )
3069  {
3070  result = mbImpl->add_adjacencies( new_quad, hexes, 2, false );RR;
3071  }
3072 
3073  if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
3074 
3075  // now update the dual
3076  result = construct_hex_dual( hexes, 2 );
3077  if( MB_SUCCESS != result ) return result;
3078 
3079  return MB_SUCCESS;
3080 }
3081 
3082 ErrorCode DualTool::fsr_get_fourth_quad( std::vector< EntityHandle >* connects,
3083  std::vector< EntityHandle >* side_quads )
3084 {
3085  // given the first three quad connectivities in ordered vectors, get the fourth,
3086  // where the fourth is really the 4 vertices originally shared by the 2 hexes
3087  // before the face shrink on them
3088 
3089  // vertex on 4th quad is in quad adj to other 3 verts
3090  for( int i = 0; i < 4; i++ )
3091  {
3092  Range start_verts, tmp_verts, quads;
3093  for( int j = 0; j < 3; j++ )
3094  start_verts.insert( connects[j][i] );
3095  ErrorCode result = mbImpl->get_adjacencies( start_verts, 2, false, quads );
3096  if( MB_SUCCESS != result ) return result;
3097  assert( quads.size() == 1 );
3098  result = mbImpl->get_adjacencies( &( *quads.begin() ), 1, 0, false, tmp_verts );RR;
3099  tmp_verts = subtract( tmp_verts, start_verts );
3100  assert( 1 == tmp_verts.size() );
3101  connects[3].push_back( *tmp_verts.begin() );
3102  }
3103 
3104  // now get the side quads
3105  for( int i = 0; i < 4; i++ )
3106  {
3107  Range dum_ents, hexes;
3108  dum_ents.insert( connects[1][i] );
3109  dum_ents.insert( connects[1][( i + 1 ) % 4] );
3110  dum_ents.insert( connects[3][i] );
3111  ErrorCode result = mbImpl->get_adjacencies( dum_ents, 3, false, hexes );
3112  if( MB_SUCCESS != result ) return result;
3113  assert( 1 == hexes.size() );
3114 
3115  hexes.insert( connects[0][i] );
3116  hexes.insert( connects[0][( i + 1 ) % 4] );
3117  hexes.insert( connects[3][i] );
3118  hexes.insert( connects[3][( i + 1 ) % 4] );
3119  dum_ents.clear();
3120  result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents );
3121  if( MB_SUCCESS != result ) return result;
3122  assert( dum_ents.size() == 1 );
3123  side_quads[0].push_back( *dum_ents.begin() );
3124 
3125  hexes.erase( connects[0][i] );
3126  hexes.erase( connects[0][( i + 1 ) % 4] );
3127  hexes.insert( connects[2][i] );
3128  hexes.insert( connects[2][( i + 1 ) % 4] );
3129  dum_ents.clear();
3130  result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents );
3131  if( MB_SUCCESS != result ) return result;
3132  side_quads[1].push_back( *dum_ents.begin() );
3133  }
3134 
3135  return MB_SUCCESS;
3136 }
3137 
3138 ErrorCode DualTool::fs_get_quads( EntityHandle odedge,
3139  EntityHandle* quads,
3140  EntityHandle* hexes,
3141  std::vector< EntityHandle >* connects )
3142 {
3143  // need to get the three quads along the chord
3144  EntityHandle chord = get_dual_hyperplane( odedge );
3145  if( 0 == chord ) return MB_FAILURE;
3146 
3147  std::vector< EntityHandle > edges;
3148  ErrorCode result = mbImpl->get_entities_by_handle( chord, edges );
3149  if( MB_FAILURE == result ) return result;
3150 
3151  std::vector< EntityHandle >::iterator vit = std::find( edges.begin(), edges.end(), odedge );
3152  // shouldn't be first or last edge on chord
3153  if( vit == edges.end() || *edges.begin() == *vit || *edges.rbegin() == *vit ) return MB_FAILURE;
3154 
3155  // get quads/connectivity for first 3 quads
3156  quads[0] = get_dual_entity( *( vit - 1 ) );
3157  quads[1] = get_dual_entity( *vit );
3158  quads[2] = get_dual_entity( *( vit + 1 ) );
3159  for( int i = 0; i < 3; i++ )
3160  {
3161  result = mbImpl->get_connectivity( &quads[i], 1, connects[i], true );
3162  if( MB_SUCCESS != result ) return result;
3163  }
3164 
3165  Range tmph;
3166  result = mbImpl->get_adjacencies( quads, 2, 3, false, tmph );
3167  if( MB_SUCCESS != result ) return result;
3168  assert( tmph.size() == 1 );
3169  hexes[0] = *tmph.begin();
3170 
3171  tmph.clear();
3172  result = mbImpl->get_adjacencies( &quads[1], 2, 3, false, tmph );
3173  if( MB_SUCCESS != result ) return result;
3174  assert( tmph.size() == 1 );
3175  hexes[1] = *tmph.begin();
3176 
3177  return MB_SUCCESS;
3178 }
3179 
3180 ErrorCode DualTool::delete_whole_dual()
3181 {
3182  // delete dual hyperplanes
3183  Range dual_surfs, dual_curves;
3184  ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs );RR;
3185  result = mbImpl->delete_entities( dual_surfs );RR;
3186  result = this->get_dual_hyperplanes( mbImpl, 1, dual_curves );RR;
3187  result = mbImpl->delete_entities( dual_curves );RR;
3188 
3189  // gather up all dual entities
3190  Range dual_ents;
3191  result = mbImpl->get_entities_by_type_and_tag( 0, MBVERTEX, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
3192  result = mbImpl->get_entities_by_type_and_tag( 0, MBEDGE, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
3193  result = mbImpl->get_entities_by_type_and_tag( 0, MBPOLYGON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
3194  result =
3195  mbImpl->get_entities_by_type_and_tag( 0, MBPOLYHEDRON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
3196 
3197  // delete them, in reverse order of dimension
3198  ErrorCode tmp_result;
3199  for( Range::reverse_iterator rit = dual_ents.rbegin(); rit != dual_ents.rend(); ++rit )
3200  {
3201  tmp_result = mbImpl->delete_entities( &( *rit ), 1 );
3202  if( MB_SUCCESS != tmp_result ) result = tmp_result;
3203  }
3204  RR;
3205 
3206  // delete dual-related tags
3207  if( 0 != dualSurfaceTag )
3208  {
3209  tmp_result = mbImpl->tag_delete( dualSurfaceTag );
3210  if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3211  }
3212  if( 0 != dualCurveTag )
3213  {
3214  tmp_result = mbImpl->tag_delete( dualCurveTag );
3215  if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3216  }
3217  if( 0 != dualEntityTag )
3218  {
3219  tmp_result = mbImpl->tag_delete( dualEntityTag );
3220  if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3221  }
3222  if( 0 != extraDualEntityTag )
3223  {
3224  tmp_result = mbImpl->tag_delete( extraDualEntityTag );
3225  if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3226  }
3227  if( 0 != dualGraphicsPointTag )
3228  {
3229  tmp_result = mbImpl->tag_delete( dualGraphicsPointTag );
3230  if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3231  }
3232 
3233  return MB_SUCCESS;
3234 }
3235 
3236 ErrorCode DualTool::check_dual_adjs()
3237 {
3238  // check primal-dual correspondence
3239 
3240  // get the primal entities
3241  Range pents[4];
3242  ErrorCode result = mbImpl->get_entities_by_type( 0, MBHEX, pents[3] );RR;
3243  for( int i = 2; i >= 0; i-- )
3244  {
3245  result = mbImpl->get_adjacencies( pents[3], 2, false, pents[2], Interface::UNION );RR;
3246  }
3247 
3248  // for each primal entity of dimension pd
3249 #define PRENT( ent ) CN::EntityTypeName( TYPE_FROM_HANDLE( ent ) ) << " " << ID_FROM_HANDLE( ent )
3250  ErrorCode overall_result = MB_SUCCESS;
3251  for( int pd = 1; pd <= 3; pd++ )
3252  {
3253  for( Range::iterator prit = pents[pd].begin(); prit != pents[pd].end(); ++prit )
3254  {
3255  // get corresponding dual entity of dimension dd = 3-pd
3256  EntityHandle dual_ent = get_dual_entity( *prit );
3257  if( 0 == dual_ent ) std::cerr << "Problem getting dual entity for " << PRENT( *prit ) << std::endl;
3258 
3259  // for each sub dimension sd = 0..pd-1
3260  for( int sd = 0; sd < pd; sd++ )
3261  {
3262  Range R1, R2, R3;
3263  // R1 = entities bounding primal entity of dim sd
3264  result = mbImpl->get_adjacencies( &( *prit ), 1, sd, false, R1 );RR;
3265 
3266  // R2 = entities bounded by dual entity, of dim 3-sd
3267  result = mbImpl->get_adjacencies( &dual_ent, 1, 3 - sd, false, R2 );RR;
3268 
3269  if( R1.size() != R2.size() )
3270  {
3271  std::cerr << PRENT( *prit ) << ": number of adj ents in "
3272  << "primal/dual don't agree for dimension " << sd << "." << std::endl;
3273  overall_result = MB_FAILURE;
3274  }
3275 
3276  // for each entity in R1, get its dual and look for it in R2
3277  for( Range::iterator r1it = R1.begin(); r1it != R1.end(); ++r1it )
3278  {
3279  EntityHandle tmp_dual = get_dual_entity( *r1it );
3280  if( R2.find( tmp_dual ) == R2.end() )
3281  {
3282  std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r1it ) << " isn't adjacent in dual."
3283  << std::endl;
3284  overall_result = MB_FAILURE;
3285  }
3286  }
3287  // ditto for R2
3288  for( Range::iterator r2it = R2.begin(); r2it != R2.end(); ++r2it )
3289  {
3290  EntityHandle tmp_prim = get_dual_entity( *r2it );
3291  if( R1.find( tmp_prim ) == R1.end() )
3292  {
3293  std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r2it ) << " isn't adjacent in primal."
3294  << std::endl;
3295  overall_result = MB_FAILURE;
3296  }
3297  }
3298  }
3299  }
3300  }
3301 
3302  return overall_result;
3303 }
3304 
3305 } // namespace moab