MOAB: Mesh Oriented datABase  (version 5.3.0)
GeomTopoTool.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/GeomTopoTool.hpp"
00017 #include "moab/GeomQueryTool.hpp"
00018 #include "moab/OrientedBoxTreeTool.hpp"
00019 #include "moab/Range.hpp"
00020 #include "MBTagConventions.hpp"
00021 #include "moab/Interface.hpp"
00022 #include "moab/CN.hpp"
00023 #include "moab/Skinner.hpp"
00024 #include "Internals.hpp"
00025 #include <iostream>
00026 
00027 namespace moab
00028 {
00029 
00030 // Tag name used for saving sense of faces in volumes.
00031 // We assume that the surface occurs in at most two volumes.
00032 // Code will error out if more than two volumes per surface.
00033 // The tag data is a pair of tag handles, representing the
00034 // forward and reverse volumes, respectively.  If a surface
00035 // is non-manifold in a single volume, the same volume will
00036 // be listed for both the forward and reverse slots.
00037 const char GEOM_SENSE_2_TAG_NAME[] = "GEOM_SENSE_2";
00038 
00039 const char GEOM_SENSE_N_ENTS_TAG_NAME[]   = "GEOM_SENSE_N_ENTS";
00040 const char GEOM_SENSE_N_SENSES_TAG_NAME[] = "GEOM_SENSE_N_SENSES";
00041 const char OBB_ROOT_TAG_NAME[]            = "OBB_ROOT";
00042 const char OBB_GSET_TAG_NAME[]            = "OBB_GSET";
00043 
00044 const char IMPLICIT_COMPLEMENT_NAME[] = "impl_complement";
00045 
00046 GeomTopoTool::GeomTopoTool( Interface* impl, bool find_geoments, EntityHandle modelRootSet, bool p_rootSets_vector,
00047                             bool restore_rootSets )
00048     : mdbImpl( impl ), sense2Tag( 0 ), senseNEntsTag( 0 ), senseNSensesTag( 0 ), geomTag( 0 ), gidTag( 0 ),
00049       obbRootTag( 0 ), obbGsetTag( 0 ), modelSet( modelRootSet ), updated( false ), setOffset( 0 ),
00050       m_rootSets_vector( p_rootSets_vector ), oneVolRootSet( 0 )
00051 {
00052 
00053     obbTree = new OrientedBoxTreeTool( impl, NULL, true );
00054 
00055     ErrorCode rval =
00056         mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_CONT( rval, "Error: Failed to create geometry dimension tag" );
00057 
00058     // global id tag is not really needed, but mbsize complains if we do not set it for
00059     // geometry entities
00060     gidTag = mdbImpl->globalId_tag();
00061 
00062     rval =
00063         mdbImpl->tag_get_handle( NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, nameTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_CONT( rval, "Error: Failed to create name tag" );
00064 
00065     rval = mdbImpl->tag_get_handle( OBB_ROOT_TAG_NAME, 1, MB_TYPE_HANDLE, obbRootTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_CONT( rval, "Error: Failed to create obb root tag" );
00066 
00067     rval = mdbImpl->tag_get_handle( OBB_GSET_TAG_NAME, 1, MB_TYPE_HANDLE, obbGsetTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_CONT( rval, "Error: Failed to create obb gset tag" );
00068 
00069     // set this value to zero for comparisons
00070     impl_compl_handle = 0;
00071 
00072     maxGlobalId[0] = maxGlobalId[1] = maxGlobalId[2] = maxGlobalId[3] = maxGlobalId[4] = 0;
00073     if( find_geoments )
00074     {
00075         find_geomsets();
00076         if( restore_rootSets )
00077         {
00078             rval = restore_obb_index();
00079             if( MB_SUCCESS != rval )
00080             {
00081                 rval = delete_all_obb_trees();MB_CHK_SET_ERR_CONT( rval, "Error: Failed to delete existing obb trees" );
00082                 rval = construct_obb_trees();MB_CHK_SET_ERR_CONT( rval, "Error: Failed to rebuild obb trees" );
00083             }
00084         }
00085     }
00086 }
00087 
00088 GeomTopoTool::~GeomTopoTool()
00089 {
00090     delete obbTree;
00091 }
00092 
00093 int GeomTopoTool::dimension( EntityHandle this_set )
00094 {
00095     ErrorCode result;
00096     if( 0 == geomTag )
00097     {
00098         result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );MB_CHK_SET_ERR( result, "Failed to get the geometry dimension tag" );
00099     }
00100 
00101     // check if the geo set belongs to this model
00102     if( modelSet )
00103     {
00104         if( !mdbImpl->contains_entities( modelSet, &this_set, 1 ) )
00105         {
00106             // this g set does not belong to the current model
00107             return -1;
00108         }
00109     }
00110     // get the data for those tags
00111     int dim;
00112     result = mdbImpl->tag_get_data( geomTag, &this_set, 1, &dim );
00113     if( MB_SUCCESS != result ) return -1;
00114     return dim;
00115 }
00116 
00117 int GeomTopoTool::global_id( EntityHandle this_set )
00118 {
00119     ErrorCode result;
00120     if( 0 == gidTag ) { gidTag = mdbImpl->globalId_tag(); }
00121 
00122     // check if the geo set belongs to this model
00123     if( modelSet )
00124     {
00125         if( !mdbImpl->contains_entities( modelSet, &this_set, 1 ) )
00126         {
00127             // this g set does not belong to the current model
00128             return -1;
00129         }
00130     }
00131 
00132     // get the data for those tags
00133     int id;
00134     result = mdbImpl->tag_get_data( gidTag, &this_set, 1, &id );
00135     if( MB_SUCCESS != result ) return -1;
00136     return id;
00137 }
00138 
00139 EntityHandle GeomTopoTool::entity_by_id( int dimension1, int id )
00140 {
00141     if( 0 > dimension1 && 3 < dimension1 ) { MB_CHK_SET_ERR_CONT( MB_FAILURE, "Incorrect dimension provided" ); };
00142     const Tag tags[]         = { gidTag, geomTag };
00143     const void* const vals[] = { &id, &dimension1 };
00144     ErrorCode rval;
00145 
00146     Range results;
00147     rval = mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, results );
00148 
00149     if( MB_SUCCESS != rval ) return 0;
00150 
00151     return results.front();
00152 }
00153 
00154 ErrorCode GeomTopoTool::other_entity( EntityHandle bounded, EntityHandle not_this, EntityHandle across,
00155                                       EntityHandle& other )
00156 {
00157     other = 0;
00158 
00159     // get all children of bounded
00160     Range bdy, tmpr;
00161     ErrorCode rval = mdbImpl->get_child_meshsets( bounded, bdy );MB_CHK_SET_ERR( rval, "Failed to get the bounded entity's child meshsets" );
00162 
00163     // get all the parents of across
00164     rval = mdbImpl->get_parent_meshsets( across, tmpr );
00165 
00166     // possible candidates is the intersection
00167     bdy = intersect( bdy, tmpr );
00168 
00169     // if only two, choose the other
00170     if( 1 == bdy.size() && *bdy.begin() == not_this ) { return MB_SUCCESS; }
00171     else if( 2 == bdy.size() )
00172     {
00173         if( *bdy.begin() == not_this ) other = *bdy.rbegin();
00174         if( *bdy.rbegin() == not_this )
00175             other = *bdy.begin();
00176         else
00177             return MB_FAILURE;
00178     }
00179     else
00180     {
00181         return MB_FAILURE;
00182     }
00183 
00184     return MB_SUCCESS;
00185 }
00186 
00187 ErrorCode GeomTopoTool::restore_obb_index()
00188 {
00189 
00190     if( m_rootSets_vector ) resize_rootSets();
00191 
00192     ErrorCode rval;
00193     EntityHandle root;
00194 
00195     for( int dim = 2; dim <= 3; dim++ )
00196         for( Range::iterator rit = geomRanges[dim].begin(); rit != geomRanges[dim].end(); ++rit )
00197         {
00198             rval = mdbImpl->tag_get_data( obbRootTag, &( *rit ), 1, &root );
00199 
00200             if( MB_SUCCESS == rval )
00201                 set_root_set( *rit, root );
00202             else
00203             {
00204                 return MB_TAG_NOT_FOUND;
00205             }
00206         }
00207 
00208     return MB_SUCCESS;
00209 }
00210 
00211 ErrorCode GeomTopoTool::find_geomsets( Range* ranges )
00212 {
00213     ErrorCode rval;
00214     // get all sets with this tag
00215     Range geom_sets;
00216 
00217     if( 0 == geomTag )
00218     {
00219         rval = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );MB_CHK_SET_ERR( rval, "Failed to get geom dimension tag handle" );
00220     }
00221 
00222     rval = mdbImpl->get_entities_by_type_and_tag( modelSet, MBENTITYSET, &geomTag, NULL, 1, geom_sets );MB_CHK_SET_ERR( rval, "Failed to get the geometry entities" );
00223 
00224     rval = separate_by_dimension( geom_sets );MB_CHK_SET_ERR( rval, "Failed to separate geometry sets by dimension" );
00225 
00226     if( ranges )
00227     {
00228         for( int i = 0; i < 5; i++ )
00229         {
00230             ranges[i] = geomRanges[i];
00231         }
00232     }
00233 
00234     return MB_SUCCESS;
00235 }
00236 
00237 ErrorCode GeomTopoTool::get_gsets_by_dimension( int dim, Range& gset )
00238 {
00239     ErrorCode rval;
00240 
00241     const int val               = dim;
00242     const void* const dim_val[] = { &val };
00243     rval = mdbImpl->get_entities_by_type_and_tag( modelSet, MBENTITYSET, &geomTag, dim_val, 1, gset );MB_CHK_SET_ERR( rval, "Failed to get entity set by type and tag" );
00244 
00245     return MB_SUCCESS;
00246 }
00247 
00248 ErrorCode GeomTopoTool::resize_rootSets()
00249 {
00250 
00251     ErrorCode rval;
00252 
00253     // store original offset for later
00254     EntityHandle orig_offset = setOffset;
00255 
00256     // get all surfaces and volumes
00257     Range surfs, vols;
00258     rval = get_gsets_by_dimension( 2, surfs );MB_CHK_SET_ERR( rval, "Could not get surface sets" );
00259     rval = get_gsets_by_dimension( 3, vols );MB_CHK_SET_ERR( rval, "Could not get volume sets" );
00260 
00261     // check the vector size
00262     Range surfs_and_vols;
00263     surfs_and_vols = vols;
00264     surfs_and_vols.merge( surfs );
00265 
00266     // update the setOffset
00267     setOffset = surfs_and_vols.front();
00268 
00269     EntityHandle exp_size = surfs_and_vols.back() - setOffset + 1;
00270 
00271     // if new EnitytHandle(s) are lower than the original offset
00272     if( setOffset < orig_offset )
00273     {
00274         // insert empty values at the beginning of the vector
00275         rootSets.insert( rootSets.begin(), orig_offset - setOffset, 0 );
00276     }
00277 
00278     if( exp_size != rootSets.size() )
00279     {
00280         // resize rootSets vector if necessary (new space will be added at the back)
00281         rootSets.resize( exp_size );
00282     }
00283 
00284     return MB_SUCCESS;
00285 }
00286 
00287 ErrorCode GeomTopoTool::is_owned_set( EntityHandle eh )
00288 {
00289     // make sure entity set is part of the model
00290     Range model_ents;
00291     ErrorCode rval = mdbImpl->get_entities_by_handle( modelSet, model_ents );MB_CHK_SET_ERR( rval, "Failed to get entities" );
00292     if( model_ents.find( eh ) == model_ents.end() ) { MB_SET_ERR( MB_FAILURE, "Entity handle not in model set" ); }
00293     return MB_SUCCESS;
00294 }
00295 
00296 ErrorCode GeomTopoTool::delete_obb_tree( EntityHandle gset, bool vol_only )
00297 {
00298 
00299     ErrorCode rval;
00300 
00301     // Make sure this set is part of the model
00302     rval = is_owned_set( gset );MB_CHK_SET_ERR( rval, "Entity set is not part of this model" );
00303 
00304     // Find the dimension of the entity
00305     int dim;
00306     rval = mdbImpl->tag_get_data( geomTag, &gset, 1, &dim );MB_CHK_SET_ERR( rval, "Failed to get dimension" );
00307 
00308     // Attempt to find a root for this set
00309     EntityHandle root;
00310     rval = get_root( gset, root );MB_CHK_SET_ERR( rval, "Failed to find an obb tree root for the entity set" );
00311 
00312     // Create range of tree nodes to delete
00313     Range nodes_to_delete;
00314     nodes_to_delete.insert( root );
00315 
00316     // If passed ent is a vol and 'vol_only' is true, delete vol root and all nodes between vol and
00317     // surf root
00318     if( dim == 3 && vol_only )
00319     {
00320         // Range of child nodes to check before adding to delete list
00321         Range child_tree_nodes;
00322         rval = mdbImpl->get_child_meshsets( root, child_tree_nodes );MB_CHK_SET_ERR( rval, "Problem getting child tree nodes" );
00323 
00324         // Traverse the tree, checking each child node until
00325         // a surface root node is reached
00326         while( child_tree_nodes.size() != 0 )
00327         {
00328             EntityHandle child = *child_tree_nodes.begin();
00329             EntityHandle surf;
00330             rval = mdbImpl->tag_get_data( obbGsetTag, &child, 1, &surf );
00331             // If the node has a gset tag, it is a surf root. Stop here.
00332             // If not, it is a tree node that needs to 1) have its children checked and
00333             //  2) be added to delete range
00334             if( MB_TAG_NOT_FOUND == rval )
00335             {
00336                 Range new_child_tree_nodes;
00337                 rval = mdbImpl->get_child_meshsets( child, new_child_tree_nodes );MB_CHK_SET_ERR( rval, "Problem getting child nodes" );
00338                 child_tree_nodes.insert_list( new_child_tree_nodes.begin(), new_child_tree_nodes.end() );
00339                 nodes_to_delete.insert( child );
00340             }
00341             // We're done checking this node, so can erase from child list
00342             child_tree_nodes.erase( child );
00343         }
00344     }
00345     // If passed ent is a surf or a vol and 'vol_only' is false, recursively gather all child nodes
00346     // and add them to delete list
00347     else
00348     {
00349         Range all_tree_nodes;
00350         rval = mdbImpl->get_child_meshsets( root, all_tree_nodes, 0 );MB_CHK_SET_ERR( rval, "Failed to get child tree node sets" );
00351         nodes_to_delete.insert_list( all_tree_nodes.begin(), all_tree_nodes.end() );
00352     }
00353 
00354     // Remove the root nodes from the GTT data structures
00355     for( Range::iterator it = nodes_to_delete.begin(); it != nodes_to_delete.end(); ++it )
00356     {
00357         // Check to see if node is a root
00358         EntityHandle vol_or_surf;
00359         rval = mdbImpl->tag_get_data( obbGsetTag, &( *it ), 1, &vol_or_surf );
00360         if( MB_SUCCESS == rval )
00361         {
00362             // Remove from set of all roots
00363             rval = remove_root( vol_or_surf );MB_CHK_SET_ERR( rval, "Failed to remove node from GTT data structure" );
00364         }
00365     }
00366 
00367     // Delete the tree nodes from the database
00368     rval = mdbImpl->delete_entities( nodes_to_delete );MB_CHK_SET_ERR( rval, "Failed to delete node set" );
00369 
00370     return MB_SUCCESS;
00371 }
00372 
00373 ErrorCode GeomTopoTool::delete_all_obb_trees()
00374 {
00375 
00376     ErrorCode rval;
00377 
00378     for( Range::iterator rit = geomRanges[3].begin(); rit != geomRanges[3].end(); ++rit )
00379     {
00380         EntityHandle root;
00381         rval = mdbImpl->tag_get_data( obbRootTag, &( *rit ), 1, &root );
00382         if( MB_SUCCESS == rval )
00383         {
00384             rval = delete_obb_tree( *rit, false );MB_CHK_SET_ERR( rval, "Failed to delete obb tree" );
00385         }
00386     }
00387 
00388     return MB_SUCCESS;
00389 }
00390 
00391 ErrorCode GeomTopoTool::construct_obb_tree( EntityHandle eh )
00392 {
00393     ErrorCode rval;
00394     int dim;
00395 
00396     rval = is_owned_set( eh );MB_CHK_SET_ERR( rval, "Entity set is not part of this model" );
00397 
00398     // get the type
00399     EntityType type = mdbImpl->type_from_handle( eh );
00400 
00401     // find the dimension of the entity
00402     rval = mdbImpl->tag_get_data( geomTag, &eh, 1, &dim );MB_CHK_SET_ERR( rval, "Failed to get dimension" );
00403 
00404     // ensure that the rootSets vector is of the correct size
00405     if( m_rootSets_vector && ( eh < setOffset || eh >= setOffset + rootSets.size() ) )
00406     {
00407         rval = resize_rootSets();MB_CHK_SET_ERR( rval, "Error setting offset and sizing rootSets vector." );
00408     }
00409 
00410     EntityHandle root;
00411     // if it's a surface
00412     if( dim == 2 && type == MBENTITYSET )
00413     {
00414         rval = get_root( eh, root );
00415         if( MB_SUCCESS == rval )
00416         {
00417             std::cerr << "Surface obb tree already exists" << std::endl;
00418             return MB_SUCCESS;
00419         }
00420         else if( MB_INDEX_OUT_OF_RANGE != rval )
00421         {
00422             MB_CHK_SET_ERR( rval, "Failed to get surface obb tree root" );
00423         }
00424 
00425         Range tris;
00426         rval = mdbImpl->get_entities_by_dimension( eh, 2, tris );MB_CHK_SET_ERR( rval, "Failed to get entities by dimension" );
00427 
00428         if( tris.empty() ) { std::cerr << "WARNING: Surface has no facets" << std::endl; }
00429 
00430         rval = obbTree->build( tris, root );MB_CHK_SET_ERR( rval, "Failed to build obb Tree for surface" );
00431 
00432         rval = mdbImpl->add_entities( root, &eh, 1 );MB_CHK_SET_ERR( rval, "Failed to add entities to root set" );
00433 
00434         // add this root to the GeomTopoTool tree root indexing
00435         set_root_set( eh, root );
00436 
00437         // if just building tree for surface, return here
00438         return MB_SUCCESS;
00439     }
00440     // if it's a volume
00441     else if( dim == 3 && type == MBENTITYSET )
00442     {
00443         // get its surfaces
00444         Range tmp_surfs, surf_trees;
00445         rval = mdbImpl->get_child_meshsets( eh, tmp_surfs );MB_CHK_SET_ERR( rval, "Failed to get surface meshsets" );
00446 
00447         // get OBB trees or create for each surface
00448         for( Range::iterator j = tmp_surfs.begin(); j != tmp_surfs.end(); ++j )
00449         {
00450             rval = get_root( *j, root );
00451             // if root doesn't exist, create obb tree
00452             if( MB_INDEX_OUT_OF_RANGE == rval )
00453             {
00454                 rval = construct_obb_tree( *j );MB_CHK_SET_ERR( rval, "Failed to get create surface obb tree" );
00455                 rval = get_root( *j, root );MB_CHK_SET_ERR( rval, "Failed to get surface obb tree root" );
00456             }
00457             else
00458             {
00459                 MB_CHK_SET_ERR( rval, "Failed to get surface obb tree root" );
00460             }
00461 
00462             surf_trees.insert( root );
00463         }
00464 
00465         // build OBB tree for volume
00466         rval = obbTree->join_trees( surf_trees, root );MB_CHK_SET_ERR( rval, "Failed to join the obb trees" );
00467 
00468         // add this root to the GeomTopoTool tree root indexing
00469         set_root_set( eh, root );
00470 
00471         return MB_SUCCESS;
00472     }
00473     else
00474     {
00475         MB_SET_ERR( MB_FAILURE, "Improper dimension or type for constructing obb tree" );
00476     }
00477 }
00478 
00479 ErrorCode GeomTopoTool::set_root_set( EntityHandle vol_or_surf, EntityHandle root )
00480 {
00481 
00482     // Tag the vol or surf with its obb root (obbRootTag)
00483     ErrorCode rval;
00484     rval = mdbImpl->tag_set_data( obbRootTag, &vol_or_surf, 1, &root );MB_CHK_SET_ERR( rval, "Failed to set the obb root tag" );
00485 
00486     // Tag obb root with corresponding gset (obbGsetTag)
00487     rval = mdbImpl->tag_set_data( obbGsetTag, &root, 1, &vol_or_surf );MB_CHK_SET_ERR( rval, "Failed to set the obb gset tag" );
00488 
00489     // Add to the set of all roots
00490     if( m_rootSets_vector )
00491         rootSets[vol_or_surf - setOffset] = root;
00492     else
00493         mapRootSets[vol_or_surf] = root;
00494 
00495     return MB_SUCCESS;
00496 }
00497 
00498 ErrorCode GeomTopoTool::remove_root( EntityHandle vol_or_surf )
00499 {
00500 
00501     // Find the root of the vol or surf
00502     ErrorCode rval;
00503     EntityHandle root;
00504     rval = mdbImpl->tag_get_data( obbRootTag, &( vol_or_surf ), 1, &root );MB_CHK_SET_ERR( rval, "Failed to get obb root tag" );
00505 
00506     // If the ent is a vol, remove its root from obbtreetool
00507     int dim;
00508     rval = mdbImpl->tag_get_data( geomTag, &vol_or_surf, 1, &dim );MB_CHK_SET_ERR( rval, "Failed to get dimension" );
00509     if( dim == 3 )
00510     {
00511         rval = obbTree->remove_root( root );MB_CHK_SET_ERR( rval, "Failed to remove root from obbTreeTool" );
00512     }
00513 
00514     // Delete the obbGsetTag data from the root
00515     rval = mdbImpl->tag_delete_data( obbGsetTag, &root, 1 );MB_CHK_SET_ERR( rval, "Failed to delete obb root tag" );
00516 
00517     // Delete the obbRootTag data from the vol or surf
00518     rval = mdbImpl->tag_delete_data( obbRootTag, &vol_or_surf, 1 );MB_CHK_SET_ERR( rval, "Failed to delete obb root tag" );
00519 
00520     // Remove the root from set of all roots
00521     if( m_rootSets_vector )
00522     {
00523         unsigned int index = vol_or_surf - setOffset;
00524         if( index < rootSets.size() ) { rootSets[index] = 0; }
00525         else
00526         {
00527             return MB_INDEX_OUT_OF_RANGE;
00528         }
00529     }
00530     else
00531     {
00532         mapRootSets[vol_or_surf] = 0;
00533     }
00534 
00535     return MB_SUCCESS;
00536 }
00537 
00538 ErrorCode GeomTopoTool::construct_obb_trees( bool make_one_vol )
00539 {
00540     ErrorCode rval;
00541     EntityHandle root;
00542 
00543     // get all surfaces and volumes
00544     Range surfs, vols, vol_trees;
00545     rval = get_gsets_by_dimension( 2, surfs );MB_CHK_SET_ERR( rval, "Could not get surface sets" );
00546     rval = get_gsets_by_dimension( 3, vols );MB_CHK_SET_ERR( rval, "Could not get volume sets" );
00547 
00548     // for surface
00549     Range one_vol_trees;
00550     for( Range::iterator i = surfs.begin(); i != surfs.end(); ++i )
00551     {
00552         rval = construct_obb_tree( *i );MB_CHK_SET_ERR( rval, "Failed to construct obb tree for surface" );
00553         // get the root set of this volume
00554         rval = get_root( *i, root );MB_CHK_SET_ERR( rval, "Failed to get obb tree root for surface" );
00555         // add to the Range of volume root sets
00556         one_vol_trees.insert( root );
00557     }
00558 
00559     // for volumes
00560     for( Range::iterator i = vols.begin(); i != vols.end(); ++i )
00561     {
00562         // create tree for this volume
00563         rval = construct_obb_tree( *i );MB_CHK_SET_ERR( rval, "Failed to construct obb tree for volume" );
00564     }
00565 
00566     // build OBB tree for volume
00567     if( make_one_vol )
00568     {
00569         rval = obbTree->join_trees( one_vol_trees, root );MB_CHK_SET_ERR( rval, "Failed to join surface trees into one volume" );
00570         oneVolRootSet = root;
00571     }
00572 
00573     return rval;
00574 }
00575 
00576 //! Restore parent/child links between GEOM_TOPO mesh sets
00577 ErrorCode GeomTopoTool::restore_topology_from_adjacency()
00578 {
00579 
00580     // look for geometric topology sets and restore parent/child links between them
00581     // algorithm:
00582     // - for each entity of dimension d=D-1..0:
00583     //   . get d-dimensional entity in entity
00584     //   . get all (d+1)-dim adjs to that entity
00585     //   . for each geom entity if dim d+1, if it contains any of the ents,
00586     //     add it to list of parents
00587     //   . make parent/child links with parents
00588 
00589     // the  geomRanges are already known, separated by dimension
00590 
00591     std::vector< EntityHandle > parents;
00592     Range tmp_parents;
00593     ErrorCode result;
00594 
00595     // loop over dimensions
00596     for( int dim = 2; dim >= 0; dim-- )
00597     {
00598         // mark entities of next higher dimension with their owners; regenerate tag
00599         // each dimension so prev dim's tag data goes away
00600         Tag owner_tag;
00601         EntityHandle dum_val = 0;
00602         result = mdbImpl->tag_get_handle( "__owner_tag", 1, MB_TYPE_HANDLE, owner_tag, MB_TAG_DENSE | MB_TAG_CREAT,
00603                                           &dum_val );
00604         if( MB_SUCCESS != result ) continue;
00605         Range dp1ents;
00606         std::vector< EntityHandle > owners;
00607         for( Range::iterator rit = geomRanges[dim + 1].begin(); rit != geomRanges[dim + 1].end(); ++rit )
00608         {
00609             dp1ents.clear();
00610             result = mdbImpl->get_entities_by_dimension( *rit, dim + 1, dp1ents );
00611             if( MB_SUCCESS != result ) continue;
00612             owners.resize( dp1ents.size() );
00613             std::fill( owners.begin(), owners.end(), *rit );
00614             result = mdbImpl->tag_set_data( owner_tag, dp1ents, &owners[0] );
00615             if( MB_SUCCESS != result ) continue;
00616         }
00617 
00618         for( Range::iterator d_it = geomRanges[dim].begin(); d_it != geomRanges[dim].end(); ++d_it )
00619         {
00620             Range dents;
00621             result = mdbImpl->get_entities_by_dimension( *d_it, dim, dents );
00622             if( MB_SUCCESS != result ) continue;
00623             if( dents.empty() ) continue;
00624 
00625             // get (d+1)-dimensional adjs
00626             dp1ents.clear();
00627             result = mdbImpl->get_adjacencies( &( *dents.begin() ), 1, dim + 1, false, dp1ents );
00628             if( MB_SUCCESS != result || dp1ents.empty() ) continue;
00629 
00630             // get owner tags
00631             parents.resize( dp1ents.size() );
00632             result = mdbImpl->tag_get_data( owner_tag, dp1ents, &parents[0] );
00633             if( MB_TAG_NOT_FOUND == result ) { MB_CHK_SET_ERR( result, "Could not find owner tag" ); }
00634             if( MB_SUCCESS != result ) continue;
00635 
00636             // compress to a range to remove duplicates
00637             tmp_parents.clear();
00638             std::copy( parents.begin(), parents.end(), range_inserter( tmp_parents ) );
00639             for( Range::iterator pit = tmp_parents.begin(); pit != tmp_parents.end(); ++pit )
00640             {
00641                 result = mdbImpl->add_parent_child( *pit, *d_it );MB_CHK_SET_ERR( result, "Failed to create parent child relationship" );
00642             }
00643 
00644             // store surface senses within regions, and edge senses within surfaces
00645             if( dim == 0 ) continue;
00646             const EntityHandle *conn3 = NULL, *conn2 = NULL;
00647             int len3 = 0, len2 = 0, err = 0, num = 0, sense = 0, offset = 0;
00648             for( size_t i = 0; i < parents.size(); ++i )
00649             {
00650                 result = mdbImpl->get_connectivity( dp1ents[i], conn3, len3, true );MB_CHK_SET_ERR( result, "Failed to get the connectivity of the element" );
00651                 result = mdbImpl->get_connectivity( dents.front(), conn2, len2, true );MB_CHK_SET_ERR( result, "Failed to get the connectivity of the first element" );
00652                 if( len2 > 4 ) { MB_SET_ERR( MB_FAILURE, "Connectivity of incorrect length" ); }
00653                 err = CN::SideNumber( TYPE_FROM_HANDLE( dp1ents[i] ), conn3, conn2, len2, dim, num, sense, offset );
00654                 if( err ) return MB_FAILURE;
00655 
00656                 result = set_sense( *d_it, parents[i], sense );
00657                 if( MB_MULTIPLE_ENTITIES_FOUND == result )
00658                 {
00659                     if( 2 == dim )
00660                         std::cerr << "Warning: Multiple volumes use surface with same sense." << std::endl
00661                                   << "         Some geometric sense data lost." << std::endl;
00662                 }
00663                 else if( MB_SUCCESS != result )
00664                 {
00665                     return result;
00666                 }
00667             }
00668         }
00669 
00670         // now delete owner tag on this dimension, automatically removes tag data
00671         result = mdbImpl->tag_delete( owner_tag );MB_CHK_SET_ERR( result, "Failed to delete the owner tag" );
00672 
00673     }  // dim
00674 
00675     return result;
00676 }
00677 
00678 ErrorCode GeomTopoTool::separate_by_dimension( const Range& geom_sets )
00679 {
00680     ErrorCode result;
00681 
00682     result = check_geom_tag();MB_CHK_SET_ERR( result, "Could not verify geometry dimension tag" );
00683 
00684     // get the data for those tags
00685     std::vector< int > tag_vals( geom_sets.size() );
00686     result = mdbImpl->tag_get_data( geomTag, geom_sets, &tag_vals[0] );MB_CHK_SET_ERR( result, "Failed to get the geometry dimension tag" );
00687 
00688     Range::const_iterator git;
00689     std::vector< int >::iterator iit;
00690 
00691     for( int i = 0; i < 5; i++ )
00692         this->geomRanges[i].clear();
00693 
00694     for( git = geom_sets.begin(), iit = tag_vals.begin(); git != geom_sets.end(); ++git, ++iit )
00695     {
00696         if( 0 <= *iit && 4 >= *iit )
00697             geomRanges[*iit].insert( *git );
00698         else
00699         {
00700             // assert(false);
00701             // do nothing for now
00702         }
00703     }
00704 
00705     // establish the max global ids so far, per dimension
00706     if( 0 == gidTag ) { gidTag = mdbImpl->globalId_tag(); }
00707 
00708     for( int i = 0; i <= 4; i++ )
00709     {
00710         maxGlobalId[i] = 0;
00711         for( Range::iterator it = geomRanges[i].begin(); it != geomRanges[i].end(); ++it )
00712         {
00713             EntityHandle set = *it;
00714             int gid;
00715 
00716             result = mdbImpl->tag_get_data( gidTag, &set, 1, &gid );
00717             if( MB_SUCCESS == result )
00718             {
00719                 if( gid > maxGlobalId[i] ) maxGlobalId[i] = gid;
00720             }
00721         }
00722     }
00723 
00724     return MB_SUCCESS;
00725 }
00726 
00727 ErrorCode GeomTopoTool::construct_vertex_ranges( const Range& geom_sets, const Tag verts_tag )
00728 {
00729     // construct the vertex range for each entity and put on that tag
00730     Range *temp_verts, temp_elems;
00731     ErrorCode result = MB_SUCCESS;
00732     for( Range::const_iterator it = geom_sets.begin(); it != geom_sets.end(); ++it )
00733     {
00734         temp_elems.clear();
00735 
00736         // get all the elements in the set, recursively
00737         result = mdbImpl->get_entities_by_handle( *it, temp_elems, true );MB_CHK_SET_ERR( result, "Failed to get the geometry set entities" );
00738 
00739         // make the new range
00740         temp_verts = new( std::nothrow ) Range();
00741         if( NULL == temp_verts ) { MB_SET_ERR( MB_FAILURE, "Could not construct Range object" ); }
00742 
00743         // get all the verts of those elements; use get_adjacencies 'cuz it handles ranges better
00744         result = mdbImpl->get_adjacencies( temp_elems, 0, false, *temp_verts, Interface::UNION );
00745         if( MB_SUCCESS != result ) { delete temp_verts; }
00746         MB_CHK_SET_ERR( result, "Failed to get the element's adjacent vertices" );
00747 
00748         // store this range as a tag on the entity
00749         result = mdbImpl->tag_set_data( verts_tag, &( *it ), 1, &temp_verts );
00750         if( MB_SUCCESS != result ) { delete temp_verts; }
00751         MB_CHK_SET_ERR( result, "Failed to get the adjacent vertex data" );
00752 
00753         delete temp_verts;
00754         temp_verts = NULL;
00755     }
00756 
00757     return result;
00758 }
00759 
00760 //! Store sense of entity relative to wrt_entity.
00761 //!\return MB_MULTIPLE_ENTITIES_FOUND if surface already has a forward volume.
00762 //!        MB_SUCCESS if successful
00763 //!        otherwise whatever internal error code occurred.
00764 ErrorCode GeomTopoTool::set_sense( EntityHandle entity, EntityHandle wrt_entity, int sense )
00765 {
00766     // entity is lower dim (edge or face), wrt_entity is face or volume
00767     int edim   = dimension( entity );
00768     int wrtdim = dimension( wrt_entity );
00769     if( -1 == edim || -1 == wrtdim ) MB_SET_ERR( MB_FAILURE, "Non-geometric entity provided" );
00770     if( wrtdim - edim != 1 ) MB_SET_ERR( MB_FAILURE, "Entity dimension mismatch" );
00771     if( sense < -1 || sense > 1 ) MB_SET_ERR( MB_FAILURE, "Invalid sense data provided" );
00772 
00773     ErrorCode rval;
00774 
00775     if( 1 == edim )
00776     {
00777         // this case is about setting the sense of an edge in a face
00778         // it could be -1, 0 (rare, non manifold), or 1
00779         rval = check_edge_sense_tags( true );MB_CHK_SET_ERR( rval, "Failed to check the curve to surface sense tag handles" );
00780         std::vector< EntityHandle > higher_ents;
00781         std::vector< int > senses;
00782         rval = get_senses( entity, higher_ents, senses );  // the tags should be defined here
00783         // if there are no higher_ents, we are fine, we will just set them
00784         // if wrt_entity is not among higher_ents, we will add it to the list
00785         // it is possible the entity (edge set) has no prior faces adjancent; in that case, the
00786         // tag would not be set, and rval could be MB_TAG_NOT_FOUND; it is not a fatal error
00787         if( MB_SUCCESS != rval && MB_TAG_NOT_FOUND != rval )
00788             MB_CHK_SET_ERR( rval, "cannot determine sense tags for edge" );
00789         bool append = true;
00790         if( !higher_ents.empty() )
00791         {
00792             std::vector< EntityHandle >::iterator it = std::find( higher_ents.begin(), higher_ents.end(), wrt_entity );
00793             if( it != higher_ents.end() )
00794             {
00795                 // we should not reset the sense, if the sense is the same
00796                 // if the sense is different, put BOTH
00797                 unsigned int idx = it - higher_ents.begin();
00798                 int oldSense     = senses[idx];
00799                 if( oldSense == sense ) return MB_SUCCESS;  // sense already set fine, do not reset
00800                 if( 0 != oldSense && oldSense + sense != 0 ) return MB_MULTIPLE_ENTITIES_FOUND;
00801                 senses[idx] = SENSE_BOTH;  // allow double senses
00802                 // do not need to add a new sense, but still need to reset the tag
00803                 // because of a new value
00804                 append = false;
00805             }
00806         }
00807         if( append )
00808         {
00809             // what happens if a var tag data was already set before, and now it is
00810             // reset with a different size??
00811             higher_ents.push_back( wrt_entity );
00812             senses.push_back( sense );
00813         }
00814         // finally, set the senses :
00815         int dum_size  = higher_ents.size();
00816         void* dum_ptr = &higher_ents[0];
00817         rval          = mdbImpl->tag_set_by_ptr( senseNEntsTag, &entity, 1, &dum_ptr, &dum_size );MB_CHK_SET_ERR( rval, "Failed to set the sense data" );
00818 
00819         dum_ptr  = &senses[0];
00820         dum_size = higher_ents.size();
00821         rval     = mdbImpl->tag_set_by_ptr( senseNSensesTag, &entity, 1, &dum_ptr, &dum_size );MB_CHK_SET_ERR( rval, "Failed to set the sense data by pointer" );
00822     }
00823     else
00824     {
00825         // this case is about a face in the volume
00826         // there could be only 2 volumes
00827 
00828         rval = check_face_sense_tag( true );MB_CHK_SET_ERR( rval, "Failed to verify the face sense tag" );
00829 
00830         EntityHandle sense_data[2] = { 0, 0 };
00831         rval                       = mdbImpl->tag_get_data( sense2Tag, &entity, 1, sense_data );
00832         if( MB_TAG_NOT_FOUND != rval && MB_SUCCESS != rval ) MB_CHK_SET_ERR( rval, "Failed to get the sense2Tag data" );
00833 
00834         if( 0 == sense )
00835         {
00836             if( 0 != sense_data[0] && wrt_entity != sense_data[0] ) return MB_MULTIPLE_ENTITIES_FOUND;
00837             if( 0 != sense_data[1] && wrt_entity != sense_data[1] ) return MB_MULTIPLE_ENTITIES_FOUND;
00838             sense_data[0] = sense_data[1] = wrt_entity;
00839         }
00840         else if( -1 == sense )
00841         {
00842             if( 0 != sense_data[1] && wrt_entity != sense_data[1] ) return MB_MULTIPLE_ENTITIES_FOUND;
00843             if( sense_data[1] == wrt_entity ) return MB_SUCCESS;  // already set as we want
00844             sense_data[1] = wrt_entity;
00845         }
00846         else if( 1 == sense )
00847         {
00848             if( 0 != sense_data[0] && wrt_entity != sense_data[0] ) return MB_MULTIPLE_ENTITIES_FOUND;
00849             if( sense_data[0] == wrt_entity ) return MB_SUCCESS;  // already set as we want
00850             sense_data[0] = wrt_entity;
00851         }
00852         return mdbImpl->tag_set_data( sense2Tag, &entity, 1, sense_data );
00853     }
00854     return MB_SUCCESS;
00855 }
00856 
00857 //! Get the sense of entity with respect to wrt_entity
00858 //! Returns MB_ENTITY_NOT_FOUND if no relationship found
00859 ErrorCode GeomTopoTool::get_sense( EntityHandle entity, EntityHandle wrt_entity, int& sense )
00860 {
00861     // entity is lower dim (edge or face), wrt_entity is face or volume
00862     int edim   = dimension( entity );
00863     int wrtdim = dimension( wrt_entity );
00864     if( -1 == edim || -1 == wrtdim ) MB_SET_ERR( MB_FAILURE, "Non-geometric entity provided" );
00865     if( wrtdim - edim != 1 ) MB_SET_ERR( MB_FAILURE, "Entity dimension mismatch" );
00866 
00867     ErrorCode rval;
00868 
00869     if( 1 == edim )
00870     {
00871         // edge in face
00872         rval = check_edge_sense_tags( false );MB_CHK_SET_ERR( rval, "Failed to check the curve to surface sense tag handles" );
00873 
00874         std::vector< EntityHandle > faces;
00875         std::vector< int > senses;
00876         rval = get_senses( entity, faces, senses );  // the tags should be defined here
00877         MB_CHK_SET_ERR( rval, "Failed to get the curve to surface sense data" );
00878 
00879         std::vector< EntityHandle >::iterator it = std::find( faces.begin(), faces.end(), wrt_entity );
00880         if( it == faces.end() ) return MB_ENTITY_NOT_FOUND;
00881         unsigned int index = it - faces.begin();
00882         sense              = senses[index];
00883     }
00884     else
00885     {
00886         // face in volume
00887         rval = check_face_sense_tag( false );MB_CHK_SET_ERR( rval, "Failed to check the surface to volume sense tag handle" );
00888         EntityHandle sense_data[2] = { 0, 0 };
00889         rval                       = mdbImpl->tag_get_data( sense2Tag, &entity, 1, sense_data );
00890         if( MB_TAG_NOT_FOUND != rval && MB_SUCCESS != rval )
00891             MB_CHK_SET_ERR( rval, "Failed to get the surface to volume sense data" );
00892         if( ( wrt_entity == sense_data[0] ) && ( wrt_entity == sense_data[1] ) )
00893             sense = 0;
00894         else if( wrt_entity == sense_data[0] )
00895             sense = 1;
00896         else if( wrt_entity == sense_data[1] )
00897             sense = -1;
00898         else
00899             return MB_ENTITY_NOT_FOUND;
00900     }
00901     return MB_SUCCESS;
00902 }
00903 
00904 ErrorCode GeomTopoTool::get_surface_senses( EntityHandle surface_ent, EntityHandle& forward_vol,
00905                                             EntityHandle& reverse_vol )
00906 {
00907     ErrorCode rval;
00908     // this method should only be called to retrieve surface to volume
00909     // sense relationships
00910     int ent_dim = dimension( surface_ent );
00911     // verify the incoming entity dimensions for this call
00912     if( ent_dim != 2 ) { MB_SET_ERR( MB_FAILURE, "Entity dimension is incorrect for surface meshset" ); }
00913 
00914     // get the sense information for this surface
00915     EntityHandle parent_vols[2] = { 0, 0 };
00916     rval                        = mdbImpl->tag_get_data( sense2Tag, &surface_ent, 1, parent_vols );MB_CHK_SET_ERR( rval, "Failed to get surface sense data" );
00917 
00918     // set the outgoing values
00919     forward_vol = parent_vols[0];
00920     reverse_vol = parent_vols[1];
00921 
00922     return MB_SUCCESS;
00923 }
00924 
00925 ErrorCode GeomTopoTool::set_surface_senses( EntityHandle surface_ent, EntityHandle forward_vol,
00926                                             EntityHandle reverse_vol )
00927 {
00928     ErrorCode rval;
00929     // this mthod should only be called to retrieve surface to volume
00930     // sense relationships
00931     int ent_dim = dimension( surface_ent );
00932     // verify the incoming entity dimensions for this call
00933     if( ent_dim != 2 ) { MB_SET_ERR( MB_FAILURE, "Entity dimension is incorrect for surface meshset" ); }
00934 
00935     // set the sense information for this surface
00936     EntityHandle parent_vols[2] = { forward_vol, reverse_vol };
00937     rval                        = mdbImpl->tag_set_data( sense2Tag, &surface_ent, 1, parent_vols );MB_CHK_SET_ERR( rval, "Failed to set surface sense data" );
00938 
00939     return MB_SUCCESS;
00940 }
00941 
00942 // get sense of surface(s) wrt volume
00943 ErrorCode GeomTopoTool::get_surface_senses( EntityHandle volume, int num_surfaces, const EntityHandle* surfaces,
00944                                             int* senses_out )
00945 {
00946 
00947     /* The sense tags do not reference the implicit complement handle.
00948        All surfaces that interact with the implicit complement should have
00949        a null handle in the direction of the implicit complement. */
00950     // if (volume == impl_compl_handle)
00951     //  volume = (EntityHandle) 0;
00952 
00953     for( int surf_num = 0; surf_num < num_surfaces; surf_num++ )
00954     {
00955         get_sense( surfaces[surf_num], volume, senses_out[surf_num] );
00956     }
00957 
00958     return MB_SUCCESS;
00959 }
00960 
00961 ErrorCode GeomTopoTool::get_senses( EntityHandle entity, std::vector< EntityHandle >& wrt_entities,
00962                                     std::vector< int >& senses )
00963 {
00964     // the question here is: the wrt_entities is supplied or not?
00965     // I assume not, we will obtain it !!
00966     int edim = dimension( entity );
00967 
00968     if( -1 == edim ) MB_SET_ERR( MB_FAILURE, "Non-geometric entity provided" );
00969 
00970     ErrorCode rval;
00971     wrt_entities.clear();
00972     senses.clear();
00973 
00974     if( 1 == edim )  // edge
00975     {
00976         rval = check_edge_sense_tags( false );MB_CHK_SET_ERR( rval, "Failed to check the curve to surface sense tag handles" );
00977         const void* dum_ptr;
00978         int num_ents;
00979         rval = mdbImpl->tag_get_by_ptr( senseNEntsTag, &entity, 1, &dum_ptr, &num_ents );MB_CHK_ERR( rval );
00980 
00981         const EntityHandle* ents_data = static_cast< const EntityHandle* >( dum_ptr );
00982         std::copy( ents_data, ents_data + num_ents, std::back_inserter( wrt_entities ) );
00983 
00984         rval = mdbImpl->tag_get_by_ptr( senseNSensesTag, &entity, 1, &dum_ptr, &num_ents );MB_CHK_ERR( rval );
00985 
00986         const int* senses_data = static_cast< const int* >( dum_ptr );
00987         std::copy( senses_data, senses_data + num_ents, std::back_inserter( senses ) );
00988     }
00989     else  // face in volume, edim == 2
00990     {
00991         rval = check_face_sense_tag( false );MB_CHK_SET_ERR( rval, "Failed to check the surface to volume sense tag handle" );
00992         EntityHandle sense_data[2] = { 0, 0 };
00993         rval                       = mdbImpl->tag_get_data( sense2Tag, &entity, 1, sense_data );MB_CHK_SET_ERR( rval, "Failed to get the surface to volume sense data" );
00994         if( sense_data[0] != 0 && sense_data[1] == sense_data[0] )
00995         {
00996             wrt_entities.push_back( sense_data[0] );
00997             senses.push_back( 0 );  // both
00998         }
00999         else
01000         {
01001             if( sense_data[0] != 0 )
01002             {
01003                 wrt_entities.push_back( sense_data[0] );
01004                 senses.push_back( 1 );
01005             }
01006             if( sense_data[1] != 0 )
01007             {
01008                 wrt_entities.push_back( sense_data[1] );
01009                 senses.push_back( -1 );
01010             }
01011         }
01012     }
01013     // filter the results with the sets that are in the model at this time
01014     // this was introduced because extracting some sets (e.g. neumann set, with mbconvert)
01015     //   from a model would leave some sense tags not defined correctly
01016     // also, the geom ent set really needs to be part of the current model set
01017     unsigned int currentSize = 0;
01018 
01019     for( unsigned int index = 0; index < wrt_entities.size(); index++ )
01020     {
01021         EntityHandle wrt_ent = wrt_entities[index];
01022         if( wrt_ent )
01023         {
01024             if( mdbImpl->contains_entities( modelSet, &wrt_ent, 1 ) )
01025             {
01026                 wrt_entities[currentSize] = wrt_entities[index];
01027                 senses[currentSize]       = senses[index];
01028                 currentSize++;
01029             }
01030         }
01031     }
01032     wrt_entities.resize( currentSize );
01033     senses.resize( currentSize );
01034     //
01035     return MB_SUCCESS;
01036 }
01037 
01038 ErrorCode GeomTopoTool::set_senses( EntityHandle entity, std::vector< EntityHandle >& wrt_entities,
01039                                     std::vector< int >& senses )
01040 {
01041     // not efficient, and maybe wrong
01042     for( unsigned int i = 0; i < wrt_entities.size(); i++ )
01043     {
01044         ErrorCode rval = set_sense( entity, wrt_entities[i], senses[i] );MB_CHK_SET_ERR( rval, "Failed to set the sense" );
01045     }
01046 
01047     return MB_SUCCESS;
01048 }
01049 
01050 ErrorCode GeomTopoTool::next_vol( EntityHandle surface, EntityHandle old_volume, EntityHandle& new_volume )
01051 {
01052     std::vector< EntityHandle > parents;
01053     ErrorCode rval = mdbImpl->get_parent_meshsets( surface, parents );
01054 
01055     if( MB_SUCCESS == rval )
01056     {
01057         if( parents.size() != 2 )
01058             rval = MB_FAILURE;
01059         else if( parents.front() == old_volume )
01060             new_volume = parents.back();
01061         else if( parents.back() == old_volume )
01062             new_volume = parents.front();
01063         else
01064             rval = MB_FAILURE;
01065     }
01066 
01067     if( rval != MB_SUCCESS )
01068     {
01069         std::cerr << "mesh error in next_vol for surf " << surface
01070                   << std::endl;  // todo: use geomtopotool to get id by entity handle
01071     }
01072 
01073     return rval;
01074 }
01075 
01076 ErrorCode GeomTopoTool::check_geom_tag( bool create )
01077 {
01078     ErrorCode rval;
01079     unsigned flags = create ? MB_TAG_DENSE | MB_TAG_CREAT : MB_TAG_DENSE;
01080     if( !geomTag )
01081     {
01082         // get any kind of tag that already exists
01083         rval = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag, flags );MB_CHK_SET_ERR( rval, "Could not get/create the geometry dimension tag" );
01084     }
01085     return MB_SUCCESS;
01086 }
01087 
01088 ErrorCode GeomTopoTool::check_gid_tag( bool create )
01089 {
01090     ErrorCode rval;
01091     unsigned flags = create ? MB_TAG_DENSE | MB_TAG_CREAT : MB_TAG_DENSE;
01092     if( !gidTag )
01093     {
01094         // get any kind of tag that already exists
01095         rval = mdbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gidTag, flags );MB_CHK_SET_ERR( rval, "Could not get/create the global id tag" );
01096     }
01097     return MB_SUCCESS;
01098 }
01099 
01100 // move the sense tag existence creation in private methods
01101 // verify sense face tag
01102 ErrorCode GeomTopoTool::check_face_sense_tag( bool create )
01103 {
01104     ErrorCode rval;
01105     unsigned flags = create ? MB_TAG_SPARSE | MB_TAG_CREAT | MB_TAG_ANY : MB_TAG_SPARSE | MB_TAG_ANY;
01106     if( !sense2Tag )
01107     {
01108         EntityHandle def_val[2] = { 0, 0 };
01109         rval = mdbImpl->tag_get_handle( GEOM_SENSE_2_TAG_NAME, 2, MB_TYPE_HANDLE, sense2Tag, flags, def_val );MB_CHK_SET_ERR( rval, "Could not get/create the sense2Tag" );
01110     }
01111     return MB_SUCCESS;
01112 }
01113 
01114 // verify sense edge tags
01115 ErrorCode GeomTopoTool::check_edge_sense_tags( bool create )
01116 {
01117     ErrorCode rval;
01118     unsigned flags = MB_TAG_VARLEN | MB_TAG_SPARSE;
01119     if( create ) flags |= MB_TAG_CREAT;
01120     if( !senseNEntsTag )
01121     {
01122         rval = mdbImpl->tag_get_handle( GEOM_SENSE_N_ENTS_TAG_NAME, 0, MB_TYPE_HANDLE, senseNEntsTag, flags );MB_CHK_SET_ERR( rval, "Failed to get the curve to surface entity tag handle" );
01123         rval = mdbImpl->tag_get_handle( GEOM_SENSE_N_SENSES_TAG_NAME, 0, MB_TYPE_INTEGER, senseNSensesTag, flags );MB_CHK_SET_ERR( rval, "Failed to get the curve to surface sense tag handle" );
01124     }
01125     return MB_SUCCESS;
01126 }
01127 
01128 ErrorCode GeomTopoTool::add_geo_set( EntityHandle set, int dim, int gid )
01129 {
01130     if( dim < 0 || dim > 4 ) MB_SET_ERR( MB_FAILURE, "Invalid geometric dimension provided" );
01131 
01132     // see if it is not already set
01133     if( geomRanges[dim].find( set ) != geomRanges[dim].end() )
01134     {
01135         return MB_SUCCESS;  // nothing to do, we already have it as a geo set of proper dimension
01136     }
01137     updated = false;  // this should trigger at least an obb recomputation
01138     // get the geom topology tag
01139     ErrorCode result;
01140     if( 0 == geomTag )
01141     {
01142         result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );MB_CHK_SET_ERR( result, "Failed to get the geometry dimension tag handle" );
01143     }
01144 
01145     if( 0 == gidTag ) { gidTag = mdbImpl->globalId_tag(); }
01146 
01147     // make sure the added set has the geom tag properly set
01148     result = mdbImpl->tag_set_data( geomTag, &set, 1, &dim );MB_CHK_SET_ERR( result, "Failed set the geometry dimension tag value" );
01149 
01150     geomRanges[dim].insert( set );
01151     // not only that, but also add it to the root model set
01152     if( modelSet )
01153     {
01154         result = mdbImpl->add_entities( modelSet, &set, 1 );MB_CHK_SET_ERR( result, "Failed to add new geometry set to the tool's modelSet" );
01155     }
01156 
01157     // set the global ID value
01158     // if passed 0, just increase the max id for the dimension
01159     if( 0 == gid ) { gid = ++maxGlobalId[dim]; }
01160 
01161     result = mdbImpl->tag_set_data( gidTag, &set, 1, &gid );MB_CHK_SET_ERR( result, "Failed to get the global id tag value for the geom entity" );
01162 
01163     return MB_SUCCESS;
01164 }
01165 
01166 // will assume no geo sets are defined for this surface
01167 // will output a mesh_set that contains everything (all sets of interest), for proper output
01168 ErrorCode GeomTopoTool::geometrize_surface_set( EntityHandle surface, EntityHandle& output )
01169 {
01170     // usual scenario is to read a surface smf file, and "geometrize" it, and output it as a
01171     //  h5m file with proper sets, tags defined for mesh-based geometry
01172 
01173     // get all triangles/quads from the surface, then build loops
01174     // we may even have to
01175     // proper care has to be given to the orientation, material to the left!!!
01176     // at some point we may have to reorient triangles, not only edges, for proper definition
01177     bool debugFlag = false;
01178 
01179     Range surface_ents, edge_ents, loop_range;
01180 
01181     // most of these should be triangles and quads
01182     ErrorCode rval = mdbImpl->get_entities_by_dimension( surface, 2, surface_ents );MB_CHK_SET_ERR( rval, "Failed to get the surface entities" );
01183 
01184     EntityHandle face = surface;
01185     if( !surface )  // in the case it is root set, create another set
01186     {
01187         rval = mdbImpl->create_meshset( MESHSET_SET, face );MB_CHK_SET_ERR( rval, "Failed to create a the new surface meshset" );
01188     }
01189     // set the geo tag
01190     rval = add_geo_set( face, 2 );MB_CHK_SET_ERR( rval, "Failed to add the geometry set to the tool" );
01191 
01192     // this will be our output set, will contain all our new geo sets
01193     rval = mdbImpl->create_meshset( MESHSET_SET, output );MB_CHK_SET_ERR( rval, "Failed to create the output meshset" );
01194 
01195     // add first geo set (face) to the output set
01196     rval = mdbImpl->add_entities( output, &face, 1 );MB_CHK_SET_ERR( rval, "Failed to add the new meshset to the output meshset" );
01197 
01198     // how many edges do we need to create?
01199     // depends on how many loops we have
01200     // also, we should avoid non-manifold topology
01201     if( !surface )
01202     {  // in this case, surface is root, so we have to add entities
01203         rval = mdbImpl->add_entities( face, surface_ents );MB_CHK_SET_ERR( rval, "Failed to add surface entities to the surface meshset" );
01204     }
01205 
01206     Skinner tool( mdbImpl );
01207     rval = tool.find_skin( 0, surface_ents, 1, edge_ents );MB_CHK_SET_ERR( rval, "Failed to skin the surface entities" );
01208     if( debugFlag )
01209     {
01210         std::cout << "skinning edges: " << edge_ents.size() << "\n";
01211         for( Range::iterator it = edge_ents.begin(); it != edge_ents.end(); ++it )
01212         {
01213             EntityHandle ed = *it;
01214             std::cout << "edge: " << mdbImpl->id_from_handle( ed ) << " type:" << mdbImpl->type_from_handle( ed )
01215                       << "\n";
01216             std::cout << mdbImpl->list_entity( ed );
01217         }
01218     }
01219 
01220     std::vector< EntityHandle > edges_loop;
01221 
01222     Range pool_of_edges = edge_ents;
01223     Range used_edges;  // these edges are already used for some loops
01224     // get a new one
01225 
01226     while( !pool_of_edges.empty() )
01227     {
01228         // get the first edge, and start a loop with it
01229         EntityHandle current_edge = pool_of_edges[0];
01230         if( debugFlag )
01231         {
01232             std::cout << "Start current edge: " << mdbImpl->id_from_handle( current_edge ) << "\n ";
01233             std::cout << mdbImpl->list_entity( current_edge );
01234         }
01235         // get its triangle / quad and see its orientation
01236         std::vector< EntityHandle > tris;
01237         rval = mdbImpl->get_adjacencies( &current_edge, 1, 2, false, tris );MB_CHK_SET_ERR( rval, "Failed to get the adjacent triangles to the current edge" );
01238         if( tris.size() != 1 ) MB_SET_ERR( MB_FAILURE, "Edge not on boundary" );
01239 
01240         int side_n, sense, offset;
01241         rval = mdbImpl->side_number( tris[0], current_edge, side_n, sense, offset );MB_CHK_SET_ERR( rval, "Failed to get the current edge's side number" );
01242 
01243         const EntityHandle* conn2;
01244         int nnodes2;
01245         rval = mdbImpl->get_connectivity( current_edge, conn2, nnodes2 );MB_CHK_SET_ERR( rval, "Failed to get the current edge's connectivity" );
01246 
01247         if( nnodes2 != 2 ) MB_SET_ERR( MB_FAILURE, "Incorrect number of nodes found." );
01248 
01249         EntityHandle start_node = conn2[0];
01250         EntityHandle next_node  = conn2[1];
01251 
01252         if( sense == -1 )
01253         {
01254             // revert the edge, and start well
01255             EntityHandle nn2[2] = { conn2[1], conn2[0] };
01256             rval                = mdbImpl->set_connectivity( current_edge, nn2, 2 );MB_CHK_SET_ERR( rval, "Failed to set the connectivity of the current edge" );
01257 
01258             start_node = nn2[0];  // or conn2[0] !!! beware: conn2 is modified
01259             next_node  = nn2[1];  // or conn2[1]   !!!
01260             // reset connectivity of edge
01261             if( debugFlag ) std::cout << " current edge needs reversed\n";
01262         }
01263         // start a new loop of edges
01264         edges_loop.clear();  // every edge loop starts fresh
01265         edges_loop.push_back( current_edge );
01266         used_edges.insert( current_edge );
01267         pool_of_edges.erase( current_edge );
01268 
01269         if( debugFlag )
01270         {
01271             std::cout << " start node: " << start_node << "\n";
01272             std::cout << mdbImpl->list_entity( start_node );
01273             std::cout << " next node: " << next_node << "\n";
01274             std::cout << mdbImpl->list_entity( next_node );
01275         }
01276         while( next_node != start_node )
01277         {
01278             // find the next edge in the skin
01279             std::vector< EntityHandle > candidate_edges;
01280             rval = mdbImpl->get_adjacencies( &next_node, 1, 1, false, candidate_edges );MB_CHK_SET_ERR( rval, "Failed to get the adjacent edges to the next node" );
01281             // filter the edges that are used, or the edges not in the skin
01282             std::vector< EntityHandle > good_edges;
01283             for( int k = 0; k < (int)candidate_edges.size(); k++ )
01284             {
01285                 EntityHandle edg = candidate_edges[k];
01286                 if( used_edges.find( edg ) != used_edges.end() ) continue;
01287                 if( pool_of_edges.find( edg ) != pool_of_edges.end() ) good_edges.push_back( edg );
01288             }
01289             if( good_edges.size() != 1 )
01290             {
01291                 std::cout << " good_edges.size()=" << good_edges.size() << " STOP\n";
01292                 MB_SET_ERR( MB_FAILURE, "Number of good edges is not one. Could not complete the loop" );
01293             }
01294             // see if the orientation is good; if not, revert it
01295 
01296             current_edge = good_edges[0];
01297             rval         = mdbImpl->get_connectivity( current_edge, conn2, nnodes2 );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the current edge" );
01298             if( nnodes2 != 2 ) MB_SET_ERR( MB_FAILURE, "Incorrect number of nodes found" );
01299 
01300             if( conn2[0] != next_node )
01301             {
01302                 if( conn2[1] != next_node )
01303                 {
01304                     // the edge is not connected then to current edge
01305                     // bail out
01306                     std::cout << "edge " << mdbImpl->id_from_handle( current_edge ) << " not connected to node "
01307                               << next_node << "\n";
01308                     MB_SET_ERR( MB_FAILURE, "Current edge is not connected to node" );
01309                     ;
01310                 }
01311                 if( debugFlag )
01312                 {
01313                     std::cout << " revert edge " << mdbImpl->id_from_handle( current_edge ) << "\n";
01314                     std::cout << mdbImpl->list_entity( current_edge );
01315                 }
01316                 // orientation should be reversed
01317                 EntityHandle nn2[2] = { conn2[1], conn2[0] };
01318                 rval                = mdbImpl->set_connectivity( current_edge, nn2, 2 );MB_CHK_SET_ERR( rval, "Failed to set the connectivity of the current edge" );
01319 
01320                 {
01321                     std::cout << "after revert edge " << mdbImpl->id_from_handle( current_edge ) << "\n";
01322                     std::cout << mdbImpl->list_entity( current_edge );
01323                     std::cout << " conn2: " << conn2[0] << " " << conn2[1] << "\n";
01324                 }
01325             }
01326             // before reversion, conn2 was something { n1, next_node}
01327             // after reversion, conn2 became {next_node, n1}, so the
01328             // new next node will be still conn2[1]; big surprise, as
01329             //  I didn't expect the conn2 to change.
01330             // it seems that const EntityHandle * conn2 means conn2 cannot be
01331             // changed, but what is pointed to by it will change when we reset connectivity for edge
01332             next_node = conn2[1];
01333 
01334             if( debugFlag )
01335             {
01336                 std::cout << " current edge: " << mdbImpl->id_from_handle( current_edge ) << "\n ";
01337                 std::cout << mdbImpl->list_entity( current_edge );
01338                 std::cout << "next node: " << next_node << "\n ";
01339                 std::cout << mdbImpl->list_entity( next_node );
01340             }
01341             edges_loop.push_back( current_edge );
01342             used_edges.insert( current_edge );
01343             pool_of_edges.erase( current_edge );
01344         }
01345         // at this point, we have a loop formed;
01346         // create a geo edge, a vertex set, and add it to our sets
01347 
01348         EntityHandle edge;
01349         rval = mdbImpl->create_meshset( MESHSET_ORDERED, edge );MB_CHK_SET_ERR( rval, "Failed to create the edge meshset" );
01350 
01351         rval = add_geo_set( edge, 1 );MB_CHK_SET_ERR( rval, "Failed to add the edge meshset to the tool's model set" );
01352         // add the mesh edges:
01353         // add loops edges to the edge set
01354         rval = mdbImpl->add_entities( edge, &edges_loop[0], edges_loop.size() );  //
01355         MB_CHK_SET_ERR( rval, "Failed to add entities to the edge meshset" );
01356         // create a vertex set
01357         EntityHandle vertex;
01358         rval = mdbImpl->create_meshset( MESHSET_SET, vertex );MB_CHK_SET_ERR( rval, "Failed to create the vertex meshset" );
01359         rval = add_geo_set( vertex, 0 );MB_CHK_SET_ERR( rval, "Failed to add the vertex meshset to the tool's model set" );
01360         // add one node to the vertex set
01361 
01362         rval = mdbImpl->add_entities( vertex, &start_node, 1 );  //
01363         MB_CHK_SET_ERR( rval, "Failed to add entities to the vertex meshset" );
01364 
01365         rval = mdbImpl->add_parent_child( face, edge );MB_CHK_SET_ERR( rval, "Failed to create the edge to face parent child relationship" );
01366 
01367         rval = mdbImpl->add_parent_child( edge, vertex );MB_CHK_SET_ERR( rval, "Failed to create the vertex to edge parent child relationship" );
01368 
01369         // the sense of the edge in face is for sure positive (forward)
01370         rval = set_sense( edge, face, 1 );  //
01371         MB_CHK_SET_ERR( rval, "Failed to set the edge to face sense" );
01372         // also add our sets to the output set, to be sure to be exported
01373 
01374         rval = mdbImpl->add_entities( output, &edge, 1 );MB_CHK_SET_ERR( rval, "Failed to add the edge meshset to the output set" );
01375         rval = mdbImpl->add_entities( output, &vertex, 1 );MB_CHK_SET_ERR( rval, "Failed to add the vertex meshset to the output set" );
01376 
01377         if( debugFlag )
01378         {
01379             std::cout << "add edge with start node " << start_node << " with " << edges_loop.size() << " edges\n";
01380         }
01381     }
01382 
01383     return MB_SUCCESS;
01384 }
01385 
01386 /*
01387  *  This would create a deep copy of the geom topo model, into a new geom topo tool
01388  * sets will be duplicated, but entities not
01389  * modelSet will be a new one
01390  * if the original set was null (root), a new model set will be created for
01391  * original model, and its entities will be the original g sets
01392  * Scenario: split a face along a ground line, then write only one surface
01393  *   the common line has 2 faces in var len tag for sense edge; if we write only one
01394  *   surface to a new database, the var len tag of the edge will be extracted with 2 values, but
01395  *   only one value will make sense, the other will be zero.
01396  *
01397  *   There is no workaround; we need to create a duplicate model that has only that surface
01398  *   and its children (and grand-children). Then the var len sense edge-face tags will have
01399  *   the right size.
01400  *
01401  */
01402 ErrorCode GeomTopoTool::duplicate_model( GeomTopoTool*& duplicate, std::vector< EntityHandle >* pvGEnts )
01403 {
01404     // will
01405     EntityHandle rootModelSet;
01406     ErrorCode rval = mdbImpl->create_meshset( MESHSET_SET, rootModelSet );MB_CHK_SET_ERR( rval, "Failed to create the rootModelSet" );
01407 
01408     if( 0 == geomTag )
01409     {
01410         rval = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );MB_CHK_SET_ERR( rval, "Failed to get the geometry dimension tag handle" );
01411     }
01412     if( 0 == gidTag ) { gidTag = mdbImpl->globalId_tag(); }
01413     // extract from the geomSet the dimension, children, and grand-children
01414     Range depSets;  // dependents of the geomSet, including the geomSet
01415     // add in this range all the dependents of this, to filter the ones that need to be deep copied
01416 
01417     if( pvGEnts != NULL )
01418     {
01419         unsigned int numGents = pvGEnts->size();
01420         for( unsigned int k = 0; k < numGents; k++ )
01421         {
01422             EntityHandle geomSet = ( *pvGEnts )[k];
01423             // will keep accumulating to the depSets range
01424             rval = mdbImpl->get_child_meshsets( geomSet, depSets, 0 );  // 0 for numHops means that all
01425             // dependents are returned, not only the direct children.
01426             MB_CHK_SET_ERR( rval, "Failed to get the geometry set's child meshsets" );
01427 
01428             depSets.insert( geomSet );
01429         }
01430     }
01431 
01432     // add to the root model set copies of the gsets, with correct sets
01433     // keep a map between sets to help in copying parent/child relations
01434     std::map< EntityHandle, EntityHandle > relate;
01435     // each set will get the same entities as the original
01436     for( int dim = 0; dim < 5; dim++ )
01437     {
01438         int gid                  = 0;
01439         unsigned int set_options = ( ( 1 != dim ) ? MESHSET_SET : MESHSET_ORDERED );
01440         for( Range::iterator it = geomRanges[dim].begin(); it != geomRanges[dim].end(); ++it )
01441         {
01442             EntityHandle set = *it;
01443             if( pvGEnts != NULL && depSets.find( set ) == depSets.end() )
01444                 continue;  // this means that this set is not of interest, skip it
01445             EntityHandle newSet;
01446             rval = mdbImpl->create_meshset( set_options, newSet );MB_CHK_SET_ERR( rval, "Failed to create new meshset" );
01447 
01448             relate[set] = newSet;
01449             rval        = mdbImpl->add_entities( rootModelSet, &newSet, 1 );MB_CHK_SET_ERR( rval, "Failed to add the new meshset to the tool's modelSet" );
01450 
01451             // make it a geo set, and give also global id in order
01452             rval = mdbImpl->tag_set_data( geomTag, &newSet, 1, &dim );MB_CHK_SET_ERR( rval, "Failed to set the new meshset's geometry dimension data" );
01453 
01454             gid++;  // increment global id, everything starts with 1 in the new model!
01455             rval = mdbImpl->tag_set_data( gidTag, &newSet, 1, &gid );MB_CHK_SET_ERR( rval, "Failed to get the new meshset's global id data" );
01456 
01457             if( dim == 1 )
01458             {
01459                 // the entities are ordered, we need to retrieve them ordered, and set them ordered
01460                 std::vector< EntityHandle > mesh_edges;
01461                 rval = mdbImpl->get_entities_by_handle( set, mesh_edges );MB_CHK_SET_ERR( rval, "Failed to get the meshset entities by handle" );
01462 
01463                 rval = mdbImpl->add_entities( newSet, &( mesh_edges[0] ), (int)mesh_edges.size() );MB_CHK_SET_ERR( rval, "Failed to add the new entities to the new meshset" );
01464             }
01465             else
01466             {
01467                 Range ents;
01468                 rval = mdbImpl->get_entities_by_handle( set, ents );MB_CHK_SET_ERR( rval, "Failed to add the entities to the existing meshset" );
01469 
01470                 rval = mdbImpl->add_entities( newSet, ents );MB_CHK_SET_ERR( rval, "Failed to add the entities to the new meshset" );
01471             }
01472             // set parent/child relations if dim>=1
01473             if( dim >= 1 )
01474             {
01475                 Range children;
01476                 // the children of geo sets are only g sets
01477                 rval = mdbImpl->get_child_meshsets( set, children );  // num_hops = 1 by default
01478                 MB_CHK_SET_ERR( rval, "Failed to get the child meshsets of the existing set" );
01479 
01480                 for( Range::iterator it2 = children.begin(); it2 != children.end(); ++it2 )
01481                 {
01482                     EntityHandle newChildSet = relate[*it2];
01483                     rval                     = mdbImpl->add_parent_child( newSet, newChildSet );MB_CHK_SET_ERR( rval, "Failed to create parent child relationship to the new meshset" );
01484                 }
01485             }
01486         }
01487     }
01488 
01489     duplicate = new GeomTopoTool( mdbImpl, true, rootModelSet );  // will retrieve the
01490     // sets and put them in ranges
01491 
01492     // this is the lazy way to it:
01493     // newgtt->restore_topology_from_adjacency(); // will reset the sense entities, and with this,
01494     // the model represented by this new gtt will be complete set senses by peeking at the old model
01495     // make sure we have the sense tags defined
01496     rval = check_face_sense_tag( true );MB_CHK_SET_ERR( rval, "Failed to check the face to volume sense tag handle" );
01497 
01498     rval = check_edge_sense_tags( true );MB_CHK_SET_ERR( rval, "Failed to check the curve to surface sense tag handles" );
01499 
01500     for( int dd = 1; dd <= 2; dd++ )  // do it for surfaces and edges
01501     {
01502         for( Range::iterator it = geomRanges[dd].begin(); it != geomRanges[dd].end(); ++it )
01503         {
01504             EntityHandle surf = *it;
01505             if( pvGEnts != NULL && depSets.find( surf ) == depSets.end() )
01506                 continue;  // this means that this set is not of interest, skip it
01507             EntityHandle newSurf = relate[surf];
01508             // we can actually look at the tag data, to be more efficient
01509             // or use the
01510             std::vector< EntityHandle > solids;
01511             std::vector< int > senses;
01512             rval = this->get_senses( surf, solids, senses );MB_CHK_SET_ERR( rval, "Failed to get the sense data for the surface with respect to its volumes" );
01513 
01514             std::vector< EntityHandle > newSolids;
01515             std::vector< int > newSenses;
01516             for( unsigned int i = 0; i < solids.size(); i++ )
01517             {
01518                 if( pvGEnts != NULL && depSets.find( solids[i] ) == depSets.end() )
01519                     continue;  // this means that this set is not of interest, skip it
01520                 EntityHandle newSolid = relate[solids[i]];
01521                 // see which "solids" are in the new model
01522                 newSolids.push_back( newSolid );
01523                 newSenses.push_back( senses[i] );
01524             }
01525             rval = duplicate->set_senses( newSurf, newSolids, newSenses );MB_CHK_SET_ERR( rval, "Failed to set the sense data for the surface with respect to the new volumes" );
01526         }
01527     }
01528     // if the original root model set for this model is 0 (root set), then create
01529     // a new set and put all the old sets in the new model set
01530     // in this way, the existing gtt remains valid (otherwise, the modelSet would contain all the
01531     // gsets, the old ones and the new ones; the root set contains everything)
01532     if( modelSet == 0 )
01533     {
01534         rval = mdbImpl->create_meshset( MESHSET_SET, modelSet );MB_CHK_SET_ERR( rval, "Failed to create the modelSet meshset" );
01535 
01536         // add to this new set all previous sets (which are still in ranges)
01537         for( int dim = 0; dim < 5; dim++ )
01538         {
01539             rval = mdbImpl->add_entities( modelSet, geomRanges[dim] );MB_CHK_SET_ERR( rval, "Failed to add the geometric meshsets to the tool's modelSet" );
01540         }
01541     }
01542     return MB_SUCCESS;
01543 }
01544 
01545 ErrorCode GeomTopoTool::get_implicit_complement( EntityHandle& implicit_complement )
01546 {
01547     if( impl_compl_handle )
01548     {
01549         implicit_complement = impl_compl_handle;
01550         return MB_SUCCESS;
01551     }
01552     else
01553     {
01554         return MB_ENTITY_NOT_FOUND;
01555     }
01556 }
01557 
01558 ErrorCode GeomTopoTool::setup_implicit_complement()
01559 {
01560 
01561     // if the implicit complement is already setup,
01562     // we're done
01563     if( impl_compl_handle != 0 )
01564     {
01565         std::cout << "IPC already exists!" << std::endl;
01566         return MB_SUCCESS;
01567     }
01568 
01569     // if not, then query for a set with it's name
01570     Range entities;
01571     const void* const tagdata[] = { IMPLICIT_COMPLEMENT_NAME };
01572     ErrorCode rval = mdbImpl->get_entities_by_type_and_tag( modelSet, MBENTITYSET, &nameTag, tagdata, 1, entities );
01573     // query error
01574     MB_CHK_SET_ERR( rval, "Unable to query for implicit complement" );
01575 
01576     // if we found exactly one, set it as the implicit complement
01577     if( entities.size() == 1 )
01578     {
01579         impl_compl_handle = entities.front();
01580         return MB_SUCCESS;
01581     }
01582 
01583     // found too many
01584     if( entities.size() > 1 ) MB_CHK_SET_ERR( MB_MULTIPLE_ENTITIES_FOUND, "Too many implicit complement sets" );
01585 
01586     // found none
01587     if( entities.empty() )
01588     {
01589         // create implicit complement if requested
01590         rval = generate_implicit_complement( impl_compl_handle );MB_CHK_SET_ERR( rval, "Could not create implicit complement" );
01591 
01592         rval = mdbImpl->tag_set_data( nameTag, &impl_compl_handle, 1, &IMPLICIT_COMPLEMENT_NAME );MB_CHK_SET_ERR( rval, "Could not set the name tag for the implicit complement" );
01593 
01594         rval = add_geo_set( impl_compl_handle, 3 );MB_CHK_SET_ERR( rval, "Failed to add implicit complement to model" );
01595 
01596         // assign category tag - this is presumably for consistency so that the
01597         // implicit complement has all the appearance of being the same as any
01598         // other volume
01599         Tag category_tag;
01600         rval = mdbImpl->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, category_tag,
01601                                         MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Could not get the category tag" );
01602 
01603         static const char volume_category[CATEGORY_TAG_SIZE] = "Volume\0";
01604         rval = mdbImpl->tag_set_data( category_tag, &impl_compl_handle, 1, volume_category );MB_CHK_SET_ERR( rval, "Could not set the category tag for the implicit complement" );
01605 
01606         return MB_SUCCESS;
01607     }
01608 
01609     return MB_FAILURE;
01610 }
01611 
01612 ErrorCode GeomTopoTool::generate_implicit_complement( EntityHandle& implicit_complement_set )
01613 {
01614 
01615     ErrorCode rval;
01616     rval = mdbImpl->create_meshset( MESHSET_SET, implicit_complement_set );MB_CHK_SET_ERR( rval, "Failed to create mesh set for implicit complement" );
01617 
01618     // make sure the sense2Tag is set
01619     if( !sense2Tag ) { check_face_sense_tag( true ); }
01620 
01621     // get all geometric surface sets
01622     Range surfs;
01623     rval = get_gsets_by_dimension( 2, surfs );MB_CHK_SET_ERR( rval, "Could not get surface sets" );
01624 
01625     // search through all surfaces
01626     std::vector< EntityHandle > parent_vols;
01627     for( Range::iterator surf_i = surfs.begin(); surf_i != surfs.end(); ++surf_i )
01628     {
01629 
01630         parent_vols.clear();
01631         // get parents of each surface
01632         rval = mdbImpl->get_parent_meshsets( *surf_i, parent_vols );MB_CHK_SET_ERR( rval, "Failed to get volume meshsets" );
01633 
01634         // if only one parent, get the OBB root for this surface
01635         if( parent_vols.size() == 1 )
01636         {
01637 
01638             // add this surf to the topology of the implicit complement volume
01639             rval = mdbImpl->add_parent_child( implicit_complement_set, *surf_i );MB_CHK_SET_ERR( rval, "Could not add surface to implicit complement set" );
01640 
01641             // get the surface sense wrt original volume
01642             EntityHandle sense_data[2] = { 0, 0 };
01643             rval                       = get_surface_senses( *surf_i, sense_data[0], sense_data[1] );MB_CHK_SET_ERR( rval, "Could not get surface sense data" );
01644 
01645             // set the surface sense wrt implicit complement volume
01646             if( 0 == sense_data[0] && 0 == sense_data[1] )
01647                 MB_SET_ERR( MB_FAILURE, "No sense data for current surface" );
01648             if( 0 == sense_data[0] )
01649                 sense_data[0] = implicit_complement_set;
01650             else if( 0 == sense_data[1] )
01651                 sense_data[1] = implicit_complement_set;
01652             else
01653                 MB_SET_ERR( MB_FAILURE, "Could not insert implicit complement into surface sense data" );
01654 
01655             // set the new sense data for this surface
01656             rval = set_surface_senses( *surf_i, sense_data[0], sense_data[1] );MB_CHK_SET_ERR( rval, "Failed to set sense tag data" );
01657         }
01658     }  // end surface loop
01659 
01660     return MB_SUCCESS;
01661 }
01662 
01663 #define RETFALSE( a, b )            \
01664     {                               \
01665         std::cout << ( a ) << "\n"; \
01666         mdbImpl->list_entity( b );  \
01667         return false;               \
01668     }
01669 bool GeomTopoTool::check_model()
01670 {
01671     // vertex sets should have one node
01672     Range::iterator rit;
01673     ErrorCode rval;
01674     for( rit = geomRanges[0].begin(); rit != geomRanges[0].end(); ++rit )
01675     {
01676         EntityHandle vSet = *rit;
01677         Range nodes;
01678         rval = mdbImpl->get_entities_by_handle( vSet, nodes );
01679         if( MB_SUCCESS != rval ) RETFALSE( " failed to get nodes from vertex set ", vSet );
01680         if( nodes.size() != 1 ) RETFALSE( " number of nodes is different from 1 ", vSet )
01681         EntityType type = mdbImpl->type_from_handle( nodes[0] );
01682         if( type != MBVERTEX ) RETFALSE( " entity in vertex set is not a node ", nodes[0] )
01683         // get all parents, and see if they belong to geomRanges[1]
01684         Range edges;
01685         rval = mdbImpl->get_parent_meshsets( vSet, edges );
01686         if( MB_SUCCESS != rval ) RETFALSE( " can't get parent edges for a node set ", vSet )
01687         Range notEdges = subtract( edges, geomRanges[1] );
01688         if( !notEdges.empty() ) RETFALSE( " some parents of a node set are not geo edges ", notEdges[0] )
01689     }
01690 
01691     // edges to be formed by continuous chain of mesh edges, oriented correctly
01692     for( rit = geomRanges[1].begin(); rit != geomRanges[1].end(); ++rit )
01693     {
01694         EntityHandle edge = *rit;
01695         std::vector< EntityHandle > mesh_edges;
01696         rval = mdbImpl->get_entities_by_type( edge, MBEDGE, mesh_edges );
01697         if( MB_SUCCESS != rval ) RETFALSE( " can't get mesh edges from edge set", edge )
01698         int num_edges = (int)mesh_edges.size();
01699         if( num_edges == 0 ) RETFALSE( " no mesh edges in edge set ", edge )
01700         EntityHandle firstNode;
01701         EntityHandle currentNode;  // will also hold the last node in chain of edges
01702         const EntityHandle* conn2;
01703         int nnodes2;
01704         // get all parents, and see if they belong to geomRanges[1]
01705         for( int i = 0; i < num_edges; i++ )
01706         {
01707             rval = mdbImpl->get_connectivity( mesh_edges[i], conn2, nnodes2 );
01708             if( MB_SUCCESS != rval || nnodes2 != 2 ) RETFALSE( " mesh edge connectivity is wrong ", mesh_edges[i] )
01709             if( i == 0 )
01710             {
01711                 firstNode   = conn2[0];
01712                 currentNode = conn2[1];
01713             }
01714 
01715             else  // if ( (i>0) )
01716             {
01717                 // check the current node is conn[0]
01718                 if( conn2[0] != currentNode )
01719                 {
01720                     std::cout << "i=" << i << " conn2:" << conn2[0] << " " << conn2[1] << " currentNode:" << currentNode
01721                               << "\n";
01722                     mdbImpl->list_entity( mesh_edges[i] );
01723                     RETFALSE( " edges are not contiguous in edge set ", edge )
01724                 }
01725                 currentNode = conn2[1];
01726             }
01727         }
01728         // check the children of the edge set; do they have the first and last node?
01729         Range vertSets;
01730         rval = mdbImpl->get_child_meshsets( edge, vertSets );
01731         if( MB_SUCCESS != rval ) RETFALSE( " can't get vertex children ", edge )
01732         Range notVertices = subtract( vertSets, geomRanges[0] );
01733         if( !notVertices.empty() ) RETFALSE( " children sets that are not vertices ", notVertices[0] )
01734         for( Range::iterator it = vertSets.begin(); it != vertSets.end(); ++it )
01735         {
01736             if( !mdbImpl->contains_entities( *it, &firstNode, 1 ) &&
01737                 !mdbImpl->contains_entities( *it, &currentNode, 1 ) )
01738                 RETFALSE( " a vertex set is not containing the first and last nodes ", *it )
01739         }
01740         // check now the faces / parents
01741         Range faceSets;
01742         rval = mdbImpl->get_parent_meshsets( edge, faceSets );
01743         if( MB_SUCCESS != rval ) RETFALSE( " can't get edge parents ", edge )
01744         Range notFaces = subtract( faceSets, geomRanges[2] );
01745         if( !notFaces.empty() ) RETFALSE( " parent sets that are not faces ", notFaces[0] )
01746 
01747         // for a geo edge, check the sense tags with respect to the adjacent faces
01748         // in general, it is sufficient to check one mesh edge (the first one)
01749         // edge/face  senses
01750         EntityHandle firstMeshEdge = mesh_edges[0];
01751         // get all entities/elements adjacent to it
01752         Range adjElem;
01753         rval = mdbImpl->get_adjacencies( &firstMeshEdge, 1, 2, false, adjElem );
01754         if( MB_SUCCESS != rval ) RETFALSE( " can't get adjacent elements to the edge ", firstMeshEdge )
01755         for( Range::iterator it2 = adjElem.begin(); it2 != adjElem.end(); ++it2 )
01756         {
01757             EntityHandle elem = *it2;
01758             // find the geo face it belongs to
01759             EntityHandle gFace = 0;
01760             for( Range::iterator fit = faceSets.begin(); fit != faceSets.end(); ++fit )
01761             {
01762                 EntityHandle possibleFace = *fit;
01763                 if( mdbImpl->contains_entities( possibleFace, &elem, 1 ) )
01764                 {
01765                     gFace = possibleFace;
01766                     break;
01767                 }
01768             }
01769             if( 0 == gFace )
01770                 RETFALSE( " can't find adjacent surface that contains the adjacent element to the edge ",
01771                           firstMeshEdge )
01772 
01773             // now, check the sense of mesh_edge in element, and the sense of gedge in gface
01774             // side_number
01775             int side_n, sense, offset;
01776             rval = mdbImpl->side_number( elem, firstMeshEdge, side_n, sense, offset );
01777             if( MB_SUCCESS != rval ) RETFALSE( " can't get sense and side number of an element ", elem )
01778             // now get the sense
01779             int topoSense;
01780             rval = this->get_sense( edge, gFace, topoSense );
01781             if( topoSense != sense ) RETFALSE( " geometric topo sense and element sense do not agree ", edge )
01782         }
01783     }
01784     // surfaces to be true meshes
01785     // for surfaces, check that the skinner will find the correct boundary
01786 
01787     // use the skinner for boundary check
01788     Skinner tool( mdbImpl );
01789 
01790     for( rit = geomRanges[2].begin(); rit != geomRanges[2].end(); ++rit )
01791     {
01792         EntityHandle faceSet = *rit;
01793         // get all boundary edges (adjacent edges)
01794 
01795         Range edges;
01796         rval = mdbImpl->get_child_meshsets( faceSet, edges );
01797         if( MB_SUCCESS != rval ) RETFALSE( " can't get children edges for a face set ", faceSet )
01798         Range notEdges = subtract( edges, geomRanges[1] );
01799         if( !notEdges.empty() ) RETFALSE( " some children of a face set are not geo edges ", notEdges[0] )
01800 
01801         Range boundary_mesh_edges;
01802         for( Range::iterator it = edges.begin(); it != edges.end(); ++it )
01803         {
01804             rval = mdbImpl->get_entities_by_type( *it, MBEDGE, boundary_mesh_edges );
01805             if( MB_SUCCESS != rval ) RETFALSE( " can't get edge elements from the edge set ", *it )
01806         }
01807         // skin the elements of the surface
01808         // most of these should be triangles and quads
01809         Range surface_ents, edge_ents;
01810         rval = mdbImpl->get_entities_by_dimension( faceSet, 2, surface_ents );
01811         if( MB_SUCCESS != rval ) RETFALSE( " can't get surface elements from the face set ", faceSet )
01812 
01813         rval = tool.find_skin( 0, surface_ents, 1, edge_ents );
01814         if( MB_SUCCESS != rval ) RETFALSE( "can't skin a surface ", surface_ents[0] )
01815 
01816         // those 2 ranges for boundary edges now must be equal
01817         if( boundary_mesh_edges != edge_ents ) RETFALSE( "boundary ranges are different", boundary_mesh_edges[0] )
01818     }
01819 
01820     // solids to be filled correctly, maybe a skin procedure too.
01821     // (maybe the solids are empty)
01822 
01823     return true;
01824 }
01825 
01826 bool GeomTopoTool::have_obb_tree()
01827 {
01828     return rootSets.size() != 0 || mapRootSets.size() != 0;
01829 }
01830 
01831 // This function gets coordinates of the minimum and maxmiumum points
01832 // from an OBB/AABB, ie. such that these points represent
01833 // the maximum and minium extents of an AABB
01834 ErrorCode GeomTopoTool::get_bounding_coords( EntityHandle volume, double minPt[3], double maxPt[3] )
01835 {
01836     double center[3], axis1[3], axis2[3], axis3[3];
01837 
01838     // get center point and vectors to OBB faces
01839     ErrorCode rval = get_obb( volume, center, axis1, axis2, axis3 );MB_CHK_SET_ERR( rval, "Failed to get the oriented bounding box of the volume" );
01840 
01841     // compute min and max vertices
01842     for( int i = 0; i < 3; i++ )
01843     {
01844         double sum = fabs( axis1[i] ) + fabs( axis2[i] ) + fabs( axis3[i] );
01845         minPt[i]   = center[i] - sum;
01846         maxPt[i]   = center[i] + sum;
01847     }
01848     return MB_SUCCESS;
01849 }
01850 
01851 ErrorCode GeomTopoTool::get_obb( EntityHandle volume, double center[3], double axis1[3], double axis2[3],
01852                                  double axis3[3] )
01853 {
01854     // find EntityHandle node_set for use in box
01855     EntityHandle root;
01856     ErrorCode rval = get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get volume's obb tree root" );
01857 
01858     // call box to get center and vectors to faces
01859     return obbTree->box( root, center, axis1, axis2, axis3 );
01860 }
01861 
01862 Range GeomTopoTool::get_ct_children_by_dimension( EntityHandle parent, int desired_dimension )
01863 {
01864     Range all_children, desired_children;
01865     Range::iterator it;
01866     int actual_dimension;
01867 
01868     desired_children.clear();
01869     all_children.clear();
01870     mdbImpl->get_child_meshsets( parent, all_children );
01871 
01872     for( it = all_children.begin(); it != all_children.end(); ++it )
01873     {
01874         mdbImpl->tag_get_data( geomTag, &( *it ), 1, &actual_dimension );
01875         if( actual_dimension == desired_dimension ) desired_children.insert( *it );
01876     }
01877 
01878     return desired_children;
01879 }
01880 
01881 // runs GeomQueryTool point_in_vol and to test if vol A is inside vol B
01882 //  returns true or false
01883 bool GeomTopoTool::A_is_in_B( EntityHandle volume_A, EntityHandle volume_B, GeomQueryTool* GQT )
01884 {
01885     ErrorCode rval;
01886 
01887     Range child_surfaces, triangles, vertices;
01888     double coord[3];  // coord[0] = x, etc.
01889     int result;       // point in vol result; 0=F, 1=T
01890 
01891     // find coordinates of any point on surface of A
01892     // get surface corresponding to volume, then get the triangles
01893     child_surfaces = get_ct_children_by_dimension( volume_A, 2 );
01894     rval           = mdbImpl->get_entities_by_type( *child_surfaces.begin(), MBTRI, triangles );MB_CHK_ERR( rval );
01895 
01896     // now get 1st triangle vertices
01897     rval = mdbImpl->get_connectivity( &( *triangles.begin() ), 1, vertices );MB_CHK_ERR( rval );
01898 
01899     // now get coordinates of first vertex
01900     rval = mdbImpl->get_coords( &( *vertices.begin() ), 1, &( coord[0] ) );MB_CHK_ERR( rval );
01901 
01902     // if point on A is inside vol B, return T; o.w. return F
01903     rval = GQT->point_in_volume( volume_B, coord, result );MB_CHK_SET_ERR( rval, "Failed to complete point in volume query." );
01904 
01905     return ( result != 0 );
01906 }
01907 
01908 ErrorCode GeomTopoTool::insert_in_tree( EntityHandle ct_root, EntityHandle volume, GeomQueryTool* GQT )
01909 {
01910     ErrorCode rval;
01911 
01912     bool inserted               = false;
01913     EntityHandle current_volume = volume;   // volume to be inserted
01914     EntityHandle tree_volume    = ct_root;  // volume already existing in the tree
01915     EntityHandle parent         = ct_root;
01916     Range child_volumes;
01917 
01918     // while not inserted in tree
01919     while( !inserted )
01920     {
01921         // if current volume is insde of tree volume-- always true if tree volume
01922         // is the root of the tree
01923         if( tree_volume == ct_root || ( tree_volume != ct_root && A_is_in_B( current_volume, tree_volume, GQT ) ) )
01924         {
01925 
01926             parent = tree_volume;
01927 
01928             // if tree_volume has children then we must test them,
01929             // (tree_volume will change)
01930             child_volumes = get_ct_children_by_dimension( tree_volume, 3 );
01931             if( child_volumes.size() > 0 ) tree_volume = child_volumes.pop_front();
01932             // otherwise current_volume is the only child of the tree volume
01933             else
01934             {
01935                 rval = mdbImpl->add_parent_child( parent, current_volume );MB_CHK_SET_ERR( rval, "Failed to add parent-child relationship." );
01936 
01937                 inserted = true;
01938             }
01939             // if current volume is not in the tree volume, the converse may be true
01940         }
01941         else
01942         {
01943             // if the tree volume is inside the current volume
01944             if( A_is_in_B( tree_volume, current_volume, GQT ) )
01945             {
01946                 // reverse their parentage
01947                 rval = mdbImpl->remove_parent_child( parent, tree_volume );MB_CHK_SET_ERR( rval, "Failed to remove parent-child relationship." );
01948                 rval = mdbImpl->add_parent_child( current_volume, tree_volume );MB_CHK_SET_ERR( rval, "Failed to add parent-child relationship." );
01949             }
01950 
01951             if( child_volumes.size() == 0 )
01952             {
01953                 rval = mdbImpl->add_parent_child( parent, current_volume );MB_CHK_SET_ERR( rval, "Failed to add parent-child relationship." );
01954                 inserted = true;
01955             }
01956             else
01957                 tree_volume = child_volumes.pop_front();
01958         }
01959     }
01960     return MB_SUCCESS;
01961 }
01962 
01963 ErrorCode GeomTopoTool::restore_topology_from_geometric_inclusion( const Range& flat_volumes )
01964 {
01965 
01966     ErrorCode rval;
01967     // local var will go out of scope if errors appear, no need to free it also
01968     GeomQueryTool GQT( this );
01969     std::map< EntityHandle, EntityHandle > volume_surface;  // map of volume
01970                                                             // to its surface
01971 
01972     EntityHandle ct_root;
01973     // create root meshset-- this will be top of tree
01974     std::string meshset_name = "build_hierarchy_root";
01975     rval                     = mdbImpl->create_meshset( MESHSET_SET, ct_root );MB_CHK_ERR( rval );
01976     rval = mdbImpl->tag_set_data( nameTag, &ct_root, 1, meshset_name.c_str() );MB_CHK_ERR( rval );
01977 
01978     for( Range::iterator vol = flat_volumes.begin(); vol != flat_volumes.end(); vol++ )
01979     {
01980         // get the surface corresponding to each volume
01981         // at this point, each volume meshset only has one 'child' surface
01982         // which exactly corresponds to that volume
01983         Range child_surfaces = get_ct_children_by_dimension( *vol, 2 );
01984         volume_surface[*vol] = *child_surfaces.begin();
01985 
01986         rval = insert_in_tree( ct_root, *vol, &GQT );MB_CHK_SET_ERR( rval, "Failed to insert volume into tree." );
01987     }
01988 
01989     // for each original volume, get its child volumes
01990     for( Range::iterator parent_it = flat_volumes.begin(); parent_it != flat_volumes.end(); parent_it++ )
01991     {
01992         Range volume_children = get_ct_children_by_dimension( *parent_it, 3 );
01993 
01994         if( volume_children.size() != 0 )
01995         {
01996             // loop over all of original volume's child volumes
01997             for( Range::iterator child_it = volume_children.begin(); child_it != volume_children.end(); ++child_it )
01998             {
01999                 // set the sense of the surface mapped to the child volume to REVERSE
02000                 // wrt the parent volume
02001                 rval = set_sense( volume_surface[*child_it], *parent_it, SENSE_REVERSE );MB_CHK_SET_ERR( rval, "Failed to set sense." );
02002 
02003                 // add the child volume's surface as a child of the original volume
02004                 // and delete the child volume as a child of original volume
02005                 rval = mdbImpl->add_parent_child( *parent_it, volume_surface[*child_it] );MB_CHK_SET_ERR( rval, "Failed to add parent-child relationship." );
02006                 rval = mdbImpl->remove_parent_child( *parent_it, *child_it );MB_CHK_SET_ERR( rval, "Failed to remove parent-child relationship." );
02007             }
02008         }
02009     }
02010 
02011     return MB_SUCCESS;
02012 }
02013 
02014 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines