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