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