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