MOAB: Mesh Oriented datABase  (version 5.2.1)
DualTool.cpp
Go to the documentation of this file.
00001 /**
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 #include "moab/DualTool.hpp"
00017 #include "moab/Range.hpp"
00018 // using Core for call to check_adjacencies
00019 #include "moab/Core.hpp"
00020 #include "Internals.hpp"
00021 #include "MBTagConventions.hpp"
00022 #include "moab/Skinner.hpp"
00023 #include "moab/Core.hpp"
00024 #include "moab/MeshTopoUtil.hpp"
00025 #include "AEntityFactory.hpp"
00026 #include "moab/CN.hpp"
00027 #include <string>
00028 #include <algorithm>
00029 #include <iostream>
00030 #include <sstream>
00031 #include <assert.h>
00032 
00033 #define RR \
00034     if( MB_SUCCESS != result ) return result
00035 #define SWAP( a, b )                    \
00036     {                                   \
00037         EntityHandle tmp_ent = a;       \
00038         a                    = b;       \
00039         b                    = tmp_ent; \
00040     }
00041 
00042 namespace moab
00043 {
00044 
00045 bool debug    = false;
00046 bool debug_ap = false;
00047 
00048 //! tag name for dual surfaces
00049 const char* DualTool::DUAL_SURFACE_TAG_NAME = "DUAL_SURFACE";
00050 
00051 //! tag name for dual curves
00052 const char* DualTool::DUAL_CURVE_TAG_NAME = "DUAL_CURVE";
00053 
00054 //! tag name for dual cells
00055 const char* DualTool::IS_DUAL_CELL_TAG_NAME = "__IS_DUAL_CELL";
00056 
00057 //! tag name for dual entities
00058 const char* DualTool::DUAL_ENTITY_TAG_NAME = "__DUAL_ENTITY";
00059 
00060 //! tag name for extra dual entities
00061 const char* DualTool::EXTRA_DUAL_ENTITY_TAG_NAME = "__EXTRA_DUAL_ENTITY";
00062 
00063 //! tag name for graphics point
00064 const char* DualTool::DUAL_GRAPHICS_POINT_TAG_NAME = "__DUAL_GRAPHICS_POINT";
00065 
00066 // const int DualTool::GP_SIZE = 20;
00067 
00068 DualTool::DualTool( Interface* impl ) : mbImpl( impl )
00069 {
00070     EntityHandle dum_handle = 0;
00071     ErrorCode result;
00072 
00073     result = mbImpl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dualSurfaceTag,
00074                                      MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle );
00075     assert( MB_SUCCESS == result );
00076 
00077     result = mbImpl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dualCurveTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00078                                      &dum_handle );
00079     assert( MB_SUCCESS == result );
00080 
00081     unsigned int dummy = 0;
00082     result             = mbImpl->tag_get_handle( IS_DUAL_CELL_TAG_NAME, 1, MB_TYPE_INTEGER, isDualCellTag,
00083                                      MB_TAG_SPARSE | MB_TAG_CREAT, &dummy );
00084     assert( MB_SUCCESS == result );
00085 
00086     result = mbImpl->tag_get_handle( DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, dualEntityTag,
00087                                      MB_TAG_DENSE | MB_TAG_CREAT, &dum_handle );
00088     assert( MB_SUCCESS == result );
00089 
00090     result = mbImpl->tag_get_handle( EXTRA_DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, extraDualEntityTag,
00091                                      MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle );
00092     assert( MB_SUCCESS == result );
00093 
00094     static const char dum_name[CATEGORY_TAG_SIZE] = { 0 };
00095     result = mbImpl->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, categoryTag,
00096                                      MB_TAG_SPARSE | MB_TAG_CREAT, dum_name );
00097     assert( MB_SUCCESS == result );
00098 
00099     DualTool::GraphicsPoint dum_pt( 0.0, 0.0, 0.0, -1 );
00100     result = mbImpl->tag_get_handle( DUAL_GRAPHICS_POINT_TAG_NAME, sizeof( DualTool::GraphicsPoint ), MB_TYPE_DOUBLE,
00101                                      dualGraphicsPointTag, MB_TAG_DENSE | MB_TAG_CREAT | MB_TAG_BYTES, &dum_pt );
00102     assert( MB_SUCCESS == result );
00103 
00104     globalIdTag = mbImpl->globalId_tag();
00105 
00106     if( MB_SUCCESS == result ) {}  // empty statement to get rid of warning.
00107 
00108     maxHexId = -1;
00109 }
00110 
00111 DualTool::~DualTool() {}
00112 
00113 //! construct the dual entities for the entire mesh
00114 ErrorCode DualTool::construct_dual( EntityHandle* entities, const int num_entities )
00115 {
00116     // allocate a dual entity for each primal entity in the mesh, starting
00117     // with highest dimension and working downward; do each dimension in a separate code
00118     // block, since they're all handled slightly differently
00119 
00120     Range regions, faces, edges, vertices;
00121     ErrorCode result;
00122 
00123     if( NULL == entities || 0 == num_entities )
00124     {
00125 
00126         // first, construct all the aentities, since they're currently needed to
00127         // compute the dual
00128         result = mbImpl->get_entities_by_dimension( 0, 0, vertices );
00129         if( MB_SUCCESS != result ) return result;
00130 
00131         result = MeshTopoUtil( mbImpl ).construct_aentities( vertices );
00132         if( MB_SUCCESS != result ) return result;
00133 
00134         // get all edges, faces and regions now, so we don't need to filter out dual
00135         // entities later
00136 
00137         result = mbImpl->get_entities_by_dimension( 0, 1, edges );
00138         if( MB_SUCCESS != result ) return result;
00139         result = mbImpl->get_entities_by_dimension( 0, 2, faces );
00140         if( MB_SUCCESS != result ) return result;
00141         result = mbImpl->get_entities_by_dimension( 0, 3, regions );
00142         if( MB_SUCCESS != result ) return result;
00143 
00144         // get the max global id for hexes, we'll need for modification ops
00145         std::vector< int > gid_vec( regions.size() );
00146         result = mbImpl->tag_get_data( globalId_tag(), regions, &gid_vec[0] );
00147         if( MB_SUCCESS != result ) return result;
00148         maxHexId = -1;
00149         Range::iterator rit;
00150         unsigned int i;
00151         for( rit = regions.begin(), i = 0; rit != regions.end(); ++rit, i++ )
00152         {
00153             if( gid_vec[i] > maxHexId && mbImpl->type_from_handle( *rit ) == MBHEX ) maxHexId = gid_vec[i];
00154         }
00155     }
00156     else
00157     {
00158         // get entities of various dimensions adjacent to these
00159         result = mbImpl->get_adjacencies( entities, num_entities, 0, true, vertices, Interface::UNION );
00160         if( MB_SUCCESS != result ) return result;
00161         result = mbImpl->get_adjacencies( entities, num_entities, 1, true, edges, Interface::UNION );
00162         if( MB_SUCCESS != result ) return result;
00163         result = mbImpl->get_adjacencies( entities, num_entities, 2, true, faces, Interface::UNION );
00164         if( MB_SUCCESS != result ) return result;
00165         result = mbImpl->get_adjacencies( entities, num_entities, 3, true, regions, Interface::UNION );
00166         if( MB_SUCCESS != result ) return result;
00167     }
00168 
00169     Range dual_verts;
00170     result = construct_dual_vertices( regions, dual_verts );
00171     if( MB_SUCCESS != result || dual_verts.size() != regions.size() ) return result;
00172     if( debug ) std::cout << "Constructed " << dual_verts.size() << " dual vertices." << std::endl;
00173 
00174     // don't really need dual edges, but construct 'em anyway
00175     Range dual_edges;
00176     result = construct_dual_edges( faces, dual_edges );
00177     if( MB_SUCCESS != result || dual_edges.size() != faces.size() ) return result;
00178     if( debug ) std::cout << "Constructed " << dual_edges.size() << " dual edges." << std::endl;
00179 
00180     // construct dual faces
00181     Range dual_faces;
00182     result = construct_dual_faces( edges, dual_faces );
00183     if( MB_SUCCESS != result || dual_faces.size() != edges.size() ) return result;
00184     if( debug ) std::cout << "Constructed " << dual_faces.size() << " dual faces." << std::endl;
00185 
00186     // construct dual cells
00187     Range dual_cells;
00188     result = construct_dual_cells( vertices, dual_cells );
00189     if( MB_SUCCESS != result || dual_cells.size() != vertices.size() ) return result;
00190     if( debug ) std::cout << "Constructed " << dual_cells.size() << " dual cells." << std::endl;
00191 
00192     return MB_SUCCESS;
00193 }
00194 
00195 ErrorCode DualTool::construct_dual_vertices( const Range& all_regions, Range& dual_ents )
00196 {
00197     if( all_regions.empty() ) return MB_SUCCESS;
00198 
00199     // make sure they're all regions
00200     assert( 3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.begin() ) ) &&
00201             3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.rbegin() ) ) );
00202 
00203     Range::const_iterator rit;
00204     EntityHandle dual_ent;
00205     ErrorCode tmp_result = MB_SUCCESS;
00206     ErrorCode result     = MB_SUCCESS;
00207 
00208     for( rit = all_regions.begin(); rit != all_regions.end(); ++rit )
00209     {
00210         if( tmp_result != MB_SUCCESS ) result = tmp_result;
00211 
00212         tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
00213         if( MB_SUCCESS == tmp_result && 0 != dual_ent )
00214         {
00215             dual_ents.insert( dual_ent );
00216             continue;
00217         }
00218         else if( MB_SUCCESS != tmp_result )
00219             continue;
00220 
00221         tmp_result = construct_dual_vertex( *rit, dual_ent, false, true );
00222         if( MB_SUCCESS != tmp_result ) continue;
00223 
00224         // save it in the list of new dual ents
00225         dual_ents.insert( dual_ent );
00226     }
00227 
00228     return result;
00229 }
00230 
00231 ErrorCode DualTool::construct_dual_vertex( EntityHandle entity, EntityHandle& dual_ent, const bool extra,
00232                                            const bool add_graphics_pt )
00233 {
00234     // no dual entity; construct one; first need the avg coordinates
00235     unsigned int is_dual = 0x1;
00236     double avg_pos[3];
00237     ErrorCode result = MeshTopoUtil( mbImpl ).get_average_position( entity, avg_pos );
00238     if( MB_SUCCESS != result ) return result;
00239 
00240     // now construct the new dual entity
00241     result = mbImpl->create_vertex( avg_pos, dual_ent );
00242     if( MB_SUCCESS != result ) return result;
00243 
00244     // tag it indicating it's a dual entity
00245     result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
00246     if( MB_SUCCESS != result ) return result;
00247 
00248     // tag the primal entity with its dual entity and vica versa
00249     if( extra )
00250         result = mbImpl->tag_set_data( extraDualEntity_tag(), &( entity ), 1, &dual_ent );
00251     else
00252         result = mbImpl->tag_set_data( dualEntity_tag(), &( entity ), 1, &dual_ent );
00253     if( MB_SUCCESS != result ) return result;
00254 
00255     result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( entity ) );
00256     if( MB_SUCCESS != result ) return result;
00257 
00258     if( add_graphics_pt )
00259         // put a graphics point on that vertex too
00260         result = add_graphics_point( dual_ent, avg_pos );
00261 
00262     return result;
00263 }
00264 
00265 ErrorCode DualTool::add_graphics_point( EntityHandle entity, double* avg_pos )
00266 {
00267     // add a graphics pt, placed at the same position as the vertex
00268     double my_pos[3];
00269     ErrorCode result;
00270 
00271     if( NULL == avg_pos )
00272     {
00273         result = MeshTopoUtil( mbImpl ).get_average_position( entity, my_pos );
00274         if( MB_SUCCESS != result ) return result;
00275     }
00276     else
00277         for( int i = 0; i < 3; i++ )
00278             my_pos[i] = avg_pos[i];
00279 
00280     DualTool::GraphicsPoint dum_pt( my_pos, -1 );
00281     result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), &entity, 1, &dum_pt );
00282     return result;
00283 }
00284 
00285 ErrorCode DualTool::construct_dual_edges( const Range& all_faces, Range& dual_ents )
00286 {
00287     if( all_faces.empty() ) return MB_SUCCESS;
00288 
00289     // make sure they're all faces
00290     assert( 2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.begin() ) ) &&
00291             2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.rbegin() ) ) );
00292 
00293     Range::const_iterator rit;
00294     EntityHandle dual_ent;
00295     unsigned int is_dual = 0x1;
00296     ErrorCode tmp_result = MB_SUCCESS;
00297     ErrorCode result     = MB_SUCCESS;
00298 
00299     for( rit = all_faces.begin(); rit != all_faces.end(); ++rit )
00300     {
00301         if( tmp_result != MB_SUCCESS ) result = tmp_result;
00302 
00303         tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
00304         if( MB_SUCCESS == tmp_result && 0 != dual_ent )
00305         {
00306             dual_ents.insert( dual_ent );
00307             continue;
00308         }
00309 
00310         // no dual entity; construct one; get the bounding regions
00311         std::vector< EntityHandle > out_ents;
00312         tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, out_ents );
00313         if( MB_SUCCESS != tmp_result || out_ents.empty() ) continue;
00314 
00315         // get the dual vertices
00316         std::vector< EntityHandle > dual_verts( out_ents.size() );
00317         tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &out_ents[0], out_ents.size(), &dual_verts[0] );
00318         if( MB_SUCCESS != tmp_result ) continue;
00319         assert( dual_verts.size() <= 2 );
00320 
00321         double avg_pos[3];
00322         bool bdy_face = ( dual_verts.size() == 1 ? true : false );
00323         if( bdy_face )
00324         {
00325             // boundary face - make a dual vertex at the face center and put in list
00326             tmp_result = construct_dual_vertex( *rit, dual_ent, true, true );
00327 
00328             // put it on vertex list
00329             dual_verts.push_back( dual_ent );
00330         }
00331 
00332         assert( dual_verts.size() == 2 );
00333 
00334         // now create the dual edge
00335         tmp_result = mbImpl->create_element( MBEDGE, &dual_verts[0], 2, dual_ent );
00336         if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue;
00337 
00338         // save it in the list of new dual ents
00339         dual_ents.insert( dual_ent );
00340 
00341         // tag the primal entity with its dual entity and vica versa
00342         tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
00343         if( MB_SUCCESS != tmp_result ) continue;
00344 
00345         tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
00346         if( MB_SUCCESS != tmp_result ) continue;
00347 
00348         // tag the edge indicating it's a dual entity
00349         tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
00350         if( MB_SUCCESS != tmp_result ) continue;
00351 
00352         // add a graphics point to the edge; position depends on whether it's a
00353         // bdy face (mid-pt of dual edge) or not (mid-pt of primal face)
00354         if( bdy_face )
00355             tmp_result = add_graphics_point( dual_ent );
00356         else
00357         {
00358             // get the face's position
00359             tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos );
00360             if( MB_SUCCESS != tmp_result ) continue;
00361             tmp_result = add_graphics_point( dual_ent, avg_pos );
00362         }
00363         if( MB_SUCCESS != tmp_result ) continue;
00364     }
00365 
00366     return result;
00367 }
00368 
00369 ErrorCode DualTool::construct_dual_faces( const Range& all_edges, Range& dual_ents )
00370 {
00371     if( all_edges.empty() ) return MB_SUCCESS;
00372 
00373     // make sure they're all edges
00374     assert( 1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.begin() ) ) &&
00375             1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.rbegin() ) ) );
00376 
00377     Range::const_iterator rit;
00378     EntityHandle dual_ent;
00379     unsigned int is_dual = 0x1;
00380     ErrorCode tmp_result = MB_SUCCESS;
00381     ErrorCode result     = MB_SUCCESS;
00382     Range equiv_edges;
00383 #define TRC                        \
00384     if( MB_SUCCESS != tmp_result ) \
00385     {                              \
00386         result = tmp_result;       \
00387         continue;                  \
00388     }
00389     for( rit = all_edges.begin(); rit != all_edges.end(); ++rit )
00390     {
00391 
00392         tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
00393         if( MB_SUCCESS == tmp_result && 0 != dual_ent )
00394         {
00395             dual_ents.insert( dual_ent );
00396             continue;
00397         }
00398 
00399         // no dual entity; construct one; get the dual vertices bounding the edge in radial order,
00400         // then construct the dual face
00401         std::vector< EntityHandle > rad_dverts;
00402         bool bdy_edge;
00403         tmp_result = get_radial_dverts( *rit, rad_dverts, bdy_edge );
00404         TRC if( rad_dverts.empty() ) continue;
00405 
00406         tmp_result = mbImpl->create_element( MBPOLYGON, &rad_dverts[0], rad_dverts.size(), dual_ent );
00407         TRC
00408 
00409             // tag it indicating it's a dual entity, and tag primal/dual with dual/primal
00410             tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
00411         TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
00412         TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
00413         TRC
00414 
00415             // save it in the list of new dual ents
00416             dual_ents.insert( dual_ent );
00417 
00418         // add a graphics point to the cell; position depends on whether it's a
00419         // bdy cell (mid-pt of cell's vertices) or not (mid-pt of primal edge)
00420         double avg_pos[3];
00421         tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos );
00422         TRC if( bdy_edge )
00423         {
00424 
00425             // add a new dual edge betw last 2 verts
00426             EntityHandle new_edge;
00427             tmp_result     = mbImpl->create_element( MBEDGE, &rad_dverts[rad_dverts.size() - 2], 2, new_edge );
00428             TRC tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &new_edge, 1, &is_dual );
00429             TRC
00430 
00431                 // tag the new dual edge with the primal edge as it's dual entity; primal
00432                 // edge IS NOT likewise tagged, since it's already tagged with the 2cell
00433                 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &new_edge, 1, &( *rit ) );
00434             TRC
00435 
00436                 // add a graphics pt, position is center of primal edge
00437                 tmp_result = add_graphics_point( dual_ent );
00438             TRC tmp_result = add_graphics_point( new_edge, avg_pos );
00439             TRC
00440         }
00441 
00442         else
00443         {
00444             // if inside, point goes on the 2cell, at primal edge mid-pt
00445             tmp_result = add_graphics_point( dual_ent, avg_pos );
00446             TRC
00447         }
00448 
00449         // check to see whether we have equiv entities; if we find any, save for later fixup
00450         Range dum_edges, dum_poly( dual_ent, dual_ent );
00451         tmp_result = mbImpl->get_adjacencies( dum_poly, 1, false, dum_edges );
00452         if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result )
00453         {
00454             // we do - need to add adjacencies to disambiguate; use the primal
00455             equiv_edges.merge( dum_edges );
00456         }
00457     }
00458 
00459     if( !equiv_edges.empty() ) result = check_dual_equiv_edges( equiv_edges );
00460 
00461     return result;
00462 }
00463 
00464 ErrorCode DualTool::check_dual_equiv_edges( Range& dual_edges )
00465 {
00466     // fix equivalent dual edges (i.e. edges whose vertices define multiple edges)
00467     // by explicitly adding adjacencies to containing polygons; adjacent polygons
00468     // found by going through primal
00469     ErrorCode tmp_result, result = MB_SUCCESS;
00470 
00471     Range all_dedges( dual_edges );
00472     // first, go through all dual edges and find equivalent edges (by looking for
00473     // up-adjacent edges on the vertices of each edge)
00474     for( Range::iterator rit = dual_edges.begin(); rit != dual_edges.end(); ++rit )
00475     {
00476         Range connect, dum_range( *rit, *rit );
00477         tmp_result = mbImpl->get_adjacencies( dum_range, 0, false, connect );
00478         if( MB_SUCCESS != tmp_result ) continue;
00479         tmp_result = mbImpl->get_adjacencies( connect, 1, false, all_dedges, Interface::UNION );
00480         if( MB_SUCCESS != tmp_result ) continue;
00481     }
00482 
00483     // save a copy for checking later
00484     Range save_all_2cells;
00485 
00486     // go through each edge
00487     while( !all_dedges.empty() )
00488     {
00489         EntityHandle this_edge = *all_dedges.begin();
00490         all_dedges.erase( all_dedges.begin() );
00491 
00492         const EntityHandle* connect;
00493         int num_connect;
00494         result = mbImpl->get_connectivity( this_edge, connect, num_connect );
00495         if( MB_SUCCESS != result ) continue;
00496 
00497         Range dum_edges, verts;
00498         verts.insert( connect[0] );
00499         verts.insert( connect[1] );
00500         tmp_result = mbImpl->get_adjacencies( verts, 1, false, dum_edges );
00501         if( MB_SUCCESS != tmp_result )
00502         {
00503             result = tmp_result;
00504             continue;
00505         }
00506         if( dum_edges.size() == 1 )
00507         {
00508             // not an equiv edge - already removed from list, so just continue
00509             continue;
00510         }
00511 
00512         // ok, have an equiv entity - fix by looking through primal
00513         // pre-get the primal of these
00514         EntityHandle dedge_quad;
00515         tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &this_edge, 1, &dedge_quad );
00516         if( MB_SUCCESS != tmp_result )
00517         {
00518             result = tmp_result;
00519             continue;
00520         }
00521 
00522         if( MBQUAD == mbImpl->type_from_handle( dedge_quad ) )
00523         {
00524 
00525             // get the primal edges adj to quad
00526             Range dum_quad_range( dedge_quad, dedge_quad ), adj_pedges;
00527             tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_pedges );
00528             if( MB_SUCCESS != tmp_result )
00529             {
00530                 result = tmp_result;
00531                 continue;
00532             }
00533             // get the dual 2cells corresponding to those pedges
00534             std::vector< EntityHandle > dcells;
00535             dcells.resize( adj_pedges.size() );
00536             tmp_result = mbImpl->tag_get_data( dualEntity_tag(), adj_pedges, &dcells[0] );
00537             if( MB_SUCCESS != tmp_result )
00538             {
00539                 result = tmp_result;
00540                 continue;
00541             }
00542             // now add explicit adjacencies from the dedge to those dcells
00543             std::vector< EntityHandle >::iterator vit;
00544             for( vit = dcells.begin(); vit != dcells.end(); ++vit )
00545             {
00546                 save_all_2cells.insert( *vit );
00547 
00548                 assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) );
00549                 tmp_result = mbImpl->add_adjacencies( this_edge, &( *vit ), 1, false );
00550                 if( MB_SUCCESS != tmp_result )
00551                 {
00552                     result = tmp_result;
00553                     continue;
00554                 }
00555                 // check that there are really adjacencies and *vit is in them
00556                 const EntityHandle* adjs;
00557                 int num_adjs;
00558                 tmp_result = reinterpret_cast< Core* >( mbImpl )->a_entity_factory()->get_adjacencies( this_edge, adjs,
00559                                                                                                        num_adjs );
00560                 if( NULL == adjs || std::find( adjs, adjs + num_adjs, *vit ) == adjs + num_adjs )
00561                     std::cout << "Add_adjacencies failed in construct_dual_faces." << std::endl;
00562             }
00563         }
00564         else
00565         {
00566             // else, have a dual edge representing a bdy edge - tie directly to
00567             // dual entity if its dual entity
00568             EntityHandle bdy_dcell;
00569             tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &dedge_quad, 1, &bdy_dcell );
00570             TRC assert( MBPOLYGON == mbImpl->type_from_handle( bdy_dcell ) );
00571 
00572             tmp_result = mbImpl->add_adjacencies( this_edge, &bdy_dcell, 1, false );
00573             if( MB_SUCCESS != tmp_result )
00574             {
00575                 result = tmp_result;
00576                 continue;
00577             }
00578         }
00579     }
00580 
00581     // sanity check - look for adj edges again, and check for equiv entities
00582     for( Range::iterator vit = save_all_2cells.begin(); vit != save_all_2cells.end(); ++vit )
00583     {
00584         Range adj_edges, dum_quad_range;
00585         dum_quad_range.insert( *vit );
00586         assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) );
00587         tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_edges );
00588         if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result )
00589         {
00590             std::cout << "Multiple entities returned for polygon " << mbImpl->id_from_handle( *vit ) << "."
00591                       << std::endl;
00592             continue;
00593         }
00594     }
00595     // success!
00596     return result;
00597 }
00598 
00599 ErrorCode DualTool::construct_dual_cells( const Range& all_verts, Range& dual_ents )
00600 {
00601     if( all_verts.empty() ) return MB_SUCCESS;
00602 
00603     // make sure they're all edges
00604     assert( 0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.begin() ) ) &&
00605             0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.rbegin() ) ) );
00606 
00607     Range::const_iterator rit;
00608     EntityHandle dual_ent;
00609     unsigned int is_dual = 0x1;
00610     ErrorCode tmp_result = MB_SUCCESS;
00611     ErrorCode result     = MB_SUCCESS;
00612     std::vector< EntityHandle > edges, dfaces;
00613 
00614     for( rit = all_verts.begin(); rit != all_verts.end(); ++rit )
00615     {
00616         if( tmp_result != MB_SUCCESS ) result = tmp_result;
00617 
00618         tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
00619         if( MB_SUCCESS == tmp_result && 0 != dual_ent )
00620         {
00621             dual_ents.insert( dual_ent );
00622             continue;
00623         }
00624 
00625         // no dual entity; construct one; get the edges bounding the vertex
00626         edges.clear();
00627         dfaces.clear();
00628         tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 1, false, edges );
00629         if( MB_SUCCESS != tmp_result ) continue;
00630 
00631         // get the dual faces corresponding to the edges
00632         dfaces.resize( edges.size() );
00633         tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &edges[0], edges.size(), &dfaces[0] );
00634         if( MB_SUCCESS != tmp_result ) continue;
00635 
00636         // create the dual cell from those faces
00637         tmp_result = mbImpl->create_element( MBPOLYHEDRON, &dfaces[0], dfaces.size(), dual_ent );
00638         if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue;
00639 
00640         // save it in the list of new dual ents
00641         dual_ents.insert( dual_ent );
00642 
00643         // tag it indicating it's a dual entity
00644         tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
00645         if( MB_SUCCESS != tmp_result ) continue;
00646 
00647         // tag the primal entity with its dual entity and vica versa
00648         tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
00649         if( MB_SUCCESS != tmp_result ) continue;
00650         tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
00651         if( MB_SUCCESS != tmp_result ) continue;
00652     }
00653 
00654     return result;
00655 }
00656 
00657 //! given an edge handle, return a list of dual vertices in radial order
00658 //! around the edge
00659 ErrorCode DualTool::get_radial_dverts( const EntityHandle edge, std::vector< EntityHandle >& rad_dverts,
00660                                        bool& bdy_edge )
00661 {
00662     rad_dverts.clear();
00663 
00664     std::vector< EntityHandle > rad_faces, rad_ents;
00665     ErrorCode result = MeshTopoUtil( mbImpl ).star_entities( edge, rad_faces, bdy_edge, 0, &rad_ents );
00666     if( MB_SUCCESS != result ) return result;
00667 
00668     if( bdy_edge )
00669     {
00670         // if we're a bdy edge, change the order back to what DualTool expects
00671         rad_ents.push_back( *rad_faces.rbegin() );
00672         rad_ents.push_back( *rad_faces.begin() );
00673     }
00674 
00675     rad_dverts.resize( rad_ents.size() );
00676     for( unsigned int i = 0; i < rad_ents.size(); i++ )
00677     {
00678         EntityHandle dual_ent;
00679         result = mbImpl->tag_get_data( dualEntity_tag(), &rad_ents[i], 1, &dual_ent );
00680         if( !bdy_edge || i < rad_ents.size() - 2 )
00681             rad_dverts[i] = dual_ent;
00682         else
00683         {
00684             // fix up this entry
00685             assert( mbImpl->type_from_handle( dual_ent ) == MBEDGE );
00686 
00687             // get connectivity of that edge
00688             const EntityHandle* connect;
00689             int num_connect;
00690             result = mbImpl->get_connectivity( dual_ent, connect, num_connect );
00691             if( MB_SUCCESS != result ) return result;
00692 
00693             // we want the one that's not already on the list; reuse last_face
00694             int last_hex           = ( i == rad_ents.size() - 1 ? 0 : i - 1 );
00695             EntityHandle last_face = ( connect[0] == rad_dverts[last_hex] ? connect[1] : connect[0] );
00696             rad_dverts[i]          = last_face;
00697         }
00698     }
00699 
00700     return result;
00701 }
00702 
00703 //! construct the dual entities for a hex mesh, including dual surfaces & curves
00704 ErrorCode DualTool::construct_hex_dual( Range& entities )
00705 {
00706     std::vector< EntityHandle > evec;
00707     std::copy( entities.begin(), entities.end(), std::back_inserter( evec ) );
00708     return construct_hex_dual( &evec[0], evec.size() );
00709 }
00710 
00711 //! construct the dual entities for a hex mesh, including dual surfaces & curves
00712 ErrorCode DualTool::construct_hex_dual( EntityHandle* entities, const int num_entities )
00713 {
00714     // really quite simple:
00715 
00716     // construct the dual...
00717     ErrorCode result = construct_dual( entities, num_entities );
00718     if( MB_SUCCESS != result )
00719     {
00720         std::cerr << "Error constructing dual entities for primal entities." << std::endl;
00721         return result;
00722     }
00723 
00724     // now traverse to build 1d and 2d hyperplanes
00725     result = construct_dual_hyperplanes( 1, entities, num_entities );
00726     if( MB_SUCCESS != result )
00727     {
00728         std::cerr << "Problem traversing 1d hyperplanes." << std::endl;
00729         return result;
00730     }
00731 
00732     result = construct_dual_hyperplanes( 2, entities, num_entities );
00733     if( MB_SUCCESS != result )
00734     {
00735         std::cerr << "Problem traversing 2d hyperplanes." << std::endl;
00736         return result;
00737     }
00738 
00739     result = construct_hp_parent_child();
00740     if( MB_SUCCESS != result )
00741     {
00742         std::cerr << "Problem constructing parent/child relations between hyperplanes." << std::endl;
00743         return result;
00744     }
00745 
00746     // see?  simple, just like I said
00747     return MB_SUCCESS;
00748 }
00749 
00750 //! get the cells of the dual
00751 ErrorCode DualTool::get_dual_entities( const int dim, EntityHandle* entities, const int num_entities, Range& dual_ents )
00752 {
00753     if( 0 == isDualCell_tag() ) return MB_SUCCESS;
00754     if( 0 > dim || 3 < dim ) return MB_INDEX_OUT_OF_RANGE;
00755 
00756     unsigned int dum              = 0x1;
00757     const void* dum_ptr           = &dum;
00758     static EntityType dual_type[] = { MBVERTEX, MBEDGE, MBPOLYGON, MBPOLYHEDRON };
00759 
00760     Range dim_ents;
00761 
00762     ErrorCode result;
00763 
00764     if( 0 == entities || 0 == num_entities )
00765     {
00766         // just get all the dual entities of this dimension
00767         result = mbImpl->get_entities_by_type_and_tag( 0, dual_type[dim], &isDualCellTag, &dum_ptr, 1, dual_ents );
00768     }
00769     else
00770     {
00771         // else look for specific dual entities
00772         result = mbImpl->get_adjacencies( entities, num_entities, 3 - dim, false, dim_ents, Interface::UNION );
00773         if( MB_SUCCESS != result ) return result;
00774         std::vector< EntityHandle > dual_ents_vec( dim_ents.size() );
00775         result = mbImpl->tag_get_data( dualEntity_tag(), dim_ents, &dual_ents_vec[0] );
00776         if( MB_SUCCESS != result ) return result;
00777         std::copy( dual_ents_vec.begin(), dual_ents_vec.end(), range_inserter( dual_ents ) );
00778     }
00779 
00780     return result;
00781 }
00782 
00783 //! get the faces of the dual
00784 ErrorCode DualTool::get_dual_entities( const int dim, EntityHandle* entities, const int num_entities,
00785                                        std::vector< EntityHandle >& dual_ents )
00786 {
00787     Range tmp_range;
00788     ErrorCode result = get_dual_entities( dim, entities, num_entities, tmp_range );
00789     if( MB_SUCCESS != result ) return result;
00790 
00791     // dual_ents.insert(dual_ents.end(), tmp_range.begin(), tmp_range.end());
00792     dual_ents.reserve( dual_ents.size() + tmp_range.size() );
00793     for( Range::const_iterator it = tmp_range.begin(); it != tmp_range.end(); ++it )
00794     {
00795         dual_ents.push_back( *it );
00796     }
00797     return MB_SUCCESS;
00798 }
00799 
00800 ErrorCode DualTool::get_dual_hyperplanes( const Interface* impl, const int dim, Range& dual_ents )
00801 {
00802     if( dim != 1 && dim != 2 ) return MB_INDEX_OUT_OF_RANGE;
00803 
00804     Tag dual_tag;
00805     ErrorCode result;
00806 
00807     if( dim == 1 )
00808         result = impl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag );
00809     else
00810         result = impl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag );
00811 
00812     if( MB_SUCCESS == result )
00813         result = impl->get_entities_by_type_and_tag( 0, MBENTITYSET, &dual_tag, NULL, 1, dual_ents, Interface::UNION );
00814 
00815     return result;
00816 }
00817 
00818 ErrorCode DualTool::construct_dual_hyperplanes( const int dim, EntityHandle* entities, const int num_entities )
00819 {
00820     // this function traverses dual faces of input dimension, constructing
00821     // dual hyperplanes of them in sets as it goes
00822 
00823     // check various inputs
00824     int num_quads, num_hexes;
00825     if(
00826         // this function only makes sense for dim == 1 or dim == 2
00827         ( dim != 1 && dim != 2 ) ||
00828         // should either be quads or hexes around
00829         mbImpl->get_number_entities_by_type( 0, MBQUAD, num_quads ) != MB_SUCCESS ||
00830         mbImpl->get_number_entities_by_type( 0, MBHEX, num_hexes ) != MB_SUCCESS ||
00831         // if we're asking for 1d dual ents, should be quads around
00832         ( num_quads == 0 && dim == 1 ) ||
00833         // if we're asking for 2d dual ents, should be hexes around
00834         ( num_hexes == 0 && dim == 2 ) )
00835         return MB_FAILURE;
00836 
00837     Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() );
00838 
00839     // two stacks: one completely untreated entities, and the other untreated
00840     // entities on the current dual hyperplane
00841     std::vector< EntityHandle > tot_untreated;
00842 
00843     // put dual entities of this dimension on the untreated list
00844     ErrorCode result = get_dual_entities( dim, entities, num_entities, tot_untreated );
00845     if( MB_SUCCESS != result ) return result;
00846 
00847     // main part of traversal loop
00848     EntityHandle this_ent;
00849     EntityHandle this_hp;
00850 
00851     while( !tot_untreated.empty() )
00852     {
00853         if( debug && dim == 2 /*(tot_untreated.size()%report == 0)*/ )
00854             std::cout << "Untreated list size " << tot_untreated.size() << "." << std::endl;
00855 
00856         this_ent = tot_untreated.back();
00857         tot_untreated.pop_back();
00858         result = mbImpl->tag_get_data( hp_tag, &this_ent, 1, &this_hp );
00859         if( MB_SUCCESS != result && MB_TAG_NOT_FOUND != result ) return result;
00860 
00861         // d for this entity having a hyperplane assignment already
00862         else if( this_hp != 0 )
00863             continue;
00864 
00865         if( 1 == dim && check_1d_loop_edge( this_ent ) ) continue;
00866 
00867         // inner loop: traverse the hyperplane 'till we don't have any more
00868         result = traverse_hyperplane( hp_tag, this_hp, this_ent );
00869         if( MB_SUCCESS != result )
00870         {
00871             std::cout << "Failed to traverse hyperplane ";
00872             if( this_hp )
00873                 std::cout << mbImpl->id_from_handle( this_hp ) << "." << std::endl;
00874             else
00875                 std::cout << "0." << std::endl;
00876             return result;
00877         }
00878 
00879         // ok, now order the edges if it's a chord
00880         if( 1 == dim ) order_chord( this_hp );
00881     }
00882 
00883     return MB_SUCCESS;
00884 }
00885 
00886 ErrorCode DualTool::traverse_hyperplane( const Tag hp_tag, EntityHandle& this_hp, EntityHandle this_ent )
00887 {
00888     Range tmp_star, star, tmp_range, new_hyperplane_ents;
00889     std::vector< EntityHandle > hp_untreated;
00890     int dim = mbImpl->dimension_from_handle( this_ent );
00891     MeshTopoUtil mtu( mbImpl );
00892     this_hp = 0;
00893     ErrorCode result;
00894 
00895     unsigned short mark_val = 0x0;
00896     Tag mark_tag;
00897     result = mbImpl->tag_get_handle( "__hyperplane_mark", 1, MB_TYPE_BIT, mark_tag, MB_TAG_CREAT | MB_TAG_BIT );
00898     if( MB_SUCCESS != result ) return result;
00899     mark_val = 0x1;
00900 
00901     while( 0 != this_ent )
00902     {
00903         EntityHandle tmp_hp = get_dual_hyperplane( this_ent );
00904         if( 0 == this_hp && 0 != tmp_hp ) this_hp = tmp_hp;
00905 
00906         if( 0 == tmp_hp ) new_hyperplane_ents.insert( this_ent );
00907 
00908         if( debug && hp_untreated.size() % 10 == 0 )
00909             std::cout << "Dual surface " << this_hp << ", hp_untreated list size = " << hp_untreated.size() << "."
00910                       << std::endl;
00911 
00912         // get the 2nd order adjacencies through lower dimension
00913         tmp_range.clear();
00914         tmp_star.clear();
00915         star.clear();
00916         result = mtu.get_bridge_adjacencies( this_ent, dim - 1, dim, star );RR;
00917 
00918         // get the bridge adjacencies through higher dimension
00919         result = mtu.get_bridge_adjacencies( this_ent, dim + 1, dim, tmp_star );RR;
00920         tmp_range = subtract( star, tmp_star );
00921 
00922         for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit )
00923         {
00924             if( new_hyperplane_ents.find( *rit ) != new_hyperplane_ents.end() ) continue;
00925 
00926             // check for tag first, 'cuz it's probably faster than checking adjacencies
00927             // assign to avoid valgrind warning
00928             unsigned short tmp_mark = 0x0;
00929             result                  = mbImpl->tag_get_data( mark_tag, &( *rit ), 1, &tmp_mark );
00930             if( MB_SUCCESS == result && mark_val == tmp_mark ) continue;
00931 
00932             // if it's on the loop, it's not eligible
00933             if( 1 == dim && check_1d_loop_edge( *rit ) ) continue;
00934 
00935             // have one on this hp; just put it on the hp_untreated list for now,
00936             // will get tagged and put in the hp set later
00937             hp_untreated.push_back( *rit );
00938             result = mbImpl->tag_set_data( mark_tag, &( *rit ), 1, &mark_val );
00939             if( MB_SUCCESS != result ) return result;
00940         }
00941 
00942         // end of inner loop; get the next this_ent, or set to zero
00943         if( hp_untreated.empty() )
00944             this_ent = 0;
00945         else
00946         {
00947             this_ent = hp_untreated.back();
00948             hp_untreated.pop_back();
00949         }
00950     }
00951 
00952     if( debug_ap )
00953     {
00954         std::string hp_name;
00955         if( 2 == dim )
00956             hp_name = "sheet";
00957         else
00958             hp_name = "chord";
00959 
00960         if( 0 == this_hp )
00961             std::cout << "Constructed new " << hp_name << " with ";
00962         else
00963         {
00964             int this_id;
00965             result = mbImpl->tag_get_data( globalId_tag(), &this_hp, 1, &this_id );RR;
00966             std::cout << "Added to " << hp_name << " " << this_id << " ";
00967         }
00968         if( dim == 2 )
00969             std::cout << "edges:" << std::endl;
00970         else
00971             std::cout << "quads:" << std::endl;
00972         std::vector< EntityHandle > pents( new_hyperplane_ents.size() );
00973         result = mbImpl->tag_get_data( dualEntity_tag(), new_hyperplane_ents, &pents[0] );RR;
00974         for( std::vector< EntityHandle >::iterator vit = pents.begin(); vit != pents.end(); ++vit )
00975         {
00976             if( vit != pents.begin() ) std::cout << ", ";
00977             std::cout << mbImpl->id_from_handle( *vit );
00978         }
00979         std::cout << std::endl;
00980     }
00981 
00982     if( 0 == this_hp )
00983     {
00984         // ok, doesn't have one; make a new hyperplane
00985         int new_id = -1;
00986         result     = construct_new_hyperplane( dim, this_hp, new_id );
00987         if( MB_SUCCESS != result ) return result;
00988 
00989         if( debug_ap )
00990         {
00991             std::cout << "New ";
00992             if( 2 == dim )
00993                 std::cout << " sheet ";
00994             else
00995                 std::cout << " chord ";
00996             std::cout << new_id << " constructed." << std::endl;
00997         }
00998     }
00999 
01000     // set the hp_val for entities which didn't have one before
01001     std::vector< EntityHandle > hp_tags( new_hyperplane_ents.size() );
01002     std::fill( hp_tags.begin(), hp_tags.end(), this_hp );
01003     result = mbImpl->tag_set_data( hp_tag, new_hyperplane_ents, &hp_tags[0] );
01004     if( MB_SUCCESS != result ) return result;
01005     result = mbImpl->add_entities( this_hp, new_hyperplane_ents );
01006     if( MB_SUCCESS != result ) return result;
01007 
01008     // unmark the entities by removing the tag
01009     result = mbImpl->tag_delete( mark_tag );
01010     if( MB_SUCCESS != result ) return result;
01011 
01012     return MB_SUCCESS;
01013 }
01014 
01015 ErrorCode DualTool::order_chord( EntityHandle chord_set )
01016 {
01017     // re-order the 1cells in the set so they are in order along the chord
01018     // start by finding the vertex dual to a quad
01019     Range verts, one_cells;
01020     ErrorCode result = mbImpl->get_entities_by_dimension( chord_set, 1, one_cells );
01021     if( MB_SUCCESS != result || one_cells.empty() ) return MB_FAILURE;
01022 
01023     result = mbImpl->get_adjacencies( one_cells, 0, false, verts, Interface::UNION );
01024     if( MB_SUCCESS != result || verts.empty() ) return MB_FAILURE;
01025 
01026     EntityHandle last_vert = 0;
01027     for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit )
01028     {
01029         if( TYPE_FROM_HANDLE( get_dual_entity( *rit ) ) == MBQUAD )
01030         {
01031             last_vert = *rit;
01032             break;
01033         }
01034     }
01035     // if there's no vertex owned by a quad, just start with 1st one
01036     if( 0 == last_vert ) last_vert = *verts.begin();
01037 
01038     // now, skip from vertex to vertex, building a list of 1cells
01039     std::vector< EntityHandle > ordered_1cells;
01040     EntityHandle last_1cell = 0;
01041     Range dum1, dum2;
01042     const EntityHandle* connect;
01043     int num_connect;
01044     ErrorCode tmp_result = MB_SUCCESS;
01045     while( ordered_1cells.size() != one_cells.size() )
01046     {
01047         dum1   = one_cells;
01048         result = mbImpl->get_adjacencies( &last_vert, 1, 1, false, dum1 );
01049         if( 0 != last_1cell ) dum1.erase( last_1cell );
01050         // assert(1 == dum1.size());
01051         if( 0 != last_1cell && 1 != dum1.size() )
01052         {
01053             std::cerr << "unexpected size traversing chord." << std::endl;
01054             tmp_result = MB_FAILURE;
01055         }
01056 
01057         last_1cell = *dum1.begin();
01058         ordered_1cells.push_back( last_1cell );
01059         result = mbImpl->get_connectivity( last_1cell, connect, num_connect );RR;
01060         if( last_vert == connect[0] )
01061             last_vert = connect[1];
01062         else
01063             last_vert = connect[0];
01064     }
01065 
01066     // now have the 1cells in order, replace them in the set
01067     if( MB_SUCCESS == tmp_result )
01068     {
01069         result = mbImpl->remove_entities( chord_set, one_cells );RR;
01070         result = mbImpl->add_entities( chord_set, &ordered_1cells[0], ordered_1cells.size() );RR;
01071     }
01072 
01073     return MB_SUCCESS;
01074 }
01075 
01076 ErrorCode DualTool::construct_new_hyperplane( const int dim, EntityHandle& new_hyperplane, int& id )
01077 {
01078     ErrorCode result;
01079     if( 1 == dim )
01080         result = mbImpl->create_meshset( ( MESHSET_ORDERED | MESHSET_TRACK_OWNER ), new_hyperplane );
01081     else
01082         result = mbImpl->create_meshset( ( MESHSET_SET | MESHSET_TRACK_OWNER ), new_hyperplane );
01083     if( MB_SUCCESS != result ) return result;
01084 
01085     if( -1 == id )
01086     {
01087         Range all_hyperplanes;
01088         result = get_dual_hyperplanes( mbImpl, dim, all_hyperplanes );RR;
01089         std::vector< int > gids( all_hyperplanes.size() );
01090         result = mbImpl->tag_get_data( globalIdTag, all_hyperplanes, ( gids.empty() ) ? NULL : &gids[0] );RR;
01091         for( unsigned int i = 0; i < gids.size(); i++ )
01092             if( gids[i] > id ) id = gids[i];
01093         id++;
01094         if( 0 == id ) id++;
01095     }
01096 
01097     result = mbImpl->tag_set_data( globalId_tag(), &new_hyperplane, 1, &id );RR;
01098     Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() );
01099     result     = mbImpl->tag_set_data( hp_tag, &new_hyperplane, 1, &new_hyperplane );RR;
01100 
01101     // assign a category name to these sets
01102     static const char dual_category_names[2][CATEGORY_TAG_SIZE] = { "Chord\0", "Sheet\0" };
01103 
01104     result = mbImpl->tag_set_data( categoryTag, &new_hyperplane, 1, dual_category_names[dim - 1] );
01105 
01106     return result;
01107 }
01108 
01109 bool DualTool::check_1d_loop_edge( EntityHandle this_ent )
01110 {
01111     // make sure it's an edge
01112     if( MBEDGE != mbImpl->type_from_handle( this_ent ) ) return false;
01113 
01114     // also has to be a dual entity
01115     unsigned int dum;
01116     ErrorCode result = mbImpl->tag_get_data( isDualCell_tag(), &this_ent, 1, &dum );
01117     if( MB_SUCCESS != result || dum != 0x1 ) return false;
01118 
01119     const EntityHandle* verts;
01120     EntityHandle vert_tags[2];
01121     int num_verts;
01122     result = mbImpl->get_connectivity( this_ent, verts, num_verts );
01123     if( MB_SUCCESS != result ) return false;
01124 
01125     result = mbImpl->tag_get_data( dualEntity_tag(), verts, 2, vert_tags );
01126     if( MB_SUCCESS != result || mbImpl->type_from_handle( vert_tags[0] ) != MBQUAD ||
01127         mbImpl->type_from_handle( vert_tags[1] ) != MBQUAD )
01128         return false;
01129 
01130     else
01131         return true;
01132 }
01133 
01134 ErrorCode DualTool::construct_hp_parent_child()
01135 {
01136     Range dual_surfs, dual_cells, dual_edges;
01137     ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs );
01138     if( MB_SUCCESS != result || dual_surfs.empty() ) return result;
01139     std::vector< EntityHandle > dual_curve_sets;
01140 
01141     for( Range::iterator surf_it = dual_surfs.begin(); surf_it != dual_surfs.end(); ++surf_it )
01142     {
01143         // get all the cells, edges in those cells, and chords for those edges
01144         dual_cells.clear();
01145         result = mbImpl->get_entities_by_handle( *surf_it, dual_cells );
01146         if( MB_SUCCESS != result ) return result;
01147         dual_edges.clear();
01148         result = mbImpl->get_adjacencies( dual_cells, 1, false, dual_edges, Interface::UNION );
01149         if( MB_SUCCESS != result ) return result;
01150         dual_curve_sets.resize( dual_edges.size() );
01151         result = mbImpl->tag_get_data( dualCurve_tag(), dual_edges, &dual_curve_sets[0] );
01152         if( MB_SUCCESS != result ) return result;
01153 
01154         // reuse dual_cells to get unique list of chord sets
01155         dual_cells.clear();
01156         for( unsigned int i = 0; i < dual_edges.size(); i++ )
01157             if( dual_curve_sets[i] != 0 ) dual_cells.insert( dual_curve_sets[i] );
01158 
01159         // now connect up this dual surf with all the 1d ones
01160         for( Range::iterator rit = dual_cells.begin(); rit != dual_cells.end(); ++rit )
01161         {
01162             result = mbImpl->add_parent_child( *surf_it, *rit );
01163             if( MB_SUCCESS != result ) return result;
01164         }
01165     }
01166 
01167     return MB_SUCCESS;
01168 }
01169 
01170 ErrorCode DualTool::get_graphics_points( EntityHandle dual_ent, std::vector< int >& npts,
01171                                          std::vector< GraphicsPoint >& points )
01172 {
01173     // shouldn't be a set
01174     assert( MBENTITYSET != mbImpl->type_from_handle( dual_ent ) );
01175 
01176     // get the graphics points comprising the given entity
01177     GraphicsPoint gp_array[DualTool::GP_SIZE];
01178 
01179     ErrorCode result = MB_SUCCESS;
01180 
01181     // switch based on topological dimension
01182     switch( mbImpl->dimension_from_handle( dual_ent ) )
01183     {
01184         case 0:
01185             // just return the vertex point
01186             result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array );
01187             if( MB_SUCCESS == result ) points.push_back( gp_array[0] );
01188 
01189             break;
01190 
01191         case 1:
01192             // get my graphics point then those of my vertices
01193             const EntityHandle* connect;
01194             int num_connect;
01195             result = mbImpl->get_connectivity( dual_ent, connect, num_connect );
01196             if( MB_SUCCESS != result ) break;
01197 
01198             result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, gp_array );
01199             if( MB_SUCCESS == result )
01200             {
01201                 points.push_back( gp_array[0] );
01202                 points.push_back( gp_array[0] );
01203                 points.push_back( gp_array[1] );
01204                 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array );
01205                 if( MB_SUCCESS == result ) points[1] = gp_array[0];
01206             }
01207 
01208             npts.push_back( 3 );
01209 
01210             break;
01211 
01212         case 2:
01213             result = get_cell_points( dual_ent, npts, points );
01214             break;
01215     }
01216 
01217     return result;
01218 }
01219 
01220 ErrorCode DualTool::get_cell_points( EntityHandle dual_ent, std::vector< int >& npts,
01221                                      std::vector< GraphicsPoint >& points )
01222 {
01223     assert( MBPOLYGON == mbImpl->type_from_handle( dual_ent ) );
01224 
01225     // get the 1cells in this 2cell
01226     Range one_cells;
01227 
01228     Range tc_range;
01229     tc_range.insert( dual_ent );
01230     ErrorCode result = mbImpl->get_adjacencies( tc_range, 1, false, one_cells, Interface::UNION );RR;
01231 
01232     int num_edges = one_cells.size();
01233     std::vector< GraphicsPoint > dum_gps( num_edges + 1 );
01234 
01235     // get graphics points for 0cells and for this cell
01236     result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), one_cells, &dum_gps[0] );RR;
01237     result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, &( dum_gps[num_edges] ) );RR;
01238 
01239     Range::iterator eit;
01240     const EntityHandle* connect;
01241     int num_connect;
01242     GraphicsPoint vert_gps[2];
01243     int i;
01244     for( i = 0, eit = one_cells.begin(); i < num_edges; i++, ++eit )
01245     {
01246         // get the vertices and the graphics points for them
01247         result = mbImpl->get_connectivity( *eit, connect, num_connect );RR;
01248         result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, vert_gps );RR;
01249 
01250         // make the 2 tris corresponding to this edge; don't worry about order
01251         // for now
01252         npts.push_back( 3 );
01253         points.push_back( dum_gps[num_edges] );
01254         points.push_back( vert_gps[0] );
01255         points.push_back( dum_gps[i] );
01256 
01257         npts.push_back( 3 );
01258         points.push_back( dum_gps[num_edges] );
01259         points.push_back( dum_gps[i] );
01260         points.push_back( vert_gps[1] );
01261     }
01262 
01263     return result;
01264 }
01265 
01266 ErrorCode DualTool::get_graphics_points( const Range& in_range, std::vector< GraphicsPoint >& points,
01267                                          const bool assign_ids, const int start_id )
01268 {
01269     // return graphics points on dual entities in in_range or in entities
01270     // in sets in in_range
01271     ErrorCode result;
01272 
01273     // for each dual hyperplane set:
01274     Range::const_iterator rit;
01275 
01276     Range two_cells, all_cells;
01277     for( rit = in_range.begin(); rit != in_range.end(); ++rit )
01278     {
01279         // for each entity:
01280         two_cells.clear();
01281         EntityType this_type = mbImpl->type_from_handle( *rit );
01282         if( MBENTITYSET == this_type )
01283         {
01284             result = mbImpl->get_entities_by_handle( *rit, two_cells );RR;
01285 
01286             std::copy( two_cells.begin(), two_cells.end(), range_inserter( all_cells ) );
01287         }
01288 
01289         else
01290         {
01291             two_cells.insert( *rit );
01292             assert( this_type == MBVERTEX || this_type == MBEDGE || this_type == MBPOLYGON ||
01293                     this_type == MBPOLYHEDRON );
01294         }
01295 
01296         result = mbImpl->get_adjacencies( two_cells, 0, false, all_cells, Interface::UNION );RR;
01297         result = mbImpl->get_adjacencies( two_cells, 1, false, all_cells, Interface::UNION );RR;
01298     }
01299 
01300     // get graphics points
01301     points.resize( all_cells.size() );
01302 
01303     result = mbImpl->tag_get_data( dualGraphicsPointTag, all_cells, &points[0] );RR;
01304 
01305     if( assign_ids )
01306     {
01307         int i = start_id;
01308 
01309         for( std::vector< GraphicsPoint >::iterator vit = points.begin(); vit != points.end(); ++vit )
01310             vit->id = i++;
01311 
01312         result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), all_cells, &points[0] );RR;
01313     }
01314 
01315     return result;
01316 }
01317 
01318 EntityHandle DualTool::next_loop_vertex( const EntityHandle last_v, const EntityHandle this_v,
01319                                          const EntityHandle dual_surf )
01320 {
01321     // given two vertices, find the next one on the loop; if one is a dual
01322     // surface, then just choose either one for that surface
01323     assert( ( 0 == last_v || mbImpl->type_from_handle( last_v ) == MBVERTEX ) &&
01324             mbImpl->type_from_handle( this_v ) == MBVERTEX && mbImpl->type_from_handle( dual_surf ) == MBENTITYSET );
01325 
01326     // get the connected vertices
01327     MeshTopoUtil tpu( mbImpl );
01328     Range other_verts;
01329     ErrorCode result = tpu.get_bridge_adjacencies( this_v, 1, 0, other_verts );
01330     if( MB_SUCCESS != result || other_verts.empty() ) return 0;
01331 
01332     // if (mbImpl->type_from_handle(last_v) == MBENTITYSET) {
01333     // dual surface, choose either; first get a 2cell on this surface
01334     Range tcells, tcells2, verts;
01335     result = mbImpl->get_entities_by_type( dual_surf, MBPOLYGON, tcells );
01336     if( MB_SUCCESS != result || tcells.empty() ) return 0;
01337 
01338     // ok, pay attention here: first get 2cells common to dual surface and this_v
01339     verts.insert( this_v );
01340     result = mbImpl->get_adjacencies( verts, 2, false, tcells );
01341     if( MB_SUCCESS != result || tcells.empty() ) return 0;
01342 
01343     // next get vertices common to both 2cells and subtract from other_verts; also
01344     // remove last_v if it's non-zero
01345     verts.clear();
01346     result = mbImpl->get_adjacencies( tcells, 0, false, verts );
01347     if( MB_SUCCESS != result || verts.empty() ) return 0;
01348 
01349     Range tmp_verts = subtract( other_verts, verts );
01350     other_verts.swap( tmp_verts );
01351     if( 0 != last_v ) other_verts.erase( last_v );
01352 
01353     // now get intersection of remaining vertices and 2 2cells vertices
01354     // look at each one successively; should find one, maybe not on both
01355     tmp_verts = other_verts;
01356     Range tmp_faces( *tcells.begin(), *tcells.begin() );
01357     result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts );
01358     if( MB_SUCCESS == result && !tmp_verts.empty() ) return *tmp_verts.begin();
01359     tmp_faces.clear();
01360     tmp_faces.insert( *tcells.rbegin() );
01361     result = mbImpl->get_adjacencies( tmp_faces, 0, false, other_verts );
01362     if( MB_SUCCESS == result && !other_verts.empty() ) return *other_verts.begin();
01363 
01364     // if we got here, there isn't any
01365     return 0;
01366 }
01367 
01368 EntityHandle DualTool::get_dual_hyperplane( const EntityHandle ncell )
01369 {
01370     // get the sheet or chord it's in
01371     std::vector< EntityHandle > adj_sets;
01372     ErrorCode result = mbImpl->get_adjacencies( &ncell, 1, 4, false, adj_sets );
01373     if( MB_SUCCESS != result ) return 0;
01374 
01375     EntityHandle dum_set;
01376     for( std::vector< EntityHandle >::iterator vit = adj_sets.begin(); vit != adj_sets.end(); ++vit )
01377     {
01378         if( mbImpl->tag_get_data( dualCurve_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND ||
01379             mbImpl->tag_get_data( dualSurface_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND )
01380             return *vit;
01381     }
01382 
01383     return 0;
01384 }
01385 
01386 //! set the dual surface or curve for an entity
01387 ErrorCode DualTool::set_dual_surface_or_curve( EntityHandle entity, const EntityHandle dual_hyperplane,
01388                                                const int dual_entity_dimension )
01389 {
01390     if( 1 == dual_entity_dimension )
01391         mbImpl->tag_set_data( dualCurve_tag(), &entity, 1, &dual_hyperplane );
01392     else if( 2 == dual_entity_dimension )
01393         mbImpl->tag_set_data( dualSurface_tag(), &entity, 1, &dual_hyperplane );
01394     else
01395         return MB_INDEX_OUT_OF_RANGE;
01396 
01397     return MB_SUCCESS;
01398 }
01399 
01400 //! return the corresponding dual entity
01401 EntityHandle DualTool::get_dual_entity( const EntityHandle this_ent ) const
01402 {
01403     EntityHandle dual_ent;
01404     ErrorCode result = mbImpl->tag_get_data( dualEntity_tag(), &this_ent, 1, &dual_ent );
01405     if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result )
01406         return 0;
01407     else
01408         return dual_ent;
01409 }
01410 
01411 //! return the corresponding dual entity
01412 EntityHandle DualTool::get_extra_dual_entity( const EntityHandle this_ent )
01413 {
01414     EntityHandle dual_ent;
01415     ErrorCode result = mbImpl->tag_get_data( extraDualEntity_tag(), &this_ent, 1, &dual_ent );
01416     if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result )
01417         return 0;
01418     else
01419         return dual_ent;
01420 }
01421 
01422 ErrorCode DualTool::atomic_pillow( EntityHandle odedge, EntityHandle& quad1, EntityHandle& quad2 )
01423 {
01424     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
01425 
01426     if( debug_ap )
01427     {
01428         Range sets;
01429         Tag ms_tag;
01430 
01431         ErrorCode result = mbImpl->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, ms_tag );
01432         if( MB_SUCCESS == result )
01433         {
01434             result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &ms_tag, NULL, 1, sets );
01435             if( MB_SUCCESS == result ) result = mbImpl->delete_entities( sets );
01436         }
01437     }
01438 
01439     std::cout << "-AP(";
01440     print_cell( odedge );
01441     std::cout << ")" << std::endl;
01442 
01443     // perform an atomic pillow operation around dedge
01444 
01445     // grab the quad before deleting the odedge
01446     quad1 = get_dual_entity( odedge );
01447     assert( 0 != quad1 );
01448 
01449     // 0. get star 2cells around odedge (before odedge changes) and 3cells around
01450     // those 2cells (going to delete the 2cells, therefore need to delete the 3cells
01451     // that depend on those too)
01452     MeshTopoUtil mtu( mbImpl );
01453     Range star_cells, tmp_cells;
01454     ErrorCode result = mbImpl->get_adjacencies( &odedge, 1, 2, false, star_cells );RR;
01455     result = mbImpl->get_adjacencies( star_cells, 3, false, tmp_cells, Interface::UNION );RR;
01456     star_cells.merge( tmp_cells );
01457     star_cells.insert( odedge );
01458 
01459     // tear down the dual entities which will be modified by the ap first
01460     result = delete_dual_entities( star_cells );RR;
01461 
01462     // now change the quad to an ap
01463     std::vector< EntityHandle > verts;
01464     result = mbImpl->get_connectivity( &quad1, 1, verts );RR;
01465 
01466     // get average position of vertices
01467     double coords[12], avg[3] = { 0.0, 0.0, 0.0 };
01468     result = mbImpl->get_coords( &verts[0], verts.size(), coords );RR;
01469     for( int i = 0; i < 4; i++ )
01470     {
01471         avg[0] += coords[3 * i];
01472         avg[1] += coords[3 * i + 1];
01473         avg[2] += coords[3 * i + 2];
01474     }
01475     for( int i = 0; i < 3; i++ )
01476         avg[i] *= 0.25;
01477 
01478     // for each position, get a corresponding position 1/2 way to avg
01479     double new_coords[12];
01480     for( int i = 0; i < 4; i++ )
01481     {
01482         new_coords[3 * i]     = avg[0] + .5 * ( coords[3 * i] - avg[0] );
01483         new_coords[3 * i + 1] = avg[1] + .5 * ( coords[3 * i + 1] - avg[1] );
01484         new_coords[3 * i + 2] = avg[2] + .5 * ( coords[3 * i + 2] - avg[2] );
01485     }
01486 
01487     // make the 4 new vertices; store in vector long enough for hex connectivity
01488     for( int i = 0; i < 4; i++ )
01489     {
01490         verts.push_back( 0 );
01491         result = mbImpl->create_vertex( &new_coords[3 * i], verts[4 + i] );RR;
01492     }
01493 
01494     // get the hexes connected to the quad
01495     Range hexes;
01496     result = mbImpl->get_adjacencies( &quad1, 1, 3, false, hexes );RR;
01497     assert( hexes.size() <= 2 );
01498 
01499     // remove any explicit adjacency from the first hex, since that'll get connected
01500     // to the new quad; add adjacency between quad and second hex, if there is a 2nd
01501     result = mbImpl->remove_adjacencies( quad1, &( *hexes.begin() ), 1 );RR;
01502     if( hexes.size() == 2 )
01503     {
01504         result = mbImpl->add_adjacencies( quad1, &( *hexes.rbegin() ), 1, false );RR;
01505     }
01506 
01507     // create the new, inner quad, and make it explicitly adjacent to 1st hex;
01508     // make the connectivity of this quad same as the original one
01509     std::vector< EntityHandle > tmp_verts;
01510     std::copy( verts.begin(), verts.end(), std::back_inserter( tmp_verts ) );
01511 
01512     result = mbImpl->create_element( MBQUAD, &tmp_verts[0], 4, quad2 );RR;
01513     result = mbImpl->add_adjacencies( quad2, &( *hexes.begin() ), 1, false );RR;
01514 
01515     // reverse the connectivity of the 1st hex
01516     std::reverse( verts.begin(), verts.begin() + 4 );
01517     std::reverse( verts.begin() + 4, verts.end() );
01518 
01519     // now make two inner hexes; note connectivity array is flipped for the two hexes
01520     EntityHandle new_hexes[2];
01521     result = mbImpl->create_element( MBHEX, &verts[0], 8, new_hexes[0] );RR;
01522     result = mbImpl->create_element( MBHEX, &tmp_verts[0], 8, new_hexes[1] );RR;
01523 
01524     // set the global id tag on the new hexes
01525     int new_hex_ids[2] = { maxHexId + 1, maxHexId + 2 };
01526     maxHexId += 2;
01527     result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 2, new_hex_ids );
01528     if( MB_SUCCESS != result ) return result;
01529 
01530     // by definition, quad1 is adj to new_hexes[0]
01531     result = mbImpl->add_adjacencies( quad1, &new_hexes[0], 1, false );RR;
01532     result = mbImpl->add_adjacencies( quad2, &new_hexes[1], 1, false );RR;
01533 
01534     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
01535 
01536     // now update the dual
01537     result = construct_hex_dual( &new_hexes[0], 2 );RR;
01538 
01539     // get the new dual surface, by getting one of the edges between the center
01540     // and outer vertex rings
01541     Range new_edge;
01542     verts[1] = verts[4];
01543     result   = mbImpl->get_adjacencies( &verts[0], 2, 1, false, new_edge );
01544     if( MB_SUCCESS != result || new_edge.size() != 1 ) return result;
01545 
01546     return MB_SUCCESS;
01547 }
01548 
01549 //! effect reverse atomic pillow operation
01550 ErrorCode DualTool::rev_atomic_pillow( EntityHandle pillow, Range& chords )
01551 {
01552     // get the dual entities associated with elements in the pillow; go through
01553     // the elements instead of the pillow sheet so you get all of them, not just
01554     // the ones on the sheet
01555     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
01556 
01557     std::cout << "-AP(";
01558     print_cell( pillow );
01559     std::cout << ")" << std::endl;
01560 
01561     Range dverts;
01562     ErrorCode result = get_dual_entities( pillow, NULL, NULL, &dverts, NULL, NULL );
01563     if( MB_SUCCESS != result ) return result;
01564     assert( 2 == dverts.size() );
01565 
01566     EntityHandle hexes[2];
01567     result = mbImpl->tag_get_data( dualEntity_tag(), dverts, hexes );RR;
01568     assert( hexes[0] != 0 && hexes[1] != 0 );
01569 
01570     std::vector< EntityHandle > dcells[4];
01571     Range pcells[4];
01572     std::copy( hexes, hexes + 2, range_inserter( pcells[3] ) );
01573     std::copy( dverts.begin(), dverts.end(), std::back_inserter( dcells[0] ) );
01574     for( int dim = 0; dim <= 2; dim++ )
01575     {
01576         result = mbImpl->get_adjacencies( hexes, 2, dim, false, pcells[dim], Interface::UNION );RR;
01577         dcells[3 - dim].resize( pcells[dim].size() );
01578         result = mbImpl->tag_get_data( dualEntity_tag(), pcells[dim], &dcells[3 - dim][0] );RR;
01579     }
01580 
01581     // delete the dual entities which are part of the original pillow
01582     result = mbImpl->delete_entities( &pillow, 1 );
01583     if( MB_SUCCESS != result ) return result;
01584 
01585     result = mbImpl->delete_entities( chords );
01586     if( MB_SUCCESS != result ) return result;
01587 
01588     for( int i = 3; i >= 0; i-- )
01589     {
01590         result = delete_dual_entities( &dcells[i][0], dcells[i].size() );RR;
01591     }
01592 
01593     // delete the primal entities inserted by the ap; be careful to get the right
01594     // faces, edges and vertices
01595     Range del_faces, del_edges, del_verts, tmp_faces, tmp_verts;
01596     // faces are the shared 5 and the 1 other one with greater handle (which
01597     // we know will be later in the range)
01598     result = mbImpl->get_adjacencies( hexes, 2, 2, false, del_faces );RR;
01599     assert( 5 == del_faces.size() );
01600     std::copy( pcells[2].begin(), pcells[2].end(), range_inserter( tmp_faces ) );
01601     tmp_faces = subtract( tmp_faces, del_faces );
01602     del_faces.insert( *tmp_faces.rbegin() );
01603     result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts );RR;
01604     std::copy( pcells[0].begin(), pcells[0].end(), range_inserter( del_verts ) );
01605     del_verts = subtract( del_verts, tmp_verts );
01606     assert( 4 == del_verts.size() );
01607     result = mbImpl->get_adjacencies( del_verts, 1, false, del_edges, Interface::UNION );RR;
01608     assert( 8 == del_edges.size() );
01609 
01610     result = mbImpl->delete_entities( hexes, 2 );RR;
01611     result = mbImpl->delete_entities( del_faces );RR;
01612     result = mbImpl->delete_entities( del_edges );RR;
01613     result = mbImpl->delete_entities( del_verts );RR;
01614 
01615     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
01616 
01617     // recompute the dual for the hexes on either side of the quad affected
01618     // by the ap removal
01619     Range tmp_hexes;
01620     result = mbImpl->get_adjacencies( tmp_verts, 3, false, tmp_hexes, Interface::UNION );RR;
01621     result = construct_hex_dual( tmp_hexes );RR;
01622 
01623     return MB_SUCCESS;
01624 }
01625 
01626 ErrorCode DualTool::delete_dual_entities( EntityHandle* entities, int num_entities )
01627 {
01628     Range tmp_ents;
01629     std::copy( entities, entities + num_entities, range_inserter( tmp_ents ) );
01630     return delete_dual_entities( tmp_ents );
01631 }
01632 
01633 ErrorCode DualTool::delete_dual_entities( Range& entities )
01634 {
01635     if( entities.empty() ) return delete_whole_dual();
01636 
01637     EntityHandle null_entity = 0;
01638     ErrorCode result;
01639     Range ents_to_delete;
01640 
01641     while( !entities.empty() )
01642     {
01643         EntityHandle this_entity = entities.pop_back();
01644 
01645         // reset the primal's dual entity
01646         EntityHandle primal = get_dual_entity( this_entity );
01647         if( get_dual_entity( primal ) == this_entity )
01648         {
01649             result = mbImpl->tag_set_data( dualEntity_tag(), &primal, 1, &null_entity );RR;
01650         }
01651         EntityHandle extra = get_extra_dual_entity( primal );
01652         if( 0 != extra )
01653         {
01654             result = mbImpl->tag_set_data( extraDualEntity_tag(), &primal, 1, &null_entity );RR;
01655         }
01656 
01657         ents_to_delete.insert( this_entity );
01658 
01659         // check for extra dual entities
01660         if( mbImpl->type_from_handle( this_entity ) == MBPOLYGON )
01661         {
01662             // for 2cell, might be a loop edge
01663             Range loop_edges;
01664             result = mbImpl->get_adjacencies( &this_entity, 1, 1, false, loop_edges );
01665             for( Range::iterator rit = loop_edges.begin(); rit != loop_edges.end(); ++rit )
01666                 if( check_1d_loop_edge( *rit ) ) entities.insert( *rit );
01667         }
01668         else if( extra && extra != this_entity )
01669             // just put it on the list; primal for which we're extra has already been
01670             // reset to not point to extra entity
01671             ents_to_delete.insert( extra );
01672     }
01673 
01674     // now delete the entities (sheets and chords will be updated automatically)
01675     return mbImpl->delete_entities( ents_to_delete );
01676 }
01677 
01678 void DualTool::print_cell( EntityHandle cell )
01679 {
01680     const EntityHandle* connect;
01681     int num_connect;
01682     ErrorCode result = mbImpl->get_connectivity( cell, connect, num_connect );
01683     if( MB_SUCCESS != result ) return;
01684     bool first = true;
01685     EntityHandle primals[20];
01686     std::vector< int > ids;
01687 
01688     assert( num_connect < 20 );
01689     result = mbImpl->tag_get_data( dualEntityTag, connect, num_connect, primals );
01690     if( MB_SUCCESS != result ) return;
01691     ids.resize( num_connect );
01692     result = mbImpl->tag_get_data( globalIdTag, primals, num_connect, &ids[0] );
01693     if( MB_SUCCESS != result ) return;
01694     for( int i = 0; i < num_connect; i++ )
01695     {
01696         if( !first ) std::cout << "-";
01697         EntityType this_type = mbImpl->type_from_handle( primals[i] );
01698         if( this_type == MBHEX )
01699             std::cout << "h";
01700         else if( this_type == MBQUAD )
01701             std::cout << "f";
01702         else
01703             std::cout << "u";
01704 
01705         if( ids[i] != 0 )
01706             std::cout << ids[i];
01707         else
01708             std::cout << mbImpl->id_from_handle( primals[i] );
01709 
01710         first = false;
01711     }
01712 }
01713 
01714 ErrorCode DualTool::face_open_collapse( EntityHandle ocl, EntityHandle ocr )
01715 {
01716     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
01717 
01718     MeshTopoUtil mtu( mbImpl );
01719 
01720     std::cout << "OC(";
01721     print_cell( ocl );
01722     std::cout << ")-(";
01723     print_cell( ocr );
01724     std::cout << ")" << std::endl;
01725 
01726     // get the primal entities we're dealing with
01727     EntityHandle split_quads[2] = { 0 }, split_edges[3] = { 0 }, split_nodes[2] = { 0 }, other_edges[6] = { 0 },
01728                  other_nodes[6] = { 0 };
01729     Range hexes;
01730     ErrorCode result = foc_get_ents( ocl, ocr, split_quads, split_edges, split_nodes, hexes, other_edges, other_nodes );RR;
01731 
01732     // get star entities around edges, separated into halves
01733     std::vector< EntityHandle > star_dp1[2], star_dp2[2];
01734     result = foc_get_stars( split_quads, split_edges, star_dp1, star_dp2 );RR;
01735 
01736     if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) )
01737         return MB_TYPE_OUT_OF_RANGE;
01738 
01739     result = foc_delete_dual( split_quads, split_edges, hexes );
01740     if( MB_SUCCESS != result ) return result;
01741 
01742     EntityHandle new_quads[2], new_edges[3], new_nodes[2];
01743     result = split_pair_nonmanifold( split_quads, split_edges, split_nodes, star_dp1, star_dp2, other_edges,
01744                                      other_nodes, new_quads, new_edges, new_nodes );
01745     if( MB_SUCCESS != result ) return result;
01746 
01747     // now merge entities, the C of foc
01748     EntityHandle keepit, deleteit;
01749 #define MIN( a, b ) ( a < b ? a : b )
01750 #define MAX( a, b ) ( a > b ? a : b )
01751 #define KEEP_DELETE( a, b, c, d ) \
01752     {                             \
01753         c = MIN( a, b );          \
01754         d = MAX( a, b );          \
01755     }
01756 
01757     // find how many shared edges there were
01758     int num_shared_edges = ( split_edges[2] ? 3 : ( split_edges[1] ? 2 : 1 ) );
01759 
01760     // first the node(s)
01761     for( int i = 0; i < 3 - num_shared_edges; i++ )
01762     {
01763         KEEP_DELETE( other_nodes[2 + 2 * i], other_nodes[3 + 2 * i], keepit, deleteit );
01764         result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
01765     }
01766 
01767     // now the edges
01768     for( int i = 0; i < 4 - num_shared_edges; i++ )
01769     {
01770         KEEP_DELETE( other_edges[2 * i], other_edges[2 * i + 1], keepit, deleteit );
01771         result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
01772     }
01773 
01774     // now the faces
01775     KEEP_DELETE( split_quads[0], split_quads[1], keepit, deleteit );
01776     result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
01777 
01778     result = mbImpl->merge_entities( new_quads[0], new_quads[1], false, true );RR;
01779 
01780     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
01781 
01782     // reconstruct dual
01783     result = construct_hex_dual( hexes );
01784     if( MB_SUCCESS != result ) return result;
01785 
01786     return check_dual_adjs();
01787 }
01788 
01789 ErrorCode DualTool::foc_get_ents( EntityHandle ocl, EntityHandle ocr, EntityHandle* split_quads,
01790                                   EntityHandle* split_edges, EntityHandle* split_nodes, Range& hexes,
01791                                   EntityHandle* other_edges, EntityHandle* other_nodes )
01792 {
01793     // get the entities used for foc; ocl and ocr are dual 1-cells
01794     // representing quads to be split; returned from this function:
01795     // quads[2] - 2 quads to be split
01796     // split_edges[2] - edge(s) to be split (2nd is 0 if only one)
01797     // split_node - node to be split, if any (otherwise 0)
01798     // hexes - connected hexes to split_edges
01799     // other_edges[0], [1] - edges in quads[0] and [1] sharing node with
01800     //        one end of split_edges[0]
01801     // other_edges[2], [3] - other end of split_edges[0] (or [1] if 2
01802     //        split_edges)
01803     // other_edges[4], [5] - edges in quads[0], [1] opposite to split_edges[0]
01804     // other_nodes[0], [1] - nodes on other_edges[0], [1] not shared with
01805     //        split_edges[0]
01806     // other_nodes[2], [3] - nodes on other_edges[2], [3] not shared with
01807     //        split_edges[0] (if 2 split_edges, there's only 1 opposite node
01808     //        in each split quad)
01809     // (for diagram, see Tim's notes from 11/12/07)
01810 
01811     split_quads[0] = get_dual_entity( ocl );
01812     split_quads[1] = get_dual_entity( ocr );
01813     if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) )
01814         return MB_TYPE_OUT_OF_RANGE;
01815 
01816     Range common_edges;
01817     ErrorCode result = mbImpl->get_adjacencies( split_quads, 2, 1, false, common_edges );
01818     if( MB_SUCCESS != result ) return result;
01819 
01820     if( common_edges.empty() ) return MB_FAILURE;
01821     for( unsigned int i = 0; i < common_edges.size(); i++ )
01822         split_edges[i] = common_edges[i];
01823 
01824     MeshTopoUtil mtu( mbImpl );
01825 
01826     if( common_edges.size() == 3 )
01827     {
01828         // find other (non-shared) edges
01829         for( int i = 0; i < 2; i++ )
01830         {
01831             Range tmp_edges;
01832             result = mbImpl->get_adjacencies( &split_quads[i], 1, 1, false, tmp_edges );
01833             if( MB_SUCCESS != result ) return result;
01834             tmp_edges = subtract( tmp_edges, common_edges );
01835             assert( tmp_edges.size() == 1 );
01836             other_edges[i] = *tmp_edges.begin();
01837         }
01838         assert( other_edges[0] && other_edges[1] && other_edges[0] != other_edges[1] );
01839 
01840         // arrange common edges so middle is in middle
01841         result = mtu.opposite_entity( split_quads[0], other_edges[0], split_edges[1] );RR;
01842         common_edges.erase( split_edges[1] );
01843         split_edges[0] = *common_edges.begin();
01844         split_edges[2] = *common_edges.rbegin();
01845         common_edges.insert( split_edges[1] );
01846 
01847         // get split nodes and other nodes
01848         split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 );
01849         split_nodes[1] = mtu.common_entity( split_edges[2], split_edges[1], 0 );
01850         other_nodes[0] = mtu.common_entity( split_edges[0], other_edges[0], 0 );
01851         other_nodes[1] = mtu.common_entity( split_edges[2], other_edges[1], 0 );
01852 
01853         assert( other_nodes[0] && other_nodes[1] && split_nodes[0] && split_nodes[1] );
01854         assert( split_edges[0] && split_edges[1] && split_edges[2] && split_edges[0] != split_edges[1] &&
01855                 split_edges[1] != split_edges[2] && split_edges[0] != split_edges[2] );
01856     }
01857     else if( common_edges.size() == 2 )
01858     {
01859         // split node is shared by split edges
01860         split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 );
01861         if( 0 == split_nodes[0] ) return MB_FAILURE;
01862         // first two other nodes are on split edges opposite split node
01863         result = mtu.opposite_entity( split_edges[0], split_nodes[0], other_nodes[0] );RR;
01864         result = mtu.opposite_entity( split_edges[1], split_nodes[0], other_nodes[1] );RR;
01865         // over split quads:
01866         for( int i = 0; i < 2; i++ )
01867         {
01868             // 1st other edge is opposite second split edge
01869             result = mtu.opposite_entity( split_quads[i], split_edges[1], other_edges[i] );RR;
01870             // 2nd other edge is opposite first split edge
01871             result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR;
01872             // last other node is opposite split node on split quad
01873             result = mtu.opposite_entity( split_quads[i], split_nodes[0], other_nodes[2 + i] );RR;
01874         }
01875     }
01876     else
01877     {
01878         const EntityHandle* connect;
01879         int num_connect;
01880         result = mbImpl->get_connectivity( split_edges[0], connect, num_connect );
01881         if( MB_SUCCESS != result ) return result;
01882         // other_nodes[0], [1] are on split edge
01883         other_nodes[0] = connect[0];
01884         other_nodes[1] = connect[1];
01885 
01886         // for each of the split quads
01887         for( int i = 0; i < 2; i++ )
01888         {
01889             // get the other edge on the split quad adj to node 0 on the split edge, by getting
01890             // edges adj to split quad and node and removing split edge; that's other_edge[i]
01891             Range tmp_range1, tmp_range2;
01892             tmp_range1.insert( connect[0] );
01893             tmp_range1.insert( split_quads[i] );
01894             result = mbImpl->get_adjacencies( tmp_range1, 1, false, tmp_range2 );
01895             if( MB_SUCCESS != result ) return result;
01896             tmp_range2.erase( split_edges[0] );
01897             assert( tmp_range2.size() == 1 );
01898             other_edges[i] = *tmp_range2.begin();
01899             // get edge connected to other node on split edge & split quad; that's
01900             // opposite prev other_edges on the split quad; that's other_edges[4+i]
01901             result = mtu.opposite_entity( split_quads[i], other_edges[i], other_edges[4 + i] );RR;
01902             // get the edge on the split quad opposite the split edge; that's other_edges[2+i]
01903             result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR;
01904             // get nodes on other side of split quad from split edge, by getting common
01905             // node between top/bottom edge and opposite edge
01906             other_nodes[2 + i] = mtu.common_entity( other_edges[i], other_edges[2 + i], 0 );
01907             other_nodes[4 + i] = mtu.common_entity( other_edges[4 + i], other_edges[2 + i], 0 );
01908             if( 0 == other_nodes[2 + i] || 0 == other_nodes[4 + i] ) return MB_FAILURE;
01909         }
01910     }
01911 
01912     result = mbImpl->get_adjacencies( split_edges, common_edges.size(), 3, false, hexes, Interface::UNION );
01913     if( MB_SUCCESS != result ) return result;
01914 
01915     assert( "split node not in other_nodes" && other_nodes[0] != split_nodes[0] && other_nodes[0] != split_nodes[1] &&
01916             other_nodes[1] != split_nodes[0] && other_nodes[1] != split_nodes[1] );
01917     assert( "each split node on an end of a split edge" && mtu.common_entity( other_nodes[0], split_edges[0], 0 ) &&
01918             ( ( ( split_edges[2] && mtu.common_entity( other_nodes[1], split_edges[2], 0 ) ) ||
01919                 ( split_edges[1] && mtu.common_entity( other_nodes[1], split_edges[1], 0 ) ) ||
01920                 mtu.common_entity( other_nodes[1], split_edges[0], 0 ) ) ) );
01921     assert( "opposite other edges meet at an other node" &&
01922             ( mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[0] ||
01923               ( split_edges[2] && mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[1] ) ) &&
01924             ( split_edges[2] ||
01925               ( split_edges[1] && mtu.common_entity( other_edges[2], other_edges[3], 0 ) == other_nodes[1] ) ||
01926               mtu.common_entity( other_edges[4], other_edges[5], 0 ) == other_nodes[1] ) );
01927 
01928     return MB_SUCCESS;
01929 }
01930 
01931 ErrorCode DualTool::split_pair_nonmanifold( EntityHandle* split_quads, EntityHandle* split_edges,
01932                                             EntityHandle* split_nodes, std::vector< EntityHandle >* star_dp1,
01933                                             std::vector< EntityHandle >* star_dp2, EntityHandle* /*other_edges*/,
01934                                             EntityHandle* /*other_nodes*/, EntityHandle* new_quads,
01935                                             EntityHandle* new_edges, EntityHandle* new_nodes )
01936 {
01937 
01938     // if there's a bdy in the star around the shared edge(s), get the quads on that
01939     // bdy so we know which edges to merge after the split-nonmanifold
01940     MeshTopoUtil mtu( mbImpl );
01941     ErrorCode result;
01942 
01943     // get which star the split faces are in, and choose the other one
01944     int new_side = -1;
01945     if( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() )
01946         new_side = 1;
01947     else if( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() )
01948         new_side = 0;
01949     assert( -1 != new_side );
01950     if( -1 == new_side ) return MB_FAILURE;
01951 
01952     //=============== split faces
01953 
01954     for( int i = 0; i < 2; i++ )
01955     {
01956         // get a hex in star_dp2[new_side] that's adj to this split quad, to tell
01957         // mtu which one the new quad should go with; there should be at least one,
01958         // if we have any hexes connected to the split quad
01959         EntityHandle gowith_hex = 0;
01960         for( std::vector< EntityHandle >::iterator vit = star_dp2[new_side].begin(); vit != star_dp2[new_side].end();
01961              ++vit )
01962         {
01963             if( mtu.common_entity( *vit, split_quads[i], 2 ) )
01964             {
01965                 gowith_hex = *vit;
01966                 break;
01967             }
01968         }
01969         assert( 0 != gowith_hex );
01970 
01971         // split manifold each of the split_quads, and put the results on the merge list
01972         result =
01973             mtu.split_entities_manifold( split_quads + i, 1, new_quads + i, NULL, ( gowith_hex ? &gowith_hex : NULL ) );RR;
01974     }
01975 
01976     // make ranges of faces which need to be explicitly adj to old, new
01977     // edge; faces come from stars and new_quads (which weren't in the stars);
01978     // new_quads go with side j, which does not have split quads
01979     Range tmp_addl_faces[2], addl_faces[2];
01980     for( int i = 0; i < 2; i++ )
01981     {
01982         std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_faces[i] ) );
01983         tmp_addl_faces[new_side].insert( new_quads[i] );
01984     }
01985 #ifndef NDEBUG
01986     bool cond1 = ( "split_quads on 1, new_quads on 0" &&
01987                    tmp_addl_faces[0].find( split_quads[0] ) == tmp_addl_faces[0].end() &&
01988                    tmp_addl_faces[0].find( split_quads[1] ) == tmp_addl_faces[0].end() &&
01989                    tmp_addl_faces[1].find( split_quads[0] ) != tmp_addl_faces[1].end() &&
01990                    tmp_addl_faces[1].find( split_quads[1] ) != tmp_addl_faces[1].end() &&
01991                    tmp_addl_faces[0].find( new_quads[0] ) != tmp_addl_faces[0].end() &&
01992                    tmp_addl_faces[0].find( new_quads[1] ) != tmp_addl_faces[0].end() &&
01993                    tmp_addl_faces[1].find( new_quads[0] ) == tmp_addl_faces[1].end() &&
01994                    tmp_addl_faces[1].find( new_quads[1] ) == tmp_addl_faces[1].end() ),
01995          cond2 = ( "split_quads on 0, new_quads on 1" &&
01996                    tmp_addl_faces[0].find( split_quads[0] ) != tmp_addl_faces[0].end() &&
01997                    tmp_addl_faces[0].find( split_quads[1] ) != tmp_addl_faces[0].end() &&
01998                    tmp_addl_faces[1].find( split_quads[0] ) == tmp_addl_faces[1].end() &&
01999                    tmp_addl_faces[1].find( split_quads[1] ) == tmp_addl_faces[1].end() &&
02000                    tmp_addl_faces[0].find( new_quads[0] ) == tmp_addl_faces[0].end() &&
02001                    tmp_addl_faces[0].find( new_quads[1] ) == tmp_addl_faces[0].end() &&
02002                    tmp_addl_faces[1].find( new_quads[0] ) != tmp_addl_faces[1].end() &&
02003                    tmp_addl_faces[1].find( new_quads[1] ) != tmp_addl_faces[1].end() );
02004 
02005     assert( cond1 || cond2 );
02006 #endif
02007 
02008     //=============== split edge(s)
02009     for( int j = 0; j < 3; j++ )
02010     {
02011         if( !split_edges[j] ) break;
02012 
02013         // filter add'l faces to only those adj to split_edges[j]
02014         addl_faces[0] = tmp_addl_faces[0];
02015         addl_faces[1] = tmp_addl_faces[1];
02016         for( int i = 0; i < 2; i++ )
02017         {
02018             result = mbImpl->get_adjacencies( &split_edges[j], 1, 2, false, addl_faces[i] );RR;
02019         }
02020 
02021         // split edge
02022         result = mtu.split_entity_nonmanifold( split_edges[j], addl_faces[1 - new_side], addl_faces[new_side],
02023                                                new_edges[j] );RR;
02024     }
02025 
02026     //=============== split node(s)
02027 
02028     for( int j = 0; j < 2; j++ )
02029     {
02030         if( !split_nodes[j] ) break;
02031 
02032         // if we're splitting multiple edges, there might be other edges that have the split
02033         // node; also need to know which side they're on
02034         Range tmp_addl_edges[2];
02035         result = foc_get_addl_ents( star_dp1, star_dp2, split_edges, split_nodes[j], tmp_addl_edges );RR;
02036 
02037         // also, we need to know which of the split/new edges go
02038         // with the split/new node; new edges go with side 0, split with 1
02039         for( int i = 0; i < 3; i++ )
02040         {
02041             if( !split_edges[i] ) break;
02042             tmp_addl_edges[new_side].insert( new_edges[i] );
02043             tmp_addl_edges[1 - new_side].insert( split_edges[i] );
02044         }
02045 
02046         // same for star faces and hexes
02047         for( int i = 0; i < 2; i++ )
02048         {
02049             std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_edges[i] ) );
02050             std::copy( star_dp2[i].begin(), star_dp2[i].end(), range_inserter( tmp_addl_edges[i] ) );
02051         }
02052 
02053         // finally, new quads
02054         for( int i = 0; i < 2; i++ )
02055             tmp_addl_edges[new_side].insert( new_quads[i] );
02056 
02057         // filter the entities, keeping only the ones adjacent to this node
02058         Range addl_edges[2];
02059         for( int i = 0; i < 2; i++ )
02060         {
02061             for( Range::reverse_iterator rit = tmp_addl_edges[i].rbegin(); rit != tmp_addl_edges[i].rend(); ++rit )
02062             {
02063                 if( mtu.common_entity( *rit, split_nodes[j], 0 ) ) addl_edges[i].insert( *rit );
02064             }
02065         }
02066 
02067         // now split the node too
02068         result = mtu.split_entity_nonmanifold( split_nodes[j], addl_edges[1 - new_side], addl_edges[new_side],
02069                                                new_nodes[j] );RR;
02070     }
02071 
02072     return MB_SUCCESS;
02073 }
02074 
02075 ErrorCode DualTool::foc_get_addl_ents( std::vector< EntityHandle >* star_dp1, std::vector< EntityHandle >* /*star_dp2*/,
02076                                        EntityHandle* split_edges, EntityHandle split_node, Range* addl_ents )
02077 {
02078     // if we're splitting 2 edges, there might be other edges that have the split
02079     // node; also need to know which side they're on
02080 
02081     // algorithm: for a given star_dp1 (faces) on a side:
02082     // - get all edges adj to all faces -> R1
02083     // - get all edges adj to split_node -> R2
02084     // - R3 = R1 & R2 (edges connected to split_node & adj to a star face)
02085     // - R3 -= split_edges (take split edges off addl_ents)
02086 
02087     Range R2;
02088     MeshTopoUtil mtu( mbImpl );
02089     ErrorCode result = mbImpl->get_adjacencies( &split_node, 1, 1, false, R2 );RR;
02090     Range::iterator rit;
02091 
02092     for( int i = 0; i < 2; i++ )
02093     {
02094         Range R1, R3;
02095         result = mbImpl->get_adjacencies( &star_dp1[i][0], star_dp1[i].size(), 1, false, R1, Interface::UNION );RR;
02096         R3 = intersect( R1, R2 );
02097         for( int j = 0; j < 3; j++ )
02098             if( split_edges[j] ) R3.erase( split_edges[j] );
02099         addl_ents[i].merge( R3 );
02100     }
02101 
02102     return MB_SUCCESS;
02103 }
02104 
02105 ErrorCode DualTool::foc_get_stars( EntityHandle* split_quads, EntityHandle* split_edges,
02106                                    std::vector< EntityHandle >* star_dp1, std::vector< EntityHandle >* star_dp2 )
02107 {
02108     bool on_bdy = false, on_bdy_tmp;
02109     ErrorCode result;
02110     MeshTopoUtil mtu( mbImpl );
02111 
02112     // get the star around the split_edge
02113     std::vector< EntityHandle > qstar, hstar;
02114     unsigned int qpos = 0;
02115 
02116     for( int i = 0; i < 3; i++ )
02117     {
02118         if( !split_edges[i] ) break;
02119 
02120         // get the star around this split edge
02121         unsigned int qpos_tmp = 0;
02122         std::vector< EntityHandle > qstar_tmp, hstar_tmp;
02123         result = mtu.star_entities( split_edges[i], qstar_tmp, on_bdy_tmp, 0, &hstar_tmp );RR;
02124         // if we're on the bdy, add a null to the hex star too
02125         if( on_bdy_tmp )
02126         {
02127             assert( hstar_tmp.size() == qstar_tmp.size() - 1 );
02128             hstar_tmp.push_back( 0 );
02129             on_bdy = true;
02130         }
02131 
02132         // get the position of first split quad in star
02133         while( qpos_tmp < qstar_tmp.size() && qstar_tmp[qpos_tmp] != split_quads[0] )
02134             qpos_tmp++;
02135         if( qpos_tmp == qstar_tmp.size() ) return MB_FAILURE;
02136 
02137         bool forward;
02138         // 1st iteration is forward by definition
02139         if( 0 == i ) forward = true;
02140 
02141         // need to be careful about direction on later iters
02142         else if( hstar[qpos] == hstar_tmp[qpos_tmp] )
02143             forward = true;
02144         else if( hstar[qpos] == hstar_tmp[( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size()] &&
02145                  hstar_tmp[qpos_tmp] == hstar[( qpos + qstar.size() - 1 ) % qstar.size()] )
02146             forward = false;
02147         else
02148             return MB_FAILURE;
02149 
02150         if( forward )
02151         {
02152             // 1st half of star
02153             // save hex right after split_quad[0] first
02154             star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
02155             qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
02156             while( qstar_tmp[qpos_tmp] != split_quads[1] )
02157             {
02158                 star_dp1[0].push_back( qstar_tmp[qpos_tmp] );
02159                 star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
02160                 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
02161             }
02162             // 2nd half of star
02163             // save hex right after split_quad[1] first
02164             star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
02165             qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
02166             while( qstar_tmp[qpos_tmp] != split_quads[0] )
02167             {
02168                 star_dp1[1].push_back( qstar_tmp[qpos_tmp] );
02169                 star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
02170                 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
02171             }
02172         }
02173         else
02174         {
02175             // go in reverse - take prev hex instead of current
02176             // one, and step in reverse
02177 
02178             // save hex right after split_quad[0] first
02179             qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
02180             star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
02181             while( qstar_tmp[qpos_tmp] != split_quads[1] )
02182             {
02183                 star_dp1[0].push_back( qstar_tmp[qpos_tmp] );
02184                 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
02185                 star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
02186             }
02187             // 2nd half of star
02188             // save hex right after split_quad[1] first
02189             qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
02190             star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
02191             while( qstar_tmp[qpos_tmp] != split_quads[0] )
02192             {
02193                 star_dp1[1].push_back( qstar_tmp[qpos_tmp] );
02194                 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
02195                 star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
02196             }
02197         }
02198 
02199         if( 0 == i )
02200         {
02201             // if we're on the first iteration, save results and continue, other iters
02202             // get compared to this one
02203             qstar.swap( qstar_tmp );
02204             hstar.swap( hstar_tmp );
02205             on_bdy = on_bdy_tmp;
02206             qpos   = qpos_tmp;
02207         }
02208     }
02209 
02210     // split quads go on list with NULLs, if any, otherwise on 2nd
02211     if( on_bdy )
02212     {
02213         if( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) != star_dp2[0].end() )
02214         {
02215             // remove *all* the zeros
02216             star_dp2[0].erase( std::remove( star_dp2[0].begin(), star_dp2[0].end(), 0 ), star_dp2[0].end() );
02217             // put the split quads on this half
02218             star_dp1[0].push_back( split_quads[0] );
02219             star_dp1[0].push_back( split_quads[1] );
02220         }
02221         else
02222         {
02223             star_dp2[1].erase( std::remove( star_dp2[1].begin(), star_dp2[1].end(), 0 ), star_dp2[1].end() );
02224             // put the split quads on this half
02225             star_dp1[1].push_back( split_quads[0] );
02226             star_dp1[1].push_back( split_quads[1] );
02227         }
02228     }
02229     else
02230     {
02231         star_dp1[1].push_back( split_quads[0] );
02232         star_dp1[1].push_back( split_quads[1] );
02233     }
02234 
02235     // some error checking
02236     if( !( ( ( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) == star_dp1[0].end() &&
02237                std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) == star_dp1[0].end() &&
02238                std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() &&
02239                std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) != star_dp1[1].end() ) ||
02240 
02241              ( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) == star_dp1[1].end() &&
02242                std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) == star_dp1[1].end() &&
02243                std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() &&
02244                std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) != star_dp1[0].end() ) ) ) )
02245     {
02246         std::cerr << "foc_get_stars: both split quads should be on the same star list half and not "
02247                   << "on the other, failed" << std::endl;
02248         return MB_FAILURE;
02249     }
02250 
02251     if( !( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) == star_dp2[0].end() &&
02252            std::find( star_dp2[1].begin(), star_dp2[1].end(), 0 ) == star_dp2[1].end() ) )
02253     {
02254         std::cerr << "foc_get_stars: no NULLs on the hstar lists, failed";
02255         return MB_FAILURE;
02256     }
02257 
02258     return MB_SUCCESS;
02259 }
02260 
02261 ErrorCode DualTool::foc_delete_dual( EntityHandle* split_quads, EntityHandle* split_edges, Range& hexes )
02262 {
02263     // special delete dual procedure, because in some cases we need to delete
02264     // a sheet too since it'll get merged into another
02265 
02266     // figure out whether we'll need to delete a sheet
02267     EntityHandle sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[0] ) );
02268     if( split_edges[1] ) sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[1] ) );
02269     EntityHandle chordl = get_dual_hyperplane( get_dual_entity( split_quads[0] ) );
02270     EntityHandle chordr = get_dual_hyperplane( get_dual_entity( split_quads[1] ) );
02271     assert( 0 != sheet1 && 0 != chordl && 0 != chordr );
02272     Range parentsl, parentsr;
02273     ErrorCode result = mbImpl->get_parent_meshsets( chordl, parentsl );
02274     if( MB_SUCCESS != result ) return result;
02275     result = mbImpl->get_parent_meshsets( chordr, parentsr );
02276     if( MB_SUCCESS != result ) return result;
02277     parentsl.erase( sheet1 );
02278     parentsr.erase( sheet1 );
02279 
02280     // before deciding which one to delete, collect the other cells which must
02281     // be deleted, and all the chords/sheets they're on
02282     Range adj_ents, dual_ents, cells1or2;
02283     for( int i = 0; i < 3; i++ )
02284     {
02285         result = mbImpl->get_adjacencies( hexes, i, false, adj_ents, Interface::UNION );
02286         if( MB_SUCCESS != result ) return result;
02287     }
02288 
02289     // cache any adjacent hexes, for rebuilding the dual later
02290     result = mbImpl->get_adjacencies( adj_ents, 3, false, hexes, Interface::UNION );
02291     if( MB_SUCCESS != result ) return result;
02292 
02293     for( Range::iterator rit = adj_ents.begin(); rit != adj_ents.end(); ++rit )
02294     {
02295         EntityHandle this_ent = get_dual_entity( *rit );
02296         dual_ents.insert( this_ent );
02297         int dim = mbImpl->dimension_from_handle( this_ent );
02298         if( 1 == dim || 2 == dim ) cells1or2.insert( this_ent );
02299     }
02300 
02301     Range dual_hps;
02302     for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit )
02303         dual_hps.insert( get_dual_hyperplane( *rit ) );
02304 
02305     result = delete_dual_entities( dual_ents );
02306     if( MB_SUCCESS != result ) return result;
02307 
02308     // now decide which sheet to delete (to be merged into the other)
02309     EntityHandle sheet_delete = 0;
02310     if( is_blind( *parentsl.begin() ) )
02311         sheet_delete = *parentsl.begin();
02312     else if( is_blind( *parentsr.begin() ) )
02313         sheet_delete = *parentsr.begin();
02314     else
02315     {
02316         // neither is blind, take the one with fewer cells
02317         Range tmp_ents;
02318         int numl, numr;
02319         result = mbImpl->get_number_entities_by_handle( *parentsl.begin(), numl );
02320         if( MB_SUCCESS != result ) return result;
02321         result = mbImpl->get_number_entities_by_handle( *parentsr.begin(), numr );
02322         if( MB_SUCCESS != result ) return result;
02323         sheet_delete = ( numl > numr ? *parentsr.begin() : *parentsl.begin() );
02324     }
02325     assert( 0 != sheet_delete );
02326 
02327     // after deleting cells, check for empty chords & sheets, and delete those too
02328     for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
02329     {
02330         Range tmp_ents;
02331         result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
02332         if( MB_SUCCESS != result ) return result;
02333         if( tmp_ents.empty() )
02334         {
02335             result = mbImpl->delete_entities( &( *rit ), 1 );
02336             if( MB_SUCCESS != result ) return result;
02337         }
02338         else if( *rit == sheet_delete )
02339         {
02340             // delete the sheet
02341             result = mbImpl->delete_entities( &( *rit ), 1 );
02342             if( MB_SUCCESS != result ) return result;
02343         }
02344     }
02345 
02346     // now just to be safe, add the hexes bridge-adjacent across vertices
02347     // to the hexes we already have
02348     Range tmp_hexes;
02349     MeshTopoUtil mtu( mbImpl );
02350     for( Range::iterator rit = hexes.begin(); rit != hexes.end(); ++rit )
02351     {
02352         result = mtu.get_bridge_adjacencies( *rit, 0, 3, tmp_hexes );
02353         if( MB_SUCCESS != result ) return result;
02354     }
02355     hexes.merge( tmp_hexes );
02356 
02357     return MB_SUCCESS;
02358 }
02359 
02360 //! returns true if all vertices are dual to hexes (not faces)
02361 bool DualTool::is_blind( const EntityHandle chord_or_sheet )
02362 {
02363     // must be an entity set
02364     if( TYPE_FROM_HANDLE( chord_or_sheet ) != MBENTITYSET ) return false;
02365 
02366     // get the vertices
02367     Range verts, ents;
02368     ErrorCode result = mbImpl->get_entities_by_handle( chord_or_sheet, ents );
02369     if( MB_SUCCESS != result || ents.empty() ) return false;
02370 
02371     result = mbImpl->get_adjacencies( ents, 0, false, verts, Interface::UNION );
02372     if( MB_SUCCESS != result || verts.empty() ) return false;
02373 
02374     for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit )
02375     {
02376         // get dual entity for this vertex
02377         EntityHandle dual_ent = get_dual_entity( *rit );
02378         if( 0 == dual_ent ) continue;
02379         if( TYPE_FROM_HANDLE( dual_ent ) == MBQUAD ) return false;
02380     }
02381 
02382     // if none of the vertices' duals were quads, chord_or_sheet must be blind
02383     return true;
02384 }
02385 
02386 //! given a 1-cell and a chord, return the neighboring vertices on the
02387 //! chord, in the same order as the 1-cell's vertices
02388 ErrorCode DualTool::get_opposite_verts( const EntityHandle middle_edge, const EntityHandle chord, EntityHandle* verts )
02389 {
02390     // get the edges on the chord, in order, and move to middle_edge
02391     std::vector< EntityHandle > chord_edges;
02392     const EntityHandle* connect;
02393     int num_connect;
02394 
02395     ErrorCode result = mbImpl->get_entities_by_handle( chord, chord_edges );RR;
02396     std::vector< EntityHandle >::iterator vit = std::find( chord_edges.begin(), chord_edges.end(), middle_edge );
02397     result                                    = mbImpl->get_connectivity( middle_edge, connect, num_connect );RR;
02398 
02399     if(
02400         // middle_edge isn't on this chord
02401         vit == chord_edges.end() ||
02402         // chord only has 1 edge
02403         chord_edges.size() == 1 ||
02404         // middle_edge is at beginning or end and chord isn't blind
02405         ( ( vit == chord_edges.begin() || vit == chord_edges.end() - 1 ) && !is_blind( chord ) ) )
02406         return MB_FAILURE;
02407 
02408     else if( chord_edges.size() == 2 )
02409     {
02410         // easier if it's a 2-edge blind chord, just get vertices in opposite order
02411         verts[0] = connect[1];
02412         verts[1] = connect[0];
02413         return MB_SUCCESS;
02414     }
02415 
02416     // get vertices with the prev edge & subtract vertices of 1-cell
02417     if( vit == chord_edges.begin() )
02418         vit = chord_edges.end() - 1;
02419     else
02420         --vit;
02421     Range dum_connect, middle_connect;
02422     result = mbImpl->get_connectivity( &middle_edge, 1, middle_connect );RR;
02423     result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR;
02424     dum_connect = subtract( dum_connect, middle_connect );
02425     if( dum_connect.size() != 1 )
02426     {
02427         std::cerr << "Trouble traversing chord." << std::endl;
02428         return MB_FAILURE;
02429     }
02430 
02431     // put in verts[0]
02432     verts[0] = *dum_connect.begin();
02433 
02434     // same with prev edge
02435     ++vit;
02436     if( vit == chord_edges.end() ) vit = chord_edges.begin();
02437     ++vit;
02438     dum_connect.clear();
02439     result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR;
02440     dum_connect = subtract( dum_connect, middle_connect );
02441     if( dum_connect.size() != 1 )
02442     {
02443         std::cerr << "Trouble traversing chord." << std::endl;
02444         return MB_FAILURE;
02445     }
02446 
02447     // put in verts[1]
02448     verts[1] = *dum_connect.begin();
02449 
02450     // if verts[0] and 1st vertex of 1cell don't have common edge, switch verts
02451     MeshTopoUtil mtu( mbImpl );
02452     if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) )
02453     {
02454         EntityHandle dum_h = verts[0];
02455         verts[0]           = verts[1];
02456         verts[1]           = dum_h;
02457     }
02458 
02459     if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) )
02460     {
02461         std::cerr << "Trouble traversing chord." << std::endl;
02462         return MB_FAILURE;
02463     }
02464 
02465     return MB_SUCCESS;
02466 }
02467 
02468 ErrorCode DualTool::get_dual_entities( const EntityHandle dual_ent, Range* dcells, Range* dedges, Range* dverts,
02469                                        Range* dverts_loop, Range* dedges_loop )
02470 {
02471     ErrorCode result = MB_SUCCESS;
02472 
02473     if( NULL != dcells )
02474     {
02475         result = mbImpl->get_entities_by_type( dual_ent, MBPOLYGON, *dcells );
02476         if( MB_SUCCESS != result ) return result;
02477     }
02478 
02479     if( NULL != dedges )
02480     {
02481         if( NULL != dcells )
02482             result = mbImpl->get_adjacencies( *dcells, 1, false, *dedges, Interface::UNION );
02483         else
02484             result = mbImpl->get_entities_by_type( dual_ent, MBEDGE, *dedges );
02485 
02486         if( MB_SUCCESS != result ) return result;
02487     }
02488 
02489     if( NULL != dverts )
02490     {
02491         if( NULL != dcells )
02492             result = mbImpl->get_adjacencies( *dcells, 0, false, *dverts, Interface::UNION );
02493         else if( NULL != dedges )
02494             result = mbImpl->get_adjacencies( *dedges, 0, false, *dverts, Interface::UNION );
02495         else
02496         {
02497             Range all_ents;
02498             result = mbImpl->get_entities_by_handle( dual_ent, all_ents );RR;
02499             result = mbImpl->get_adjacencies( all_ents, 0, false, *dverts, Interface::UNION );
02500         }
02501 
02502         if( MB_SUCCESS != result ) return result;
02503     }
02504 
02505     if( NULL != dverts_loop && NULL != dverts )
02506     {
02507         static std::vector< EntityHandle > dual_ents;
02508         dual_ents.resize( dverts->size() );
02509         result = mbImpl->tag_get_data( dualEntity_tag(), *dverts, &dual_ents[0] );
02510         if( MB_SUCCESS != result ) return result;
02511         Range::iterator rit;
02512         unsigned int i;
02513         for( rit = dverts->begin(), i = 0; rit != dverts->end(); ++rit, i++ )
02514             if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBQUAD ) dverts_loop->insert( *rit );
02515     }
02516 
02517     if( NULL != dedges_loop && NULL != dedges )
02518     {
02519         static std::vector< EntityHandle > dual_ents;
02520         dual_ents.resize( dedges->size() );
02521         result = mbImpl->tag_get_data( dualEntity_tag(), *dedges, &dual_ents[0] );
02522         if( MB_SUCCESS != result ) return result;
02523         Range::iterator rit;
02524         unsigned int i;
02525         for( rit = dedges->begin(), i = 0; rit != dedges->end(); ++rit, i++ )
02526             if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBEDGE ) dedges_loop->insert( *rit );
02527     }
02528 
02529     return result;
02530 }
02531 
02532 ErrorCode DualTool::list_entities( const EntityHandle* entities, const int num_entities ) const
02533 {
02534     Range temp_range;
02535     ErrorCode result;
02536     if( NULL == entities && 0 == num_entities )
02537         return mbImpl->list_entities( entities, num_entities );
02538 
02539     else if( NULL == entities && 0 < num_entities )
02540     {
02541 
02542         // list all entities of all types
02543         std::cout << std::endl;
02544         for( EntityType this_type = MBVERTEX; this_type < MBMAXTYPE; this_type++ )
02545         {
02546             result = mbImpl->get_entities_by_type( 0, this_type, temp_range );
02547             if( MB_SUCCESS != result ) return result;
02548         }
02549     }
02550 
02551     else
02552     {
02553         std::copy( entities, entities + num_entities, range_inserter( temp_range ) );
02554     }
02555 
02556     return list_entities( temp_range );
02557 }
02558 
02559 ErrorCode DualTool::list_entities( const Range& entities ) const
02560 {
02561     // now print each entity, listing the dual information first then calling Interface to do
02562     // the rest
02563     ErrorCode result = MB_SUCCESS, tmp_result;
02564     for( Range::const_iterator iter = entities.begin(); iter != entities.end(); ++iter )
02565     {
02566         EntityType this_type = TYPE_FROM_HANDLE( *iter );
02567         std::cout << CN::EntityTypeName( this_type ) << " " << ID_FROM_HANDLE( *iter ) << ":" << std::endl;
02568 
02569         EntityHandle dual_ent = get_dual_entity( *iter );
02570         if( 0 != dual_ent )
02571         {
02572             std::cout << "Dual to " << CN::EntityTypeName( mbImpl->type_from_handle( dual_ent ) ) << " "
02573                       << mbImpl->id_from_handle( dual_ent ) << std::endl;
02574         }
02575 
02576         if( TYPE_FROM_HANDLE( *iter ) == MBENTITYSET )
02577         {
02578             EntityHandle chord = 0, sheet = 0;
02579             int id;
02580             result = mbImpl->tag_get_data( dualCurve_tag(), &( *iter ), 1, &chord );
02581             if( MB_SUCCESS != result ) return result;
02582             result = mbImpl->tag_get_data( dualSurface_tag(), &( *iter ), 1, &sheet );
02583             if( MB_SUCCESS != result ) return result;
02584             result = mbImpl->tag_get_data( globalId_tag(), &( *iter ), 1, &id );
02585             if( MB_SUCCESS != result ) return result;
02586 
02587             if( 0 != chord ) std::cout << "(Dual chord " << id << ")" << std::endl;
02588             if( 0 != sheet ) std::cout << "(Dual sheet " << id << ")" << std::endl;
02589         }
02590 
02591         tmp_result = mbImpl->list_entity( *iter );
02592         if( MB_SUCCESS != tmp_result ) result = tmp_result;
02593     }
02594 
02595     return result;
02596 }
02597 
02598 ErrorCode DualTool::face_shrink( EntityHandle odedge )
02599 {
02600     // some preliminary checking
02601     if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE;
02602 
02603     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
02604 
02605     std::cout << "FS(";
02606     print_cell( odedge );
02607     std::cout << ")" << std::endl;
02608 
02609     EntityHandle quads[4], hexes[2];
02610     std::vector< EntityHandle > connects[4], side_quads[2];
02611 
02612     // get the quads along the chord through the 2 hexes, and the vertices
02613     // for those quads
02614     ErrorCode result = fs_get_quads( odedge, quads, hexes, connects );
02615     if( MB_SUCCESS != result ) return result;
02616 
02617     // flip/rotate connect arrays so they align & are same sense
02618     result = fs_check_quad_sense( hexes[0], quads[0], connects );
02619     if( MB_SUCCESS != result ) return result;
02620 
02621     // get the quad loops along the "side" surfaces
02622     result = fs_get_quad_loops( hexes, connects, side_quads );
02623     if( MB_SUCCESS != result ) return result;
02624 
02625     // ok, done with setup; now delete dual entities affected by this operation,
02626     // which is all the entities adjacent to vertices of dual edge
02627     Range adj_verts, adj_edges, dual_ents, cells1or2;
02628     MeshTopoUtil mtu( mbImpl );
02629     result = mtu.get_bridge_adjacencies( odedge, 0, 1, adj_edges );
02630     if( MB_SUCCESS != result ) return result;
02631     result = mbImpl->get_adjacencies( adj_edges, 0, false, adj_verts, Interface::UNION );
02632     if( MB_SUCCESS != result ) return result;
02633     for( int i = 1; i <= 3; i++ )
02634     {
02635         result = mbImpl->get_adjacencies( adj_verts, i, false, dual_ents, Interface::UNION );
02636         if( MB_SUCCESS != result ) return result;
02637     }
02638 
02639     // before deleting dual, grab the 1- and 2-cells
02640     for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit )
02641     {
02642         int dim = mbImpl->dimension_from_handle( *rit );
02643         if( 1 == dim || 2 == dim ) cells1or2.insert( *rit );
02644     }
02645     Range dual_hps;
02646     for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit )
02647         dual_hps.insert( get_dual_hyperplane( *rit ) );
02648 
02649     dual_ents.insert( odedge );
02650     result = delete_dual_entities( dual_ents );
02651     if( MB_SUCCESS != result ) return result;
02652 
02653     // after deleting cells, check for empty chords & sheets, and delete those too
02654     for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
02655     {
02656         Range tmp_ents;
02657         result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
02658         if( MB_SUCCESS != result ) return result;
02659         if( tmp_ents.empty() )
02660         {
02661             result = mbImpl->delete_entities( &( *rit ), 1 );
02662             if( MB_SUCCESS != result ) return result;
02663         }
02664     }
02665 
02666     // remove any explicit adjacencies between side quads and hexes; don't check
02667     // for error, since there might not be adjacencies
02668     for( int i = 0; i < 4; i++ )
02669     {
02670         for( int j = 0; j < 2; j++ )
02671         {
02672             result = mbImpl->remove_adjacencies( side_quads[j][i], &hexes[j], 1 );
02673         }
02674     }
02675 
02676     // make inner ring of vertices
02677     // get centroid of quad2
02678     double q2coords[12], avg[3] = { 0.0, 0.0, 0.0 };
02679     result = mbImpl->get_coords( &connects[1][0], 4, q2coords );
02680     if( MB_SUCCESS != result ) return result;
02681     for( int i = 0; i < 4; i++ )
02682     {
02683         avg[0] += q2coords[3 * i];
02684         avg[1] += q2coords[3 * i + 1];
02685         avg[2] += q2coords[3 * i + 2];
02686     }
02687     avg[0] *= .25;
02688     avg[1] *= .25;
02689     avg[2] *= .25;
02690     // position new vertices
02691     connects[3].resize( 4 );
02692     for( int i = 0; i < 4; i++ )
02693     {
02694         q2coords[3 * i]     = avg[0] + .25 * ( q2coords[3 * i] - avg[0] );
02695         q2coords[3 * i + 1] = avg[1] + .25 * ( q2coords[3 * i + 1] - avg[1] );
02696         q2coords[3 * i + 2] = avg[2] + .25 * ( q2coords[3 * i + 2] - avg[2] );
02697         result              = mbImpl->create_vertex( &q2coords[3 * i], connects[3][i] );
02698         if( MB_SUCCESS != result ) return result;
02699     }
02700 
02701     // ok, now have the 4 connectivity arrays for 4 quads; construct hexes
02702     EntityHandle hconnect[8], new_hexes[4];
02703     int new_hex_ids[4];
02704 
02705     for( int i = 0; i < 4; i++ )
02706     {
02707         int i1 = i, i2 = ( i + 1 ) % 4;
02708         hconnect[0] = connects[0][i1];
02709         hconnect[1] = connects[0][i2];
02710         hconnect[2] = connects[3][i2];
02711         hconnect[3] = connects[3][i1];
02712 
02713         hconnect[4] = connects[1][i1];
02714         hconnect[5] = connects[1][i2];
02715         hconnect[6] = connects[2][i2];
02716         hconnect[7] = connects[2][i1];
02717 
02718         result = mbImpl->create_element( MBHEX, hconnect, 8, new_hexes[i] );
02719         if( MB_SUCCESS != result ) return result;
02720 
02721         // test for equiv entities from the side quads, and make explicit adjacencies
02722         // if there are any
02723         for( int j = 0; j < 2; j++ )
02724         {
02725             if( mtu.equivalent_entities( side_quads[j][i] ) )
02726             {
02727                 result = mbImpl->add_adjacencies( side_quads[j][i], &new_hexes[i], 1, false );
02728                 if( MB_SUCCESS != result ) return result;
02729             }
02730         }
02731 
02732         new_hex_ids[i] = ++maxHexId;
02733     }
02734 
02735     // set the global id tag on the new hexes
02736     result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 4, new_hex_ids );
02737     if( MB_SUCCESS != result ) return result;
02738 
02739     // now fixup other two hexes; start by getting hex through quads 0, 1
02740     // make this first hex switch to the other side, to make the dual look like
02741     // a hex push
02742     int tmp_ids[2];
02743     result = mbImpl->tag_get_data( globalId_tag(), hexes, 2, tmp_ids );
02744     if( MB_SUCCESS != result ) return result;
02745 
02746     result = mbImpl->delete_entities( hexes, 2 );
02747     if( MB_SUCCESS != result ) return result;
02748     result = mbImpl->delete_entities( &quads[1], 1 );
02749     if( MB_SUCCESS != result ) return result;
02750     for( int i = 0; i < 4; i++ )
02751     {
02752         hconnect[i]     = connects[3][i];
02753         hconnect[4 + i] = connects[2][i];
02754     }
02755     result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[0] );
02756     if( MB_SUCCESS != result ) return result;
02757 
02758     for( int i = 0; i < 4; i++ )
02759     {
02760         hconnect[i]     = connects[0][i];
02761         hconnect[4 + i] = connects[3][i];
02762     }
02763     result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[1] );
02764     if( MB_SUCCESS != result ) return result;
02765 
02766     // check for and fix any explicit adjacencies on either end quad
02767     if( mtu.equivalent_entities( quads[0] ) ) mbImpl->add_adjacencies( quads[0], &hexes[1], 1, false );
02768     if( mtu.equivalent_entities( quads[2] ) ) mbImpl->add_adjacencies( quads[2], &hexes[0], 1, false );
02769 
02770     // re-set the global ids for the hexes to what they were
02771     result = mbImpl->tag_set_data( globalId_tag(), hexes, 2, tmp_ids );
02772     if( MB_SUCCESS != result ) return result;
02773 
02774     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
02775 
02776     // now update the dual
02777     Range tmph;
02778     result = mtu.get_bridge_adjacencies( hexes[0], 0, 3, tmph );
02779     if( MB_SUCCESS != result ) return result;
02780     result = mtu.get_bridge_adjacencies( hexes[1], 0, 3, tmph );
02781     if( MB_SUCCESS != result ) return result;
02782     tmph.insert( hexes[1] );
02783     result = construct_hex_dual( tmph );
02784     if( MB_SUCCESS != result ) return result;
02785 
02786     return result;
02787 }
02788 
02789 ErrorCode DualTool::fs_get_quad_loops( EntityHandle* hexes, std::vector< EntityHandle >* connects,
02790                                        std::vector< EntityHandle >* side_quads )
02791 {
02792     for( int i = 0; i < 4; i++ )
02793     {
02794         for( int j = 0; j < 2; j++ )
02795         {
02796             Range adj_ents, dum_quads;
02797             adj_ents.insert( hexes[j] );
02798             adj_ents.insert( connects[j][i] );
02799             adj_ents.insert( connects[j][( i + 1 ) % 4] );
02800             adj_ents.insert( connects[j + 1][i] );
02801             adj_ents.insert( connects[j + 1][( i + 1 ) % 4] );
02802 
02803             ErrorCode result = mbImpl->get_adjacencies( adj_ents, 2, false, dum_quads );
02804             if( MB_SUCCESS != result ) return result;
02805             assert( 1 == dum_quads.size() );
02806             side_quads[j].push_back( *dum_quads.begin() );
02807         }
02808     }
02809 
02810     return MB_SUCCESS;
02811 }
02812 
02813 ErrorCode DualTool::fs_check_quad_sense( EntityHandle hex0, EntityHandle quad0, std::vector< EntityHandle >* connects )
02814 {
02815     // check sense of 0th quad wrt hex; since sense is out of element,
02816     // switch if quad is NOT reversed wrt hex
02817     int dum1, dum2, sense = 0;
02818     ErrorCode result = mbImpl->side_number( hex0, quad0, dum1, sense, dum2 );
02819     if( MB_SUCCESS != result ) return result;
02820     assert( 0 != sense );
02821     if( 1 == sense )
02822     {
02823         // just switch sense of this one; others will get switched next
02824         EntityHandle dum = connects[0][0];
02825         connects[0][0]   = connects[0][2];
02826         connects[0][2]   = dum;
02827     }
02828 
02829     // check sense of 1st, 2nd quads, rotate if necessary to align connect arrays
02830     int index0 = -1, index2 = -1, sense0 = 0, sense2 = 0;
02831     MeshTopoUtil mtu( mbImpl );
02832     for( int i = 0; i < 4; i++ )
02833     {
02834         if( 0 != mtu.common_entity( connects[0][0], connects[1][i], 1 ) )
02835         {
02836             index0 = i;
02837             if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 1 ) % 4], 1 ) )
02838                 sense0 = 1;
02839             else if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 4 - 1 ) % 4], 1 ) )
02840                 sense0 = -1;
02841             break;
02842         }
02843     }
02844 
02845     assert( index0 != -1 && sense0 != 0 );
02846 
02847     if( sense0 == -1 )
02848     {
02849         EntityHandle dumh = connects[1][0];
02850         connects[1][0]    = connects[1][2];
02851         connects[1][2]    = dumh;
02852         if( index0 % 2 == 0 ) index0 = ( index0 + 2 ) % 4;
02853     }
02854 
02855     if( index0 != 0 )
02856     {
02857         std::vector< EntityHandle > tmpc;
02858         for( int i = 0; i < 4; i++ )
02859             tmpc.push_back( connects[1][( index0 + i ) % 4] );
02860         connects[1].swap( tmpc );
02861     }
02862 
02863     for( int i = 0; i < 4; i++ )
02864     {
02865         if( 0 != mtu.common_entity( connects[1][0], connects[2][i], 1 ) )
02866         {
02867             index2 = i;
02868             if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 1 ) % 4], 1 ) )
02869                 sense2 = 1;
02870             else if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 4 - 1 ) % 4], 1 ) )
02871                 sense2 = -1;
02872             break;
02873         }
02874     }
02875 
02876     assert( index2 != -1 && sense2 != 0 );
02877 
02878     if( sense2 == -1 )
02879     {
02880         EntityHandle dumh = connects[2][0];
02881         connects[2][0]    = connects[2][2];
02882         connects[2][2]    = dumh;
02883         if( index2 % 2 == 0 ) index2 = ( index2 + 2 ) % 4;
02884     }
02885 
02886     if( index2 != 0 )
02887     {
02888         std::vector< EntityHandle > tmpc;
02889         for( int i = 0; i < 4; i++ )
02890             tmpc.push_back( connects[2][( index2 + i ) % 4] );
02891         connects[2].swap( tmpc );
02892     }
02893 
02894     return MB_SUCCESS;
02895 }
02896 
02897 //! effect reverse face shrink operation
02898 ErrorCode DualTool::rev_face_shrink( EntityHandle odedge )
02899 {
02900     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
02901 
02902     // some preliminary checking
02903     if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE;
02904 
02905     std::cout << "-FS(";
02906     print_cell( odedge );
02907     std::cout << ")" << std::endl;
02908 
02909     EntityHandle quads[4], hexes[2];
02910     std::vector< EntityHandle > connects[4], side_quads[2];
02911 
02912     // get three quads (shared quad & 2 end quads), hexes, and quad
02913     // connects
02914     ErrorCode result = fs_get_quads( odedge, quads, hexes, connects );
02915     if( MB_SUCCESS != result ) return result;
02916 
02917     // adjust sense & rotation so they're aligned, together & wrt first
02918     // hex
02919     result = fs_check_quad_sense( hexes[0], quads[0], connects );
02920     if( MB_SUCCESS != result ) return result;
02921 
02922     result = fsr_get_fourth_quad( connects, side_quads );
02923     if( MB_SUCCESS != result )
02924     {
02925         std::cout << "Can't do -FS here, two hexes must be adjacent to ring of 4 hexes." << std::endl;
02926         return result;
02927     }
02928 
02929     Range adj_ents, outer_hexes, all_adjs;
02930 
02931     // first get the entities connected to interior 4 verts
02932     for( int i = 1; i <= 3; i++ )
02933     {
02934         result = mbImpl->get_adjacencies( &connects[1][0], 4, i, false, adj_ents, Interface::UNION );
02935         if( MB_SUCCESS != result ) return result;
02936     }
02937 
02938     // next get all entities adjacent to those; these will have their dual
02939     // entities deleted
02940     for( int i = 0; i < 3; i++ )
02941     {
02942         result = mbImpl->get_adjacencies( adj_ents, i, false, all_adjs, Interface::UNION );
02943         if( MB_SUCCESS != result ) return result;
02944     }
02945 
02946     // get the dual entities and delete them
02947     Range dual_ents, dual_hps;
02948     for( Range::iterator rit = all_adjs.begin(); rit != all_adjs.end(); ++rit )
02949     {
02950         EntityHandle this_ent = get_dual_entity( *rit );
02951         dual_ents.insert( this_ent );
02952     }
02953 
02954     // before deleting dual, grab the 1- and 2-cells
02955     for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit )
02956     {
02957         int dim = mbImpl->dimension_from_handle( *rit );
02958         if( 1 == dim || 2 == dim ) dual_hps.insert( get_dual_hyperplane( *rit ) );
02959     }
02960 
02961     result = delete_dual_entities( dual_ents );
02962     if( MB_SUCCESS != result ) return result;
02963 
02964     // after deleting cells, check for empty chords & sheets, and delete those too
02965     for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
02966     {
02967         Range tmp_ents;
02968         result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
02969         if( MB_SUCCESS != result ) return result;
02970         if( tmp_ents.empty() )
02971         {
02972             result = mbImpl->delete_entities( &( *rit ), 1 );
02973             if( MB_SUCCESS != result ) return result;
02974         }
02975     }
02976 
02977     // before re-connecting two hexes, check for existing quad on 4th quad vertices;
02978     // if there is a quad there, need to add explicit adjs to any adj hexes, since
02979     // by definition there'll be another quad on those vertices
02980     bool need_explicit = false;
02981     Range adj_quads;
02982     result = mbImpl->get_adjacencies( &connects[3][0], 4, 2, false, adj_quads );
02983     if( MB_MULTIPLE_ENTITIES_FOUND == result || !adj_quads.empty() )
02984     {
02985         // there's already a quad for these 4 vertices; by definition,
02986         // we'll be creating equivalent entities, so that original quad
02987         // needs explicit adj's to its bounding elements
02988         need_explicit = true;
02989         for( Range::iterator rit = adj_quads.begin(); rit != adj_quads.end(); ++rit )
02990         {
02991             Range adj_hexes;
02992             result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, adj_hexes );RR;
02993             result = mbImpl->add_adjacencies( *rit, adj_hexes, false );RR;
02994         }
02995     }
02996 
02997     // re-connect the two hexes
02998     std::vector< EntityHandle > new_connect;
02999     std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) );
03000     std::copy( connects[2].begin(), connects[2].end(), std::back_inserter( new_connect ) );
03001     result = mbImpl->set_connectivity( hexes[0], &new_connect[0], 8 );
03002     if( MB_SUCCESS != result ) return result;
03003 
03004     new_connect.clear();
03005     std::copy( connects[0].begin(), connects[0].end(), std::back_inserter( new_connect ) );
03006     std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) );
03007     result = mbImpl->set_connectivity( hexes[1], &new_connect[0], 8 );
03008     if( MB_SUCCESS != result ) return result;
03009 
03010     // test for equiv entities from the side quads, and make explicit
03011     // adjacencies if there are any
03012     MeshTopoUtil mtu( mbImpl );
03013     for( int j = 0; j < 2; j++ )
03014     {
03015         for( int i = 0; i < 4; i++ )
03016         {
03017             if( mtu.equivalent_entities( side_quads[j][i] ) )
03018             {
03019                 result = mbImpl->add_adjacencies( side_quads[j][i], &hexes[j], 1, false );
03020                 if( MB_SUCCESS != result ) return result;
03021             }
03022         }
03023     }
03024 
03025     // remove hexes we want to keep
03026     adj_ents.erase( hexes[0] );
03027     adj_ents.erase( hexes[1] );
03028 
03029     // delete the other interior entities
03030     result = mbImpl->delete_entities( adj_ents );
03031     if( MB_SUCCESS != result ) return result;
03032 
03033     EntityHandle new_quad;
03034     result = mbImpl->create_element( MBQUAD, &connects[3][0], 4, new_quad );RR;
03035     if( need_explicit )
03036     {
03037         result = mbImpl->add_adjacencies( new_quad, hexes, 2, false );RR;
03038     }
03039 
03040     if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
03041 
03042     // now update the dual
03043     result = construct_hex_dual( hexes, 2 );
03044     if( MB_SUCCESS != result ) return result;
03045 
03046     return MB_SUCCESS;
03047 }
03048 
03049 ErrorCode DualTool::fsr_get_fourth_quad( std::vector< EntityHandle >* connects,
03050                                          std::vector< EntityHandle >* side_quads )
03051 {
03052     // given the first three quad connectivities in ordered vectors, get the fourth,
03053     // where the fourth is really the 4 vertices originally shared by the 2 hexes
03054     // before the face shrink on them
03055 
03056     // vertex on 4th quad is in quad adj to other 3 verts
03057     for( int i = 0; i < 4; i++ )
03058     {
03059         Range start_verts, tmp_verts, quads;
03060         for( int j = 0; j < 3; j++ )
03061             start_verts.insert( connects[j][i] );
03062         ErrorCode result = mbImpl->get_adjacencies( start_verts, 2, false, quads );
03063         if( MB_SUCCESS != result ) return result;
03064         assert( quads.size() == 1 );
03065         result = mbImpl->get_adjacencies( &( *quads.begin() ), 1, 0, false, tmp_verts );RR;
03066         tmp_verts = subtract( tmp_verts, start_verts );
03067         assert( 1 == tmp_verts.size() );
03068         connects[3].push_back( *tmp_verts.begin() );
03069     }
03070 
03071     // now get the side quads
03072     for( int i = 0; i < 4; i++ )
03073     {
03074         Range dum_ents, hexes;
03075         dum_ents.insert( connects[1][i] );
03076         dum_ents.insert( connects[1][( i + 1 ) % 4] );
03077         dum_ents.insert( connects[3][i] );
03078         ErrorCode result = mbImpl->get_adjacencies( dum_ents, 3, false, hexes );
03079         if( MB_SUCCESS != result ) return result;
03080         assert( 1 == hexes.size() );
03081 
03082         hexes.insert( connects[0][i] );
03083         hexes.insert( connects[0][( i + 1 ) % 4] );
03084         hexes.insert( connects[3][i] );
03085         hexes.insert( connects[3][( i + 1 ) % 4] );
03086         dum_ents.clear();
03087         result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents );
03088         if( MB_SUCCESS != result ) return result;
03089         assert( dum_ents.size() == 1 );
03090         side_quads[0].push_back( *dum_ents.begin() );
03091 
03092         hexes.erase( connects[0][i] );
03093         hexes.erase( connects[0][( i + 1 ) % 4] );
03094         hexes.insert( connects[2][i] );
03095         hexes.insert( connects[2][( i + 1 ) % 4] );
03096         dum_ents.clear();
03097         result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents );
03098         if( MB_SUCCESS != result ) return result;
03099         side_quads[1].push_back( *dum_ents.begin() );
03100     }
03101 
03102     return MB_SUCCESS;
03103 }
03104 
03105 ErrorCode DualTool::fs_get_quads( EntityHandle odedge, EntityHandle* quads, EntityHandle* hexes,
03106                                   std::vector< EntityHandle >* connects )
03107 {
03108     // need to get the three quads along the chord
03109     EntityHandle chord = get_dual_hyperplane( odedge );
03110     if( 0 == chord ) return MB_FAILURE;
03111 
03112     std::vector< EntityHandle > edges;
03113     ErrorCode result = mbImpl->get_entities_by_handle( chord, edges );
03114     if( MB_FAILURE == result ) return result;
03115 
03116     std::vector< EntityHandle >::iterator vit = std::find( edges.begin(), edges.end(), odedge );
03117     // shouldn't be first or last edge on chord
03118     if( vit == edges.end() || *edges.begin() == *vit || *edges.rbegin() == *vit ) return MB_FAILURE;
03119 
03120     // get quads/connectivity for first 3 quads
03121     quads[0] = get_dual_entity( *( vit - 1 ) );
03122     quads[1] = get_dual_entity( *vit );
03123     quads[2] = get_dual_entity( *( vit + 1 ) );
03124     for( int i = 0; i < 3; i++ )
03125     {
03126         result = mbImpl->get_connectivity( &quads[i], 1, connects[i], true );
03127         if( MB_SUCCESS != result ) return result;
03128     }
03129 
03130     Range tmph;
03131     result = mbImpl->get_adjacencies( quads, 2, 3, false, tmph );
03132     if( MB_SUCCESS != result ) return result;
03133     assert( tmph.size() == 1 );
03134     hexes[0] = *tmph.begin();
03135 
03136     tmph.clear();
03137     result = mbImpl->get_adjacencies( &quads[1], 2, 3, false, tmph );
03138     if( MB_SUCCESS != result ) return result;
03139     assert( tmph.size() == 1 );
03140     hexes[1] = *tmph.begin();
03141 
03142     return MB_SUCCESS;
03143 }
03144 
03145 ErrorCode DualTool::delete_whole_dual()
03146 {
03147     // delete dual hyperplanes
03148     Range dual_surfs, dual_curves;
03149     ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs );RR;
03150     result = mbImpl->delete_entities( dual_surfs );RR;
03151     result = this->get_dual_hyperplanes( mbImpl, 1, dual_curves );RR;
03152     result = mbImpl->delete_entities( dual_curves );RR;
03153 
03154     // gather up all dual entities
03155     Range dual_ents;
03156     result = mbImpl->get_entities_by_type_and_tag( 0, MBVERTEX, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
03157     result = mbImpl->get_entities_by_type_and_tag( 0, MBEDGE, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
03158     result = mbImpl->get_entities_by_type_and_tag( 0, MBPOLYGON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
03159     result =
03160         mbImpl->get_entities_by_type_and_tag( 0, MBPOLYHEDRON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
03161 
03162     // delete them, in reverse order of dimension
03163     ErrorCode tmp_result;
03164     for( Range::reverse_iterator rit = dual_ents.rbegin(); rit != dual_ents.rend(); ++rit )
03165     {
03166         tmp_result = mbImpl->delete_entities( &( *rit ), 1 );
03167         if( MB_SUCCESS != tmp_result ) result = tmp_result;
03168     }
03169     RR;
03170 
03171     // delete dual-related tags
03172     if( 0 != dualSurfaceTag )
03173     {
03174         tmp_result = mbImpl->tag_delete( dualSurfaceTag );
03175         if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
03176     }
03177     if( 0 != dualCurveTag )
03178     {
03179         tmp_result = mbImpl->tag_delete( dualCurveTag );
03180         if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
03181     }
03182     if( 0 != dualEntityTag )
03183     {
03184         tmp_result = mbImpl->tag_delete( dualEntityTag );
03185         if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
03186     }
03187     if( 0 != extraDualEntityTag )
03188     {
03189         tmp_result = mbImpl->tag_delete( extraDualEntityTag );
03190         if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
03191     }
03192     if( 0 != dualGraphicsPointTag )
03193     {
03194         tmp_result = mbImpl->tag_delete( dualGraphicsPointTag );
03195         if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
03196     }
03197 
03198     return MB_SUCCESS;
03199 }
03200 
03201 ErrorCode DualTool::check_dual_adjs()
03202 {
03203     // check primal-dual correspondence
03204 
03205     // get the primal entities
03206     Range pents[4];
03207     ErrorCode result = mbImpl->get_entities_by_type( 0, MBHEX, pents[3] );RR;
03208     for( int i = 2; i >= 0; i-- )
03209     {
03210         result = mbImpl->get_adjacencies( pents[3], 2, false, pents[2], Interface::UNION );RR;
03211     }
03212 
03213     // for each primal entity of dimension pd
03214 #define PRENT( ent ) CN::EntityTypeName( TYPE_FROM_HANDLE( ent ) ) << " " << ID_FROM_HANDLE( ent )
03215     ErrorCode overall_result = MB_SUCCESS;
03216     for( int pd = 1; pd <= 3; pd++ )
03217     {
03218         for( Range::iterator prit = pents[pd].begin(); prit != pents[pd].end(); ++prit )
03219         {
03220             // get corresponding dual entity of dimension dd = 3-pd
03221             EntityHandle dual_ent = get_dual_entity( *prit );
03222             if( 0 == dual_ent ) std::cerr << "Problem getting dual entity for " << PRENT( *prit ) << std::endl;
03223 
03224             // for each sub dimension sd = 0..pd-1
03225             for( int sd = 0; sd < pd; sd++ )
03226             {
03227                 Range R1, R2, R3;
03228                 //   R1 = entities bounding primal entity of dim sd
03229                 result = mbImpl->get_adjacencies( &( *prit ), 1, sd, false, R1 );RR;
03230 
03231                 //   R2 = entities bounded by dual entity, of dim 3-sd
03232                 result = mbImpl->get_adjacencies( &dual_ent, 1, 3 - sd, false, R2 );RR;
03233 
03234                 if( R1.size() != R2.size() )
03235                 {
03236                     std::cerr << PRENT( *prit ) << ": number of adj ents in "
03237                               << "primal/dual don't agree for dimension " << sd << "." << std::endl;
03238                     overall_result = MB_FAILURE;
03239                 }
03240 
03241                 // for each entity in R1, get its dual and look for it in R2
03242                 for( Range::iterator r1it = R1.begin(); r1it != R1.end(); ++r1it )
03243                 {
03244                     EntityHandle tmp_dual = get_dual_entity( *r1it );
03245                     if( R2.find( tmp_dual ) == R2.end() )
03246                     {
03247                         std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r1it ) << " isn't adjacent in dual."
03248                                   << std::endl;
03249                         overall_result = MB_FAILURE;
03250                     }
03251                 }
03252                 // ditto for R2
03253                 for( Range::iterator r2it = R2.begin(); r2it != R2.end(); ++r2it )
03254                 {
03255                     EntityHandle tmp_prim = get_dual_entity( *r2it );
03256                     if( R1.find( tmp_prim ) == R1.end() )
03257                     {
03258                         std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r2it ) << " isn't adjacent in primal."
03259                                   << std::endl;
03260                         overall_result = MB_FAILURE;
03261                     }
03262                 }
03263             }
03264         }
03265     }
03266 
03267     return overall_result;
03268 }
03269 
03270 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines