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