Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
Skinner.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 __MFC_VER
00017 #pragma warning( disable : 4786 )
00018 #endif
00019 
00020 #ifdef WIN32               /* windows */
00021 #define _USE_MATH_DEFINES  // For M_PI
00022 #endif
00023 #include "moab/Skinner.hpp"
00024 #include "moab/Range.hpp"
00025 #include "moab/CN.hpp"
00026 #include <vector>
00027 #include <set>
00028 #include <algorithm>
00029 #include <cmath>
00030 #include <cassert>
00031 #include <iostream>
00032 #include "moab/Util.hpp"
00033 #include "Internals.hpp"
00034 #include "MBTagConventions.hpp"
00035 #include "moab/Core.hpp"
00036 #include "AEntityFactory.hpp"
00037 #include "moab/ScdInterface.hpp"
00038 
00039 #ifdef M_PI
00040 #define SKINNER_PI M_PI
00041 #else
00042 #define SKINNER_PI 3.1415926535897932384626
00043 #endif
00044 
00045 namespace moab
00046 {
00047 
00048 Skinner::~Skinner()
00049 {
00050     // delete the adjacency tag
00051 }
00052 
00053 ErrorCode Skinner::initialize()
00054 {
00055     // go through and mark all the target dimension entities
00056     // that already exist as not deleteable
00057     // also get the connectivity tags for each type
00058     // also populate adjacency information
00059     EntityType type;
00060     DimensionPair target_ent_types = CN::TypeDimensionMap[mTargetDim];
00061 
00062     void* null_ptr = NULL;
00063 
00064     ErrorCode result = thisMB->tag_get_handle( "skinner adj", sizeof( void* ), MB_TYPE_OPAQUE, mAdjTag,
00065                                                MB_TAG_DENSE | MB_TAG_CREAT, &null_ptr );MB_CHK_ERR( result );
00066 
00067     if( mDeletableMBTag == 0 )
00068     {
00069         result =
00070             thisMB->tag_get_handle( "skinner deletable", 1, MB_TYPE_BIT, mDeletableMBTag, MB_TAG_BIT | MB_TAG_CREAT );MB_CHK_ERR( result );
00071     }
00072 
00073     Range entities;
00074 
00075     // go through each type at this dimension
00076     for( type = target_ent_types.first; type <= target_ent_types.second; ++type )
00077     {
00078         // get the entities of this type in the MB
00079         thisMB->get_entities_by_type( 0, type, entities );
00080 
00081         // go through each entity of this type in the MB
00082         // and set its deletable tag to NO
00083         Range::iterator iter, end_iter;
00084         end_iter = entities.end();
00085         for( iter = entities.begin(); iter != end_iter; ++iter )
00086         {
00087             unsigned char bit = 0x1;
00088             result            = thisMB->tag_set_data( mDeletableMBTag, &( *iter ), 1, &bit );
00089             assert( MB_SUCCESS == result );
00090             // add adjacency information too
00091             if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) add_adjacency( *iter );
00092         }
00093     }
00094 
00095     return MB_SUCCESS;
00096 }
00097 
00098 ErrorCode Skinner::deinitialize()
00099 {
00100     ErrorCode result;
00101     if( 0 != mDeletableMBTag )
00102     {
00103         result          = thisMB->tag_delete( mDeletableMBTag );
00104         mDeletableMBTag = 0;MB_CHK_ERR( result );
00105     }
00106 
00107     // remove the adjacency tag
00108     std::vector< std::vector< EntityHandle >* > adj_arr;
00109     std::vector< std::vector< EntityHandle >* >::iterator i;
00110     if( 0 != mAdjTag )
00111     {
00112         for( EntityType t = MBVERTEX; t != MBMAXTYPE; ++t )
00113         {
00114             Range entities;
00115             result = thisMB->get_entities_by_type_and_tag( 0, t, &mAdjTag, 0, 1, entities );MB_CHK_ERR( result );
00116             adj_arr.resize( entities.size() );
00117             result = thisMB->tag_get_data( mAdjTag, entities, &adj_arr[0] );MB_CHK_ERR( result );
00118             for( i = adj_arr.begin(); i != adj_arr.end(); ++i )
00119                 delete *i;
00120         }
00121 
00122         result  = thisMB->tag_delete( mAdjTag );
00123         mAdjTag = 0;MB_CHK_ERR( result );
00124     }
00125 
00126     return MB_SUCCESS;
00127 }
00128 
00129 ErrorCode Skinner::add_adjacency( EntityHandle entity )
00130 {
00131     std::vector< EntityHandle >* adj = NULL;
00132     const EntityHandle* nodes;
00133     int num_nodes;
00134     ErrorCode result = thisMB->get_connectivity( entity, nodes, num_nodes, true );MB_CHK_ERR( result );
00135     const EntityHandle* iter = std::min_element( nodes, nodes + num_nodes );
00136 
00137     if( iter == nodes + num_nodes ) return MB_SUCCESS;
00138 
00139     // add this entity to the node
00140     if( thisMB->tag_get_data( mAdjTag, iter, 1, &adj ) == MB_SUCCESS && adj != NULL )
00141     {
00142         adj->push_back( entity );
00143     }
00144     // create a new vector and add it
00145     else
00146     {
00147         adj = new std::vector< EntityHandle >;
00148         adj->push_back( entity );
00149         result = thisMB->tag_set_data( mAdjTag, iter, 1, &adj );MB_CHK_ERR( result );
00150     }
00151 
00152     return MB_SUCCESS;
00153 }
00154 
00155 void Skinner::add_adjacency( EntityHandle entity, const EntityHandle* nodes, const int num_nodes )
00156 {
00157     std::vector< EntityHandle >* adj = NULL;
00158     const EntityHandle* iter         = std::min_element( nodes, nodes + num_nodes );
00159 
00160     if( iter == nodes + num_nodes ) return;
00161 
00162     // should not be setting adjacency lists in ho-nodes
00163     assert( TYPE_FROM_HANDLE( entity ) == MBPOLYGON ||
00164             num_nodes == CN::VerticesPerEntity( TYPE_FROM_HANDLE( entity ) ) );
00165 
00166     // add this entity to the node
00167     if( thisMB->tag_get_data( mAdjTag, iter, 1, &adj ) == MB_SUCCESS && adj != NULL )
00168     {
00169         adj->push_back( entity );
00170     }
00171     // create a new vector and add it
00172     else
00173     {
00174         adj = new std::vector< EntityHandle >;
00175         adj->push_back( entity );
00176         thisMB->tag_set_data( mAdjTag, iter, 1, &adj );
00177     }
00178 }
00179 
00180 ErrorCode Skinner::find_geometric_skin( const EntityHandle meshset, Range& forward_target_entities )
00181 {
00182     // attempts to find whole model skin, using geom topo sets first then
00183     // normal find_skin function
00184     bool debug = true;
00185 
00186     // look for geom topo sets
00187     Tag geom_tag;
00188     ErrorCode result =
00189         thisMB->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
00190 
00191     if( MB_SUCCESS != result ) return result;
00192 
00193     // get face sets (dimension = 2)
00194     Range face_sets;
00195     int two             = 2;
00196     const void* two_ptr = &two;
00197     result = thisMB->get_entities_by_type_and_tag( meshset, MBENTITYSET, &geom_tag, &two_ptr, 1, face_sets );
00198 
00199     Range::iterator it;
00200     if( MB_SUCCESS != result )
00201         return result;
00202     else if( face_sets.empty() )
00203         return MB_ENTITY_NOT_FOUND;
00204 
00205     // ok, we have face sets; use those to determine skin
00206     Range skin_sets;
00207     if( debug ) std::cout << "Found " << face_sets.size() << " face sets total..." << std::endl;
00208 
00209     for( it = face_sets.begin(); it != face_sets.end(); ++it )
00210     {
00211         int num_parents;
00212         result = thisMB->num_parent_meshsets( *it, &num_parents );
00213         if( MB_SUCCESS != result )
00214             return result;
00215         else if( num_parents == 1 )
00216             skin_sets.insert( *it );
00217     }
00218 
00219     if( debug ) std::cout << "Found " << skin_sets.size() << " 1-parent face sets..." << std::endl;
00220 
00221     if( skin_sets.empty() ) return MB_FAILURE;
00222 
00223     // ok, we have the shell; gather up the elements, putting them all in forward for now
00224     for( it = skin_sets.begin(); it != skin_sets.end(); ++it )
00225     {
00226         result = thisMB->get_entities_by_handle( *it, forward_target_entities, true );
00227         if( MB_SUCCESS != result ) return result;
00228     }
00229 
00230     return result;
00231 }
00232 
00233 ErrorCode Skinner::find_skin( const EntityHandle meshset,
00234                               const Range& source_entities,
00235                               bool get_vertices,
00236                               Range& output_handles,
00237                               Range* output_reverse_handles,
00238                               bool create_vert_elem_adjs,
00239                               bool create_skin_elements,
00240                               bool look_for_scd )
00241 {
00242     if( source_entities.empty() ) return MB_SUCCESS;
00243 
00244     if( look_for_scd )
00245     {
00246         ErrorCode rval = find_skin_scd( source_entities, get_vertices, output_handles, create_skin_elements );
00247         // if it returns success, it's all scd, and we don't need to do anything more
00248         if( MB_SUCCESS == rval ) return rval;
00249     }
00250 
00251     Core* this_core = dynamic_cast< Core* >( thisMB );
00252     if( this_core && create_vert_elem_adjs && !this_core->a_entity_factory()->vert_elem_adjacencies() )
00253         this_core->a_entity_factory()->create_vert_elem_adjacencies();
00254 
00255     return find_skin_vertices( meshset, source_entities, get_vertices ? &output_handles : 0,
00256                                get_vertices ? 0 : &output_handles, output_reverse_handles, create_skin_elements );
00257 }
00258 
00259 ErrorCode Skinner::find_skin_scd( const Range& source_entities,
00260                                   bool get_vertices,
00261                                   Range& output_handles,
00262                                   bool create_skin_elements )
00263 {
00264     // get the scd interface and check if it's been initialized
00265     ScdInterface* scdi = NULL;
00266     ErrorCode rval     = thisMB->query_interface( scdi );
00267     if( !scdi ) return MB_FAILURE;
00268 
00269     // ok, there's scd mesh; see if the entities passed in are all in a scd box
00270     // a box needs to be wholly included in entities for this to work
00271     std::vector< ScdBox* > boxes, myboxes;
00272     Range myrange;
00273     rval = scdi->find_boxes( boxes );
00274     if( MB_SUCCESS != rval ) return rval;
00275     for( std::vector< ScdBox* >::iterator bit = boxes.begin(); bit != boxes.end(); ++bit )
00276     {
00277         Range belems( ( *bit )->start_element(), ( *bit )->start_element() + ( *bit )->num_elements() - 1 );
00278         if( source_entities.contains( belems ) )
00279         {
00280             myboxes.push_back( *bit );
00281             myrange.merge( belems );
00282         }
00283     }
00284     if( myboxes.empty() || myrange.size() != source_entities.size() ) return MB_FAILURE;
00285 
00286     // ok, we're all structured; get the skin for each box
00287     for( std::vector< ScdBox* >::iterator bit = boxes.begin(); bit != boxes.end(); ++bit )
00288     {
00289         rval = skin_box( *bit, get_vertices, output_handles, create_skin_elements );
00290         if( MB_SUCCESS != rval ) return rval;
00291     }
00292 
00293     return MB_SUCCESS;
00294 }
00295 
00296 ErrorCode Skinner::skin_box( ScdBox* box, bool get_vertices, Range& output_handles, bool create_skin_elements )
00297 {
00298     HomCoord bmin = box->box_min(), bmax = box->box_max();
00299 
00300     // don't support 1d boxes
00301     if( bmin.j() == bmax.j() && bmin.k() == bmax.k() ) return MB_FAILURE;
00302 
00303     int dim = ( bmin.k() == bmax.k() ? 1 : 2 );
00304 
00305     ErrorCode rval;
00306     EntityHandle ent;
00307 
00308     // i=min
00309     for( int k = bmin.k(); k < bmax.k(); k++ )
00310     {
00311         for( int j = bmin.j(); j < bmax.j(); j++ )
00312         {
00313             ent  = 0;
00314             rval = box->get_adj_edge_or_face( dim, bmin.i(), j, k, 0, ent, create_skin_elements );
00315             if( MB_SUCCESS != rval ) return rval;
00316             if( ent ) output_handles.insert( ent );
00317         }
00318     }
00319     // i=max
00320     for( int k = bmin.k(); k < bmax.k(); k++ )
00321     {
00322         for( int j = bmin.j(); j < bmax.j(); j++ )
00323         {
00324             ent  = 0;
00325             rval = box->get_adj_edge_or_face( dim, bmax.i(), j, k, 0, ent, create_skin_elements );
00326             if( MB_SUCCESS != rval ) return rval;
00327             if( ent ) output_handles.insert( ent );
00328         }
00329     }
00330     // j=min
00331     for( int k = bmin.k(); k < bmax.k(); k++ )
00332     {
00333         for( int i = bmin.i(); i < bmax.i(); i++ )
00334         {
00335             ent  = 0;
00336             rval = box->get_adj_edge_or_face( dim, i, bmin.j(), k, 1, ent, create_skin_elements );
00337             if( MB_SUCCESS != rval ) return rval;
00338             if( ent ) output_handles.insert( ent );
00339         }
00340     }
00341     // j=max
00342     for( int k = bmin.k(); k < bmax.k(); k++ )
00343     {
00344         for( int i = bmin.i(); i < bmax.i(); i++ )
00345         {
00346             ent  = 0;
00347             rval = box->get_adj_edge_or_face( dim, i, bmax.j(), k, 1, ent, create_skin_elements );
00348             if( MB_SUCCESS != rval ) return rval;
00349             if( ent ) output_handles.insert( ent );
00350         }
00351     }
00352     // k=min
00353     for( int j = bmin.j(); j < bmax.j(); j++ )
00354     {
00355         for( int i = bmin.i(); i < bmax.i(); i++ )
00356         {
00357             ent  = 0;
00358             rval = box->get_adj_edge_or_face( dim, i, j, bmin.k(), 2, ent, create_skin_elements );
00359             if( MB_SUCCESS != rval ) return rval;
00360             if( ent ) output_handles.insert( ent );
00361         }
00362     }
00363     // k=max
00364     for( int j = bmin.j(); j < bmax.j(); j++ )
00365     {
00366         for( int i = bmin.i(); i < bmax.i(); i++ )
00367         {
00368             ent  = 0;
00369             rval = box->get_adj_edge_or_face( dim, i, j, bmax.k(), 2, ent, create_skin_elements );
00370             if( MB_SUCCESS != rval ) return rval;
00371             if( ent ) output_handles.insert( ent );
00372         }
00373     }
00374 
00375     if( get_vertices )
00376     {
00377         Range verts;
00378         rval = thisMB->get_adjacencies( output_handles, 0, true, verts, Interface::UNION );
00379         if( MB_SUCCESS != rval ) return rval;
00380         output_handles.merge( verts );
00381     }
00382 
00383     return MB_SUCCESS;
00384 }
00385 
00386 ErrorCode Skinner::find_skin_noadj(const Range &source_entities,
00387                                  Range &forward_target_entities,
00388                                  Range &reverse_target_entities/*,
00389                                  bool create_vert_elem_adjs*/)
00390 {
00391     if( source_entities.empty() ) return MB_FAILURE;
00392 
00393     // get our working dimensions
00394     EntityType type      = thisMB->type_from_handle( *( source_entities.begin() ) );
00395     const int source_dim = CN::Dimension( type );
00396     mTargetDim           = source_dim - 1;
00397 
00398     // make sure we can handle the working dimensions
00399     if( mTargetDim < 0 || source_dim > 3 ) return MB_FAILURE;
00400 
00401     initialize();
00402 
00403     Range::const_iterator iter, end_iter;
00404     end_iter = source_entities.end();
00405     const EntityHandle* conn;
00406     EntityHandle match;
00407 
00408     direction direct;
00409     ErrorCode result;
00410     // assume we'll never have more than 32 vertices on a facet (checked
00411     // with assert later)
00412     EntityHandle sub_conn[32];
00413     std::vector< EntityHandle > tmp_conn_vec;
00414     int num_nodes, num_sub_nodes, num_sides;
00415     int sub_indices[32];  // Also, assume that no polygon has more than 32 nodes
00416     // we could increase that, but we will not display it right in visit moab h5m , anyway
00417     EntityType sub_type;
00418 
00419     // for each source entity
00420     for( iter = source_entities.begin(); iter != end_iter; ++iter )
00421     {
00422         // get the connectivity of this entity
00423         int actual_num_nodes_polygon = 0;
00424         result                       = thisMB->get_connectivity( *iter, conn, num_nodes, false, &tmp_conn_vec );
00425         if( MB_SUCCESS != result ) return result;
00426 
00427         type = thisMB->type_from_handle( *iter );
00428         Range::iterator seek_iter;
00429 
00430         // treat separately polygons (also, polyhedra will need special handling)
00431         if( MBPOLYGON == type )
00432         {
00433             // treat padded polygons, if existing; count backwards, see how many of the last nodes
00434             // are repeated assume connectivity is fine, otherwise we could be in trouble
00435             actual_num_nodes_polygon = num_nodes;
00436             while( actual_num_nodes_polygon >= 3 &&
00437                    conn[actual_num_nodes_polygon - 1] == conn[actual_num_nodes_polygon - 2] )
00438                 actual_num_nodes_polygon--;
00439             num_sides     = actual_num_nodes_polygon;
00440             sub_type      = MBEDGE;
00441             num_sub_nodes = 2;
00442         }
00443         else  // get connectivity of each n-1 dimension entity
00444             num_sides = CN::NumSubEntities( type, mTargetDim );
00445         for( int i = 0; i < num_sides; i++ )
00446         {
00447             if( MBPOLYGON == type )
00448             {
00449                 sub_conn[0] = conn[i];
00450                 sub_conn[1] = conn[i + 1];
00451                 if( i + 1 == actual_num_nodes_polygon ) sub_conn[1] = conn[0];
00452             }
00453             else
00454             {
00455                 CN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, sub_type, num_sub_nodes, sub_indices );
00456                 assert( (size_t)num_sub_nodes <= sizeof( sub_indices ) / sizeof( sub_indices[0] ) );
00457                 for( int j = 0; j < num_sub_nodes; j++ )
00458                     sub_conn[j] = conn[sub_indices[j]];
00459             }
00460 
00461             // see if we can match this connectivity with
00462             // an existing entity
00463             find_match( sub_type, sub_conn, num_sub_nodes, match, direct );
00464 
00465             // if there is no match, create a new entity
00466             if( match == 0 )
00467             {
00468                 EntityHandle tmphndl = 0;
00469                 int indices[MAX_SUB_ENTITY_VERTICES];
00470                 EntityType new_type;
00471                 int num_new_nodes;
00472                 if( MBPOLYGON == type )
00473                 {
00474                     new_type      = MBEDGE;
00475                     num_new_nodes = 2;
00476                 }
00477                 else
00478                 {
00479                     CN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, new_type, num_new_nodes, indices );
00480                     for( int j = 0; j < num_new_nodes; j++ )
00481                         sub_conn[j] = conn[indices[j]];
00482                 }
00483                 result = thisMB->create_element( new_type, sub_conn, num_new_nodes, tmphndl );
00484                 assert( MB_SUCCESS == result );
00485                 add_adjacency( tmphndl, sub_conn, CN::VerticesPerEntity( new_type ) );
00486                 forward_target_entities.insert( tmphndl );
00487             }
00488             // if there is a match, delete the matching entity
00489             // if we can.
00490             else
00491             {
00492                 if( ( seek_iter = forward_target_entities.find( match ) ) != forward_target_entities.end() )
00493                 {
00494                     forward_target_entities.erase( seek_iter );
00495                     remove_adjacency( match );
00496                     if( /*!use_adjs &&*/ entity_deletable( match ) )
00497                     {
00498                         result = thisMB->delete_entities( &match, 1 );
00499                         assert( MB_SUCCESS == result );
00500                     }
00501                 }
00502                 else if( ( seek_iter = reverse_target_entities.find( match ) ) != reverse_target_entities.end() )
00503                 {
00504                     reverse_target_entities.erase( seek_iter );
00505                     remove_adjacency( match );
00506                     if( /*!use_adjs &&*/ entity_deletable( match ) )
00507                     {
00508                         result = thisMB->delete_entities( &match, 1 );
00509                         assert( MB_SUCCESS == result );
00510                     }
00511                 }
00512                 else
00513                 {
00514                     if( direct == FORWARD )
00515                     {
00516                         forward_target_entities.insert( match );
00517                     }
00518                     else
00519                     {
00520                         reverse_target_entities.insert( match );
00521                     }
00522                 }
00523             }
00524         }
00525     }
00526 
00527     deinitialize();
00528 
00529     return MB_SUCCESS;
00530 }
00531 
00532 void Skinner::find_match( EntityType type,
00533                           const EntityHandle* conn,
00534                           const int num_nodes,
00535                           EntityHandle& match,
00536                           Skinner::direction& direct )
00537 {
00538     match = 0;
00539 
00540     if( type == MBVERTEX )
00541     {
00542         match  = *conn;
00543         direct = FORWARD;
00544         return;
00545     }
00546 
00547     const EntityHandle* iter = std::min_element( conn, conn + num_nodes );
00548 
00549     std::vector< EntityHandle >* adj = NULL;
00550 
00551     ErrorCode result = thisMB->tag_get_data( mAdjTag, iter, 1, &adj );
00552     if( result == MB_FAILURE || adj == NULL )
00553     {
00554         return;
00555     }
00556 
00557     std::vector< EntityHandle >::iterator jter, end_jter;
00558     end_jter = adj->end();
00559 
00560     const EntityHandle* tmp;
00561     int num_verts;
00562 
00563     for( jter = adj->begin(); jter != end_jter; ++jter )
00564     {
00565         EntityType tmp_type;
00566         tmp_type = thisMB->type_from_handle( *jter );
00567 
00568         if( type != tmp_type ) continue;
00569 
00570         result = thisMB->get_connectivity( *jter, tmp, num_verts, false );
00571         assert( MB_SUCCESS == result && num_verts >= CN::VerticesPerEntity( type ) );
00572         // FIXME: connectivity_match appears to work only for linear elements,
00573         //        so ignore higher-order nodes.
00574         if( connectivity_match( conn, tmp, CN::VerticesPerEntity( type ), direct ) )
00575         {
00576             match = *jter;
00577             break;
00578         }
00579     }
00580 }
00581 
00582 bool Skinner::connectivity_match( const EntityHandle* conn1,
00583                                   const EntityHandle* conn2,
00584                                   const int num_verts,
00585                                   Skinner::direction& direct )
00586 {
00587     const EntityHandle* iter = std::find( conn2, conn2 + num_verts, conn1[0] );
00588     if( iter == conn2 + num_verts ) return false;
00589 
00590     bool they_match = true;
00591 
00592     int i;
00593     unsigned int j = iter - conn2;
00594 
00595     // first compare forward
00596     for( i = 1; i < num_verts; ++i )
00597     {
00598         if( conn1[i] != conn2[( j + i ) % num_verts] )
00599         {
00600             they_match = false;
00601             break;
00602         }
00603     }
00604 
00605     if( they_match == true )
00606     {
00607         // need to check for reversed edges here
00608         direct = ( num_verts == 2 && j ) ? REVERSE : FORWARD;
00609         return true;
00610     }
00611 
00612     they_match = true;
00613 
00614     // then compare reverse
00615     j += num_verts;
00616     for( i = 1; i < num_verts; )
00617     {
00618         if( conn1[i] != conn2[( j - i ) % num_verts] )
00619         {
00620             they_match = false;
00621             break;
00622         }
00623         ++i;
00624     }
00625     if( they_match )
00626     {
00627         direct = REVERSE;
00628     }
00629     return they_match;
00630 }
00631 
00632 ErrorCode Skinner::remove_adjacency( EntityHandle entity )
00633 {
00634     std::vector< EntityHandle > nodes, *adj = NULL;
00635     ErrorCode result = thisMB->get_connectivity( &entity, 1, nodes );MB_CHK_ERR( result );
00636     std::vector< EntityHandle >::iterator iter = std::min_element( nodes.begin(), nodes.end() );
00637 
00638     if( iter == nodes.end() ) return MB_FAILURE;
00639 
00640     // remove this entity from the node
00641     if( thisMB->tag_get_data( mAdjTag, &( *iter ), 1, &adj ) == MB_SUCCESS && adj != NULL )
00642     {
00643         iter = std::find( adj->begin(), adj->end(), entity );
00644         if( iter != adj->end() ) adj->erase( iter );
00645     }
00646 
00647     return result;
00648 }
00649 
00650 bool Skinner::entity_deletable( EntityHandle entity )
00651 {
00652     unsigned char deletable = 0;
00653     ErrorCode result        = thisMB->tag_get_data( mDeletableMBTag, &entity, 1, &deletable );
00654     assert( MB_SUCCESS == result );
00655     if( MB_SUCCESS == result && deletable == 1 ) return false;
00656     return true;
00657 }
00658 
00659 ErrorCode Skinner::classify_2d_boundary( const Range& boundary,
00660                                          const Range& bar_elements,
00661                                          EntityHandle boundary_edges,
00662                                          EntityHandle inferred_edges,
00663                                          EntityHandle non_manifold_edges,
00664                                          EntityHandle other_edges,
00665                                          int& number_boundary_nodes )
00666 {
00667     Range bedges, iedges, nmedges, oedges;
00668     ErrorCode result =
00669         classify_2d_boundary( boundary, bar_elements, bedges, iedges, nmedges, oedges, number_boundary_nodes );MB_CHK_ERR( result );
00670 
00671     // now set the input meshsets to the output ranges
00672     result = thisMB->clear_meshset( &boundary_edges, 1 );MB_CHK_ERR( result );
00673     result = thisMB->add_entities( boundary_edges, bedges );MB_CHK_ERR( result );
00674 
00675     result = thisMB->clear_meshset( &inferred_edges, 1 );MB_CHK_ERR( result );
00676     result = thisMB->add_entities( inferred_edges, iedges );MB_CHK_ERR( result );
00677 
00678     result = thisMB->clear_meshset( &non_manifold_edges, 1 );MB_CHK_ERR( result );
00679     result = thisMB->add_entities( non_manifold_edges, nmedges );MB_CHK_ERR( result );
00680 
00681     result = thisMB->clear_meshset( &other_edges, 1 );MB_CHK_ERR( result );
00682     result = thisMB->add_entities( other_edges, oedges );MB_CHK_ERR( result );
00683 
00684     return MB_SUCCESS;
00685 }
00686 
00687 ErrorCode Skinner::classify_2d_boundary( const Range& boundary,
00688                                          const Range& bar_elements,
00689                                          Range& boundary_edges,
00690                                          Range& inferred_edges,
00691                                          Range& non_manifold_edges,
00692                                          Range& other_edges,
00693                                          int& number_boundary_nodes )
00694 {
00695 
00696     // clear out the edge lists
00697 
00698     boundary_edges.clear();
00699     inferred_edges.clear();
00700     non_manifold_edges.clear();
00701     other_edges.clear();
00702 
00703     number_boundary_nodes = 0;
00704 
00705     // make sure we have something to work with
00706     if( boundary.empty() )
00707     {
00708         return MB_FAILURE;
00709     }
00710 
00711     // get our working dimensions
00712     EntityType type      = thisMB->type_from_handle( *( boundary.begin() ) );
00713     const int source_dim = CN::Dimension( type );
00714 
00715     // make sure we can handle the working dimensions
00716     if( source_dim != 2 )
00717     {
00718         return MB_FAILURE;
00719     }
00720     mTargetDim = source_dim - 1;
00721 
00722     // initialize
00723     initialize();
00724 
00725     // additional initialization for this routine
00726     // define a tag for MBEDGE which counts the occurrences of the edge below
00727     // default should be 0 for existing edges, if any
00728 
00729     Tag count_tag;
00730     int default_count = 0;
00731     ErrorCode result =
00732         thisMB->tag_get_handle( 0, 1, MB_TYPE_INTEGER, count_tag, MB_TAG_DENSE | MB_TAG_CREAT, &default_count );MB_CHK_ERR( result );
00733 
00734     Range::const_iterator iter, end_iter;
00735     end_iter = boundary.end();
00736 
00737     std::vector< EntityHandle > conn;
00738     EntityHandle sub_conn[2];
00739     EntityHandle match;
00740 
00741     Range edge_list;
00742     Range boundary_nodes;
00743     Skinner::direction direct;
00744 
00745     EntityType sub_type;
00746     int num_edge, num_sub_ent_vert;
00747     const short* edge_verts;
00748 
00749     // now, process each entity in the boundary
00750 
00751     for( iter = boundary.begin(); iter != end_iter; ++iter )
00752     {
00753         // get the connectivity of this entity
00754         conn.clear();
00755         result = thisMB->get_connectivity( &( *iter ), 1, conn, false );
00756         assert( MB_SUCCESS == result );
00757 
00758         // add node handles to boundary_node range
00759         std::copy( conn.begin(), conn.begin() + CN::VerticesPerEntity( type ), range_inserter( boundary_nodes ) );
00760 
00761         type = thisMB->type_from_handle( *iter );
00762 
00763         // get connectivity of each n-1 dimension entity (edge in this case)
00764         const struct CN::ConnMap* conn_map = &( CN::mConnectivityMap[type][0] );
00765         num_edge                           = CN::NumSubEntities( type, 1 );
00766         for( int i = 0; i < num_edge; i++ )
00767         {
00768             edge_verts = CN::SubEntityVertexIndices( type, 1, i, sub_type, num_sub_ent_vert );
00769             assert( sub_type == MBEDGE && num_sub_ent_vert == 2 );
00770             sub_conn[0]       = conn[edge_verts[0]];
00771             sub_conn[1]       = conn[edge_verts[1]];
00772             int num_sub_nodes = conn_map->num_corners_per_sub_element[i];
00773 
00774             // see if we can match this connectivity with
00775             // an existing entity
00776             find_match( MBEDGE, sub_conn, num_sub_nodes, match, direct );
00777 
00778             // if there is no match, create a new entity
00779             if( match == 0 )
00780             {
00781                 EntityHandle tmphndl = 0;
00782                 int indices[MAX_SUB_ENTITY_VERTICES];
00783                 EntityType new_type;
00784                 int num_new_nodes;
00785                 CN::SubEntityNodeIndices( type, conn.size(), 1, i, new_type, num_new_nodes, indices );
00786                 for( int j = 0; j < num_new_nodes; j++ )
00787                     sub_conn[j] = conn[indices[j]];
00788 
00789                 result = thisMB->create_element( new_type, sub_conn, num_new_nodes, tmphndl );
00790                 assert( MB_SUCCESS == result );
00791                 add_adjacency( tmphndl, sub_conn, num_sub_nodes );
00792                 // target_entities.insert(tmphndl);
00793                 edge_list.insert( tmphndl );
00794                 int count;
00795                 result = thisMB->tag_get_data( count_tag, &tmphndl, 1, &count );
00796                 assert( MB_SUCCESS == result );
00797                 count++;
00798                 result = thisMB->tag_set_data( count_tag, &tmphndl, 1, &count );
00799                 assert( MB_SUCCESS == result );
00800             }
00801             else
00802             {
00803                 // We found a match, we must increment the count on the match
00804                 int count;
00805                 result = thisMB->tag_get_data( count_tag, &match, 1, &count );
00806                 assert( MB_SUCCESS == result );
00807                 count++;
00808                 result = thisMB->tag_set_data( count_tag, &match, 1, &count );
00809                 assert( MB_SUCCESS == result );
00810 
00811                 // if the entity is not deletable, it was pre-existing in
00812                 // the database.  We therefore may need to add it to the
00813                 // edge_list.  Since it will not hurt the range, we add
00814                 // whether it was added before or not
00815                 if( !entity_deletable( match ) )
00816                 {
00817                     edge_list.insert( match );
00818                 }
00819             }
00820         }
00821     }
00822 
00823     // Any bar elements in the model should be classified separately
00824     // If the element is in the skin edge_list, then it should be put in
00825     // the non-manifold edge list.  Edges not in the edge_list are stand-alone
00826     // bars, and we make them simply boundary elements
00827 
00828     if( !bar_elements.empty() )
00829     {
00830         Range::iterator bar_iter;
00831         for( iter = bar_elements.begin(); iter != bar_elements.end(); ++iter )
00832         {
00833             EntityHandle handle = *iter;
00834             bar_iter            = edge_list.find( handle );
00835             if( bar_iter != edge_list.end() )
00836             {
00837                 // it is in the list, erase it and put in non-manifold list
00838                 edge_list.erase( bar_iter );
00839                 non_manifold_edges.insert( handle );
00840             }
00841             else
00842             {
00843                 // not in the edge list, make it a boundary edge
00844                 boundary_edges.insert( handle );
00845             }
00846         }
00847     }
00848 
00849     // now all edges should be classified.  Go through the edge_list,
00850     // and put all in the appropriate lists
00851 
00852     Range::iterator edge_iter, edge_end_iter;
00853     edge_end_iter = edge_list.end();
00854     int count;
00855     for( edge_iter = edge_list.begin(); edge_iter != edge_end_iter; ++edge_iter )
00856     {
00857         // check the count_tag
00858         result = thisMB->tag_get_data( count_tag, &( *edge_iter ), 1, &count );
00859         assert( MB_SUCCESS == result );
00860         if( count == 1 )
00861         {
00862             boundary_edges.insert( *edge_iter );
00863         }
00864         else if( count == 2 )
00865         {
00866             other_edges.insert( *edge_iter );
00867         }
00868         else
00869         {
00870             non_manifold_edges.insert( *edge_iter );
00871         }
00872     }
00873 
00874     // find the inferred edges from the other_edge_list
00875 
00876     double min_angle_degrees = 20.0;
00877     find_inferred_edges( const_cast< Range& >( boundary ), other_edges, inferred_edges, min_angle_degrees );
00878 
00879     // we now want to remove the inferred_edges from the other_edges
00880 
00881     Range temp_range;
00882 
00883     std::set_difference( other_edges.begin(), other_edges.end(), inferred_edges.begin(), inferred_edges.end(),
00884                          range_inserter( temp_range ), std::less< EntityHandle >() );
00885 
00886     other_edges = temp_range;
00887 
00888     // get rid of count tag and deinitialize
00889 
00890     result = thisMB->tag_delete( count_tag );
00891     assert( MB_SUCCESS == result );
00892     deinitialize();
00893 
00894     // set the node count
00895     number_boundary_nodes = boundary_nodes.size();
00896 
00897     return MB_SUCCESS;
00898 }
00899 
00900 void Skinner::find_inferred_edges( Range& skin_boundary,
00901                                    Range& candidate_edges,
00902                                    Range& inferred_edges,
00903                                    double reference_angle_degrees )
00904 {
00905 
00906     // mark all the entities in the skin boundary
00907     Tag mark_tag;
00908     ErrorCode result = thisMB->tag_get_handle( 0, 1, MB_TYPE_BIT, mark_tag, MB_TAG_CREAT );
00909     assert( MB_SUCCESS == result );
00910     unsigned char bit = true;
00911     result            = thisMB->tag_clear_data( mark_tag, skin_boundary, &bit );
00912     assert( MB_SUCCESS == result );
00913 
00914     // find the cosine of the reference angle
00915 
00916     double reference_cosine = cos( reference_angle_degrees * SKINNER_PI / 180.0 );
00917 
00918     // check all candidate edges for an angle greater than the minimum
00919 
00920     Range::iterator iter, end_iter = candidate_edges.end();
00921     std::vector< EntityHandle > adjacencies;
00922     std::vector< EntityHandle >::iterator adj_iter;
00923     EntityHandle face[2];
00924 
00925     for( iter = candidate_edges.begin(); iter != end_iter; ++iter )
00926     {
00927 
00928         // get the 2D elements connected to this edge
00929         adjacencies.clear();
00930         result = thisMB->get_adjacencies( &( *iter ), 1, 2, false, adjacencies );
00931         if( MB_SUCCESS != result ) continue;
00932 
00933         // there should be exactly two, that is why the edge is classified as nonBoundary
00934         // and manifold
00935 
00936         int faces_found = 0;
00937         for( adj_iter = adjacencies.begin(); adj_iter != adjacencies.end() && faces_found < 2; ++adj_iter )
00938         {
00939             // we need to find two of these which are in the skin
00940             unsigned char is_marked = 0;
00941             result                  = thisMB->tag_get_data( mark_tag, &( *adj_iter ), 1, &is_marked );
00942             assert( MB_SUCCESS == result );
00943             if( is_marked )
00944             {
00945                 face[faces_found] = *adj_iter;
00946                 faces_found++;
00947             }
00948         }
00949 
00950         //    assert(faces_found == 2 || faces_found == 0);
00951         if( 2 != faces_found ) continue;
00952 
00953         // see if the two entities have a sufficient angle
00954 
00955         if( has_larger_angle( face[0], face[1], reference_cosine ) )
00956         {
00957             inferred_edges.insert( *iter );
00958         }
00959     }
00960 
00961     result = thisMB->tag_delete( mark_tag );
00962     assert( MB_SUCCESS == result );
00963 }
00964 
00965 bool Skinner::has_larger_angle( EntityHandle& entity1, EntityHandle& entity2, double reference_angle_cosine )
00966 {
00967     // compare normals to get angle.  We assume that the surface quads
00968     // which we test here will be approximately planar
00969 
00970     double norm[2][3];
00971     Util::normal( thisMB, entity1, norm[0][0], norm[0][1], norm[0][2] );
00972     Util::normal( thisMB, entity2, norm[1][0], norm[1][1], norm[1][2] );
00973 
00974     double cosine = norm[0][0] * norm[1][0] + norm[0][1] * norm[1][1] + norm[0][2] * norm[1][2];
00975 
00976     if( cosine < reference_angle_cosine )
00977     {
00978         return true;
00979     }
00980 
00981     return false;
00982 }
00983 
00984 // get skin entities of prescribed dimension
00985 ErrorCode Skinner::find_skin( const EntityHandle this_set,
00986                               const Range& entities,
00987                               int dim,
00988                               Range& skin_entities,
00989                               bool create_vert_elem_adjs,
00990                               bool create_skin_elements )
00991 {
00992     Range tmp_skin;
00993     ErrorCode result =
00994         find_skin( this_set, entities, ( dim == 0 ), tmp_skin, 0, create_vert_elem_adjs, create_skin_elements );
00995     if( MB_SUCCESS != result || tmp_skin.empty() ) return result;
00996 
00997     if( tmp_skin.all_of_dimension( dim ) )
00998     {
00999         if( skin_entities.empty() )
01000             skin_entities.swap( tmp_skin );
01001         else
01002             skin_entities.merge( tmp_skin );
01003     }
01004     else
01005     {
01006         result = thisMB->get_adjacencies( tmp_skin, dim, create_skin_elements, skin_entities, Interface::UNION );MB_CHK_ERR( result );
01007         if( this_set ) result = thisMB->add_entities( this_set, skin_entities );
01008     }
01009 
01010     return result;
01011 }
01012 
01013 ErrorCode Skinner::find_skin_vertices( const EntityHandle this_set,
01014                                        const Range& entities,
01015                                        Range* skin_verts,
01016                                        Range* skin_elems,
01017                                        Range* skin_rev_elems,
01018                                        bool create_skin_elems,
01019                                        bool corners_only )
01020 {
01021     ErrorCode rval;
01022     if( entities.empty() ) return MB_SUCCESS;
01023 
01024     const int dim = CN::Dimension( TYPE_FROM_HANDLE( entities.front() ) );
01025     if( dim < 1 || dim > 3 || !entities.all_of_dimension( dim ) ) return MB_TYPE_OUT_OF_RANGE;
01026 
01027     // are we skinning all entities
01028     size_t count = entities.size();
01029     int num_total;
01030     rval = thisMB->get_number_entities_by_dimension( this_set, dim, num_total );
01031     if( MB_SUCCESS != rval ) return rval;
01032     bool all = ( count == (size_t)num_total );
01033 
01034     // Create a bit tag for fast intersection with input entities range.
01035     // If we're skinning all the entities in the mesh, we really don't
01036     // need the tag.  To save memory, just create it with a default value
01037     // of one and don't set it.  That way MOAB will return 1 for all
01038     // entities.
01039     Tag tag;
01040     char bit = all ? 1 : 0;
01041     rval     = thisMB->tag_get_handle( NULL, 1, MB_TYPE_BIT, tag, MB_TAG_CREAT, &bit );
01042     if( MB_SUCCESS != rval ) return rval;
01043 
01044     // tag all entities in input range
01045     if( !all )
01046     {
01047         std::vector< unsigned char > vect( count, 1 );
01048         rval = thisMB->tag_set_data( tag, entities, &vect[0] );
01049         if( MB_SUCCESS != rval )
01050         {
01051             thisMB->tag_delete( tag );
01052             return rval;
01053         }
01054     }
01055 
01056     switch( dim )
01057     {
01058         case 1:
01059             if( skin_verts )
01060                 rval = find_skin_vertices_1D( tag, entities, *skin_verts );
01061             else if( skin_elems )
01062                 rval = find_skin_vertices_1D( tag, entities, *skin_elems );
01063             else
01064                 rval = MB_SUCCESS;
01065             break;
01066         case 2:
01067             rval = find_skin_vertices_2D( this_set, tag, entities, skin_verts, skin_elems, skin_rev_elems,
01068                                           create_skin_elems, corners_only );
01069             break;
01070         case 3:
01071             rval = find_skin_vertices_3D( this_set, tag, entities, skin_verts, skin_elems, skin_rev_elems,
01072                                           create_skin_elems, corners_only );
01073             break;
01074         default:
01075             rval = MB_TYPE_OUT_OF_RANGE;
01076             break;
01077     }
01078 
01079     thisMB->tag_delete( tag );
01080     return rval;
01081 }
01082 
01083 ErrorCode Skinner::find_skin_vertices_1D( Tag tag, const Range& edges, Range& skin_verts )
01084 {
01085     // This rather simple algorithm is provided for completeness
01086     // (not sure how often one really wants the 'skin' of a chain
01087     // or tangle of edges.)
01088     //
01089     // A vertex is on the skin of the edges if it is contained in exactly
01090     // one of the edges *in the input range*.
01091     //
01092     // This function expects the caller to have tagged all edges in the
01093     // input range with a value of one for the passed bit tag, and all
01094     // other edges with a value of zero.  This allows us to do a faster
01095     // intersection with the input range and the edges adjacent to a vertex.
01096 
01097     ErrorCode rval;
01098     Range::iterator hint = skin_verts.begin();
01099 
01100     // All input entities must be edges.
01101     if( !edges.all_of_dimension( 1 ) ) return MB_TYPE_OUT_OF_RANGE;
01102 
01103     // get all the vertices
01104     Range verts;
01105     rval = thisMB->get_connectivity( edges, verts, true );
01106     if( MB_SUCCESS != rval ) return rval;
01107 
01108     // Test how many edges each input vertex is adjacent to.
01109     std::vector< char > tag_vals;
01110     std::vector< EntityHandle > adj;
01111     int n;
01112     for( Range::const_iterator it = verts.begin(); it != verts.end(); ++it )
01113     {
01114         // get edges adjacent to vertex
01115         adj.clear();
01116         rval = thisMB->get_adjacencies( &*it, 1, 1, false, adj );
01117         if( MB_SUCCESS != rval ) return rval;
01118         if( adj.empty() ) continue;
01119 
01120         // Intersect adjacent edges with the input list of edges
01121         tag_vals.resize( adj.size() );
01122         rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
01123         if( MB_SUCCESS != rval ) return rval;
01124 #ifdef MOAB_OLD_STD_COUNT
01125         n = 0;
01126         std::count( tag_vals.begin(), tag_vals.end(), '\001' );
01127 #else
01128         n = std::count( tag_vals.begin(), tag_vals.end(), '\001' );
01129 #endif
01130         // If adjacent to only one input edge, then vertex is on skin
01131         if( n == 1 )
01132         {
01133             hint = skin_verts.insert( hint, *it );
01134         }
01135     }
01136 
01137     return MB_SUCCESS;
01138 }
01139 
01140 // A Container for storing a list of element sides adjacent
01141 // to a vertex.  The template parameter is the number of
01142 // corners for the side.
01143 template < unsigned CORNERS >
01144 class AdjSides
01145 {
01146   public:
01147     /**
01148      * This struct is used to for a reduced representation of element
01149      * "sides" adjacent to a give vertex.  As such, it
01150      * a) does not store the initial vertex that all sides are adjacent to
01151      * b) orders the remaining vertices in a specific way for fast comparison.
01152      *
01153      * For edge elements, only the opposite vertex is stored.
01154      * For triangle elements, only the other two vertices are stored,
01155      *   and they are stored with the smaller of those two handles first.
01156      * For quad elements, only the other three vertices are stored.
01157      *  They are stored such that the vertex opposite the implicit (not
01158      *  stored) vertex is always in slot 1.  The other two vertices
01159      *  (in slots 0 and 2) are ordered such that the handle of the one in
01160      *  slot 0 is smaller than the handle in slot 2.
01161      *
01162      * For each side, the adj_elem field is used to store the element that
01163      * it is a side of as long as the element is considered to be on the skin.
01164      * The adj_elem field is cleared (set to zero) to indicate that this
01165      * side is no longer considered to be on the skin (and is the side of
01166      * more than one element.)
01167      */
01168     struct Side
01169     {
01170         EntityHandle handles[CORNERS - 1];  //!< side vertices, except for implicit one
01171         EntityHandle adj_elem;              //!< element that this is a side of, or zero
01172         bool skin() const
01173         {
01174             return 0 != adj_elem;
01175         }
01176 
01177         /** construct from connectivity of side
01178          *\param array The connectivity of the element side.
01179          *\param idx   The index of the implicit vertex (contained
01180          *             in all sides in the list.)
01181          *\param adj   The element that this is a side of.
01182          */
01183         Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short ) : adj_elem( adj )
01184         {
01185             switch( CORNERS )
01186             {
01187                 case 3:
01188                     handles[1] = array[( idx + 2 ) % CORNERS];
01189                     // fall through
01190                 case 2:
01191                     if( 3 == CORNERS ) handles[1] = array[( idx + 2 ) % CORNERS];
01192                     if( 2 <= CORNERS ) handles[0] = array[( idx + 1 ) % CORNERS];
01193                     break;
01194                 default:
01195                     assert( false );
01196                     break;
01197             }
01198             if( CORNERS == 3 && handles[1] > handles[0] ) std::swap( handles[0], handles[1] );
01199         }
01200 
01201         /** construct from connectivity of parent element
01202          *\param array The connectivity of the parent element
01203          *\param idx   The index of the implicit vertex (contained
01204          *             in all sides in the list.)  This is an index
01205          *             into 'indices', not 'array'.
01206          *\param adj   The element that this is a side of.
01207          *\param indices  The indices into 'array' at which the vertices
01208          *             representing the side occur.
01209          */
01210         Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short, const short* indices )
01211             : adj_elem( adj )
01212         {
01213             switch( CORNERS )
01214             {
01215                 case 3:
01216                     handles[1] = array[indices[( idx + 2 ) % CORNERS]];
01217                     // fall through
01218                 case 2:
01219                     if( 3 == CORNERS ) handles[1] = array[indices[( idx + 2 ) % CORNERS]];
01220                     if( 2 <= CORNERS ) handles[0] = array[indices[( idx + 1 ) % CORNERS]];
01221                     break;
01222                 default:
01223                     assert( false );
01224                     break;
01225             }
01226             if( CORNERS == 3 && handles[1] > handles[0] ) std::swap( handles[0], handles[1] );
01227         }
01228 
01229         // Compare two side instances.  Relies in the ordering of
01230         // vertex handles as described above.
01231         bool operator==( const Side& other ) const
01232         {
01233             switch( CORNERS )
01234             {
01235                 case 3:
01236                     return handles[0] == other.handles[0] && handles[1] == other.handles[1];
01237                 case 2:
01238                     return handles[0] == other.handles[0];
01239                 default:
01240                     assert( false );
01241                     return false;
01242             }
01243         }
01244     };
01245 
01246   private:
01247     std::vector< Side > data;  //!< List of sides
01248     size_t skin_count;         //!< Cached count of sides that are skin
01249 
01250   public:
01251     typedef typename std::vector< Side >::iterator iterator;
01252     typedef typename std::vector< Side >::const_iterator const_iterator;
01253     const_iterator begin() const
01254     {
01255         return data.begin();
01256     }
01257     const_iterator end() const
01258     {
01259         return data.end();
01260     }
01261 
01262     void clear()
01263     {
01264         data.clear();
01265         skin_count = 0;
01266     }
01267     bool empty() const
01268     {
01269         return data.empty();
01270     }
01271 
01272     AdjSides() : skin_count( 0 ) {}
01273 
01274     size_t num_skin() const
01275     {
01276         return skin_count;
01277     }
01278 
01279     /** \brief insert side, specifying side connectivity
01280      *
01281      * Either insert a new side, or if the side is already in the
01282      * list, mark it as not on the skin.
01283      *
01284      *\param handles The connectivity of the element side.
01285      *\param skip_idx The index of the implicit vertex (contained
01286      *             in all sides in the list.)
01287      *\param adj_elem The element that this is a side of.
01288      *\param elem_side Which side of adj_elem are we storing
01289      *             (CN side number.)
01290      */
01291     void insert( const EntityHandle* handles, int skip_idx, EntityHandle adj_elem, unsigned short elem_side )
01292     {
01293         Side side( handles, skip_idx, adj_elem, elem_side );
01294         iterator p = std::find( data.begin(), data.end(), side );
01295         if( p == data.end() )
01296         {
01297             data.push_back( side );
01298             ++skin_count;  // not in list yet, so skin side (so far)
01299         }
01300         else if( p->adj_elem )
01301         {
01302             p->adj_elem = 0;  // mark as not on skin
01303             --skin_count;     // decrement cached count of skin elements
01304         }
01305     }
01306 
01307     /** \brief insert side, specifying list of indices into parent element
01308      * connectivity.
01309      *
01310      * Either insert a new side, or if the side is already in the
01311      * list, mark it as not on the skin.
01312      *
01313      *\param handles The connectivity of the parent element
01314      *\param skip_idx The index of the implicit vertex (contained
01315      *             in all sides in the list.)  This is an index
01316      *             into 'indices', not 'handles'.
01317      *\param adj_elem The element that this is a side of (parent handle).
01318      *\param indices  The indices into 'handles' at which the vertices
01319      *             representing the side occur.
01320      *\param elem_side Which side of adj_elem are we storing
01321      *             (CN side number.)
01322      */
01323     void insert( const EntityHandle* handles,
01324                  int skip_idx,
01325                  EntityHandle adj_elem,
01326                  unsigned short elem_side,
01327                  const short* indices )
01328     {
01329         Side side( handles, skip_idx, adj_elem, elem_side, indices );
01330         iterator p = std::find( data.begin(), data.end(), side );
01331         if( p == data.end() )
01332         {
01333             data.push_back( side );
01334             ++skin_count;  // not in list yet, so skin side (so far)
01335         }
01336         else if( p->adj_elem )
01337         {
01338             p->adj_elem = 0;  // mark as not on skin
01339             --skin_count;     // decrement cached count of skin elements
01340         }
01341     }
01342 
01343     /**\brief Search list for a given side, and if found, mark as not skin.
01344      *
01345      *\param other   Connectivity of side
01346      *\param skip_index Index in 'other' at which implicit vertex occurs.
01347      *\param elem_out If return value is true, the element that the side is a
01348      *                side of.  If return value is false, not set.
01349      *\return true if found and marked not-skin, false if not found.
01350      *
01351      * Given the connectivity of some existing element, check if it occurs
01352      * in the list.  If it does, clear the "is skin" state of the side so
01353      * that we know that we don't need to later create the side element.
01354      */
01355     bool find_and_unmark( const EntityHandle* other, int skip_index, EntityHandle& elem_out )
01356     {
01357         Side s( other, skip_index, 0, 0 );
01358         iterator p = std::find( data.begin(), data.end(), s );
01359         if( p == data.end() || !p->adj_elem )
01360             return false;
01361         else
01362         {
01363             elem_out    = p->adj_elem;
01364             p->adj_elem = 0;  // clear "is skin" state for side
01365             --skin_count;     // decrement cached count of skin sides
01366             return true;
01367         }
01368     }
01369 };
01370 
01371 /** construct from connectivity of side
01372  *\param array The connectivity of the element side.
01373  *\param idx   The index of the implicit vertex (contained
01374  *             in all sides in the list.)
01375  *\param adj   The element that this is a side of.
01376  */
01377 template <>
01378 AdjSides< 4 >::Side::Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short ) : adj_elem( adj )
01379 {
01380     const unsigned int CORNERS = 4;
01381     handles[2]                 = array[( idx + 3 ) % CORNERS];
01382     handles[1]                 = array[( idx + 2 ) % CORNERS];
01383     handles[0]                 = array[( idx + 1 ) % CORNERS];
01384     if( handles[2] > handles[0] ) std::swap( handles[0], handles[2] );
01385 }
01386 
01387 /** construct from connectivity of parent element
01388  *\param array The connectivity of the parent element
01389  *\param idx   The index of the implicit vertex (contained
01390  *             in all sides in the list.)  This is an index
01391  *             into 'indices', not 'array'.
01392  *\param adj   The element that this is a side of.
01393  *\param indices  The indices into 'array' at which the vertices
01394  *             representing the side occur.
01395  */
01396 template <>
01397 AdjSides< 4 >::Side::Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short, const short* indices )
01398     : adj_elem( adj )
01399 {
01400     const unsigned int CORNERS = 4;
01401     handles[2]                 = array[indices[( idx + 3 ) % CORNERS]];
01402     handles[1]                 = array[indices[( idx + 2 ) % CORNERS]];
01403     handles[0]                 = array[indices[( idx + 1 ) % CORNERS]];
01404     if( handles[2] > handles[0] ) std::swap( handles[0], handles[2] );
01405 }
01406 
01407 // Compare two side instances.  Relies in the ordering of
01408 // vertex handles as described above.
01409 template <>
01410 bool AdjSides< 4 >::Side::operator==( const Side& other ) const
01411 {
01412     return handles[0] == other.handles[0] && handles[1] == other.handles[1] && handles[2] == other.handles[2];
01413 }
01414 
01415 // Utility function used by find_skin_vertices_2D and
01416 // find_skin_vertices_3D to create elements representing
01417 // the skin side of a higher-dimension element if one
01418 // does not already exist.
01419 //
01420 // Some arguments may seem redundant, but they are used
01421 // to create the correct order of element when the input
01422 // element contains higher-order nodes.
01423 //
01424 // This function always creates elements that have a "forward"
01425 // orientation with respect to the parent element (have
01426 // nodes ordered the same as CN returns for the "side").
01427 //
01428 // elem - The higher-dimension element for which to create
01429 //        a lower-dim element representing the side.
01430 // side_type - The EntityType of the desired side.
01431 // side_conn - The connectivity of the new side.
01432 ErrorCode Skinner::create_side( const EntityHandle this_set,
01433                                 EntityHandle elem,
01434                                 EntityType side_type,
01435                                 const EntityHandle* side_conn,
01436                                 EntityHandle& side_elem )
01437 {
01438     const int max_side = 9;
01439     const EntityHandle* conn;
01440     int len, side_len, side, sense, offset, indices[max_side];
01441     ErrorCode rval;
01442     EntityType type   = TYPE_FROM_HANDLE( elem ), tmp_type;
01443     const int ncorner = CN::VerticesPerEntity( side_type );
01444     const int d       = CN::Dimension( side_type );
01445     std::vector< EntityHandle > storage;
01446 
01447     // Get the connectivity of the parent element
01448     rval = thisMB->get_connectivity( elem, conn, len, false, &storage );
01449     if( MB_SUCCESS != rval ) return rval;
01450 
01451     // treat separately MBPOLYGON; we want to create the edge in the
01452     // forward sense always ; so figure out the sense first, then get out
01453     if( MBPOLYGON == type && 1 == d && MBEDGE == side_type )
01454     {
01455         // first find the first vertex in the conn list
01456         int i = 0;
01457         for( i = 0; i < len; i++ )
01458         {
01459             if( conn[i] == side_conn[0] ) break;
01460         }
01461         if( len == i ) return MB_FAILURE;  // not found, big error
01462         // now, what if the polygon is padded?
01463         // the previous index is fine always. but the next one could be trouble :(
01464         int prevIndex = ( i + len - 1 ) % len;
01465         int nextIndex = ( i + 1 ) % len;
01466         // if the next index actually point to the same node, as current, it means it is padded
01467         if( conn[nextIndex] == conn[i] )
01468         {
01469             // it really means we are at the end of proper nodes, the last nodes are repeated, so it
01470             // should be the first node
01471             nextIndex = 0;  // this is the first node!
01472         }
01473         EntityHandle conn2[2] = { side_conn[0], side_conn[1] };
01474         if( conn[prevIndex] == side_conn[1] )
01475         {
01476             // reverse, so the edge will be forward
01477             conn2[0] = side_conn[1];
01478             conn2[1] = side_conn[0];
01479         }
01480         else if( conn[nextIndex] != side_conn[1] )
01481             return MB_FAILURE;  // it is not adjacent to the polygon
01482 
01483         rval = thisMB->create_element( MBEDGE, conn2, 2, side_elem );MB_CHK_ERR( rval );
01484         if( this_set )
01485         {
01486             rval = thisMB->add_entities( this_set, &side_elem, 1 );MB_CHK_ERR( rval );
01487         }
01488         return MB_SUCCESS;
01489     }
01490     // Find which side we are creating and get indices of all nodes
01491     // (including higher-order, if any.)
01492     CN::SideNumber( type, conn, side_conn, ncorner, d, side, sense, offset );
01493     CN::SubEntityNodeIndices( type, len, d, side, tmp_type, side_len, indices );
01494     assert( side_len <= max_side );
01495     assert( side_type == tmp_type );
01496 
01497     // NOTE: re-create conn array even when no higher-order nodes
01498     //      because we want it to always be forward with respect
01499     //      to the side ordering.
01500     EntityHandle side_conn_full[max_side];
01501     for( int i = 0; i < side_len; ++i )
01502         side_conn_full[i] = conn[indices[i]];
01503 
01504     rval = thisMB->create_element( side_type, side_conn_full, side_len, side_elem );MB_CHK_ERR( rval );
01505     if( this_set )
01506     {
01507         rval = thisMB->add_entities( this_set, &side_elem, 1 );MB_CHK_ERR( rval );
01508     }
01509     return MB_SUCCESS;
01510     ;
01511 }
01512 
01513 // Test if an edge is reversed with respect CN's ordering
01514 // for the "side" of a face.
01515 bool Skinner::edge_reversed( EntityHandle face, const EntityHandle* edge_ends )
01516 {
01517     const EntityHandle* conn;
01518     int len, idx;
01519     ErrorCode rval = thisMB->get_connectivity( face, conn, len, true );
01520     if( MB_SUCCESS != rval )
01521     {
01522         assert( false );
01523         return false;
01524     }
01525     idx = std::find( conn, conn + len, edge_ends[0] ) - conn;
01526     if( idx == len )
01527     {
01528         assert( false );
01529         return false;
01530     }
01531     return ( edge_ends[1] == conn[( idx + len - 1 ) % len] );
01532 }
01533 
01534 // Test if a 2D element representing the side or face of a
01535 // volume element is reversed with respect to the CN node
01536 // ordering for the corresponding region element side.
01537 bool Skinner::face_reversed( EntityHandle region, const EntityHandle* face_corners, EntityType face_type )
01538 {
01539     const EntityHandle* conn;
01540     int len, side, sense, offset;
01541     ErrorCode rval = thisMB->get_connectivity( region, conn, len, true );
01542     if( MB_SUCCESS != rval )
01543     {
01544         assert( false );
01545         return false;
01546     }
01547     short r = CN::SideNumber( TYPE_FROM_HANDLE( region ), conn, face_corners, CN::VerticesPerEntity( face_type ),
01548                               CN::Dimension( face_type ), side, sense, offset );
01549     assert( 0 == r );
01550     return ( !r && sense == -1 );
01551 }
01552 
01553 ErrorCode Skinner::find_skin_vertices_2D( const EntityHandle this_set,
01554                                           Tag tag,
01555                                           const Range& faces,
01556                                           Range* skin_verts,
01557                                           Range* skin_edges,
01558                                           Range* reversed_edges,
01559                                           bool create_edges,
01560                                           bool corners_only )
01561 {
01562     // This function iterates over all the vertices contained in the
01563     // input face list.  For each such vertex, it then iterates over
01564     // all of the sides of the face elements adjacent to the vertex.
01565     // If an adjacent side is the side of only one of the input
01566     // faces, then that side is on the skin.
01567     //
01568     // This algorithm will visit each skin vertex exactly once.  It
01569     // will visit each skin side once for each vertex in the side.
01570     //
01571     // This function expects the caller to have created the passed bit
01572     // tag and set it to one only for the faces in the passed range.  This
01573     // tag is used to do a fast intersection of the faces adjacent to a
01574     // vertex with the faces in the input range (discard any for which the
01575     // tag is not set to one.)
01576 
01577     ErrorCode rval;
01578     std::vector< EntityHandle >::iterator i, j;
01579     Range::iterator hint;
01580     if( skin_verts ) hint = skin_verts->begin();
01581     std::vector< EntityHandle > storage;
01582     const EntityHandle* conn;
01583     int len;
01584     bool find_edges                      = skin_edges || create_edges;
01585     bool printed_nonconformal_ho_warning = false;
01586     EntityHandle face;
01587 
01588     if( !faces.all_of_dimension( 2 ) ) return MB_TYPE_OUT_OF_RANGE;
01589 
01590     // get all the vertices
01591     Range verts;
01592     rval = thisMB->get_connectivity( faces, verts, true );
01593     if( MB_SUCCESS != rval ) return rval;
01594 
01595     std::vector< char > tag_vals;
01596     std::vector< EntityHandle > adj;
01597     AdjSides< 2 > adj_edges;
01598     for( Range::const_iterator it = verts.begin(); it != verts.end(); ++it )
01599     {
01600         bool higher_order = false;
01601 
01602         // get all adjacent faces
01603         adj.clear();
01604         rval = thisMB->get_adjacencies( &*it, 1, 2, false, adj );
01605         if( MB_SUCCESS != rval ) return rval;
01606         if( adj.empty() ) continue;
01607 
01608         // remove those not in the input list (intersect with input list)
01609         i = j = adj.begin();
01610         tag_vals.resize( adj.size() );
01611         rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
01612         if( MB_SUCCESS != rval ) return rval;
01613         // remove non-tagged entries
01614         i = j = adj.begin();
01615         for( ; i != adj.end(); ++i )
01616             if( tag_vals[i - adj.begin()] ) *( j++ ) = *i;
01617         adj.erase( j, adj.end() );
01618 
01619         // For each adjacent face, check the edges adjacent to the current vertex
01620         adj_edges.clear();  // other vertex for adjacent edges
01621         for( i = adj.begin(); i != adj.end(); ++i )
01622         {
01623             rval = thisMB->get_connectivity( *i, conn, len, false, &storage );
01624             if( MB_SUCCESS != rval ) return rval;
01625 
01626             // For a single face element adjacent to this vertex, there
01627             // will be exactly two sides (edges) adjacent to the vertex.
01628             // Find the other vertex for each of the two edges.
01629 
01630             EntityHandle prev, next;  // vertices of two adjacent edge-sides
01631             const int idx = std::find( conn, conn + len, *it ) - conn;
01632             assert( idx != len );
01633 
01634             if( TYPE_FROM_HANDLE( *i ) == MBTRI && len > 3 )
01635             {
01636                 len          = 3;
01637                 higher_order = true;
01638                 if( idx > 2 )
01639                 {  // skip higher-order nodes for now
01640                     if( !printed_nonconformal_ho_warning )
01641                     {
01642                         printed_nonconformal_ho_warning = true;
01643                         std::cerr << "Non-conformal higher-order mesh detected in skinner: "
01644                                   << "vertex " << ID_FROM_HANDLE( *it ) << " is a corner in "
01645                                   << "some elements and a higher-order node in others" << std::endl;
01646                     }
01647                     continue;
01648                 }
01649             }
01650             else if( TYPE_FROM_HANDLE( *i ) == MBQUAD && len > 4 )
01651             {
01652                 len          = 4;
01653                 higher_order = true;
01654                 if( idx > 3 )
01655                 {  // skip higher-order nodes for now
01656                     if( !printed_nonconformal_ho_warning )
01657                     {
01658                         printed_nonconformal_ho_warning = true;
01659                         std::cerr << "Non-conformal higher-order mesh detected in skinner: "
01660                                   << "vertex " << ID_FROM_HANDLE( *it ) << " is a corner in "
01661                                   << "some elements and a higher-order node in others" << std::endl;
01662                     }
01663                     continue;
01664                 }
01665             }
01666 
01667             // so it must be a MBPOLYGON
01668             const int prev_idx = ( idx + len - 1 ) % len;  // this should be fine, always, even for padded case
01669             prev               = conn[prev_idx];
01670             next               = conn[( idx + 1 ) % len];
01671             if( next == conn[idx] )  // it must be the padded case, so roll to the beginning
01672                 next = conn[0];
01673 
01674             // Insert sides (edges) in our list of candidate skin sides
01675             adj_edges.insert( &prev, 1, *i, prev_idx );
01676             adj_edges.insert( &next, 1, *i, idx );
01677         }
01678 
01679         // If vertex is not on skin, advance to next vertex.
01680         // adj_edges handled checking for duplicates on insertion.
01681         // If every candidate skin edge occurred more than once (was
01682         // not in fact on the skin), then we're done with this vertex.
01683         if( 0 == adj_edges.num_skin() ) continue;
01684 
01685         // If user requested Range of *vertices* on the skin...
01686         if( skin_verts )
01687         {
01688             // Put skin vertex in output list
01689             hint = skin_verts->insert( hint, *it );
01690 
01691             // Add mid edge nodes to vertex list
01692             if( !corners_only && higher_order )
01693             {
01694                 for( AdjSides< 2 >::const_iterator p = adj_edges.begin(); p != adj_edges.end(); ++p )
01695                 {
01696                     if( p->skin() )
01697                     {
01698                         face            = p->adj_elem;
01699                         EntityType type = TYPE_FROM_HANDLE( face );
01700 
01701                         rval = thisMB->get_connectivity( face, conn, len, false );
01702                         if( MB_SUCCESS != rval ) return rval;
01703                         if( !CN::HasMidEdgeNodes( type, len ) ) continue;
01704 
01705                         EntityHandle ec[2] = { *it, p->handles[0] };
01706                         int side, sense, offset;
01707                         CN::SideNumber( type, conn, ec, 2, 1, side, sense, offset );
01708                         offset = CN::HONodeIndex( type, len, 1, side );
01709                         assert( offset >= 0 && offset < len );
01710                         skin_verts->insert( conn[offset] );
01711                     }
01712                 }
01713             }
01714         }
01715 
01716         // If user requested Range of *edges* on the skin...
01717         if( find_edges )
01718         {
01719             // Search list of existing adjacent edges for any that are on the skin
01720             adj.clear();
01721             rval = thisMB->get_adjacencies( &*it, 1, 1, false, adj );
01722             if( MB_SUCCESS != rval ) return rval;
01723             for( i = adj.begin(); i != adj.end(); ++i )
01724             {
01725                 rval = thisMB->get_connectivity( *i, conn, len, true );
01726                 if( MB_SUCCESS != rval ) return rval;
01727 
01728                 // bool equality expression within find_and_unmark call
01729                 // will be evaluate to the index of *it in the conn array.
01730                 //
01731                 // Note that the order of the terms in the if statement is important.
01732                 // We want to unmark any existing skin edges even if we aren't
01733                 // returning them.  Otherwise we'll end up creating duplicates
01734                 // if create_edges is true and skin_edges is not.
01735                 if( adj_edges.find_and_unmark( conn, ( conn[1] == *it ), face ) && skin_edges )
01736                 {
01737                     if( reversed_edges && edge_reversed( face, conn ) )
01738                         reversed_edges->insert( *i );
01739                     else
01740                         skin_edges->insert( *i );
01741                 }
01742             }
01743         }
01744 
01745         // If the user requested that we create new edges for sides
01746         // on the skin for which there is no existing edge, and there
01747         // are still skin sides for which no corresponding edge was found...
01748         if( create_edges && adj_edges.num_skin() )
01749         {
01750             // Create any skin edges that don't exist
01751             for( AdjSides< 2 >::const_iterator p = adj_edges.begin(); p != adj_edges.end(); ++p )
01752             {
01753                 if( p->skin() )
01754                 {
01755                     EntityHandle edge, ec[] = { *it, p->handles[0] };
01756                     rval = create_side( this_set, p->adj_elem, MBEDGE, ec, edge );
01757                     if( MB_SUCCESS != rval ) return rval;
01758                     if( skin_edges ) skin_edges->insert( edge );
01759                 }
01760             }
01761         }
01762 
01763     }  // end for each vertex
01764 
01765     return MB_SUCCESS;
01766 }
01767 
01768 ErrorCode Skinner::find_skin_vertices_3D( const EntityHandle this_set,
01769                                           Tag tag,
01770                                           const Range& entities,
01771                                           Range* skin_verts,
01772                                           Range* skin_faces,
01773                                           Range* reversed_faces,
01774                                           bool create_faces,
01775                                           bool corners_only )
01776 {
01777     // This function iterates over all the vertices contained in the
01778     // input vol elem list.  For each such vertex, it then iterates over
01779     // all of the sides of the vol elements adjacent to the vertex.
01780     // If an adjacent side is the side of only one of the input
01781     // elements, then that side is on the skin.
01782     //
01783     // This algorithm will visit each skin vertex exactly once.  It
01784     // will visit each skin side once for each vertex in the side.
01785     //
01786     // This function expects the caller to have created the passed bit
01787     // tag and set it to one only for the elements in the passed range.  This
01788     // tag is used to do a fast intersection of the elements adjacent to a
01789     // vertex with the elements in the input range (discard any for which the
01790     // tag is not set to one.)
01791     //
01792     // For each vertex, iterate over each adjacent element.  Construct
01793     // lists of the sides of each adjacent element that contain the vertex.
01794     //
01795     // A list of three-vertex sides is kept for all triangular faces,
01796     // included three-vertex faces of type MBPOLYGON.  Putting polygons
01797     // in the same list ensures that we find polyhedron and non-polyhedron
01798     // elements that are adjacent.
01799     //
01800     // A list of four-vertex sides is kept for all quadrilateral faces,
01801     // including four-vertex faces of type MBPOLYGON.
01802     //
01803     // Sides with more than four vertices must have an explicit MBPOLYGON
01804     // element representing them because MBPOLYHEDRON connectivity is a
01805     // list of faces rather than vertices.  So the third list (vertices>=5),
01806     // need contain only the handle of the face rather than the vertex handles.
01807 
01808     ErrorCode rval;
01809     std::vector< EntityHandle >::iterator i, j;
01810     Range::iterator hint;
01811     if( skin_verts ) hint = skin_verts->begin();
01812     std::vector< EntityHandle > storage, storage2;  // temp storage for conn lists
01813     const EntityHandle *conn, *conn2;
01814     int len, len2;
01815     bool find_faces = skin_faces || create_faces;
01816     int clen, side, sense, offset, indices[9];
01817     EntityType face_type;
01818     EntityHandle elem;
01819     bool printed_nonconformal_ho_warning = false;
01820 
01821     if( !entities.all_of_dimension( 3 ) ) return MB_TYPE_OUT_OF_RANGE;
01822 
01823     // get all the vertices
01824     Range verts;
01825     rval = thisMB->get_connectivity( entities, verts, true );
01826     if( MB_SUCCESS != rval ) return rval;
01827     // if there are polyhedra in the input list, need to make another
01828     // call to get vertices from faces
01829     if( !verts.all_of_dimension( 0 ) )
01830     {
01831         Range::iterator it = verts.upper_bound( MBVERTEX );
01832         Range pfaces;
01833         pfaces.merge( it, verts.end() );
01834         verts.erase( it, verts.end() );
01835         rval = thisMB->get_connectivity( pfaces, verts, true );
01836         if( MB_SUCCESS != rval ) return rval;
01837         assert( verts.all_of_dimension( 0 ) );
01838     }
01839 
01840     AdjSides< 4 > adj_quads;  // 4-node sides adjacent to a vertex
01841     AdjSides< 3 > adj_tris;   // 3-node sides adjacent to a vertex
01842     AdjSides< 2 > adj_poly;   // n-node sides (n>5) adjacent to vertex
01843                               // (must have an explicit polygon, so store
01844                               // polygon handle rather than vertices.)
01845     std::vector< char > tag_vals;
01846     std::vector< EntityHandle > adj;
01847     for( Range::const_iterator it = verts.begin(); it != verts.end(); ++it )
01848     {
01849         bool higher_order = false;
01850 
01851         // get all adjacent elements
01852         adj.clear();
01853         rval = thisMB->get_adjacencies( &*it, 1, 3, false, adj );
01854         if( MB_SUCCESS != rval ) return rval;
01855         if( adj.empty() ) continue;
01856 
01857         // remove those not tagged (intersect with input range)
01858         i = j = adj.begin();
01859         tag_vals.resize( adj.size() );
01860         rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
01861         if( MB_SUCCESS != rval ) return rval;
01862         for( ; i != adj.end(); ++i )
01863             if( tag_vals[i - adj.begin()] ) *( j++ ) = *i;
01864         adj.erase( j, adj.end() );
01865 
01866         // Build lists of sides of 3D element adjacent to the current vertex
01867         adj_quads.clear();  // store three other vertices for each adjacent quad face
01868         adj_tris.clear();   // store two other vertices for each adjacent tri face
01869         adj_poly.clear();   // store handle of each adjacent polygonal face
01870         int idx;
01871         for( i = adj.begin(); i != adj.end(); ++i )
01872         {
01873             const EntityType type = TYPE_FROM_HANDLE( *i );
01874 
01875             // Special case for POLYHEDRA
01876             if( type == MBPOLYHEDRON )
01877             {
01878                 rval = thisMB->get_connectivity( *i, conn, len );
01879                 if( MB_SUCCESS != rval ) return rval;
01880                 for( int k = 0; k < len; ++k )
01881                 {
01882                     rval = thisMB->get_connectivity( conn[k], conn2, len2, true, &storage2 );
01883                     if( MB_SUCCESS != rval ) return rval;
01884                     idx = std::find( conn2, conn2 + len2, *it ) - conn2;
01885                     if( idx == len2 )  // vertex not in this face
01886                         continue;
01887 
01888                     // Treat 3- and 4-vertex faces specially, so that
01889                     // if the mesh contains both elements and polyhedra,
01890                     // we don't miss one type adjacent to the other.
01891                     switch( len2 )
01892                     {
01893                         case 3:
01894                             adj_tris.insert( conn2, idx, *i, k );
01895                             break;
01896                         case 4:
01897                             adj_quads.insert( conn2, idx, *i, k );
01898                             break;
01899                         default:
01900                             adj_poly.insert( conn + k, 1, *i, k );
01901                             break;
01902                     }
01903                 }
01904             }
01905             else
01906             {
01907                 rval = thisMB->get_connectivity( *i, conn, len, false, &storage );
01908                 if( MB_SUCCESS != rval ) return rval;
01909 
01910                 idx = std::find( conn, conn + len, *it ) - conn;
01911                 assert( idx != len );
01912 
01913                 if( len > CN::VerticesPerEntity( type ) )
01914                 {
01915                     higher_order = true;
01916                     // skip higher-order nodes for now
01917                     if( idx >= CN::VerticesPerEntity( type ) )
01918                     {
01919                         if( !printed_nonconformal_ho_warning )
01920                         {
01921                             printed_nonconformal_ho_warning = true;
01922                             std::cerr << "Non-conformal higher-order mesh detected in skinner: "
01923                                       << "vertex " << ID_FROM_HANDLE( *it ) << " is a corner in "
01924                                       << "some elements and a higher-order node in others" << std::endl;
01925                         }
01926                         continue;
01927                     }
01928                 }
01929 
01930                 // For each side of the element...
01931                 const int num_faces = CN::NumSubEntities( type, 2 );
01932                 for( int f = 0; f < num_faces; ++f )
01933                 {
01934                     int num_vtx;
01935                     const short* face_indices = CN::SubEntityVertexIndices( type, 2, f, face_type, num_vtx );
01936                     const short face_idx = std::find( face_indices, face_indices + num_vtx, (short)idx ) - face_indices;
01937                     // skip sides that do not contain vertex from outer loop
01938                     if( face_idx == num_vtx ) continue;  // current vertex not in this face
01939 
01940                     assert( num_vtx <= 4 );  // polyhedra handled above
01941                     switch( face_type )
01942                     {
01943                         case MBTRI:
01944                             adj_tris.insert( conn, face_idx, *i, f, face_indices );
01945                             break;
01946                         case MBQUAD:
01947                             adj_quads.insert( conn, face_idx, *i, f, face_indices );
01948                             break;
01949                         default:
01950                             return MB_TYPE_OUT_OF_RANGE;
01951                     }
01952                 }
01953             }
01954         }  // end for (adj[3])
01955 
01956         // If vertex is not on skin, advance to next vertex
01957         if( 0 == ( adj_tris.num_skin() + adj_quads.num_skin() + adj_poly.num_skin() ) ) continue;
01958 
01959         // If user requested that skin *vertices* be passed back...
01960         if( skin_verts )
01961         {
01962             // Put skin vertex in output list
01963             hint = skin_verts->insert( hint, *it );
01964 
01965             // Add mid-edge and mid-face nodes to vertex list
01966             if( !corners_only && higher_order )
01967             {
01968                 for( AdjSides< 3 >::const_iterator t = adj_tris.begin(); t != adj_tris.end(); ++t )
01969                 {
01970                     if( t->skin() )
01971                     {
01972                         elem            = t->adj_elem;
01973                         EntityType type = TYPE_FROM_HANDLE( elem );
01974 
01975                         rval = thisMB->get_connectivity( elem, conn, len, false );
01976                         if( MB_SUCCESS != rval ) return rval;
01977                         if( !CN::HasMidNodes( type, len ) ) continue;
01978 
01979                         EntityHandle ec[3] = { *it, t->handles[0], t->handles[1] };
01980                         CN::SideNumber( type, conn, ec, 3, 2, side, sense, offset );
01981                         CN::SubEntityNodeIndices( type, len, 2, side, face_type, clen, indices );
01982                         assert( MBTRI == face_type );
01983                         for( int k = 3; k < clen; ++k )
01984                             skin_verts->insert( conn[indices[k]] );
01985                     }
01986                 }
01987                 for( AdjSides< 4 >::const_iterator q = adj_quads.begin(); q != adj_quads.end(); ++q )
01988                 {
01989                     if( q->skin() )
01990                     {
01991                         elem            = q->adj_elem;
01992                         EntityType type = TYPE_FROM_HANDLE( elem );
01993 
01994                         rval = thisMB->get_connectivity( elem, conn, len, false );
01995                         if( MB_SUCCESS != rval ) return rval;
01996                         if( !CN::HasMidNodes( type, len ) ) continue;
01997 
01998                         EntityHandle ec[4] = { *it, q->handles[0], q->handles[1], q->handles[2] };
01999                         CN::SideNumber( type, conn, ec, 4, 2, side, sense, offset );
02000                         CN::SubEntityNodeIndices( type, len, 2, side, face_type, clen, indices );
02001                         assert( MBQUAD == face_type );
02002                         for( int k = 4; k < clen; ++k )
02003                             skin_verts->insert( conn[indices[k]] );
02004                     }
02005                 }
02006             }
02007         }
02008 
02009         // If user requested that we pass back the list of 2D elements
02010         // representing the skin of the mesh...
02011         if( find_faces )
02012         {
02013             // Search list of adjacent faces for any that are on the skin
02014             adj.clear();
02015             rval = thisMB->get_adjacencies( &*it, 1, 2, false, adj );
02016             if( MB_SUCCESS != rval ) return rval;
02017 
02018             for( i = adj.begin(); i != adj.end(); ++i )
02019             {
02020                 rval = thisMB->get_connectivity( *i, conn, len, true );
02021                 if( MB_SUCCESS != rval ) return rval;
02022                 const int idx2 = std::find( conn, conn + len, *it ) - conn;
02023                 if( idx2 >= len )
02024                 {
02025                     assert( printed_nonconformal_ho_warning );
02026                     continue;
02027                 }
02028 
02029                 // Note that the order of the terms in the if statements below
02030                 // is important.  We want to unmark any existing skin faces even
02031                 // if we aren't returning them.  Otherwise we'll end up creating
02032                 // duplicates if create_faces is true.
02033                 if( 3 == len )
02034                 {
02035                     if( adj_tris.find_and_unmark( conn, idx2, elem ) && skin_faces )
02036                     {
02037                         if( reversed_faces && face_reversed( elem, conn, MBTRI ) )
02038                             reversed_faces->insert( *i );
02039                         else
02040                             skin_faces->insert( *i );
02041                     }
02042                 }
02043                 else if( 4 == len )
02044                 {
02045                     if( adj_quads.find_and_unmark( conn, idx2, elem ) && skin_faces )
02046                     {
02047                         if( reversed_faces && face_reversed( elem, conn, MBQUAD ) )
02048                             reversed_faces->insert( *i );
02049                         else
02050                             skin_faces->insert( *i );
02051                     }
02052                 }
02053                 else
02054                 {
02055                     if( adj_poly.find_and_unmark( &*i, 1, elem ) && skin_faces ) skin_faces->insert( *i );
02056                 }
02057             }
02058         }
02059 
02060         // If user does not want use to create new faces representing
02061         // sides for which there is currently no explicit element,
02062         // skip the remaining code and advance the outer loop to the
02063         // next vertex.
02064         if( !create_faces ) continue;
02065 
02066         // Polyhedra always have explicitly defined faces, so
02067         // there is no way we could need to create such a face.
02068         assert( 0 == adj_poly.num_skin() );
02069 
02070         // Create any skin tris that don't exist
02071         if( adj_tris.num_skin() )
02072         {
02073             for( AdjSides< 3 >::const_iterator t = adj_tris.begin(); t != adj_tris.end(); ++t )
02074             {
02075                 if( t->skin() )
02076                 {
02077                     EntityHandle tri, c[3] = { *it, t->handles[0], t->handles[1] };
02078                     rval = create_side( this_set, t->adj_elem, MBTRI, c, tri );
02079                     if( MB_SUCCESS != rval ) return rval;
02080                     if( skin_faces ) skin_faces->insert( tri );
02081                 }
02082             }
02083         }
02084 
02085         // Create any skin quads that don't exist
02086         if( adj_quads.num_skin() )
02087         {
02088             for( AdjSides< 4 >::const_iterator q = adj_quads.begin(); q != adj_quads.end(); ++q )
02089             {
02090                 if( q->skin() )
02091                 {
02092                     EntityHandle quad, c[4] = { *it, q->handles[0], q->handles[1], q->handles[2] };
02093                     rval = create_side( this_set, q->adj_elem, MBQUAD, c, quad );
02094                     if( MB_SUCCESS != rval ) return rval;
02095                     if( skin_faces ) skin_faces->insert( quad );
02096                 }
02097             }
02098         }
02099     }  // end for each vertex
02100 
02101     return MB_SUCCESS;
02102 }
02103 
02104 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines