MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /** 00002 * MOAB, a Mesh-Oriented datABase, is a software component for creating, 00003 * storing and accessing finite element mesh data. 00004 * 00005 * Copyright 2004 Sandia Corporation. Under the terms of Contract 00006 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 00007 * retains certain rights in this software. 00008 * 00009 * This library is free software; you can redistribute it and/or 00010 * modify it under the terms of the GNU Lesser General Public 00011 * License as published by the Free Software Foundation; either 00012 * version 2.1 of the License, or (at your option) any later version. 00013 * 00014 */ 00015 00016 #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