Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
MeshTopoUtil.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 #ifdef WIN32
00017 #pragma warning( disable : 4786 )
00018 #endif
00019 
00020 #include "moab/MeshTopoUtil.hpp"
00021 #include "moab/Range.hpp"
00022 #include "Internals.hpp"
00023 #include "moab/Interface.hpp"
00024 #include "moab/CN.hpp"
00025 
00026 #include <cassert>
00027 
00028 #define RR                                        \
00029     {                                             \
00030         if( MB_SUCCESS != result ) return result; \
00031     }
00032 
00033 namespace moab
00034 {
00035 
00036 //! generate all the AEntities bounding the vertices
00037 ErrorCode MeshTopoUtil::construct_aentities( const Range& vertices )
00038 {
00039     Range out_range;
00040     ErrorCode result;
00041     result = mbImpl->get_adjacencies( vertices, 1, true, out_range, Interface::UNION );
00042     if( MB_SUCCESS != result ) return result;
00043     out_range.clear();
00044     result = mbImpl->get_adjacencies( vertices, 2, true, out_range, Interface::UNION );
00045     if( MB_SUCCESS != result ) return result;
00046     out_range.clear();
00047     result = mbImpl->get_adjacencies( vertices, 3, true, out_range, Interface::UNION );
00048 
00049     return result;
00050 }
00051 
00052 //! given an entity, get its average position (avg vertex locations)
00053 ErrorCode MeshTopoUtil::get_average_position( Range& entities, double* avg_position )
00054 {
00055     std::vector< EntityHandle > ent_vec;
00056     std::copy( entities.begin(), entities.end(), std::back_inserter( ent_vec ) );
00057     return get_average_position( &ent_vec[0], ent_vec.size(), avg_position );
00058 }
00059 
00060 //! given an entity, get its average position (avg vertex locations)
00061 ErrorCode MeshTopoUtil::get_average_position( const EntityHandle* entities,
00062                                               const int num_entities,
00063                                               double* avg_position )
00064 {
00065     double dum_pos[3];
00066     avg_position[0] = avg_position[1] = avg_position[2] = 0.0;
00067 
00068     Range connect;
00069     ErrorCode result = mbImpl->get_adjacencies( entities, num_entities, 0, false, connect, Interface::UNION );
00070     if( MB_SUCCESS != result ) return result;
00071 
00072     if( connect.empty() ) return MB_FAILURE;
00073 
00074     for( Range::iterator rit = connect.begin(); rit != connect.end(); ++rit )
00075     {
00076         result = mbImpl->get_coords( &( *rit ), 1, dum_pos );
00077         if( MB_SUCCESS != result ) return result;
00078         avg_position[0] += dum_pos[0];
00079         avg_position[1] += dum_pos[1];
00080         avg_position[2] += dum_pos[2];
00081     }
00082     avg_position[0] /= (double)connect.size();
00083     avg_position[1] /= (double)connect.size();
00084     avg_position[2] /= (double)connect.size();
00085 
00086     return MB_SUCCESS;
00087 }
00088 
00089 //! given an entity, get its average position (avg vertex locations)
00090 ErrorCode MeshTopoUtil::get_average_position( const EntityHandle entity, double* avg_position )
00091 {
00092     const EntityHandle* connect = NULL;
00093     int num_connect             = 0;
00094     if( MBVERTEX == mbImpl->type_from_handle( entity ) ) return mbImpl->get_coords( &entity, 1, avg_position );
00095 
00096     ErrorCode result = mbImpl->get_connectivity( entity, connect, num_connect );
00097     if( MB_SUCCESS != result ) return result;
00098 
00099     return get_average_position( connect, num_connect, avg_position );
00100 }
00101 
00102 // given an entity, find the entities of next higher dimension around
00103 // that entity, ordered by connection through next higher dimension entities;
00104 // if any of the star entities is in only one entity of next higher dimension,
00105 // on_boundary is returned true
00106 ErrorCode MeshTopoUtil::star_entities( const EntityHandle star_center,
00107                                        std::vector< EntityHandle >& star_ents,
00108                                        bool& bdy_entity,
00109                                        const EntityHandle starting_star_entity,
00110                                        std::vector< EntityHandle >* star_entities_dp2,
00111                                        Range* star_candidates_dp2 )
00112 {
00113     // now start the traversal
00114     bdy_entity               = false;
00115     EntityHandle last_entity = starting_star_entity, last_dp2 = 0, next_entity, next_dp2;
00116     std::vector< EntityHandle > star_dp2;
00117     ErrorCode result;
00118     int center_dim = mbImpl->dimension_from_handle( star_center );
00119 
00120     Range tmp_candidates_dp2;
00121     if( NULL != star_candidates_dp2 )
00122         tmp_candidates_dp2 = *star_candidates_dp2;
00123     else
00124     {
00125         result = mbImpl->get_adjacencies( &star_center, 1, center_dim + 2, false, tmp_candidates_dp2 );
00126         if( MB_SUCCESS != result ) return result;
00127     }
00128 
00129     do
00130     {
00131         // get the next star entity
00132         result = star_next_entity( star_center, last_entity, last_dp2, &tmp_candidates_dp2, next_entity, next_dp2 );
00133         if( MB_SUCCESS != result ) return result;
00134 
00135         // special case: if starting_star_entity isn't connected to any entities of next
00136         // higher dimension, it's the only entity in the star; put it on the list and return
00137         if( star_ents.empty() && next_entity == 0 && next_dp2 == 0 )
00138         {
00139             star_ents.push_back( last_entity );
00140             bdy_entity = true;
00141             return MB_SUCCESS;
00142         }
00143 
00144         // if we're at a bdy and bdy_entity hasn't been set yet, we're at the
00145         // first bdy; reverse the lists and start traversing in the other direction; but,
00146         // pop the last star entity off the list and find it again, so that we properly
00147         // check for next_dp2
00148         if( 0 == next_dp2 && !bdy_entity )
00149         {
00150             star_ents.push_back( next_entity );
00151             bdy_entity = true;
00152             std::reverse( star_ents.begin(), star_ents.end() );
00153             star_ents.pop_back();
00154             last_entity = star_ents.back();
00155             if( !star_dp2.empty() )
00156             {
00157                 std::reverse( star_dp2.begin(), star_dp2.end() );
00158                 last_dp2 = star_dp2.back();
00159             }
00160         }
00161         // else if we're not on the bdy and next_entity is already in star, that means
00162         // we've come all the way around; don't put next_entity on list again, and
00163         // zero out last_dp2 to terminate while loop
00164         else if( !bdy_entity && std::find( star_ents.begin(), star_ents.end(), next_entity ) != star_ents.end() &&
00165                  ( std::find( star_dp2.begin(), star_dp2.end(), next_dp2 ) != star_dp2.end() || !next_dp2 ) )
00166         {
00167             last_dp2 = 0;
00168         }
00169 
00170         // else, just assign last entities seen and go on to next iteration
00171         else
00172         {
00173             if( std::find( star_ents.begin(), star_ents.end(), next_entity ) == star_ents.end() )
00174                 star_ents.push_back( next_entity );
00175             if( 0 != next_dp2 )
00176             {
00177                 star_dp2.push_back( next_dp2 );
00178                 tmp_candidates_dp2.erase( next_dp2 );
00179             }
00180             last_entity = next_entity;
00181             last_dp2    = next_dp2;
00182         }
00183     } while( 0 != last_dp2 );
00184 
00185     // copy over the star_dp2 list, if requested
00186     if( NULL != star_entities_dp2 ) ( *star_entities_dp2 ).swap( star_dp2 );
00187 
00188     return MB_SUCCESS;
00189 }
00190 
00191 ErrorCode MeshTopoUtil::star_next_entity( const EntityHandle star_center,
00192                                           const EntityHandle last_entity,
00193                                           const EntityHandle last_dp1,
00194                                           Range* star_candidates_dp1,
00195                                           EntityHandle& next_entity,
00196                                           EntityHandle& next_dp1 )
00197 {
00198     // given a star_center, a last_entity (whose dimension should be 1 greater than center)
00199     // and last_dp1 (dimension 2 higher than center), returns the next star entity across
00200     // last_dp1, and the next dp1 entity sharing next_entity; if star_candidates is non-empty,
00201     // star must come from those
00202     Range from_ents, to_ents;
00203     from_ents.insert( star_center );
00204     if( 0 != last_dp1 ) from_ents.insert( last_dp1 );
00205 
00206     int dim = mbImpl->dimension_from_handle( star_center );
00207 
00208     ErrorCode result = mbImpl->get_adjacencies( from_ents, dim + 1, true, to_ents );
00209     if( MB_SUCCESS != result ) return result;
00210 
00211     // remove last_entity from result, and should only have 1 left, if any
00212     if( 0 != last_entity ) to_ents.erase( last_entity );
00213 
00214     // if no last_dp1, contents of to_ents should share dp1-dimensional entity with last_entity
00215     if( 0 != last_entity && 0 == last_dp1 )
00216     {
00217         Range tmp_to_ents;
00218         for( Range::iterator rit = to_ents.begin(); rit != to_ents.end(); ++rit )
00219         {
00220             if( 0 != common_entity( last_entity, *rit, dim + 2 ) ) tmp_to_ents.insert( *rit );
00221         }
00222         to_ents = tmp_to_ents;
00223     }
00224 
00225     if( 0 == last_dp1 && to_ents.size() > 1 && NULL != star_candidates_dp1 && !star_candidates_dp1->empty() )
00226     {
00227         // if we have a choice of to_ents and no previous dp1 and there are dp1 candidates,
00228         // the one we choose needs to be adjacent to one of the candidates
00229         result = mbImpl->get_adjacencies( *star_candidates_dp1, dim + 1, true, from_ents, Interface::UNION );
00230         if( MB_SUCCESS != result ) return result;
00231         to_ents = intersect( to_ents, from_ents );
00232     }
00233 
00234     if( !to_ents.empty() )
00235         next_entity = *to_ents.begin();
00236     else
00237     {
00238         next_entity = 0;
00239         next_dp1    = 0;
00240         return MB_SUCCESS;
00241     }
00242 
00243     // get next_dp1
00244     if( 0 != star_candidates_dp1 )
00245         to_ents = *star_candidates_dp1;
00246     else
00247         to_ents.clear();
00248 
00249     result = mbImpl->get_adjacencies( &next_entity, 1, dim + 2, true, to_ents );
00250     if( MB_SUCCESS != result ) return result;
00251 
00252     // can't be last one
00253     if( 0 != last_dp1 ) to_ents.erase( last_dp1 );
00254 
00255     if( !to_ents.empty() ) next_dp1 = *to_ents.begin();
00256 
00257     // could be zero, means we're at bdy
00258     else
00259         next_dp1 = 0;
00260 
00261     return MB_SUCCESS;
00262 }
00263 
00264 ErrorCode MeshTopoUtil::star_entities_nonmanifold( const EntityHandle star_entity,
00265                                                    std::vector< std::vector< EntityHandle > >& stars,
00266                                                    std::vector< bool >* bdy_flags,
00267                                                    std::vector< std::vector< EntityHandle > >* dp2_stars )
00268 {
00269     // Get a series of (d+1)-dimensional stars around a d-dimensional entity, such that
00270     // each star is on a (d+2)-manifold containing the d-dimensional entity; each star
00271     // is either open or closed, and also defines a (d+2)-star whose entities are bounded by
00272     // (d+1)-entities on the star and on the (d+2)-manifold
00273     //
00274     // Algorithm:
00275     // get the (d+2)-manifold entities; for d=1 / d+2=3, just assume all connected elements, since
00276     //   we don't do 4d yet
00277     // get intersection of (d+1)-entities adjacent to star entity and union of (d+1)-entities
00278     //   adjacent to (d+2)-manifold entities; these will be the entities in the star
00279     // while (d+1)-entities
00280     //   remove (d+1)-entity from (d+1)-entities
00281     //   get the (d+1)-star and (d+2)-star around that (d+1)-entity (using star_entities)
00282     //   save that star to the star list, and the bdy flag and (d+2)-star if requested
00283     //   remove (d+2)-entities from the (d+2)-manifold entities
00284     //   remove (d+1)-entities from the (d+1)-entities
00285     // (end while)
00286 
00287     int this_dim = mbImpl->dimension_from_handle( star_entity );
00288     if( 3 <= this_dim || 0 > this_dim ) return MB_FAILURE;
00289 
00290     // get the (d+2)-manifold entities; for d=1 / d+2=3, just assume all connected elements, since
00291     //   we don't do 4d yet
00292     Range dp2_manifold;
00293     ErrorCode result = get_manifold( star_entity, this_dim + 2, dp2_manifold );
00294     if( MB_SUCCESS != result ) return result;
00295 
00296     // get intersection of (d+1)-entities adjacent to star and union of (d+1)-entities
00297     //   adjacent to (d+2)-manifold entities; also add manifold (d+1)-entities, to catch
00298     //   any not connected to (d+2)-entities
00299     Range dp1_manifold;
00300     result = mbImpl->get_adjacencies( dp2_manifold, this_dim + 1, false, dp1_manifold, Interface::UNION );
00301     if( MB_SUCCESS != result ) return result;
00302 
00303     result = mbImpl->get_adjacencies( &star_entity, 1, this_dim + 1, false, dp1_manifold );
00304     if( MB_SUCCESS != result ) return result;
00305 
00306     result = get_manifold( star_entity, this_dim + 1, dp1_manifold );
00307     if( MB_SUCCESS != result ) return result;
00308 
00309     // while (d+1)-entities
00310     while( !dp1_manifold.empty() )
00311     {
00312 
00313         //   get (d+1)-entity from (d+1)-entities (don't remove it until after star,
00314         //     since the star entities must come from dp1_manifold)
00315         EntityHandle this_ent = *dp1_manifold.begin();
00316 
00317         //   get the (d+1)-star and (d+2)-star around that (d+1)-entity (using star_entities)
00318         std::vector< EntityHandle > this_star_dp1, this_star_dp2;
00319         bool on_bdy;
00320         result = star_entities( star_entity, this_star_dp1, on_bdy, this_ent, &this_star_dp2, &dp2_manifold );
00321         if( MB_SUCCESS != result ) return result;
00322 
00323         // if there's no star entities, it must mean this_ent isn't bounded by any dp2
00324         // entities (wasn't put into star in star_entities 'cuz we're passing in non-null
00325         // dp2_manifold above); put it in
00326         if( this_star_dp1.empty() )
00327         {
00328             Range dum_range;
00329             result = mbImpl->get_adjacencies( &this_ent, 1, this_dim + 2, false, dum_range );
00330             if( MB_SUCCESS != result ) return result;
00331             if( dum_range.empty() ) this_star_dp1.push_back( this_ent );
00332         }
00333 
00334         // now we can remove it
00335         dp1_manifold.erase( dp1_manifold.begin() );
00336 
00337         //   save that star to the star list, and the bdy flag and (d+2)-star if requested
00338         if( !this_star_dp1.empty() )
00339         {
00340             stars.push_back( this_star_dp1 );
00341             if( NULL != bdy_flags ) bdy_flags->push_back( on_bdy );
00342             if( NULL != dp2_stars ) dp2_stars->push_back( this_star_dp2 );
00343         }
00344 
00345         //   remove (d+2)-entities from the (d+2)-manifold entities
00346         for( std::vector< EntityHandle >::iterator vit = this_star_dp2.begin(); vit != this_star_dp2.end(); ++vit )
00347             dp2_manifold.erase( *vit );
00348 
00349         //   remove (d+1)-entities from the (d+1)-entities
00350         for( std::vector< EntityHandle >::iterator vit = this_star_dp1.begin(); vit != this_star_dp1.end(); ++vit )
00351             dp1_manifold.erase( *vit );
00352 
00353         // (end while)
00354     }
00355 
00356     // check for leftover dp2 manifold entities, these should be in one of the
00357     // stars
00358     if( !dp2_manifold.empty() )
00359     {
00360         for( Range::iterator rit = dp2_manifold.begin(); rit != dp2_manifold.end(); ++rit )
00361         {
00362         }
00363     }
00364 
00365     return MB_SUCCESS;
00366 }
00367 
00368 //! get (target_dim)-dimensional manifold entities connected to star_entity; that is,
00369 //! the entities with <= 1 connected (target_dim+2)-dimensional adjacent entities;
00370 //! for target_dim=3, just return all of them
00371 //! just insert into the list, w/o clearing manifold list first
00372 ErrorCode MeshTopoUtil::get_manifold( const EntityHandle star_entity, const int target_dim, Range& manifold )
00373 {
00374     // get all the entities of target dimension connected to star
00375     Range tmp_range;
00376     ErrorCode result = mbImpl->get_adjacencies( &star_entity, 1, target_dim, false, tmp_range );
00377     if( MB_SUCCESS != result ) return result;
00378 
00379     // now save the ones which are (target_dim+1)-dimensional manifold;
00380     // for target_dim=3, just return whole range, since we don't do 4d
00381     if( target_dim == 3 )
00382     {
00383         manifold.merge( tmp_range );
00384         return MB_SUCCESS;
00385     }
00386 
00387     for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit )
00388     {
00389         Range dum_range;
00390         // get (target_dim+1)-dimensional entities
00391         result = mbImpl->get_adjacencies( &( *rit ), 1, target_dim + 1, false, dum_range );
00392         if( MB_SUCCESS != result ) return result;
00393 
00394         // if there are only 1 or zero, add to manifold list
00395         if( 1 >= dum_range.size() ) manifold.insert( *rit );
00396     }
00397 
00398     return MB_SUCCESS;
00399 }
00400 
00401 //! get "bridge" or "2nd order" adjacencies, going through dimension bridge_dim
00402 ErrorCode MeshTopoUtil::get_bridge_adjacencies( Range& from_entities,
00403                                                 int bridge_dim,
00404                                                 int to_dim,
00405                                                 Range& to_ents,
00406                                                 int num_layers )
00407 {
00408     Range bridge_ents, accum_layers, new_toents( from_entities );
00409     ErrorCode result;
00410     if( 0 == num_layers || from_entities.empty() ) return MB_FAILURE;
00411 
00412     // for each layer, get bridge-adj entities and accumulate
00413     for( int nl = 0; nl < num_layers; nl++ )
00414     {
00415         Range new_bridges;
00416         // get bridge ents
00417         result = mbImpl->get_adjacencies( new_toents, bridge_dim, true, new_bridges, Interface::UNION );
00418         if( MB_SUCCESS != result ) return result;
00419 
00420         // get to_dim adjacencies, merge into to_ents
00421         Range new_layer;
00422         if( -1 == to_dim )
00423         {
00424             result = mbImpl->get_adjacencies( new_bridges, 3, false, new_layer, Interface::UNION );
00425             if( MB_SUCCESS != result ) return result;
00426             for( int d = 2; d >= 1; d-- )
00427             {
00428                 result = mbImpl->get_adjacencies( to_ents, d, true, new_layer, Interface::UNION );
00429                 if( MB_SUCCESS != result ) return result;
00430             }
00431         }
00432         else
00433         {
00434             result = mbImpl->get_adjacencies( new_bridges, to_dim, false, new_layer, Interface::UNION );
00435             if( MB_SUCCESS != result ) return result;
00436         }
00437 
00438         // subtract last_toents to get new_toents
00439         accum_layers.merge( new_layer );
00440         if( nl < num_layers - 1 ) new_toents = subtract( new_layer, new_toents );
00441     }
00442     to_ents.merge( accum_layers );
00443 
00444     return MB_SUCCESS;
00445 }
00446 
00447 //! get "bridge" or "2nd order" adjacencies, going through dimension bridge_dim
00448 ErrorCode MeshTopoUtil::get_bridge_adjacencies( const EntityHandle from_entity,
00449                                                 const int bridge_dim,
00450                                                 const int to_dim,
00451                                                 Range& to_adjs )
00452 {
00453     // get pointer to connectivity for this entity
00454     const EntityHandle* connect;
00455     int num_connect;
00456     ErrorCode result     = MB_SUCCESS;
00457     EntityType from_type = TYPE_FROM_HANDLE( from_entity );
00458     if( from_type == MBVERTEX )
00459     {
00460         connect     = &from_entity;
00461         num_connect = 1;
00462     }
00463     else
00464     {
00465         result = mbImpl->get_connectivity( from_entity, connect, num_connect );
00466         if( MB_SUCCESS != result ) return result;
00467     }
00468 
00469     if( from_type >= MBENTITYSET ) return MB_FAILURE;
00470 
00471     int from_dim = CN::Dimension( from_type );
00472 
00473     Range to_ents;
00474 
00475     if( bridge_dim < from_dim )
00476     {
00477         // looping over each sub-entity of dimension bridge_dim...
00478         if( MBPOLYGON == from_type )
00479         {
00480             for( int i = 0; i < num_connect; i++ )
00481             {
00482                 // loop over edges, and get the vertices
00483                 EntityHandle verts_on_edge[2] = { connect[i], connect[( i + 1 ) % num_connect] };
00484                 to_ents.clear();
00485                 ErrorCode tmp_result =
00486                     mbImpl->get_adjacencies( verts_on_edge, 2, to_dim, false, to_ents, Interface::INTERSECT );
00487                 if( MB_SUCCESS != tmp_result ) result = tmp_result;
00488                 to_adjs.merge( to_ents );
00489             }
00490         }
00491         else
00492         {
00493             EntityHandle bridge_verts[MAX_SUB_ENTITIES];
00494             int bridge_indices[MAX_SUB_ENTITIES];
00495             for( int i = 0; i < CN::NumSubEntities( from_type, bridge_dim ); i++ )
00496             {
00497 
00498                 // get the vertices making up this sub-entity
00499                 int num_bridge_verts = CN::VerticesPerEntity( CN::SubEntityType( from_type, bridge_dim, i ) );
00500                 assert( num_bridge_verts >= 0 && num_bridge_verts <= MAX_SUB_ENTITIES );
00501                 CN::SubEntityVertexIndices( from_type, bridge_dim, i, bridge_indices );
00502                 for( int j = 0; j < num_bridge_verts; ++j )
00503                 {
00504                     if( bridge_indices[j] >= 0 && bridge_indices[j] < num_connect )
00505                         bridge_verts[j] = connect[bridge_indices[j]];
00506                     else
00507                         bridge_verts[j] = 0;
00508                 }
00509                 // CN::SubEntityConn(connect, from_type, bridge_dim, i, &bridge_verts[0],
00510                 // num_bridge_verts);
00511 
00512                 // get the to_dim entities adjacent
00513                 to_ents.clear();
00514                 ErrorCode tmp_result = mbImpl->get_adjacencies( bridge_verts, num_bridge_verts, to_dim, false, to_ents,
00515                                                                 Interface::INTERSECT );
00516                 if( MB_SUCCESS != tmp_result ) result = tmp_result;
00517 
00518                 to_adjs.merge( to_ents );
00519             }
00520         }
00521     }
00522 
00523     // now get the direct ones too, or only in the case where we're
00524     // going to higher dimension for bridge
00525     Range bridge_ents, tmp_ents;
00526     tmp_ents.insert( from_entity );
00527     ErrorCode tmp_result = mbImpl->get_adjacencies( tmp_ents, bridge_dim, false, bridge_ents, Interface::UNION );
00528     if( MB_SUCCESS != tmp_result ) return tmp_result;
00529 
00530     tmp_result = mbImpl->get_adjacencies( bridge_ents, to_dim, false, to_adjs, Interface::UNION );
00531     if( MB_SUCCESS != tmp_result ) return tmp_result;
00532 
00533     // if to_dimension is same as that of from_entity, make sure from_entity isn't
00534     // in list
00535     if( to_dim == from_dim ) to_adjs.erase( from_entity );
00536 
00537     return result;
00538 }
00539 
00540 //! return a common entity of the specified dimension, or 0 if there isn't one
00541 EntityHandle MeshTopoUtil::common_entity( const EntityHandle ent1, const EntityHandle ent2, const int dim )
00542 {
00543     Range tmp_range, tmp_range2;
00544     tmp_range.insert( ent1 );
00545     tmp_range.insert( ent2 );
00546     ErrorCode result = mbImpl->get_adjacencies( tmp_range, dim, false, tmp_range2 );
00547     if( MB_SUCCESS != result || tmp_range2.empty() )
00548         return 0;
00549     else
00550         return *tmp_range2.begin();
00551 }
00552 
00553 //! return the opposite side entity given a parent and bounding entity.
00554 //! This function is only defined for certain types of parent/child types;
00555 //! See CN.hpp::OppositeSide for details.
00556 //!
00557 //! \param parent The parent element
00558 //! \param child The child element
00559 //! \param opposite_element The index of the opposite element
00560 ErrorCode MeshTopoUtil::opposite_entity( const EntityHandle parent,
00561                                          const EntityHandle child,
00562                                          EntityHandle& opposite_element )
00563 {
00564     // get the side no.
00565     int side_no, offset, sense;
00566     ErrorCode result = mbImpl->side_number( parent, child, side_no, offset, sense );
00567     if( MB_SUCCESS != result ) return result;
00568 
00569     // get the child index from CN
00570     int opposite_index, opposite_dim;
00571     int status = CN::OppositeSide( mbImpl->type_from_handle( parent ), side_no, mbImpl->dimension_from_handle( child ),
00572                                    opposite_index, opposite_dim );
00573     if( 0 != status ) return MB_FAILURE;
00574 
00575     // now get the side element from MOAB
00576     result = mbImpl->side_element( parent, opposite_dim, opposite_index, opposite_element );
00577     if( MB_SUCCESS != result ) return result;
00578 
00579     return MB_SUCCESS;
00580 }
00581 
00582 ErrorCode MeshTopoUtil::split_entities_manifold( Range& entities, Range& new_entities, Range* fill_entities )
00583 {
00584     Range tmp_range, *tmp_ptr_fill_entity;
00585     if( NULL != fill_entities )
00586         tmp_ptr_fill_entity = &tmp_range;
00587     else
00588         tmp_ptr_fill_entity = NULL;
00589 
00590     for( Range::iterator rit = entities.begin(); rit != entities.end(); ++rit )
00591     {
00592         EntityHandle new_entity;
00593         if( NULL != tmp_ptr_fill_entity ) tmp_ptr_fill_entity->clear();
00594 
00595         EntityHandle this_ent = *rit;
00596         ErrorCode result      = split_entities_manifold( &this_ent, 1, &new_entity, tmp_ptr_fill_entity );
00597         if( MB_SUCCESS != result ) return result;
00598 
00599         new_entities.insert( new_entity );
00600         if( NULL != fill_entities ) fill_entities->merge( *tmp_ptr_fill_entity );
00601     }
00602 
00603     return MB_SUCCESS;
00604 }
00605 
00606 ErrorCode MeshTopoUtil::split_entities_manifold( EntityHandle* entities,
00607                                                  const int num_entities,
00608                                                  EntityHandle* new_entities,
00609                                                  Range* fill_entities,
00610                                                  EntityHandle* gowith_ents )
00611 {
00612     // split entities by duplicating them; splitting manifold means that there is at
00613     // most two higher-dimension entities bounded by a given entity; after split, the
00614     // new entity bounds one and the original entity bounds the other
00615 
00616 #define ITERATE_RANGE( range, it ) for( Range::iterator it = ( range ).begin(); ( it ) != ( range ).end(); ++( it ) )
00617 #define GET_CONNECT_DECL( ent, connect, num_connect )                                     \
00618     const EntityHandle* connect = NULL;                                                   \
00619     int num_connect             = 0;                                                      \
00620     {                                                                                     \
00621         ErrorCode connect_result = mbImpl->get_connectivity( ent, connect, num_connect ); \
00622         if( MB_SUCCESS != connect_result ) return connect_result;                         \
00623     }
00624 #define GET_CONNECT( ent, connect, num_connect )                                          \
00625     {                                                                                     \
00626         ErrorCode connect_result = mbImpl->get_connectivity( ent, connect, num_connect ); \
00627         if( MB_SUCCESS != connect_result ) return connect_result;                         \
00628     }
00629 #define TC                         \
00630     if( MB_SUCCESS != tmp_result ) \
00631     {                              \
00632         result = tmp_result;       \
00633         continue;                  \
00634     }
00635 
00636     ErrorCode result = MB_SUCCESS;
00637     for( int i = 0; i < num_entities; i++ )
00638     {
00639         ErrorCode tmp_result;
00640 
00641         // get original higher-dimensional bounding entities
00642         Range up_adjs[4];
00643         // can only do a split_manifold if there are at most 2 entities of each
00644         // higher dimension; otherwise it's a split non-manifold
00645         bool valid_up_adjs = true;
00646         for( int dim = 1; dim <= 3; dim++ )
00647         {
00648             tmp_result = mbImpl->get_adjacencies( entities + i, 1, dim, false, up_adjs[dim] );
00649             TC;
00650             if( dim > CN::Dimension( TYPE_FROM_HANDLE( entities[i] ) ) && up_adjs[dim].size() > 2 )
00651             {
00652                 valid_up_adjs = false;
00653                 break;
00654             }
00655         }
00656         if( !valid_up_adjs ) return MB_FAILURE;
00657 
00658         // ok to split; create the new entity, with connectivity of the original
00659         GET_CONNECT_DECL( entities[i], connect, num_connect );
00660         EntityHandle new_entity;
00661         result = mbImpl->create_element( mbImpl->type_from_handle( entities[i] ), connect, num_connect, new_entity );
00662         TC;
00663 
00664         // by definition, new entity and original will be equivalent; need to add explicit
00665         // adjs to distinguish them; don't need to check if there's already one there,
00666         // 'cuz add_adjacency does that for us
00667         for( int dim = 1; dim <= 3; dim++ )
00668         {
00669             if( up_adjs[dim].empty() || dim == CN::Dimension( TYPE_FROM_HANDLE( entities[i] ) ) ) continue;
00670 
00671             if( dim < CN::Dimension( TYPE_FROM_HANDLE( entities[i] ) ) )
00672             {
00673                 // adjacencies from other entities to this one; if any of those are equivalent
00674                 // entities, need to make explicit adjacency to new entity too
00675                 for( Range::iterator rit = up_adjs[dim].begin(); rit != up_adjs[dim].end(); ++rit )
00676                 {
00677                     if( equivalent_entities( *rit ) ) result = mbImpl->add_adjacencies( *rit, &new_entity, 1, false );
00678                 }
00679             }
00680             else
00681             {
00682 
00683                 // get the two up-elements
00684                 EntityHandle up_elem1 = *( up_adjs[dim].begin() ),
00685                              up_elem2 = ( up_adjs[dim].size() > 1 ? *( up_adjs[dim].rbegin() ) : 0 );
00686 
00687                 // if two, and a gowith entity was input, make sure the new entity goes with
00688                 // that one
00689                 if( gowith_ents && up_elem2 && gowith_ents[i] != up_elem1 && gowith_ents[i] == up_elem2 )
00690                 {
00691                     EntityHandle tmp_elem = up_elem1;
00692                     up_elem1              = up_elem2;
00693                     up_elem2              = tmp_elem;
00694                 }
00695 
00696                 mbImpl->remove_adjacencies( entities[i], &up_elem1, 1 );
00697                 // (ok if there's an error, that just means there wasn't an explicit adj)
00698 
00699                 tmp_result = mbImpl->add_adjacencies( new_entity, &up_elem1, 1, false );
00700                 TC;
00701                 if( !up_elem2 ) continue;
00702 
00703                 // add adj to other up_adj
00704                 tmp_result = mbImpl->add_adjacencies( entities[i], &up_elem2, 1, false );
00705                 TC;
00706             }
00707         }
00708 
00709         // if we're asked to build a next-higher-dimension object, do so
00710         EntityHandle fill_entity = 0;
00711         EntityHandle tmp_ents[2];
00712         if( NULL != fill_entities )
00713         {
00714             // how to do this depends on dimension
00715             switch( CN::Dimension( TYPE_FROM_HANDLE( entities[i] ) ) )
00716             {
00717                 case 0:
00718                     tmp_ents[0] = entities[i];
00719                     tmp_ents[1] = new_entity;
00720                     tmp_result  = mbImpl->create_element( MBEDGE, tmp_ents, 2, fill_entity );
00721                     TC;
00722                     break;
00723                 case 1:
00724                     tmp_result = mbImpl->create_element( MBPOLYGON, connect, 2, fill_entity );
00725                     TC;
00726                     // need to create explicit adj in this case
00727                     tmp_result = mbImpl->add_adjacencies( entities[i], &fill_entity, 1, false );
00728                     TC;
00729                     tmp_result = mbImpl->add_adjacencies( new_entity, &fill_entity, 1, false );
00730                     TC;
00731                     break;
00732                 case 2:
00733                     tmp_ents[0] = entities[i];
00734                     tmp_ents[1] = new_entity;
00735                     tmp_result  = mbImpl->create_element( MBPOLYHEDRON, tmp_ents, 2, fill_entity );
00736                     TC;
00737                     break;
00738             }
00739             if( 0 == fill_entity )
00740             {
00741                 result = MB_FAILURE;
00742                 continue;
00743             }
00744             fill_entities->insert( fill_entity );
00745         }
00746 
00747         new_entities[i] = new_entity;
00748 
00749     }  // end for over input entities
00750 
00751     return result;
00752 }
00753 
00754 ErrorCode MeshTopoUtil::split_entity_nonmanifold( EntityHandle split_ent,
00755                                                   Range& old_adjs,
00756                                                   Range& new_adjs,
00757                                                   EntityHandle& new_entity )
00758 {
00759     // split an entity into two entities; new entity gets explicit adj to new_adjs,
00760     // old to old_adjs
00761 
00762     // make new entities and add adjacencies
00763     // create the new entity
00764     EntityType split_type = mbImpl->type_from_handle( split_ent );
00765 
00766     ErrorCode result;
00767     if( MBVERTEX == split_type )
00768     {
00769         double coords[3];
00770         result = mbImpl->get_coords( &split_ent, 1, coords );RR;
00771         result = mbImpl->create_vertex( coords, new_entity );RR;
00772     }
00773     else
00774     {
00775         const EntityHandle* connect;
00776         int num_connect;
00777         result = mbImpl->get_connectivity( split_ent, connect, num_connect );RR;
00778         result = mbImpl->create_element( split_type, connect, num_connect, new_entity );RR;
00779 
00780         // remove any explicit adjacencies between new_adjs and split entity
00781         for( Range::iterator rit = new_adjs.begin(); rit != new_adjs.end(); ++rit )
00782             mbImpl->remove_adjacencies( split_ent, &( *rit ), 1 );
00783     }
00784 
00785     if( MBVERTEX != split_type )
00786     {
00787         //  add adj's between new_adjs & new entity, old_adjs & split_entity
00788         for( Range::iterator rit = new_adjs.begin(); rit != new_adjs.end(); ++rit )
00789             mbImpl->add_adjacencies( new_entity, &( *rit ), 1, true );
00790         for( Range::iterator rit = old_adjs.begin(); rit != old_adjs.end(); ++rit )
00791             mbImpl->add_adjacencies( split_ent, &( *rit ), 1, true );
00792     }
00793     else if( split_ent != new_entity )
00794     {
00795         // in addition to explicit adjs, need to check if vertex is part of any
00796         // other entities, and check those entities against ents in old and new adjs
00797         Range other_adjs;
00798         for( int i = 1; i < 4; i++ )
00799         {
00800             result = mbImpl->get_adjacencies( &split_ent, 1, i, false, other_adjs, Interface::UNION );RR;
00801         }
00802         other_adjs = subtract( other_adjs, old_adjs );
00803         other_adjs = subtract( other_adjs, new_adjs );
00804         for( Range::iterator rit1 = other_adjs.begin(); rit1 != other_adjs.end(); ++rit1 )
00805         {
00806             // find an adjacent lower-dimensional entity in old_ or new_ adjs
00807             bool found = false;
00808             for( Range::iterator rit2 = old_adjs.begin(); rit2 != old_adjs.end(); ++rit2 )
00809             {
00810                 if( mbImpl->dimension_from_handle( *rit1 ) != mbImpl->dimension_from_handle( *rit2 ) &&
00811                     common_entity( *rit1, *rit2, mbImpl->dimension_from_handle( *rit1 ) ) )
00812                 {
00813                     found = true;
00814                     old_adjs.insert( *rit1 );
00815                     break;
00816                 }
00817             }
00818             if( found ) continue;
00819             for( Range::iterator rit2 = new_adjs.begin(); rit2 != new_adjs.end(); ++rit2 )
00820             {
00821                 if( mbImpl->dimension_from_handle( *rit1 ) != mbImpl->dimension_from_handle( *rit2 ) &&
00822                     common_entity( *rit1, *rit2, mbImpl->dimension_from_handle( *rit1 ) ) )
00823                 {
00824                     found = true;
00825                     new_adjs.insert( *rit1 );
00826                     break;
00827                 }
00828             }
00829             if( !found ) return MB_FAILURE;
00830         }
00831 
00832         // instead of adjs replace in connectivity
00833         std::vector< EntityHandle > connect;
00834         for( Range::iterator rit = new_adjs.begin(); rit != new_adjs.end(); ++rit )
00835         {
00836             connect.clear();
00837             result = mbImpl->get_connectivity( &( *rit ), 1, connect );RR;
00838             std::replace( connect.begin(), connect.end(), split_ent, new_entity );
00839             result = mbImpl->set_connectivity( *rit, &connect[0], connect.size() );RR;
00840         }
00841     }
00842 
00843     return result;
00844 
00845     /*
00846 
00847     Commented out for now, because I decided to do a different implementation
00848     for the sake of brevity.  However, I still think this function is the right
00849     way to do it, if I ever get the time.  Sigh.
00850 
00851         // split entity d, producing entity nd; generates various new entities,
00852         // see algorithm description in notes from 2/25/05
00853       const EntityHandle split_types = {MBEDGE, MBPOLYGON, MBPOLYHEDRON};
00854       ErrorCode result = MB_SUCCESS;
00855       const int dim = CN::Dimension(TYPE_FROM_HANDLE(d));
00856       MeshTopoUtil mtu(this);
00857 
00858         // get all (d+2)-, (d+1)-cells connected to d
00859       Range dp2s, dp1s, dp1s_manif, dp2s_manif;
00860       result = get_adjacencies(&d, 1, dim+2, false, dp2s); RR;
00861       result = get_adjacencies(&d, 1, dim+1, false, dp1s); RR;
00862 
00863         // also get (d+1)-cells connected to d which are manifold
00864       get_manifold_dp1s(d, dp1s_manif);
00865       get_manifold_dp2s(d, dp2s_manif);
00866 
00867         // make new cell nd, then ndp1
00868       result = copy_entity(d, nd); RR;
00869       EntityHandle tmp_connect[] = {d, nd};
00870       EntityHandle ndp1;
00871       result = create_element(split_types[dim],
00872                               tmp_connect, 2, ndp1); RR;
00873 
00874         // modify (d+2)-cells, depending on what type they are
00875       ITERATE_RANGE(dp2s, dp2) {
00876           // first, get number of connected manifold (d+1)-entities
00877         Range tmp_range, tmp_range2(dp1s_manif);
00878         tmp_range.insert(*dp2);
00879         tmp_range.insert(d);
00880         tmp_result = get_adjacencies(tmp_range, 1, false, tmp_range2); TC;
00881         EntityHandle ndp2;
00882 
00883           // a. manif (d+1)-cells is zero
00884         if (tmp_range2.empty()) {
00885             // construct new (d+1)-cell
00886           EntityHandle ndp1a;
00887           EntityHandle tmp_result = create_element(split_types[dim],
00888                                                      tmp_connect, 2, ndp1a); TC;
00889             // now make new (d+2)-cell
00890           EntityHandle tmp_connect2[] = {ndp1, ndp1a};
00891           tmp_result = create_element(split_types[dim+1],
00892                                       tmp_connect2, 2, ndp2); TC;
00893             // need to add explicit adjacencies, since by definition ndp1, ndp1a will be equivalent
00894           tmp_result = add_adjacencies(ndp1a, &dp2, 1, false); TC;
00895           tmp_result = add_adjacencies(ndp1a, &ndp2, 1, false); TC;
00896           tmp_result = add_adjacencies(ndp1, &ndp2, 1, false); TC;
00897 
00898             // now insert nd into connectivity of dp2, right after d if dim < 1
00899           std::vector<EntityHandle> connect;
00900           tmp_result = get_connectivity(&dp2, 1, connect); TC;
00901           if (dim < 1) {
00902             std::vector<EntityHandle>::iterator vit = std::find(connect.begin(), connect.end(), d);
00903             if (vit == connect.end()) {
00904               result = MB_FAILURE;
00905               continue;
00906             }
00907             connect.insert(vit, nd);
00908           }
00909           else
00910             connect.push_back(nd);
00911           tmp_result = set_connectivity(dp2, connect); TC;
00912 
00913             // if dim < 1, need to add explicit adj from ndp2 to higher-dim ents, since it'll
00914             // be equiv to other dp2 entities
00915           if (dim < 1) {
00916             Range tmp_dp3s;
00917             tmp_result = get_adjacencies(&dp2, 1, dim+3, false, tmp_dp3s); TC;
00918             tmp_result = add_adjacencies(ndp2, tmp_dp3s, false); TC;
00919           }
00920         } // end if (tmp_range2.empty())
00921 
00922           // b. single manifold (d+1)-cell, which isn't adjacent to manifold (d+2)-cell
00923         else if (tmp_range2.size() == 1) {
00924             // b1. check validity, and skip if not valid
00925 
00926             // only change if not dp1-adjacent to manifold dp2cell; check that...
00927           Range tmp_adjs(dp2s_manif);
00928           tmp_result = get_adjacencies(&(*tmp_range2.begin()), 1, dim+2, false, tmp_adjs); TC;
00929           if (!tmp_adjs.empty()) continue;
00930 
00931           EntityHandle dp1 = *tmp_range2.begin();
00932 
00933             // b2. make new (d+1)- and (d+2)-cell next to dp2
00934 
00935             // get the (d+2)-cell on the other side of dp1
00936           tmp_result = get_adjacencies(&dp1, 1, dim+2, false, tmp_adjs); TC;
00937           EntityHandle odp2 = *tmp_adjs.begin();
00938           if (odp2 == dp2) odp2 = *tmp_adjs.rbegin();
00939 
00940             // get od, the d-cell on dp1_manif which isn't d
00941           tmp_result = get_adjacencies(&dp1_manif, 1, dim, false, tmp_adjs); TC;
00942           tmp_adjs.erase(d);
00943           if (tmp_adjs.size() != 1) {
00944             result = MB_FAILURE;
00945             continue;
00946           }
00947           EntityHandle od = *tmp_adjs.begin();
00948 
00949             // make a new (d+1)-cell from od and nd
00950           tmp_adjs.insert(nd);
00951           tmp_result = create_element(split_types[1], tmp_adjs, ndp1a); TC;
00952 
00953             // construct new (d+2)-cell from dp1, ndp1, ndp1a
00954           tmp_adjs.clear();
00955           tmp_adjs.insert(dp1); tmp_adjs.insert(ndp1); tmp_adjs.insert(ndp1a);
00956           tmp_result = create_element(split_types[2], tmp_adjs, ndp2); TC;
00957 
00958             // b3. replace d, dp1 in connect/adjs of odp2
00959           std::vector<EntityHandle> connect;
00960           tmp_result = get_connectivity(&odp2, 1, connect); TC;
00961           if (dim == 0) {
00962             *(std::find(connect.begin(), connect.end(), d)) = nd;
00963             remove_adjacency(dp1, odp2);
00964 
00965 
00966 
00967             // if dp1 was explicitly adj to odp2, remove it
00968           remove_adjacency
00969 
00970     ...
00971 
00972     */
00973 }
00974 
00975 //! return whether entity is equivalent to any other of same type and same vertices;
00976 //! if equivalent entity is found, it's returned in equiv_ents and return value is true,
00977 //! false otherwise.
00978 bool MeshTopoUtil::equivalent_entities( const EntityHandle entity, Range* equiv_ents )
00979 {
00980     const EntityHandle* connect = NULL;
00981     int num_connect             = 0;
00982     ErrorCode result            = mbImpl->get_connectivity( entity, connect, num_connect );
00983     if( MB_SUCCESS != result ) return false;
00984 
00985     Range dum;
00986     result = mbImpl->get_adjacencies( connect, num_connect, mbImpl->dimension_from_handle( entity ), false, dum );
00987     dum.erase( entity );
00988 
00989     if( NULL != equiv_ents )
00990     {
00991         equiv_ents->swap( dum );
00992     }
00993 
00994     if( !dum.empty() )
00995         return true;
00996     else
00997         return false;
00998 }
00999 
01000 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines