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