MOAB: Mesh Oriented datABase
(version 5.4.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 #include "moab/DualTool.hpp" 00017 #include "moab/Range.hpp" 00018 // using Core for call to check_adjacencies 00019 #include "moab/Core.hpp" 00020 #include "Internals.hpp" 00021 #include "MBTagConventions.hpp" 00022 #include "moab/Skinner.hpp" 00023 #include "moab/Core.hpp" 00024 #include "moab/MeshTopoUtil.hpp" 00025 #include "AEntityFactory.hpp" 00026 #include "moab/CN.hpp" 00027 #include <string> 00028 #include <algorithm> 00029 #include <iostream> 00030 #include <sstream> 00031 #include <cassert> 00032 00033 #define RR \ 00034 if( MB_SUCCESS != result ) return result 00035 #define SWAP( a, b ) \ 00036 { \ 00037 EntityHandle tmp_ent = a; \ 00038 ( a ) = b; \ 00039 ( b ) = tmp_ent; \ 00040 } 00041 00042 namespace moab 00043 { 00044 00045 bool debug = false; 00046 bool debug_ap = false; 00047 00048 //! tag name for dual surfaces 00049 const char* DualTool::DUAL_SURFACE_TAG_NAME = "DUAL_SURFACE"; 00050 00051 //! tag name for dual curves 00052 const char* DualTool::DUAL_CURVE_TAG_NAME = "DUAL_CURVE"; 00053 00054 //! tag name for dual cells 00055 const char* DualTool::IS_DUAL_CELL_TAG_NAME = "__IS_DUAL_CELL"; 00056 00057 //! tag name for dual entities 00058 const char* DualTool::DUAL_ENTITY_TAG_NAME = "__DUAL_ENTITY"; 00059 00060 //! tag name for extra dual entities 00061 const char* DualTool::EXTRA_DUAL_ENTITY_TAG_NAME = "__EXTRA_DUAL_ENTITY"; 00062 00063 //! tag name for graphics point 00064 const char* DualTool::DUAL_GRAPHICS_POINT_TAG_NAME = "__DUAL_GRAPHICS_POINT"; 00065 00066 // const int DualTool::GP_SIZE = 20; 00067 00068 DualTool::DualTool( Interface* impl ) : mbImpl( impl ) 00069 { 00070 EntityHandle dum_handle = 0; 00071 ErrorCode result; 00072 00073 result = mbImpl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dualSurfaceTag, 00074 MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle ); 00075 assert( MB_SUCCESS == result ); 00076 00077 result = mbImpl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dualCurveTag, MB_TAG_SPARSE | MB_TAG_CREAT, 00078 &dum_handle ); 00079 assert( MB_SUCCESS == result ); 00080 00081 unsigned int dummy = 0; 00082 result = mbImpl->tag_get_handle( IS_DUAL_CELL_TAG_NAME, 1, MB_TYPE_INTEGER, isDualCellTag, 00083 MB_TAG_SPARSE | MB_TAG_CREAT, &dummy ); 00084 assert( MB_SUCCESS == result ); 00085 00086 result = mbImpl->tag_get_handle( DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, dualEntityTag, 00087 MB_TAG_DENSE | MB_TAG_CREAT, &dum_handle ); 00088 assert( MB_SUCCESS == result ); 00089 00090 result = mbImpl->tag_get_handle( EXTRA_DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, extraDualEntityTag, 00091 MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle ); 00092 assert( MB_SUCCESS == result ); 00093 00094 static const char dum_name[CATEGORY_TAG_SIZE] = { 0 }; 00095 result = mbImpl->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, categoryTag, 00096 MB_TAG_SPARSE | MB_TAG_CREAT, dum_name ); 00097 assert( MB_SUCCESS == result ); 00098 00099 DualTool::GraphicsPoint dum_pt( 0.0, 0.0, 0.0, -1 ); 00100 result = mbImpl->tag_get_handle( DUAL_GRAPHICS_POINT_TAG_NAME, sizeof( DualTool::GraphicsPoint ), MB_TYPE_DOUBLE, 00101 dualGraphicsPointTag, MB_TAG_DENSE | MB_TAG_CREAT | MB_TAG_BYTES, &dum_pt ); 00102 assert( MB_SUCCESS == result ); 00103 00104 globalIdTag = mbImpl->globalId_tag(); 00105 00106 if( MB_SUCCESS == result ) 00107 { 00108 } // empty statement to get rid of warning. 00109 00110 maxHexId = -1; 00111 } 00112 00113 DualTool::~DualTool() {} 00114 00115 //! construct the dual entities for the entire mesh 00116 ErrorCode DualTool::construct_dual( EntityHandle* entities, const int num_entities ) 00117 { 00118 // allocate a dual entity for each primal entity in the mesh, starting 00119 // with highest dimension and working downward; do each dimension in a separate code 00120 // block, since they're all handled slightly differently 00121 00122 Range regions, faces, edges, vertices; 00123 ErrorCode result; 00124 00125 if( NULL == entities || 0 == num_entities ) 00126 { 00127 00128 // first, construct all the aentities, since they're currently needed to 00129 // compute the dual 00130 result = mbImpl->get_entities_by_dimension( 0, 0, vertices ); 00131 if( MB_SUCCESS != result ) return result; 00132 00133 result = MeshTopoUtil( mbImpl ).construct_aentities( vertices ); 00134 if( MB_SUCCESS != result ) return result; 00135 00136 // get all edges, faces and regions now, so we don't need to filter out dual 00137 // entities later 00138 00139 result = mbImpl->get_entities_by_dimension( 0, 1, edges ); 00140 if( MB_SUCCESS != result ) return result; 00141 result = mbImpl->get_entities_by_dimension( 0, 2, faces ); 00142 if( MB_SUCCESS != result ) return result; 00143 result = mbImpl->get_entities_by_dimension( 0, 3, regions ); 00144 if( MB_SUCCESS != result ) return result; 00145 00146 // get the max global id for hexes, we'll need for modification ops 00147 std::vector< int > gid_vec( regions.size() ); 00148 result = mbImpl->tag_get_data( globalId_tag(), regions, &gid_vec[0] ); 00149 if( MB_SUCCESS != result ) return result; 00150 maxHexId = -1; 00151 Range::iterator rit; 00152 unsigned int i; 00153 for( rit = regions.begin(), i = 0; rit != regions.end(); ++rit, i++ ) 00154 { 00155 if( gid_vec[i] > maxHexId && mbImpl->type_from_handle( *rit ) == MBHEX ) maxHexId = gid_vec[i]; 00156 } 00157 } 00158 else 00159 { 00160 // get entities of various dimensions adjacent to these 00161 result = mbImpl->get_adjacencies( entities, num_entities, 0, true, vertices, Interface::UNION ); 00162 if( MB_SUCCESS != result ) return result; 00163 result = mbImpl->get_adjacencies( entities, num_entities, 1, true, edges, Interface::UNION ); 00164 if( MB_SUCCESS != result ) return result; 00165 result = mbImpl->get_adjacencies( entities, num_entities, 2, true, faces, Interface::UNION ); 00166 if( MB_SUCCESS != result ) return result; 00167 result = mbImpl->get_adjacencies( entities, num_entities, 3, true, regions, Interface::UNION ); 00168 if( MB_SUCCESS != result ) return result; 00169 } 00170 00171 Range dual_verts; 00172 result = construct_dual_vertices( regions, dual_verts ); 00173 if( MB_SUCCESS != result || dual_verts.size() != regions.size() ) return result; 00174 if( debug ) std::cout << "Constructed " << dual_verts.size() << " dual vertices." << std::endl; 00175 00176 // don't really need dual edges, but construct 'em anyway 00177 Range dual_edges; 00178 result = construct_dual_edges( faces, dual_edges ); 00179 if( MB_SUCCESS != result || dual_edges.size() != faces.size() ) return result; 00180 if( debug ) std::cout << "Constructed " << dual_edges.size() << " dual edges." << std::endl; 00181 00182 // construct dual faces 00183 Range dual_faces; 00184 result = construct_dual_faces( edges, dual_faces ); 00185 if( MB_SUCCESS != result || dual_faces.size() != edges.size() ) return result; 00186 if( debug ) std::cout << "Constructed " << dual_faces.size() << " dual faces." << std::endl; 00187 00188 // construct dual cells 00189 Range dual_cells; 00190 result = construct_dual_cells( vertices, dual_cells ); 00191 if( MB_SUCCESS != result || dual_cells.size() != vertices.size() ) return result; 00192 if( debug ) std::cout << "Constructed " << dual_cells.size() << " dual cells." << std::endl; 00193 00194 return MB_SUCCESS; 00195 } 00196 00197 ErrorCode DualTool::construct_dual_vertices( const Range& all_regions, Range& dual_ents ) 00198 { 00199 if( all_regions.empty() ) return MB_SUCCESS; 00200 00201 // make sure they're all regions 00202 assert( 3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.begin() ) ) && 00203 3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.rbegin() ) ) ); 00204 00205 Range::const_iterator rit; 00206 EntityHandle dual_ent; 00207 ErrorCode tmp_result = MB_SUCCESS; 00208 ErrorCode result = MB_SUCCESS; 00209 00210 for( rit = all_regions.begin(); rit != all_regions.end(); ++rit ) 00211 { 00212 if( tmp_result != MB_SUCCESS ) result = tmp_result; 00213 00214 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent ); 00215 if( MB_SUCCESS == tmp_result && 0 != dual_ent ) 00216 { 00217 dual_ents.insert( dual_ent ); 00218 continue; 00219 } 00220 else if( MB_SUCCESS != tmp_result ) 00221 continue; 00222 00223 tmp_result = construct_dual_vertex( *rit, dual_ent, false, true ); 00224 if( MB_SUCCESS != tmp_result ) continue; 00225 00226 // save it in the list of new dual ents 00227 dual_ents.insert( dual_ent ); 00228 } 00229 00230 return result; 00231 } 00232 00233 ErrorCode DualTool::construct_dual_vertex( EntityHandle entity, 00234 EntityHandle& dual_ent, 00235 const bool extra, 00236 const bool add_graphics_pt ) 00237 { 00238 // no dual entity; construct one; first need the avg coordinates 00239 unsigned int is_dual = 0x1; 00240 double avg_pos[3]; 00241 ErrorCode result = MeshTopoUtil( mbImpl ).get_average_position( entity, avg_pos ); 00242 if( MB_SUCCESS != result ) return result; 00243 00244 // now construct the new dual entity 00245 result = mbImpl->create_vertex( avg_pos, dual_ent ); 00246 if( MB_SUCCESS != result ) return result; 00247 00248 // tag it indicating it's a dual entity 00249 result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual ); 00250 if( MB_SUCCESS != result ) return result; 00251 00252 // tag the primal entity with its dual entity and vica versa 00253 if( extra ) 00254 result = mbImpl->tag_set_data( extraDualEntity_tag(), &( entity ), 1, &dual_ent ); 00255 else 00256 result = mbImpl->tag_set_data( dualEntity_tag(), &( entity ), 1, &dual_ent ); 00257 if( MB_SUCCESS != result ) return result; 00258 00259 result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( entity ) ); 00260 if( MB_SUCCESS != result ) return result; 00261 00262 if( add_graphics_pt ) 00263 // put a graphics point on that vertex too 00264 result = add_graphics_point( dual_ent, avg_pos ); 00265 00266 return result; 00267 } 00268 00269 ErrorCode DualTool::add_graphics_point( EntityHandle entity, double* avg_pos ) 00270 { 00271 // add a graphics pt, placed at the same position as the vertex 00272 double my_pos[3]; 00273 ErrorCode result; 00274 00275 if( NULL == avg_pos ) 00276 { 00277 result = MeshTopoUtil( mbImpl ).get_average_position( entity, my_pos ); 00278 if( MB_SUCCESS != result ) return result; 00279 } 00280 else 00281 for( int i = 0; i < 3; i++ ) 00282 my_pos[i] = avg_pos[i]; 00283 00284 DualTool::GraphicsPoint dum_pt( my_pos, -1 ); 00285 result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), &entity, 1, &dum_pt ); 00286 return result; 00287 } 00288 00289 ErrorCode DualTool::construct_dual_edges( const Range& all_faces, Range& dual_ents ) 00290 { 00291 if( all_faces.empty() ) return MB_SUCCESS; 00292 00293 // make sure they're all faces 00294 assert( 2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.begin() ) ) && 00295 2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.rbegin() ) ) ); 00296 00297 Range::const_iterator rit; 00298 EntityHandle dual_ent; 00299 unsigned int is_dual = 0x1; 00300 ErrorCode tmp_result = MB_SUCCESS; 00301 ErrorCode result = MB_SUCCESS; 00302 00303 for( rit = all_faces.begin(); rit != all_faces.end(); ++rit ) 00304 { 00305 if( tmp_result != MB_SUCCESS ) result = tmp_result; 00306 00307 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent ); 00308 if( MB_SUCCESS == tmp_result && 0 != dual_ent ) 00309 { 00310 dual_ents.insert( dual_ent ); 00311 continue; 00312 } 00313 00314 // no dual entity; construct one; get the bounding regions 00315 std::vector< EntityHandle > out_ents; 00316 tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, out_ents ); 00317 if( MB_SUCCESS != tmp_result || out_ents.empty() ) continue; 00318 00319 // get the dual vertices 00320 std::vector< EntityHandle > dual_verts( out_ents.size() ); 00321 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &out_ents[0], out_ents.size(), &dual_verts[0] ); 00322 if( MB_SUCCESS != tmp_result ) continue; 00323 assert( dual_verts.size() <= 2 ); 00324 00325 double avg_pos[3]; 00326 bool bdy_face = ( dual_verts.size() == 1 ? true : false ); 00327 if( bdy_face ) 00328 { 00329 // boundary face - make a dual vertex at the face center and put in list 00330 tmp_result = construct_dual_vertex( *rit, dual_ent, true, true ); 00331 00332 // put it on vertex list 00333 dual_verts.push_back( dual_ent ); 00334 } 00335 00336 assert( dual_verts.size() == 2 ); 00337 00338 // now create the dual edge 00339 tmp_result = mbImpl->create_element( MBEDGE, &dual_verts[0], 2, dual_ent ); 00340 if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue; 00341 00342 // save it in the list of new dual ents 00343 dual_ents.insert( dual_ent ); 00344 00345 // tag the primal entity with its dual entity and vica versa 00346 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent ); 00347 if( MB_SUCCESS != tmp_result ) continue; 00348 00349 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) ); 00350 if( MB_SUCCESS != tmp_result ) continue; 00351 00352 // tag the edge indicating it's a dual entity 00353 tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual ); 00354 if( MB_SUCCESS != tmp_result ) continue; 00355 00356 // add a graphics point to the edge; position depends on whether it's a 00357 // bdy face (mid-pt of dual edge) or not (mid-pt of primal face) 00358 if( bdy_face ) 00359 tmp_result = add_graphics_point( dual_ent ); 00360 else 00361 { 00362 // get the face's position 00363 tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos ); 00364 if( MB_SUCCESS != tmp_result ) continue; 00365 tmp_result = add_graphics_point( dual_ent, avg_pos ); 00366 } 00367 if( MB_SUCCESS != tmp_result ) continue; 00368 } 00369 00370 return result; 00371 } 00372 00373 ErrorCode DualTool::construct_dual_faces( const Range& all_edges, Range& dual_ents ) 00374 { 00375 if( all_edges.empty() ) return MB_SUCCESS; 00376 00377 // make sure they're all edges 00378 assert( 1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.begin() ) ) && 00379 1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.rbegin() ) ) ); 00380 00381 Range::const_iterator rit; 00382 EntityHandle dual_ent; 00383 unsigned int is_dual = 0x1; 00384 ErrorCode tmp_result = MB_SUCCESS; 00385 ErrorCode result = MB_SUCCESS; 00386 Range equiv_edges; 00387 #define TRC \ 00388 if( MB_SUCCESS != tmp_result ) \ 00389 { \ 00390 result = tmp_result; \ 00391 continue; \ 00392 } 00393 for( rit = all_edges.begin(); rit != all_edges.end(); ++rit ) 00394 { 00395 00396 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent ); 00397 if( MB_SUCCESS == tmp_result && 0 != dual_ent ) 00398 { 00399 dual_ents.insert( dual_ent ); 00400 continue; 00401 } 00402 00403 // no dual entity; construct one; get the dual vertices bounding the edge in radial order, 00404 // then construct the dual face 00405 std::vector< EntityHandle > rad_dverts; 00406 bool bdy_edge; 00407 tmp_result = get_radial_dverts( *rit, rad_dverts, bdy_edge ); 00408 TRC if( rad_dverts.empty() ) continue; 00409 00410 tmp_result = mbImpl->create_element( MBPOLYGON, &rad_dverts[0], rad_dverts.size(), dual_ent ); 00411 TRC 00412 00413 // tag it indicating it's a dual entity, and tag primal/dual with dual/primal 00414 tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual ); 00415 TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent ); 00416 TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) ); 00417 TRC 00418 00419 // save it in the list of new dual ents 00420 dual_ents.insert( dual_ent ); 00421 00422 // add a graphics point to the cell; position depends on whether it's a 00423 // bdy cell (mid-pt of cell's vertices) or not (mid-pt of primal edge) 00424 double avg_pos[3]; 00425 tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos ); 00426 TRC if( bdy_edge ) 00427 { 00428 00429 // add a new dual edge betw last 2 verts 00430 EntityHandle new_edge; 00431 tmp_result = mbImpl->create_element( MBEDGE, &rad_dverts[rad_dverts.size() - 2], 2, new_edge ); 00432 TRC tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &new_edge, 1, &is_dual ); 00433 TRC 00434 00435 // tag the new dual edge with the primal edge as it's dual entity; primal 00436 // edge IS NOT likewise tagged, since it's already tagged with the 2cell 00437 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &new_edge, 1, &( *rit ) ); 00438 TRC 00439 00440 // add a graphics pt, position is center of primal edge 00441 tmp_result = add_graphics_point( dual_ent ); 00442 TRC tmp_result = add_graphics_point( new_edge, avg_pos ); 00443 TRC 00444 } 00445 00446 else 00447 { 00448 // if inside, point goes on the 2cell, at primal edge mid-pt 00449 tmp_result = add_graphics_point( dual_ent, avg_pos ); 00450 TRC 00451 } 00452 00453 // check to see whether we have equiv entities; if we find any, save for later fixup 00454 Range dum_edges, dum_poly( dual_ent, dual_ent ); 00455 tmp_result = mbImpl->get_adjacencies( dum_poly, 1, false, dum_edges ); 00456 if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result ) 00457 { 00458 // we do - need to add adjacencies to disambiguate; use the primal 00459 equiv_edges.merge( dum_edges ); 00460 } 00461 } 00462 00463 if( !equiv_edges.empty() ) result = check_dual_equiv_edges( equiv_edges ); 00464 00465 return result; 00466 } 00467 00468 ErrorCode DualTool::check_dual_equiv_edges( Range& dual_edges ) 00469 { 00470 // fix equivalent dual edges (i.e. edges whose vertices define multiple edges) 00471 // by explicitly adding adjacencies to containing polygons; adjacent polygons 00472 // found by going through primal 00473 ErrorCode tmp_result, result = MB_SUCCESS; 00474 00475 Range all_dedges( dual_edges ); 00476 // first, go through all dual edges and find equivalent edges (by looking for 00477 // up-adjacent edges on the vertices of each edge) 00478 for( Range::iterator rit = dual_edges.begin(); rit != dual_edges.end(); ++rit ) 00479 { 00480 Range connect, dum_range( *rit, *rit ); 00481 tmp_result = mbImpl->get_adjacencies( dum_range, 0, false, connect ); 00482 if( MB_SUCCESS != tmp_result ) continue; 00483 tmp_result = mbImpl->get_adjacencies( connect, 1, false, all_dedges, Interface::UNION ); 00484 if( MB_SUCCESS != tmp_result ) continue; 00485 } 00486 00487 // save a copy for checking later 00488 Range save_all_2cells; 00489 00490 // go through each edge 00491 while( !all_dedges.empty() ) 00492 { 00493 EntityHandle this_edge = *all_dedges.begin(); 00494 all_dedges.erase( all_dedges.begin() ); 00495 00496 const EntityHandle* connect; 00497 int num_connect; 00498 result = mbImpl->get_connectivity( this_edge, connect, num_connect ); 00499 if( MB_SUCCESS != result ) continue; 00500 00501 Range dum_edges, verts; 00502 verts.insert( connect[0] ); 00503 verts.insert( connect[1] ); 00504 tmp_result = mbImpl->get_adjacencies( verts, 1, false, dum_edges ); 00505 if( MB_SUCCESS != tmp_result ) 00506 { 00507 result = tmp_result; 00508 continue; 00509 } 00510 if( dum_edges.size() == 1 ) 00511 { 00512 // not an equiv edge - already removed from list, so just continue 00513 continue; 00514 } 00515 00516 // ok, have an equiv entity - fix by looking through primal 00517 // pre-get the primal of these 00518 EntityHandle dedge_quad; 00519 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &this_edge, 1, &dedge_quad ); 00520 if( MB_SUCCESS != tmp_result ) 00521 { 00522 result = tmp_result; 00523 continue; 00524 } 00525 00526 if( MBQUAD == mbImpl->type_from_handle( dedge_quad ) ) 00527 { 00528 00529 // get the primal edges adj to quad 00530 Range dum_quad_range( dedge_quad, dedge_quad ), adj_pedges; 00531 tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_pedges ); 00532 if( MB_SUCCESS != tmp_result ) 00533 { 00534 result = tmp_result; 00535 continue; 00536 } 00537 // get the dual 2cells corresponding to those pedges 00538 std::vector< EntityHandle > dcells; 00539 dcells.resize( adj_pedges.size() ); 00540 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), adj_pedges, &dcells[0] ); 00541 if( MB_SUCCESS != tmp_result ) 00542 { 00543 result = tmp_result; 00544 continue; 00545 } 00546 // now add explicit adjacencies from the dedge to those dcells 00547 std::vector< EntityHandle >::iterator vit; 00548 for( vit = dcells.begin(); vit != dcells.end(); ++vit ) 00549 { 00550 save_all_2cells.insert( *vit ); 00551 00552 assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) ); 00553 tmp_result = mbImpl->add_adjacencies( this_edge, &( *vit ), 1, false ); 00554 if( MB_SUCCESS != tmp_result ) 00555 { 00556 result = tmp_result; 00557 continue; 00558 } 00559 // check that there are really adjacencies and *vit is in them 00560 const EntityHandle* adjs; 00561 int num_adjs; 00562 tmp_result = reinterpret_cast< Core* >( mbImpl )->a_entity_factory()->get_adjacencies( this_edge, adjs, 00563 num_adjs ); 00564 if( NULL == adjs || std::find( adjs, adjs + num_adjs, *vit ) == adjs + num_adjs ) 00565 std::cout << "Add_adjacencies failed in construct_dual_faces." << std::endl; 00566 } 00567 } 00568 else 00569 { 00570 // else, have a dual edge representing a bdy edge - tie directly to 00571 // dual entity if its dual entity 00572 EntityHandle bdy_dcell; 00573 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &dedge_quad, 1, &bdy_dcell ); 00574 TRC assert( MBPOLYGON == mbImpl->type_from_handle( bdy_dcell ) ); 00575 00576 tmp_result = mbImpl->add_adjacencies( this_edge, &bdy_dcell, 1, false ); 00577 if( MB_SUCCESS != tmp_result ) 00578 { 00579 result = tmp_result; 00580 continue; 00581 } 00582 } 00583 } 00584 00585 // sanity check - look for adj edges again, and check for equiv entities 00586 for( Range::iterator vit = save_all_2cells.begin(); vit != save_all_2cells.end(); ++vit ) 00587 { 00588 Range adj_edges, dum_quad_range; 00589 dum_quad_range.insert( *vit ); 00590 assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) ); 00591 tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_edges ); 00592 if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result ) 00593 { 00594 std::cout << "Multiple entities returned for polygon " << mbImpl->id_from_handle( *vit ) << "." 00595 << std::endl; 00596 continue; 00597 } 00598 } 00599 // success! 00600 return result; 00601 } 00602 00603 ErrorCode DualTool::construct_dual_cells( const Range& all_verts, Range& dual_ents ) 00604 { 00605 if( all_verts.empty() ) return MB_SUCCESS; 00606 00607 // make sure they're all edges 00608 assert( 0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.begin() ) ) && 00609 0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.rbegin() ) ) ); 00610 00611 Range::const_iterator rit; 00612 EntityHandle dual_ent; 00613 unsigned int is_dual = 0x1; 00614 ErrorCode tmp_result = MB_SUCCESS; 00615 ErrorCode result = MB_SUCCESS; 00616 std::vector< EntityHandle > edges, dfaces; 00617 00618 for( rit = all_verts.begin(); rit != all_verts.end(); ++rit ) 00619 { 00620 if( tmp_result != MB_SUCCESS ) result = tmp_result; 00621 00622 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent ); 00623 if( MB_SUCCESS == tmp_result && 0 != dual_ent ) 00624 { 00625 dual_ents.insert( dual_ent ); 00626 continue; 00627 } 00628 00629 // no dual entity; construct one; get the edges bounding the vertex 00630 edges.clear(); 00631 dfaces.clear(); 00632 tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 1, false, edges ); 00633 if( MB_SUCCESS != tmp_result ) continue; 00634 00635 // get the dual faces corresponding to the edges 00636 dfaces.resize( edges.size() ); 00637 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &edges[0], edges.size(), &dfaces[0] ); 00638 if( MB_SUCCESS != tmp_result ) continue; 00639 00640 // create the dual cell from those faces 00641 tmp_result = mbImpl->create_element( MBPOLYHEDRON, &dfaces[0], dfaces.size(), dual_ent ); 00642 if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue; 00643 00644 // save it in the list of new dual ents 00645 dual_ents.insert( dual_ent ); 00646 00647 // tag it indicating it's a dual entity 00648 tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual ); 00649 if( MB_SUCCESS != tmp_result ) continue; 00650 00651 // tag the primal entity with its dual entity and vica versa 00652 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent ); 00653 if( MB_SUCCESS != tmp_result ) continue; 00654 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) ); 00655 if( MB_SUCCESS != tmp_result ) continue; 00656 } 00657 00658 return result; 00659 } 00660 00661 //! given an edge handle, return a list of dual vertices in radial order 00662 //! around the edge 00663 ErrorCode DualTool::get_radial_dverts( const EntityHandle edge, 00664 std::vector< EntityHandle >& rad_dverts, 00665 bool& bdy_edge ) 00666 { 00667 rad_dverts.clear(); 00668 00669 std::vector< EntityHandle > rad_faces, rad_ents; 00670 ErrorCode result = MeshTopoUtil( mbImpl ).star_entities( edge, rad_faces, bdy_edge, 0, &rad_ents ); 00671 if( MB_SUCCESS != result ) return result; 00672 00673 if( bdy_edge ) 00674 { 00675 // if we're a bdy edge, change the order back to what DualTool expects 00676 rad_ents.push_back( *rad_faces.rbegin() ); 00677 rad_ents.push_back( *rad_faces.begin() ); 00678 } 00679 00680 rad_dverts.resize( rad_ents.size() ); 00681 for( unsigned int i = 0; i < rad_ents.size(); i++ ) 00682 { 00683 EntityHandle dual_ent; 00684 result = mbImpl->tag_get_data( dualEntity_tag(), &rad_ents[i], 1, &dual_ent ); 00685 if( !bdy_edge || i < rad_ents.size() - 2 ) 00686 rad_dverts[i] = dual_ent; 00687 else 00688 { 00689 // fix up this entry 00690 assert( mbImpl->type_from_handle( dual_ent ) == MBEDGE ); 00691 00692 // get connectivity of that edge 00693 const EntityHandle* connect; 00694 int num_connect; 00695 result = mbImpl->get_connectivity( dual_ent, connect, num_connect ); 00696 if( MB_SUCCESS != result ) return result; 00697 00698 // we want the one that's not already on the list; reuse last_face 00699 int last_hex = ( i == rad_ents.size() - 1 ? 0 : i - 1 ); 00700 EntityHandle last_face = ( connect[0] == rad_dverts[last_hex] ? connect[1] : connect[0] ); 00701 rad_dverts[i] = last_face; 00702 } 00703 } 00704 00705 return result; 00706 } 00707 00708 //! construct the dual entities for a hex mesh, including dual surfaces & curves 00709 ErrorCode DualTool::construct_hex_dual( Range& entities ) 00710 { 00711 std::vector< EntityHandle > evec; 00712 std::copy( entities.begin(), entities.end(), std::back_inserter( evec ) ); 00713 return construct_hex_dual( &evec[0], evec.size() ); 00714 } 00715 00716 //! construct the dual entities for a hex mesh, including dual surfaces & curves 00717 ErrorCode DualTool::construct_hex_dual( EntityHandle* entities, const int num_entities ) 00718 { 00719 // really quite simple: 00720 00721 // construct the dual... 00722 ErrorCode result = construct_dual( entities, num_entities ); 00723 if( MB_SUCCESS != result ) 00724 { 00725 std::cerr << "Error constructing dual entities for primal entities." << std::endl; 00726 return result; 00727 } 00728 00729 // now traverse to build 1d and 2d hyperplanes 00730 result = construct_dual_hyperplanes( 1, entities, num_entities ); 00731 if( MB_SUCCESS != result ) 00732 { 00733 std::cerr << "Problem traversing 1d hyperplanes." << std::endl; 00734 return result; 00735 } 00736 00737 result = construct_dual_hyperplanes( 2, entities, num_entities ); 00738 if( MB_SUCCESS != result ) 00739 { 00740 std::cerr << "Problem traversing 2d hyperplanes." << std::endl; 00741 return result; 00742 } 00743 00744 result = construct_hp_parent_child(); 00745 if( MB_SUCCESS != result ) 00746 { 00747 std::cerr << "Problem constructing parent/child relations between hyperplanes." << std::endl; 00748 return result; 00749 } 00750 00751 // see? simple, just like I said 00752 return MB_SUCCESS; 00753 } 00754 00755 //! get the cells of the dual 00756 ErrorCode DualTool::get_dual_entities( const int dim, EntityHandle* entities, const int num_entities, Range& dual_ents ) 00757 { 00758 if( 0 == isDualCell_tag() ) return MB_SUCCESS; 00759 if( 0 > dim || 3 < dim ) return MB_INDEX_OUT_OF_RANGE; 00760 00761 unsigned int dum = 0x1; 00762 const void* dum_ptr = &dum; 00763 static EntityType dual_type[] = { MBVERTEX, MBEDGE, MBPOLYGON, MBPOLYHEDRON }; 00764 00765 Range dim_ents; 00766 00767 ErrorCode result; 00768 00769 if( 0 == entities || 0 == num_entities ) 00770 { 00771 // just get all the dual entities of this dimension 00772 result = mbImpl->get_entities_by_type_and_tag( 0, dual_type[dim], &isDualCellTag, &dum_ptr, 1, dual_ents ); 00773 } 00774 else 00775 { 00776 // else look for specific dual entities 00777 result = mbImpl->get_adjacencies( entities, num_entities, 3 - dim, false, dim_ents, Interface::UNION ); 00778 if( MB_SUCCESS != result ) return result; 00779 std::vector< EntityHandle > dual_ents_vec( dim_ents.size() ); 00780 result = mbImpl->tag_get_data( dualEntity_tag(), dim_ents, &dual_ents_vec[0] ); 00781 if( MB_SUCCESS != result ) return result; 00782 std::copy( dual_ents_vec.begin(), dual_ents_vec.end(), range_inserter( dual_ents ) ); 00783 } 00784 00785 return result; 00786 } 00787 00788 //! get the faces of the dual 00789 ErrorCode DualTool::get_dual_entities( const int dim, 00790 EntityHandle* entities, 00791 const int num_entities, 00792 std::vector< EntityHandle >& dual_ents ) 00793 { 00794 Range tmp_range; 00795 ErrorCode result = get_dual_entities( dim, entities, num_entities, tmp_range ); 00796 if( MB_SUCCESS != result ) return result; 00797 00798 // dual_ents.insert(dual_ents.end(), tmp_range.begin(), tmp_range.end()); 00799 dual_ents.reserve( dual_ents.size() + tmp_range.size() ); 00800 for( Range::const_iterator it = tmp_range.begin(); it != tmp_range.end(); ++it ) 00801 { 00802 dual_ents.push_back( *it ); 00803 } 00804 return MB_SUCCESS; 00805 } 00806 00807 ErrorCode DualTool::get_dual_hyperplanes( const Interface* impl, const int dim, Range& dual_ents ) 00808 { 00809 if( dim != 1 && dim != 2 ) return MB_INDEX_OUT_OF_RANGE; 00810 00811 Tag dual_tag; 00812 ErrorCode result; 00813 00814 if( dim == 1 ) 00815 result = impl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag ); 00816 else 00817 result = impl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag ); 00818 00819 if( MB_SUCCESS == result ) 00820 result = impl->get_entities_by_type_and_tag( 0, MBENTITYSET, &dual_tag, NULL, 1, dual_ents, Interface::UNION ); 00821 00822 return result; 00823 } 00824 00825 ErrorCode DualTool::construct_dual_hyperplanes( const int dim, EntityHandle* entities, const int num_entities ) 00826 { 00827 // this function traverses dual faces of input dimension, constructing 00828 // dual hyperplanes of them in sets as it goes 00829 00830 // check various inputs 00831 int num_quads, num_hexes; 00832 if( 00833 // this function only makes sense for dim == 1 or dim == 2 00834 ( dim != 1 && dim != 2 ) || 00835 // should either be quads or hexes around 00836 mbImpl->get_number_entities_by_type( 0, MBQUAD, num_quads ) != MB_SUCCESS || 00837 mbImpl->get_number_entities_by_type( 0, MBHEX, num_hexes ) != MB_SUCCESS || 00838 // if we're asking for 1d dual ents, should be quads around 00839 ( num_quads == 0 && dim == 1 ) || 00840 // if we're asking for 2d dual ents, should be hexes around 00841 ( num_hexes == 0 && dim == 2 ) ) 00842 return MB_FAILURE; 00843 00844 Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() ); 00845 00846 // two stacks: one completely untreated entities, and the other untreated 00847 // entities on the current dual hyperplane 00848 std::vector< EntityHandle > tot_untreated; 00849 00850 // put dual entities of this dimension on the untreated list 00851 ErrorCode result = get_dual_entities( dim, entities, num_entities, tot_untreated ); 00852 if( MB_SUCCESS != result ) return result; 00853 00854 // main part of traversal loop 00855 EntityHandle this_ent; 00856 EntityHandle this_hp; 00857 00858 while( !tot_untreated.empty() ) 00859 { 00860 if( debug && dim == 2 /*(tot_untreated.size()%report == 0)*/ ) 00861 std::cout << "Untreated list size " << tot_untreated.size() << "." << std::endl; 00862 00863 this_ent = tot_untreated.back(); 00864 tot_untreated.pop_back(); 00865 result = mbImpl->tag_get_data( hp_tag, &this_ent, 1, &this_hp ); 00866 if( MB_SUCCESS != result && MB_TAG_NOT_FOUND != result ) return result; 00867 00868 // d for this entity having a hyperplane assignment already 00869 else if( this_hp != 0 ) 00870 continue; 00871 00872 if( 1 == dim && check_1d_loop_edge( this_ent ) ) continue; 00873 00874 // inner loop: traverse the hyperplane 'till we don't have any more 00875 result = traverse_hyperplane( hp_tag, this_hp, this_ent ); 00876 if( MB_SUCCESS != result ) 00877 { 00878 std::cout << "Failed to traverse hyperplane "; 00879 if( this_hp ) 00880 std::cout << mbImpl->id_from_handle( this_hp ) << "." << std::endl; 00881 else 00882 std::cout << "0." << std::endl; 00883 return result; 00884 } 00885 00886 // ok, now order the edges if it's a chord 00887 if( 1 == dim ) order_chord( this_hp ); 00888 } 00889 00890 return MB_SUCCESS; 00891 } 00892 00893 ErrorCode DualTool::traverse_hyperplane( const Tag hp_tag, EntityHandle& this_hp, EntityHandle this_ent ) 00894 { 00895 Range tmp_star, star, tmp_range, new_hyperplane_ents; 00896 std::vector< EntityHandle > hp_untreated; 00897 int dim = mbImpl->dimension_from_handle( this_ent ); 00898 MeshTopoUtil mtu( mbImpl ); 00899 this_hp = 0; 00900 ErrorCode result; 00901 00902 unsigned short mark_val = 0x0; 00903 Tag mark_tag; 00904 result = mbImpl->tag_get_handle( "__hyperplane_mark", 1, MB_TYPE_BIT, mark_tag, MB_TAG_CREAT | MB_TAG_BIT ); 00905 if( MB_SUCCESS != result ) return result; 00906 mark_val = 0x1; 00907 00908 while( 0 != this_ent ) 00909 { 00910 EntityHandle tmp_hp = get_dual_hyperplane( this_ent ); 00911 if( 0 == this_hp && 0 != tmp_hp ) this_hp = tmp_hp; 00912 00913 if( 0 == tmp_hp ) new_hyperplane_ents.insert( this_ent ); 00914 00915 if( debug && hp_untreated.size() % 10 == 0 ) 00916 std::cout << "Dual surface " << this_hp << ", hp_untreated list size = " << hp_untreated.size() << "." 00917 << std::endl; 00918 00919 // get the 2nd order adjacencies through lower dimension 00920 tmp_range.clear(); 00921 tmp_star.clear(); 00922 star.clear(); 00923 result = mtu.get_bridge_adjacencies( this_ent, dim - 1, dim, star );RR; 00924 00925 // get the bridge adjacencies through higher dimension 00926 result = mtu.get_bridge_adjacencies( this_ent, dim + 1, dim, tmp_star );RR; 00927 tmp_range = subtract( star, tmp_star ); 00928 00929 for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit ) 00930 { 00931 if( new_hyperplane_ents.find( *rit ) != new_hyperplane_ents.end() ) continue; 00932 00933 // check for tag first, 'cuz it's probably faster than checking adjacencies 00934 // assign to avoid valgrind warning 00935 unsigned short tmp_mark = 0x0; 00936 result = mbImpl->tag_get_data( mark_tag, &( *rit ), 1, &tmp_mark ); 00937 if( MB_SUCCESS == result && mark_val == tmp_mark ) continue; 00938 00939 // if it's on the loop, it's not eligible 00940 if( 1 == dim && check_1d_loop_edge( *rit ) ) continue; 00941 00942 // have one on this hp; just put it on the hp_untreated list for now, 00943 // will get tagged and put in the hp set later 00944 hp_untreated.push_back( *rit ); 00945 result = mbImpl->tag_set_data( mark_tag, &( *rit ), 1, &mark_val ); 00946 if( MB_SUCCESS != result ) return result; 00947 } 00948 00949 // end of inner loop; get the next this_ent, or set to zero 00950 if( hp_untreated.empty() ) 00951 this_ent = 0; 00952 else 00953 { 00954 this_ent = hp_untreated.back(); 00955 hp_untreated.pop_back(); 00956 } 00957 } 00958 00959 if( debug_ap ) 00960 { 00961 std::string hp_name; 00962 if( 2 == dim ) 00963 hp_name = "sheet"; 00964 else 00965 hp_name = "chord"; 00966 00967 if( 0 == this_hp ) 00968 std::cout << "Constructed new " << hp_name << " with "; 00969 else 00970 { 00971 int this_id; 00972 result = mbImpl->tag_get_data( globalId_tag(), &this_hp, 1, &this_id );RR; 00973 std::cout << "Added to " << hp_name << " " << this_id << " "; 00974 } 00975 if( dim == 2 ) 00976 std::cout << "edges:" << std::endl; 00977 else 00978 std::cout << "quads:" << std::endl; 00979 std::vector< EntityHandle > pents( new_hyperplane_ents.size() ); 00980 result = mbImpl->tag_get_data( dualEntity_tag(), new_hyperplane_ents, &pents[0] );RR; 00981 for( std::vector< EntityHandle >::iterator vit = pents.begin(); vit != pents.end(); ++vit ) 00982 { 00983 if( vit != pents.begin() ) std::cout << ", "; 00984 std::cout << mbImpl->id_from_handle( *vit ); 00985 } 00986 std::cout << std::endl; 00987 } 00988 00989 if( 0 == this_hp ) 00990 { 00991 // ok, doesn't have one; make a new hyperplane 00992 int new_id = -1; 00993 result = construct_new_hyperplane( dim, this_hp, new_id ); 00994 if( MB_SUCCESS != result ) return result; 00995 00996 if( debug_ap ) 00997 { 00998 std::cout << "New "; 00999 if( 2 == dim ) 01000 std::cout << " sheet "; 01001 else 01002 std::cout << " chord "; 01003 std::cout << new_id << " constructed." << std::endl; 01004 } 01005 } 01006 01007 // set the hp_val for entities which didn't have one before 01008 std::vector< EntityHandle > hp_tags( new_hyperplane_ents.size() ); 01009 std::fill( hp_tags.begin(), hp_tags.end(), this_hp ); 01010 result = mbImpl->tag_set_data( hp_tag, new_hyperplane_ents, &hp_tags[0] ); 01011 if( MB_SUCCESS != result ) return result; 01012 result = mbImpl->add_entities( this_hp, new_hyperplane_ents ); 01013 if( MB_SUCCESS != result ) return result; 01014 01015 // unmark the entities by removing the tag 01016 result = mbImpl->tag_delete( mark_tag ); 01017 if( MB_SUCCESS != result ) return result; 01018 01019 return MB_SUCCESS; 01020 } 01021 01022 ErrorCode DualTool::order_chord( EntityHandle chord_set ) 01023 { 01024 // re-order the 1cells in the set so they are in order along the chord 01025 // start by finding the vertex dual to a quad 01026 Range verts, one_cells; 01027 ErrorCode result = mbImpl->get_entities_by_dimension( chord_set, 1, one_cells ); 01028 if( MB_SUCCESS != result || one_cells.empty() ) return MB_FAILURE; 01029 01030 result = mbImpl->get_adjacencies( one_cells, 0, false, verts, Interface::UNION ); 01031 if( MB_SUCCESS != result || verts.empty() ) return MB_FAILURE; 01032 01033 EntityHandle last_vert = 0; 01034 for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit ) 01035 { 01036 if( TYPE_FROM_HANDLE( get_dual_entity( *rit ) ) == MBQUAD ) 01037 { 01038 last_vert = *rit; 01039 break; 01040 } 01041 } 01042 // if there's no vertex owned by a quad, just start with 1st one 01043 if( 0 == last_vert ) last_vert = *verts.begin(); 01044 01045 // now, skip from vertex to vertex, building a list of 1cells 01046 std::vector< EntityHandle > ordered_1cells; 01047 EntityHandle last_1cell = 0; 01048 Range dum1, dum2; 01049 const EntityHandle* connect; 01050 int num_connect; 01051 ErrorCode tmp_result = MB_SUCCESS; 01052 while( ordered_1cells.size() != one_cells.size() ) 01053 { 01054 dum1 = one_cells; 01055 result = mbImpl->get_adjacencies( &last_vert, 1, 1, false, dum1 ); 01056 if( 0 != last_1cell ) dum1.erase( last_1cell ); 01057 // assert(1 == dum1.size()); 01058 if( 0 != last_1cell && 1 != dum1.size() ) 01059 { 01060 std::cerr << "unexpected size traversing chord." << std::endl; 01061 tmp_result = MB_FAILURE; 01062 } 01063 01064 last_1cell = *dum1.begin(); 01065 ordered_1cells.push_back( last_1cell ); 01066 result = mbImpl->get_connectivity( last_1cell, connect, num_connect );RR; 01067 if( last_vert == connect[0] ) 01068 last_vert = connect[1]; 01069 else 01070 last_vert = connect[0]; 01071 } 01072 01073 // now have the 1cells in order, replace them in the set 01074 if( MB_SUCCESS == tmp_result ) 01075 { 01076 result = mbImpl->remove_entities( chord_set, one_cells );RR; 01077 result = mbImpl->add_entities( chord_set, &ordered_1cells[0], ordered_1cells.size() );RR; 01078 } 01079 01080 return MB_SUCCESS; 01081 } 01082 01083 ErrorCode DualTool::construct_new_hyperplane( const int dim, EntityHandle& new_hyperplane, int& id ) 01084 { 01085 ErrorCode result; 01086 if( 1 == dim ) 01087 result = mbImpl->create_meshset( ( MESHSET_ORDERED | MESHSET_TRACK_OWNER ), new_hyperplane ); 01088 else 01089 result = mbImpl->create_meshset( ( MESHSET_SET | MESHSET_TRACK_OWNER ), new_hyperplane ); 01090 if( MB_SUCCESS != result ) return result; 01091 01092 if( -1 == id ) 01093 { 01094 Range all_hyperplanes; 01095 result = get_dual_hyperplanes( mbImpl, dim, all_hyperplanes );RR; 01096 std::vector< int > gids( all_hyperplanes.size() ); 01097 result = mbImpl->tag_get_data( globalIdTag, all_hyperplanes, ( gids.empty() ) ? NULL : &gids[0] );RR; 01098 for( unsigned int i = 0; i < gids.size(); i++ ) 01099 if( gids[i] > id ) id = gids[i]; 01100 id++; 01101 if( 0 == id ) id++; 01102 } 01103 01104 result = mbImpl->tag_set_data( globalId_tag(), &new_hyperplane, 1, &id );RR; 01105 Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() ); 01106 result = mbImpl->tag_set_data( hp_tag, &new_hyperplane, 1, &new_hyperplane );RR; 01107 01108 // assign a category name to these sets 01109 static const char dual_category_names[2][CATEGORY_TAG_SIZE] = { "Chord\0", "Sheet\0" }; 01110 01111 result = mbImpl->tag_set_data( categoryTag, &new_hyperplane, 1, dual_category_names[dim - 1] ); 01112 01113 return result; 01114 } 01115 01116 bool DualTool::check_1d_loop_edge( EntityHandle this_ent ) 01117 { 01118 // make sure it's an edge 01119 if( MBEDGE != mbImpl->type_from_handle( this_ent ) ) return false; 01120 01121 // also has to be a dual entity 01122 unsigned int dum; 01123 ErrorCode result = mbImpl->tag_get_data( isDualCell_tag(), &this_ent, 1, &dum ); 01124 if( MB_SUCCESS != result || dum != 0x1 ) return false; 01125 01126 const EntityHandle* verts; 01127 EntityHandle vert_tags[2]; 01128 int num_verts; 01129 result = mbImpl->get_connectivity( this_ent, verts, num_verts ); 01130 if( MB_SUCCESS != result ) return false; 01131 01132 result = mbImpl->tag_get_data( dualEntity_tag(), verts, 2, vert_tags ); 01133 if( MB_SUCCESS != result || mbImpl->type_from_handle( vert_tags[0] ) != MBQUAD || 01134 mbImpl->type_from_handle( vert_tags[1] ) != MBQUAD ) 01135 return false; 01136 01137 else 01138 return true; 01139 } 01140 01141 ErrorCode DualTool::construct_hp_parent_child() 01142 { 01143 Range dual_surfs, dual_cells, dual_edges; 01144 ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs ); 01145 if( MB_SUCCESS != result || dual_surfs.empty() ) return result; 01146 std::vector< EntityHandle > dual_curve_sets; 01147 01148 for( Range::iterator surf_it = dual_surfs.begin(); surf_it != dual_surfs.end(); ++surf_it ) 01149 { 01150 // get all the cells, edges in those cells, and chords for those edges 01151 dual_cells.clear(); 01152 result = mbImpl->get_entities_by_handle( *surf_it, dual_cells ); 01153 if( MB_SUCCESS != result ) return result; 01154 dual_edges.clear(); 01155 result = mbImpl->get_adjacencies( dual_cells, 1, false, dual_edges, Interface::UNION ); 01156 if( MB_SUCCESS != result ) return result; 01157 dual_curve_sets.resize( dual_edges.size() ); 01158 result = mbImpl->tag_get_data( dualCurve_tag(), dual_edges, &dual_curve_sets[0] ); 01159 if( MB_SUCCESS != result ) return result; 01160 01161 // reuse dual_cells to get unique list of chord sets 01162 dual_cells.clear(); 01163 for( unsigned int i = 0; i < dual_edges.size(); i++ ) 01164 if( dual_curve_sets[i] != 0 ) dual_cells.insert( dual_curve_sets[i] ); 01165 01166 // now connect up this dual surf with all the 1d ones 01167 for( Range::iterator rit = dual_cells.begin(); rit != dual_cells.end(); ++rit ) 01168 { 01169 result = mbImpl->add_parent_child( *surf_it, *rit ); 01170 if( MB_SUCCESS != result ) return result; 01171 } 01172 } 01173 01174 return MB_SUCCESS; 01175 } 01176 01177 ErrorCode DualTool::get_graphics_points( EntityHandle dual_ent, 01178 std::vector< int >& npts, 01179 std::vector< GraphicsPoint >& points ) 01180 { 01181 // shouldn't be a set 01182 assert( MBENTITYSET != mbImpl->type_from_handle( dual_ent ) ); 01183 01184 // get the graphics points comprising the given entity 01185 GraphicsPoint gp_array[DualTool::GP_SIZE]; 01186 01187 ErrorCode result = MB_SUCCESS; 01188 01189 // switch based on topological dimension 01190 switch( mbImpl->dimension_from_handle( dual_ent ) ) 01191 { 01192 case 0: 01193 // just return the vertex point 01194 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array ); 01195 if( MB_SUCCESS == result ) points.push_back( gp_array[0] ); 01196 01197 break; 01198 01199 case 1: 01200 // get my graphics point then those of my vertices 01201 const EntityHandle* connect; 01202 int num_connect; 01203 result = mbImpl->get_connectivity( dual_ent, connect, num_connect ); 01204 if( MB_SUCCESS != result ) break; 01205 01206 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, gp_array ); 01207 if( MB_SUCCESS == result ) 01208 { 01209 points.push_back( gp_array[0] ); 01210 points.push_back( gp_array[0] ); 01211 points.push_back( gp_array[1] ); 01212 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array ); 01213 if( MB_SUCCESS == result ) points[1] = gp_array[0]; 01214 } 01215 01216 npts.push_back( 3 ); 01217 01218 break; 01219 01220 case 2: 01221 result = get_cell_points( dual_ent, npts, points ); 01222 break; 01223 } 01224 01225 return result; 01226 } 01227 01228 ErrorCode DualTool::get_cell_points( EntityHandle dual_ent, 01229 std::vector< int >& npts, 01230 std::vector< GraphicsPoint >& points ) 01231 { 01232 assert( MBPOLYGON == mbImpl->type_from_handle( dual_ent ) ); 01233 01234 // get the 1cells in this 2cell 01235 Range one_cells; 01236 01237 Range tc_range; 01238 tc_range.insert( dual_ent ); 01239 ErrorCode result = mbImpl->get_adjacencies( tc_range, 1, false, one_cells, Interface::UNION );RR; 01240 01241 int num_edges = one_cells.size(); 01242 std::vector< GraphicsPoint > dum_gps( num_edges + 1 ); 01243 01244 // get graphics points for 0cells and for this cell 01245 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), one_cells, &dum_gps[0] );RR; 01246 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, &( dum_gps[num_edges] ) );RR; 01247 01248 Range::iterator eit; 01249 const EntityHandle* connect; 01250 int num_connect; 01251 GraphicsPoint vert_gps[2]; 01252 int i; 01253 for( i = 0, eit = one_cells.begin(); i < num_edges; i++, ++eit ) 01254 { 01255 // get the vertices and the graphics points for them 01256 result = mbImpl->get_connectivity( *eit, connect, num_connect );RR; 01257 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, vert_gps );RR; 01258 01259 // make the 2 tris corresponding to this edge; don't worry about order 01260 // for now 01261 npts.push_back( 3 ); 01262 points.push_back( dum_gps[num_edges] ); 01263 points.push_back( vert_gps[0] ); 01264 points.push_back( dum_gps[i] ); 01265 01266 npts.push_back( 3 ); 01267 points.push_back( dum_gps[num_edges] ); 01268 points.push_back( dum_gps[i] ); 01269 points.push_back( vert_gps[1] ); 01270 } 01271 01272 return result; 01273 } 01274 01275 ErrorCode DualTool::get_graphics_points( const Range& in_range, 01276 std::vector< GraphicsPoint >& points, 01277 const bool assign_ids, 01278 const int start_id ) 01279 { 01280 // return graphics points on dual entities in in_range or in entities 01281 // in sets in in_range 01282 ErrorCode result; 01283 01284 // for each dual hyperplane set: 01285 Range::const_iterator rit; 01286 01287 Range two_cells, all_cells; 01288 for( rit = in_range.begin(); rit != in_range.end(); ++rit ) 01289 { 01290 // for each entity: 01291 two_cells.clear(); 01292 EntityType this_type = mbImpl->type_from_handle( *rit ); 01293 if( MBENTITYSET == this_type ) 01294 { 01295 result = mbImpl->get_entities_by_handle( *rit, two_cells );RR; 01296 01297 std::copy( two_cells.begin(), two_cells.end(), range_inserter( all_cells ) ); 01298 } 01299 01300 else 01301 { 01302 two_cells.insert( *rit ); 01303 assert( this_type == MBVERTEX || this_type == MBEDGE || this_type == MBPOLYGON || 01304 this_type == MBPOLYHEDRON ); 01305 } 01306 01307 result = mbImpl->get_adjacencies( two_cells, 0, false, all_cells, Interface::UNION );RR; 01308 result = mbImpl->get_adjacencies( two_cells, 1, false, all_cells, Interface::UNION );RR; 01309 } 01310 01311 // get graphics points 01312 points.resize( all_cells.size() ); 01313 01314 result = mbImpl->tag_get_data( dualGraphicsPointTag, all_cells, &points[0] );RR; 01315 01316 if( assign_ids ) 01317 { 01318 int i = start_id; 01319 01320 for( std::vector< GraphicsPoint >::iterator vit = points.begin(); vit != points.end(); ++vit ) 01321 vit->id = i++; 01322 01323 result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), all_cells, &points[0] );RR; 01324 } 01325 01326 return result; 01327 } 01328 01329 EntityHandle DualTool::next_loop_vertex( const EntityHandle last_v, 01330 const EntityHandle this_v, 01331 const EntityHandle dual_surf ) 01332 { 01333 // given two vertices, find the next one on the loop; if one is a dual 01334 // surface, then just choose either one for that surface 01335 assert( ( 0 == last_v || mbImpl->type_from_handle( last_v ) == MBVERTEX ) && 01336 mbImpl->type_from_handle( this_v ) == MBVERTEX && mbImpl->type_from_handle( dual_surf ) == MBENTITYSET ); 01337 01338 // get the connected vertices 01339 MeshTopoUtil tpu( mbImpl ); 01340 Range other_verts; 01341 ErrorCode result = tpu.get_bridge_adjacencies( this_v, 1, 0, other_verts ); 01342 if( MB_SUCCESS != result || other_verts.empty() ) return 0; 01343 01344 // if (mbImpl->type_from_handle(last_v) == MBENTITYSET) { 01345 // dual surface, choose either; first get a 2cell on this surface 01346 Range tcells, tcells2, verts; 01347 result = mbImpl->get_entities_by_type( dual_surf, MBPOLYGON, tcells ); 01348 if( MB_SUCCESS != result || tcells.empty() ) return 0; 01349 01350 // ok, pay attention here: first get 2cells common to dual surface and this_v 01351 verts.insert( this_v ); 01352 result = mbImpl->get_adjacencies( verts, 2, false, tcells ); 01353 if( MB_SUCCESS != result || tcells.empty() ) return 0; 01354 01355 // next get vertices common to both 2cells and subtract from other_verts; also 01356 // remove last_v if it's non-zero 01357 verts.clear(); 01358 result = mbImpl->get_adjacencies( tcells, 0, false, verts ); 01359 if( MB_SUCCESS != result || verts.empty() ) return 0; 01360 01361 Range tmp_verts = subtract( other_verts, verts ); 01362 other_verts.swap( tmp_verts ); 01363 if( 0 != last_v ) other_verts.erase( last_v ); 01364 01365 // now get intersection of remaining vertices and 2 2cells vertices 01366 // look at each one successively; should find one, maybe not on both 01367 tmp_verts = other_verts; 01368 Range tmp_faces( *tcells.begin(), *tcells.begin() ); 01369 result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts ); 01370 if( MB_SUCCESS == result && !tmp_verts.empty() ) return *tmp_verts.begin(); 01371 tmp_faces.clear(); 01372 tmp_faces.insert( *tcells.rbegin() ); 01373 result = mbImpl->get_adjacencies( tmp_faces, 0, false, other_verts ); 01374 if( MB_SUCCESS == result && !other_verts.empty() ) return *other_verts.begin(); 01375 01376 // if we got here, there isn't any 01377 return 0; 01378 } 01379 01380 EntityHandle DualTool::get_dual_hyperplane( const EntityHandle ncell ) 01381 { 01382 // get the sheet or chord it's in 01383 std::vector< EntityHandle > adj_sets; 01384 ErrorCode result = mbImpl->get_adjacencies( &ncell, 1, 4, false, adj_sets ); 01385 if( MB_SUCCESS != result ) return 0; 01386 01387 EntityHandle dum_set; 01388 for( std::vector< EntityHandle >::iterator vit = adj_sets.begin(); vit != adj_sets.end(); ++vit ) 01389 { 01390 if( mbImpl->tag_get_data( dualCurve_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND || 01391 mbImpl->tag_get_data( dualSurface_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND ) 01392 return *vit; 01393 } 01394 01395 return 0; 01396 } 01397 01398 //! set the dual surface or curve for an entity 01399 ErrorCode DualTool::set_dual_surface_or_curve( EntityHandle entity, 01400 const EntityHandle dual_hyperplane, 01401 const int dual_entity_dimension ) 01402 { 01403 if( 1 == dual_entity_dimension ) 01404 mbImpl->tag_set_data( dualCurve_tag(), &entity, 1, &dual_hyperplane ); 01405 else if( 2 == dual_entity_dimension ) 01406 mbImpl->tag_set_data( dualSurface_tag(), &entity, 1, &dual_hyperplane ); 01407 else 01408 return MB_INDEX_OUT_OF_RANGE; 01409 01410 return MB_SUCCESS; 01411 } 01412 01413 //! return the corresponding dual entity 01414 EntityHandle DualTool::get_dual_entity( const EntityHandle this_ent ) const 01415 { 01416 EntityHandle dual_ent; 01417 ErrorCode result = mbImpl->tag_get_data( dualEntity_tag(), &this_ent, 1, &dual_ent ); 01418 if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result ) 01419 return 0; 01420 else 01421 return dual_ent; 01422 } 01423 01424 //! return the corresponding dual entity 01425 EntityHandle DualTool::get_extra_dual_entity( const EntityHandle this_ent ) 01426 { 01427 EntityHandle dual_ent; 01428 ErrorCode result = mbImpl->tag_get_data( extraDualEntity_tag(), &this_ent, 1, &dual_ent ); 01429 if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result ) 01430 return 0; 01431 else 01432 return dual_ent; 01433 } 01434 01435 ErrorCode DualTool::atomic_pillow( EntityHandle odedge, EntityHandle& quad1, EntityHandle& quad2 ) 01436 { 01437 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 01438 01439 if( debug_ap ) 01440 { 01441 Range sets; 01442 Tag ms_tag; 01443 01444 ErrorCode result = mbImpl->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, ms_tag ); 01445 if( MB_SUCCESS == result ) 01446 { 01447 result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &ms_tag, NULL, 1, sets ); 01448 if( MB_SUCCESS == result ) result = mbImpl->delete_entities( sets ); 01449 } 01450 } 01451 01452 std::cout << "-AP("; 01453 print_cell( odedge ); 01454 std::cout << ")" << std::endl; 01455 01456 // perform an atomic pillow operation around dedge 01457 01458 // grab the quad before deleting the odedge 01459 quad1 = get_dual_entity( odedge ); 01460 assert( 0 != quad1 ); 01461 01462 // 0. get star 2cells around odedge (before odedge changes) and 3cells around 01463 // those 2cells (going to delete the 2cells, therefore need to delete the 3cells 01464 // that depend on those too) 01465 MeshTopoUtil mtu( mbImpl ); 01466 Range star_cells, tmp_cells; 01467 ErrorCode result = mbImpl->get_adjacencies( &odedge, 1, 2, false, star_cells );RR; 01468 result = mbImpl->get_adjacencies( star_cells, 3, false, tmp_cells, Interface::UNION );RR; 01469 star_cells.merge( tmp_cells ); 01470 star_cells.insert( odedge ); 01471 01472 // tear down the dual entities which will be modified by the ap first 01473 result = delete_dual_entities( star_cells );RR; 01474 01475 // now change the quad to an ap 01476 std::vector< EntityHandle > verts; 01477 result = mbImpl->get_connectivity( &quad1, 1, verts );RR; 01478 01479 // get average position of vertices 01480 double coords[12], avg[3] = { 0.0, 0.0, 0.0 }; 01481 result = mbImpl->get_coords( &verts[0], verts.size(), coords );RR; 01482 for( int i = 0; i < 4; i++ ) 01483 { 01484 avg[0] += coords[3 * i]; 01485 avg[1] += coords[3 * i + 1]; 01486 avg[2] += coords[3 * i + 2]; 01487 } 01488 for( int i = 0; i < 3; i++ ) 01489 avg[i] *= 0.25; 01490 01491 // for each position, get a corresponding position 1/2 way to avg 01492 double new_coords[12]; 01493 for( int i = 0; i < 4; i++ ) 01494 { 01495 new_coords[3 * i] = avg[0] + .5 * ( coords[3 * i] - avg[0] ); 01496 new_coords[3 * i + 1] = avg[1] + .5 * ( coords[3 * i + 1] - avg[1] ); 01497 new_coords[3 * i + 2] = avg[2] + .5 * ( coords[3 * i + 2] - avg[2] ); 01498 } 01499 01500 // make the 4 new vertices; store in vector long enough for hex connectivity 01501 for( int i = 0; i < 4; i++ ) 01502 { 01503 verts.push_back( 0 ); 01504 result = mbImpl->create_vertex( &new_coords[3 * i], verts[4 + i] );RR; 01505 } 01506 01507 // get the hexes connected to the quad 01508 Range hexes; 01509 result = mbImpl->get_adjacencies( &quad1, 1, 3, false, hexes );RR; 01510 assert( hexes.size() <= 2 ); 01511 01512 // remove any explicit adjacency from the first hex, since that'll get connected 01513 // to the new quad; add adjacency between quad and second hex, if there is a 2nd 01514 result = mbImpl->remove_adjacencies( quad1, &( *hexes.begin() ), 1 );RR; 01515 if( hexes.size() == 2 ) 01516 { 01517 result = mbImpl->add_adjacencies( quad1, &( *hexes.rbegin() ), 1, false );RR; 01518 } 01519 01520 // create the new, inner quad, and make it explicitly adjacent to 1st hex; 01521 // make the connectivity of this quad same as the original one 01522 std::vector< EntityHandle > tmp_verts; 01523 std::copy( verts.begin(), verts.end(), std::back_inserter( tmp_verts ) ); 01524 01525 result = mbImpl->create_element( MBQUAD, &tmp_verts[0], 4, quad2 );RR; 01526 result = mbImpl->add_adjacencies( quad2, &( *hexes.begin() ), 1, false );RR; 01527 01528 // reverse the connectivity of the 1st hex 01529 std::reverse( verts.begin(), verts.begin() + 4 ); 01530 std::reverse( verts.begin() + 4, verts.end() ); 01531 01532 // now make two inner hexes; note connectivity array is flipped for the two hexes 01533 EntityHandle new_hexes[2]; 01534 result = mbImpl->create_element( MBHEX, &verts[0], 8, new_hexes[0] );RR; 01535 result = mbImpl->create_element( MBHEX, &tmp_verts[0], 8, new_hexes[1] );RR; 01536 01537 // set the global id tag on the new hexes 01538 int new_hex_ids[2] = { maxHexId + 1, maxHexId + 2 }; 01539 maxHexId += 2; 01540 result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 2, new_hex_ids ); 01541 if( MB_SUCCESS != result ) return result; 01542 01543 // by definition, quad1 is adj to new_hexes[0] 01544 result = mbImpl->add_adjacencies( quad1, &new_hexes[0], 1, false );RR; 01545 result = mbImpl->add_adjacencies( quad2, &new_hexes[1], 1, false );RR; 01546 01547 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 01548 01549 // now update the dual 01550 result = construct_hex_dual( &new_hexes[0], 2 );RR; 01551 01552 // get the new dual surface, by getting one of the edges between the center 01553 // and outer vertex rings 01554 Range new_edge; 01555 verts[1] = verts[4]; 01556 result = mbImpl->get_adjacencies( &verts[0], 2, 1, false, new_edge ); 01557 if( MB_SUCCESS != result || new_edge.size() != 1 ) return result; 01558 01559 return MB_SUCCESS; 01560 } 01561 01562 //! effect reverse atomic pillow operation 01563 ErrorCode DualTool::rev_atomic_pillow( EntityHandle pillow, Range& chords ) 01564 { 01565 // get the dual entities associated with elements in the pillow; go through 01566 // the elements instead of the pillow sheet so you get all of them, not just 01567 // the ones on the sheet 01568 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 01569 01570 std::cout << "-AP("; 01571 print_cell( pillow ); 01572 std::cout << ")" << std::endl; 01573 01574 Range dverts; 01575 ErrorCode result = get_dual_entities( pillow, NULL, NULL, &dverts, NULL, NULL ); 01576 if( MB_SUCCESS != result ) return result; 01577 assert( 2 == dverts.size() ); 01578 01579 EntityHandle hexes[2]; 01580 result = mbImpl->tag_get_data( dualEntity_tag(), dverts, hexes );RR; 01581 assert( hexes[0] != 0 && hexes[1] != 0 ); 01582 01583 std::vector< EntityHandle > dcells[4]; 01584 Range pcells[4]; 01585 std::copy( hexes, hexes + 2, range_inserter( pcells[3] ) ); 01586 std::copy( dverts.begin(), dverts.end(), std::back_inserter( dcells[0] ) ); 01587 for( int dim = 0; dim <= 2; dim++ ) 01588 { 01589 result = mbImpl->get_adjacencies( hexes, 2, dim, false, pcells[dim], Interface::UNION );RR; 01590 dcells[3 - dim].resize( pcells[dim].size() ); 01591 result = mbImpl->tag_get_data( dualEntity_tag(), pcells[dim], &dcells[3 - dim][0] );RR; 01592 } 01593 01594 // delete the dual entities which are part of the original pillow 01595 result = mbImpl->delete_entities( &pillow, 1 ); 01596 if( MB_SUCCESS != result ) return result; 01597 01598 result = mbImpl->delete_entities( chords ); 01599 if( MB_SUCCESS != result ) return result; 01600 01601 for( int i = 3; i >= 0; i-- ) 01602 { 01603 result = delete_dual_entities( &dcells[i][0], dcells[i].size() );RR; 01604 } 01605 01606 // delete the primal entities inserted by the ap; be careful to get the right 01607 // faces, edges and vertices 01608 Range del_faces, del_edges, del_verts, tmp_faces, tmp_verts; 01609 // faces are the shared 5 and the 1 other one with greater handle (which 01610 // we know will be later in the range) 01611 result = mbImpl->get_adjacencies( hexes, 2, 2, false, del_faces );RR; 01612 assert( 5 == del_faces.size() ); 01613 std::copy( pcells[2].begin(), pcells[2].end(), range_inserter( tmp_faces ) ); 01614 tmp_faces = subtract( tmp_faces, del_faces ); 01615 del_faces.insert( *tmp_faces.rbegin() ); 01616 result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts );RR; 01617 std::copy( pcells[0].begin(), pcells[0].end(), range_inserter( del_verts ) ); 01618 del_verts = subtract( del_verts, tmp_verts ); 01619 assert( 4 == del_verts.size() ); 01620 result = mbImpl->get_adjacencies( del_verts, 1, false, del_edges, Interface::UNION );RR; 01621 assert( 8 == del_edges.size() ); 01622 01623 result = mbImpl->delete_entities( hexes, 2 );RR; 01624 result = mbImpl->delete_entities( del_faces );RR; 01625 result = mbImpl->delete_entities( del_edges );RR; 01626 result = mbImpl->delete_entities( del_verts );RR; 01627 01628 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 01629 01630 // recompute the dual for the hexes on either side of the quad affected 01631 // by the ap removal 01632 Range tmp_hexes; 01633 result = mbImpl->get_adjacencies( tmp_verts, 3, false, tmp_hexes, Interface::UNION );RR; 01634 result = construct_hex_dual( tmp_hexes );RR; 01635 01636 return MB_SUCCESS; 01637 } 01638 01639 ErrorCode DualTool::delete_dual_entities( EntityHandle* entities, int num_entities ) 01640 { 01641 Range tmp_ents; 01642 std::copy( entities, entities + num_entities, range_inserter( tmp_ents ) ); 01643 return delete_dual_entities( tmp_ents ); 01644 } 01645 01646 ErrorCode DualTool::delete_dual_entities( Range& entities ) 01647 { 01648 if( entities.empty() ) return delete_whole_dual(); 01649 01650 EntityHandle null_entity = 0; 01651 ErrorCode result; 01652 Range ents_to_delete; 01653 01654 while( !entities.empty() ) 01655 { 01656 EntityHandle this_entity = entities.pop_back(); 01657 01658 // reset the primal's dual entity 01659 EntityHandle primal = get_dual_entity( this_entity ); 01660 if( get_dual_entity( primal ) == this_entity ) 01661 { 01662 result = mbImpl->tag_set_data( dualEntity_tag(), &primal, 1, &null_entity );RR; 01663 } 01664 EntityHandle extra = get_extra_dual_entity( primal ); 01665 if( 0 != extra ) 01666 { 01667 result = mbImpl->tag_set_data( extraDualEntity_tag(), &primal, 1, &null_entity );RR; 01668 } 01669 01670 ents_to_delete.insert( this_entity ); 01671 01672 // check for extra dual entities 01673 if( mbImpl->type_from_handle( this_entity ) == MBPOLYGON ) 01674 { 01675 // for 2cell, might be a loop edge 01676 Range loop_edges; 01677 result = mbImpl->get_adjacencies( &this_entity, 1, 1, false, loop_edges ); 01678 for( Range::iterator rit = loop_edges.begin(); rit != loop_edges.end(); ++rit ) 01679 if( check_1d_loop_edge( *rit ) ) entities.insert( *rit ); 01680 } 01681 else if( extra && extra != this_entity ) 01682 // just put it on the list; primal for which we're extra has already been 01683 // reset to not point to extra entity 01684 ents_to_delete.insert( extra ); 01685 } 01686 01687 // now delete the entities (sheets and chords will be updated automatically) 01688 return mbImpl->delete_entities( ents_to_delete ); 01689 } 01690 01691 void DualTool::print_cell( EntityHandle cell ) 01692 { 01693 const EntityHandle* connect; 01694 int num_connect; 01695 ErrorCode result = mbImpl->get_connectivity( cell, connect, num_connect ); 01696 if( MB_SUCCESS != result ) return; 01697 bool first = true; 01698 EntityHandle primals[20]; 01699 std::vector< int > ids; 01700 01701 assert( num_connect < 20 ); 01702 result = mbImpl->tag_get_data( dualEntityTag, connect, num_connect, primals ); 01703 if( MB_SUCCESS != result ) return; 01704 ids.resize( num_connect ); 01705 result = mbImpl->tag_get_data( globalIdTag, primals, num_connect, &ids[0] ); 01706 if( MB_SUCCESS != result ) return; 01707 for( int i = 0; i < num_connect; i++ ) 01708 { 01709 if( !first ) std::cout << "-"; 01710 EntityType this_type = mbImpl->type_from_handle( primals[i] ); 01711 if( this_type == MBHEX ) 01712 std::cout << "h"; 01713 else if( this_type == MBQUAD ) 01714 std::cout << "f"; 01715 else 01716 std::cout << "u"; 01717 01718 if( ids[i] != 0 ) 01719 std::cout << ids[i]; 01720 else 01721 std::cout << mbImpl->id_from_handle( primals[i] ); 01722 01723 first = false; 01724 } 01725 } 01726 01727 ErrorCode DualTool::face_open_collapse( EntityHandle ocl, EntityHandle ocr ) 01728 { 01729 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 01730 01731 MeshTopoUtil mtu( mbImpl ); 01732 01733 std::cout << "OC("; 01734 print_cell( ocl ); 01735 std::cout << ")-("; 01736 print_cell( ocr ); 01737 std::cout << ")" << std::endl; 01738 01739 // get the primal entities we're dealing with 01740 EntityHandle split_quads[2] = { 0 }, split_edges[3] = { 0 }, split_nodes[2] = { 0 }, other_edges[6] = { 0 }, 01741 other_nodes[6] = { 0 }; 01742 Range hexes; 01743 ErrorCode result = foc_get_ents( ocl, ocr, split_quads, split_edges, split_nodes, hexes, other_edges, other_nodes );RR; 01744 01745 // get star entities around edges, separated into halves 01746 std::vector< EntityHandle > star_dp1[2], star_dp2[2]; 01747 result = foc_get_stars( split_quads, split_edges, star_dp1, star_dp2 );RR; 01748 01749 if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) ) 01750 return MB_TYPE_OUT_OF_RANGE; 01751 01752 result = foc_delete_dual( split_quads, split_edges, hexes ); 01753 if( MB_SUCCESS != result ) return result; 01754 01755 EntityHandle new_quads[2], new_edges[3], new_nodes[2]; 01756 result = split_pair_nonmanifold( split_quads, split_edges, split_nodes, star_dp1, star_dp2, other_edges, 01757 other_nodes, new_quads, new_edges, new_nodes ); 01758 if( MB_SUCCESS != result ) return result; 01759 01760 // now merge entities, the C of foc 01761 EntityHandle keepit, deleteit; 01762 #define MIN( a, b ) ( ( a ) < ( b ) ? ( a ) : ( b ) ) 01763 #define MAX( a, b ) ( ( a ) > ( b ) ? ( a ) : ( b ) ) 01764 #define KEEP_DELETE( a, b, c, d ) \ 01765 { \ 01766 ( c ) = MIN( a, b ); \ 01767 ( d ) = MAX( a, b ); \ 01768 } 01769 01770 // find how many shared edges there were 01771 int num_shared_edges = ( split_edges[2] ? 3 : ( split_edges[1] ? 2 : 1 ) ); 01772 01773 // first the node(s) 01774 for( int i = 0; i < 3 - num_shared_edges; i++ ) 01775 { 01776 KEEP_DELETE( other_nodes[2 + 2 * i], other_nodes[3 + 2 * i], keepit, deleteit ); 01777 result = mbImpl->merge_entities( keepit, deleteit, false, true );RR; 01778 } 01779 01780 // now the edges 01781 for( int i = 0; i < 4 - num_shared_edges; i++ ) 01782 { 01783 KEEP_DELETE( other_edges[2 * i], other_edges[2 * i + 1], keepit, deleteit ); 01784 result = mbImpl->merge_entities( keepit, deleteit, false, true );RR; 01785 } 01786 01787 // now the faces 01788 KEEP_DELETE( split_quads[0], split_quads[1], keepit, deleteit ); 01789 result = mbImpl->merge_entities( keepit, deleteit, false, true );RR; 01790 01791 result = mbImpl->merge_entities( new_quads[0], new_quads[1], false, true );RR; 01792 01793 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 01794 01795 // reconstruct dual 01796 result = construct_hex_dual( hexes ); 01797 if( MB_SUCCESS != result ) return result; 01798 01799 return check_dual_adjs(); 01800 } 01801 01802 ErrorCode DualTool::foc_get_ents( EntityHandle ocl, 01803 EntityHandle ocr, 01804 EntityHandle* split_quads, 01805 EntityHandle* split_edges, 01806 EntityHandle* split_nodes, 01807 Range& hexes, 01808 EntityHandle* other_edges, 01809 EntityHandle* other_nodes ) 01810 { 01811 // get the entities used for foc; ocl and ocr are dual 1-cells 01812 // representing quads to be split; returned from this function: 01813 // quads[2] - 2 quads to be split 01814 // split_edges[2] - edge(s) to be split (2nd is 0 if only one) 01815 // split_node - node to be split, if any (otherwise 0) 01816 // hexes - connected hexes to split_edges 01817 // other_edges[0], [1] - edges in quads[0] and [1] sharing node with 01818 // one end of split_edges[0] 01819 // other_edges[2], [3] - other end of split_edges[0] (or [1] if 2 01820 // split_edges) 01821 // other_edges[4], [5] - edges in quads[0], [1] opposite to split_edges[0] 01822 // other_nodes[0], [1] - nodes on other_edges[0], [1] not shared with 01823 // split_edges[0] 01824 // other_nodes[2], [3] - nodes on other_edges[2], [3] not shared with 01825 // split_edges[0] (if 2 split_edges, there's only 1 opposite node 01826 // in each split quad) 01827 // (for diagram, see Tim's notes from 11/12/07) 01828 01829 split_quads[0] = get_dual_entity( ocl ); 01830 split_quads[1] = get_dual_entity( ocr ); 01831 if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) ) 01832 return MB_TYPE_OUT_OF_RANGE; 01833 01834 Range common_edges; 01835 ErrorCode result = mbImpl->get_adjacencies( split_quads, 2, 1, false, common_edges ); 01836 if( MB_SUCCESS != result ) return result; 01837 01838 if( common_edges.empty() ) return MB_FAILURE; 01839 for( unsigned int i = 0; i < common_edges.size(); i++ ) 01840 split_edges[i] = common_edges[i]; 01841 01842 MeshTopoUtil mtu( mbImpl ); 01843 01844 if( common_edges.size() == 3 ) 01845 { 01846 // find other (non-shared) edges 01847 for( int i = 0; i < 2; i++ ) 01848 { 01849 Range tmp_edges; 01850 result = mbImpl->get_adjacencies( &split_quads[i], 1, 1, false, tmp_edges ); 01851 if( MB_SUCCESS != result ) return result; 01852 tmp_edges = subtract( tmp_edges, common_edges ); 01853 assert( tmp_edges.size() == 1 ); 01854 other_edges[i] = *tmp_edges.begin(); 01855 } 01856 assert( other_edges[0] && other_edges[1] && other_edges[0] != other_edges[1] ); 01857 01858 // arrange common edges so middle is in middle 01859 result = mtu.opposite_entity( split_quads[0], other_edges[0], split_edges[1] );RR; 01860 common_edges.erase( split_edges[1] ); 01861 split_edges[0] = *common_edges.begin(); 01862 split_edges[2] = *common_edges.rbegin(); 01863 common_edges.insert( split_edges[1] ); 01864 01865 // get split nodes and other nodes 01866 split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 ); 01867 split_nodes[1] = mtu.common_entity( split_edges[2], split_edges[1], 0 ); 01868 other_nodes[0] = mtu.common_entity( split_edges[0], other_edges[0], 0 ); 01869 other_nodes[1] = mtu.common_entity( split_edges[2], other_edges[1], 0 ); 01870 01871 assert( other_nodes[0] && other_nodes[1] && split_nodes[0] && split_nodes[1] ); 01872 assert( split_edges[0] && split_edges[1] && split_edges[2] && split_edges[0] != split_edges[1] && 01873 split_edges[1] != split_edges[2] && split_edges[0] != split_edges[2] ); 01874 } 01875 else if( common_edges.size() == 2 ) 01876 { 01877 // split node is shared by split edges 01878 split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 ); 01879 if( 0 == split_nodes[0] ) return MB_FAILURE; 01880 // first two other nodes are on split edges opposite split node 01881 result = mtu.opposite_entity( split_edges[0], split_nodes[0], other_nodes[0] );RR; 01882 result = mtu.opposite_entity( split_edges[1], split_nodes[0], other_nodes[1] );RR; 01883 // over split quads: 01884 for( int i = 0; i < 2; i++ ) 01885 { 01886 // 1st other edge is opposite second split edge 01887 result = mtu.opposite_entity( split_quads[i], split_edges[1], other_edges[i] );RR; 01888 // 2nd other edge is opposite first split edge 01889 result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR; 01890 // last other node is opposite split node on split quad 01891 result = mtu.opposite_entity( split_quads[i], split_nodes[0], other_nodes[2 + i] );RR; 01892 } 01893 } 01894 else 01895 { 01896 const EntityHandle* connect; 01897 int num_connect; 01898 result = mbImpl->get_connectivity( split_edges[0], connect, num_connect ); 01899 if( MB_SUCCESS != result ) return result; 01900 // other_nodes[0], [1] are on split edge 01901 other_nodes[0] = connect[0]; 01902 other_nodes[1] = connect[1]; 01903 01904 // for each of the split quads 01905 for( int i = 0; i < 2; i++ ) 01906 { 01907 // get the other edge on the split quad adj to node 0 on the split edge, by getting 01908 // edges adj to split quad and node and removing split edge; that's other_edge[i] 01909 Range tmp_range1, tmp_range2; 01910 tmp_range1.insert( connect[0] ); 01911 tmp_range1.insert( split_quads[i] ); 01912 result = mbImpl->get_adjacencies( tmp_range1, 1, false, tmp_range2 ); 01913 if( MB_SUCCESS != result ) return result; 01914 tmp_range2.erase( split_edges[0] ); 01915 assert( tmp_range2.size() == 1 ); 01916 other_edges[i] = *tmp_range2.begin(); 01917 // get edge connected to other node on split edge & split quad; that's 01918 // opposite prev other_edges on the split quad; that's other_edges[4+i] 01919 result = mtu.opposite_entity( split_quads[i], other_edges[i], other_edges[4 + i] );RR; 01920 // get the edge on the split quad opposite the split edge; that's other_edges[2+i] 01921 result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR; 01922 // get nodes on other side of split quad from split edge, by getting common 01923 // node between top/bottom edge and opposite edge 01924 other_nodes[2 + i] = mtu.common_entity( other_edges[i], other_edges[2 + i], 0 ); 01925 other_nodes[4 + i] = mtu.common_entity( other_edges[4 + i], other_edges[2 + i], 0 ); 01926 if( 0 == other_nodes[2 + i] || 0 == other_nodes[4 + i] ) return MB_FAILURE; 01927 } 01928 } 01929 01930 result = mbImpl->get_adjacencies( split_edges, common_edges.size(), 3, false, hexes, Interface::UNION ); 01931 if( MB_SUCCESS != result ) return result; 01932 01933 assert( "split node not in other_nodes" && other_nodes[0] != split_nodes[0] && other_nodes[0] != split_nodes[1] && 01934 other_nodes[1] != split_nodes[0] && other_nodes[1] != split_nodes[1] ); 01935 assert( "each split node on an end of a split edge" && mtu.common_entity( other_nodes[0], split_edges[0], 0 ) && 01936 ( ( ( split_edges[2] && mtu.common_entity( other_nodes[1], split_edges[2], 0 ) ) || 01937 ( split_edges[1] && mtu.common_entity( other_nodes[1], split_edges[1], 0 ) ) || 01938 mtu.common_entity( other_nodes[1], split_edges[0], 0 ) ) ) ); 01939 assert( "opposite other edges meet at an other node" && 01940 ( mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[0] || 01941 ( split_edges[2] && mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[1] ) ) && 01942 ( split_edges[2] || 01943 ( split_edges[1] && mtu.common_entity( other_edges[2], other_edges[3], 0 ) == other_nodes[1] ) || 01944 mtu.common_entity( other_edges[4], other_edges[5], 0 ) == other_nodes[1] ) ); 01945 01946 return MB_SUCCESS; 01947 } 01948 01949 ErrorCode DualTool::split_pair_nonmanifold( EntityHandle* split_quads, 01950 EntityHandle* split_edges, 01951 EntityHandle* split_nodes, 01952 std::vector< EntityHandle >* star_dp1, 01953 std::vector< EntityHandle >* star_dp2, 01954 EntityHandle* /*other_edges*/, 01955 EntityHandle* /*other_nodes*/, 01956 EntityHandle* new_quads, 01957 EntityHandle* new_edges, 01958 EntityHandle* new_nodes ) 01959 { 01960 01961 // if there's a bdy in the star around the shared edge(s), get the quads on that 01962 // bdy so we know which edges to merge after the split-nonmanifold 01963 MeshTopoUtil mtu( mbImpl ); 01964 ErrorCode result; 01965 01966 // get which star the split faces are in, and choose the other one 01967 int new_side = -1; 01968 if( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() ) 01969 new_side = 1; 01970 else if( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() ) 01971 new_side = 0; 01972 assert( -1 != new_side ); 01973 if( -1 == new_side ) return MB_FAILURE; 01974 01975 //=============== split faces 01976 01977 for( int i = 0; i < 2; i++ ) 01978 { 01979 // get a hex in star_dp2[new_side] that's adj to this split quad, to tell 01980 // mtu which one the new quad should go with; there should be at least one, 01981 // if we have any hexes connected to the split quad 01982 EntityHandle gowith_hex = 0; 01983 for( std::vector< EntityHandle >::iterator vit = star_dp2[new_side].begin(); vit != star_dp2[new_side].end(); 01984 ++vit ) 01985 { 01986 if( mtu.common_entity( *vit, split_quads[i], 2 ) ) 01987 { 01988 gowith_hex = *vit; 01989 break; 01990 } 01991 } 01992 assert( 0 != gowith_hex ); 01993 01994 // split manifold each of the split_quads, and put the results on the merge list 01995 result = 01996 mtu.split_entities_manifold( split_quads + i, 1, new_quads + i, NULL, ( gowith_hex ? &gowith_hex : NULL ) );RR; 01997 } 01998 01999 // make ranges of faces which need to be explicitly adj to old, new 02000 // edge; faces come from stars and new_quads (which weren't in the stars); 02001 // new_quads go with side j, which does not have split quads 02002 Range tmp_addl_faces[2], addl_faces[2]; 02003 for( int i = 0; i < 2; i++ ) 02004 { 02005 std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_faces[i] ) ); 02006 tmp_addl_faces[new_side].insert( new_quads[i] ); 02007 } 02008 #ifndef NDEBUG 02009 bool cond1 = ( "split_quads on 1, new_quads on 0" && 02010 tmp_addl_faces[0].find( split_quads[0] ) == tmp_addl_faces[0].end() && 02011 tmp_addl_faces[0].find( split_quads[1] ) == tmp_addl_faces[0].end() && 02012 tmp_addl_faces[1].find( split_quads[0] ) != tmp_addl_faces[1].end() && 02013 tmp_addl_faces[1].find( split_quads[1] ) != tmp_addl_faces[1].end() && 02014 tmp_addl_faces[0].find( new_quads[0] ) != tmp_addl_faces[0].end() && 02015 tmp_addl_faces[0].find( new_quads[1] ) != tmp_addl_faces[0].end() && 02016 tmp_addl_faces[1].find( new_quads[0] ) == tmp_addl_faces[1].end() && 02017 tmp_addl_faces[1].find( new_quads[1] ) == tmp_addl_faces[1].end() ), 02018 cond2 = ( "split_quads on 0, new_quads on 1" && 02019 tmp_addl_faces[0].find( split_quads[0] ) != tmp_addl_faces[0].end() && 02020 tmp_addl_faces[0].find( split_quads[1] ) != tmp_addl_faces[0].end() && 02021 tmp_addl_faces[1].find( split_quads[0] ) == tmp_addl_faces[1].end() && 02022 tmp_addl_faces[1].find( split_quads[1] ) == tmp_addl_faces[1].end() && 02023 tmp_addl_faces[0].find( new_quads[0] ) == tmp_addl_faces[0].end() && 02024 tmp_addl_faces[0].find( new_quads[1] ) == tmp_addl_faces[0].end() && 02025 tmp_addl_faces[1].find( new_quads[0] ) != tmp_addl_faces[1].end() && 02026 tmp_addl_faces[1].find( new_quads[1] ) != tmp_addl_faces[1].end() ); 02027 02028 assert( cond1 || cond2 ); 02029 #endif 02030 02031 //=============== split edge(s) 02032 for( int j = 0; j < 3; j++ ) 02033 { 02034 if( !split_edges[j] ) break; 02035 02036 // filter add'l faces to only those adj to split_edges[j] 02037 addl_faces[0] = tmp_addl_faces[0]; 02038 addl_faces[1] = tmp_addl_faces[1]; 02039 for( int i = 0; i < 2; i++ ) 02040 { 02041 result = mbImpl->get_adjacencies( &split_edges[j], 1, 2, false, addl_faces[i] );RR; 02042 } 02043 02044 // split edge 02045 result = mtu.split_entity_nonmanifold( split_edges[j], addl_faces[1 - new_side], addl_faces[new_side], 02046 new_edges[j] );RR; 02047 } 02048 02049 //=============== split node(s) 02050 02051 for( int j = 0; j < 2; j++ ) 02052 { 02053 if( !split_nodes[j] ) break; 02054 02055 // if we're splitting multiple edges, there might be other edges that have the split 02056 // node; also need to know which side they're on 02057 Range tmp_addl_edges[2]; 02058 result = foc_get_addl_ents( star_dp1, star_dp2, split_edges, split_nodes[j], tmp_addl_edges );RR; 02059 02060 // also, we need to know which of the split/new edges go 02061 // with the split/new node; new edges go with side 0, split with 1 02062 for( int i = 0; i < 3; i++ ) 02063 { 02064 if( !split_edges[i] ) break; 02065 tmp_addl_edges[new_side].insert( new_edges[i] ); 02066 tmp_addl_edges[1 - new_side].insert( split_edges[i] ); 02067 } 02068 02069 // same for star faces and hexes 02070 for( int i = 0; i < 2; i++ ) 02071 { 02072 std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_edges[i] ) ); 02073 std::copy( star_dp2[i].begin(), star_dp2[i].end(), range_inserter( tmp_addl_edges[i] ) ); 02074 } 02075 02076 // finally, new quads 02077 for( int i = 0; i < 2; i++ ) 02078 tmp_addl_edges[new_side].insert( new_quads[i] ); 02079 02080 // filter the entities, keeping only the ones adjacent to this node 02081 Range addl_edges[2]; 02082 for( int i = 0; i < 2; i++ ) 02083 { 02084 for( Range::reverse_iterator rit = tmp_addl_edges[i].rbegin(); rit != tmp_addl_edges[i].rend(); ++rit ) 02085 { 02086 if( mtu.common_entity( *rit, split_nodes[j], 0 ) ) addl_edges[i].insert( *rit ); 02087 } 02088 } 02089 02090 // now split the node too 02091 result = mtu.split_entity_nonmanifold( split_nodes[j], addl_edges[1 - new_side], addl_edges[new_side], 02092 new_nodes[j] );RR; 02093 } 02094 02095 return MB_SUCCESS; 02096 } 02097 02098 ErrorCode DualTool::foc_get_addl_ents( std::vector< EntityHandle >* star_dp1, 02099 std::vector< EntityHandle >* /*star_dp2*/, 02100 EntityHandle* split_edges, 02101 EntityHandle split_node, 02102 Range* addl_ents ) 02103 { 02104 // if we're splitting 2 edges, there might be other edges that have the split 02105 // node; also need to know which side they're on 02106 02107 // algorithm: for a given star_dp1 (faces) on a side: 02108 // - get all edges adj to all faces -> R1 02109 // - get all edges adj to split_node -> R2 02110 // - R3 = R1 & R2 (edges connected to split_node & adj to a star face) 02111 // - R3 -= split_edges (take split edges off addl_ents) 02112 02113 Range R2; 02114 MeshTopoUtil mtu( mbImpl ); 02115 ErrorCode result = mbImpl->get_adjacencies( &split_node, 1, 1, false, R2 );RR; 02116 Range::iterator rit; 02117 02118 for( int i = 0; i < 2; i++ ) 02119 { 02120 Range R1, R3; 02121 result = mbImpl->get_adjacencies( &star_dp1[i][0], star_dp1[i].size(), 1, false, R1, Interface::UNION );RR; 02122 R3 = intersect( R1, R2 ); 02123 for( int j = 0; j < 3; j++ ) 02124 if( split_edges[j] ) R3.erase( split_edges[j] ); 02125 addl_ents[i].merge( R3 ); 02126 } 02127 02128 return MB_SUCCESS; 02129 } 02130 02131 ErrorCode DualTool::foc_get_stars( EntityHandle* split_quads, 02132 EntityHandle* split_edges, 02133 std::vector< EntityHandle >* star_dp1, 02134 std::vector< EntityHandle >* star_dp2 ) 02135 { 02136 bool on_bdy = false, on_bdy_tmp; 02137 ErrorCode result; 02138 MeshTopoUtil mtu( mbImpl ); 02139 02140 // get the star around the split_edge 02141 std::vector< EntityHandle > qstar, hstar; 02142 unsigned int qpos = 0; 02143 02144 for( int i = 0; i < 3; i++ ) 02145 { 02146 if( !split_edges[i] ) break; 02147 02148 // get the star around this split edge 02149 unsigned int qpos_tmp = 0; 02150 std::vector< EntityHandle > qstar_tmp, hstar_tmp; 02151 result = mtu.star_entities( split_edges[i], qstar_tmp, on_bdy_tmp, 0, &hstar_tmp );RR; 02152 // if we're on the bdy, add a null to the hex star too 02153 if( on_bdy_tmp ) 02154 { 02155 assert( hstar_tmp.size() == qstar_tmp.size() - 1 ); 02156 hstar_tmp.push_back( 0 ); 02157 on_bdy = true; 02158 } 02159 02160 // get the position of first split quad in star 02161 while( qpos_tmp < qstar_tmp.size() && qstar_tmp[qpos_tmp] != split_quads[0] ) 02162 qpos_tmp++; 02163 if( qpos_tmp == qstar_tmp.size() ) return MB_FAILURE; 02164 02165 bool forward; 02166 // 1st iteration is forward by definition 02167 if( 0 == i ) forward = true; 02168 02169 // need to be careful about direction on later iters 02170 else if( hstar[qpos] == hstar_tmp[qpos_tmp] ) 02171 forward = true; 02172 else if( hstar[qpos] == hstar_tmp[( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size()] && 02173 hstar_tmp[qpos_tmp] == hstar[( qpos + qstar.size() - 1 ) % qstar.size()] ) 02174 forward = false; 02175 else 02176 return MB_FAILURE; 02177 02178 if( forward ) 02179 { 02180 // 1st half of star 02181 // save hex right after split_quad[0] first 02182 star_dp2[0].push_back( hstar_tmp[qpos_tmp] ); 02183 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size(); 02184 while( qstar_tmp[qpos_tmp] != split_quads[1] ) 02185 { 02186 star_dp1[0].push_back( qstar_tmp[qpos_tmp] ); 02187 star_dp2[0].push_back( hstar_tmp[qpos_tmp] ); 02188 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size(); 02189 } 02190 // 2nd half of star 02191 // save hex right after split_quad[1] first 02192 star_dp2[1].push_back( hstar_tmp[qpos_tmp] ); 02193 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size(); 02194 while( qstar_tmp[qpos_tmp] != split_quads[0] ) 02195 { 02196 star_dp1[1].push_back( qstar_tmp[qpos_tmp] ); 02197 star_dp2[1].push_back( hstar_tmp[qpos_tmp] ); 02198 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size(); 02199 } 02200 } 02201 else 02202 { 02203 // go in reverse - take prev hex instead of current 02204 // one, and step in reverse 02205 02206 // save hex right after split_quad[0] first 02207 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size(); 02208 star_dp2[0].push_back( hstar_tmp[qpos_tmp] ); 02209 while( qstar_tmp[qpos_tmp] != split_quads[1] ) 02210 { 02211 star_dp1[0].push_back( qstar_tmp[qpos_tmp] ); 02212 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size(); 02213 star_dp2[0].push_back( hstar_tmp[qpos_tmp] ); 02214 } 02215 // 2nd half of star 02216 // save hex right after split_quad[1] first 02217 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size(); 02218 star_dp2[1].push_back( hstar_tmp[qpos_tmp] ); 02219 while( qstar_tmp[qpos_tmp] != split_quads[0] ) 02220 { 02221 star_dp1[1].push_back( qstar_tmp[qpos_tmp] ); 02222 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size(); 02223 star_dp2[1].push_back( hstar_tmp[qpos_tmp] ); 02224 } 02225 } 02226 02227 if( 0 == i ) 02228 { 02229 // if we're on the first iteration, save results and continue, other iters 02230 // get compared to this one 02231 qstar.swap( qstar_tmp ); 02232 hstar.swap( hstar_tmp ); 02233 on_bdy = on_bdy_tmp; 02234 qpos = qpos_tmp; 02235 } 02236 } 02237 02238 // split quads go on list with NULLs, if any, otherwise on 2nd 02239 if( on_bdy ) 02240 { 02241 if( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) != star_dp2[0].end() ) 02242 { 02243 // remove *all* the zeros 02244 star_dp2[0].erase( std::remove( star_dp2[0].begin(), star_dp2[0].end(), 0 ), star_dp2[0].end() ); 02245 // put the split quads on this half 02246 star_dp1[0].push_back( split_quads[0] ); 02247 star_dp1[0].push_back( split_quads[1] ); 02248 } 02249 else 02250 { 02251 star_dp2[1].erase( std::remove( star_dp2[1].begin(), star_dp2[1].end(), 0 ), star_dp2[1].end() ); 02252 // put the split quads on this half 02253 star_dp1[1].push_back( split_quads[0] ); 02254 star_dp1[1].push_back( split_quads[1] ); 02255 } 02256 } 02257 else 02258 { 02259 star_dp1[1].push_back( split_quads[0] ); 02260 star_dp1[1].push_back( split_quads[1] ); 02261 } 02262 02263 // some error checking 02264 if( !( ( ( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) == star_dp1[0].end() && 02265 std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) == star_dp1[0].end() && 02266 std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() && 02267 std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) != star_dp1[1].end() ) || 02268 02269 ( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) == star_dp1[1].end() && 02270 std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) == star_dp1[1].end() && 02271 std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() && 02272 std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) != star_dp1[0].end() ) ) ) ) 02273 { 02274 std::cerr << "foc_get_stars: both split quads should be on the same star list half and not " 02275 << "on the other, failed" << std::endl; 02276 return MB_FAILURE; 02277 } 02278 02279 if( !( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) == star_dp2[0].end() && 02280 std::find( star_dp2[1].begin(), star_dp2[1].end(), 0 ) == star_dp2[1].end() ) ) 02281 { 02282 std::cerr << "foc_get_stars: no NULLs on the hstar lists, failed"; 02283 return MB_FAILURE; 02284 } 02285 02286 return MB_SUCCESS; 02287 } 02288 02289 ErrorCode DualTool::foc_delete_dual( EntityHandle* split_quads, EntityHandle* split_edges, Range& hexes ) 02290 { 02291 // special delete dual procedure, because in some cases we need to delete 02292 // a sheet too since it'll get merged into another 02293 02294 // figure out whether we'll need to delete a sheet 02295 EntityHandle sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[0] ) ); 02296 if( split_edges[1] ) sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[1] ) ); 02297 EntityHandle chordl = get_dual_hyperplane( get_dual_entity( split_quads[0] ) ); 02298 EntityHandle chordr = get_dual_hyperplane( get_dual_entity( split_quads[1] ) ); 02299 assert( 0 != sheet1 && 0 != chordl && 0 != chordr ); 02300 Range parentsl, parentsr; 02301 ErrorCode result = mbImpl->get_parent_meshsets( chordl, parentsl ); 02302 if( MB_SUCCESS != result ) return result; 02303 result = mbImpl->get_parent_meshsets( chordr, parentsr ); 02304 if( MB_SUCCESS != result ) return result; 02305 parentsl.erase( sheet1 ); 02306 parentsr.erase( sheet1 ); 02307 02308 // before deciding which one to delete, collect the other cells which must 02309 // be deleted, and all the chords/sheets they're on 02310 Range adj_ents, dual_ents, cells1or2; 02311 for( int i = 0; i < 3; i++ ) 02312 { 02313 result = mbImpl->get_adjacencies( hexes, i, false, adj_ents, Interface::UNION ); 02314 if( MB_SUCCESS != result ) return result; 02315 } 02316 02317 // cache any adjacent hexes, for rebuilding the dual later 02318 result = mbImpl->get_adjacencies( adj_ents, 3, false, hexes, Interface::UNION ); 02319 if( MB_SUCCESS != result ) return result; 02320 02321 for( Range::iterator rit = adj_ents.begin(); rit != adj_ents.end(); ++rit ) 02322 { 02323 EntityHandle this_ent = get_dual_entity( *rit ); 02324 dual_ents.insert( this_ent ); 02325 int dim = mbImpl->dimension_from_handle( this_ent ); 02326 if( 1 == dim || 2 == dim ) cells1or2.insert( this_ent ); 02327 } 02328 02329 Range dual_hps; 02330 for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit ) 02331 dual_hps.insert( get_dual_hyperplane( *rit ) ); 02332 02333 result = delete_dual_entities( dual_ents ); 02334 if( MB_SUCCESS != result ) return result; 02335 02336 // now decide which sheet to delete (to be merged into the other) 02337 EntityHandle sheet_delete = 0; 02338 if( is_blind( *parentsl.begin() ) ) 02339 sheet_delete = *parentsl.begin(); 02340 else if( is_blind( *parentsr.begin() ) ) 02341 sheet_delete = *parentsr.begin(); 02342 else 02343 { 02344 // neither is blind, take the one with fewer cells 02345 Range tmp_ents; 02346 int numl, numr; 02347 result = mbImpl->get_number_entities_by_handle( *parentsl.begin(), numl ); 02348 if( MB_SUCCESS != result ) return result; 02349 result = mbImpl->get_number_entities_by_handle( *parentsr.begin(), numr ); 02350 if( MB_SUCCESS != result ) return result; 02351 sheet_delete = ( numl > numr ? *parentsr.begin() : *parentsl.begin() ); 02352 } 02353 assert( 0 != sheet_delete ); 02354 02355 // after deleting cells, check for empty chords & sheets, and delete those too 02356 for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit ) 02357 { 02358 Range tmp_ents; 02359 result = mbImpl->get_entities_by_handle( *rit, tmp_ents ); 02360 if( MB_SUCCESS != result ) return result; 02361 if( tmp_ents.empty() ) 02362 { 02363 result = mbImpl->delete_entities( &( *rit ), 1 ); 02364 if( MB_SUCCESS != result ) return result; 02365 } 02366 else if( *rit == sheet_delete ) 02367 { 02368 // delete the sheet 02369 result = mbImpl->delete_entities( &( *rit ), 1 ); 02370 if( MB_SUCCESS != result ) return result; 02371 } 02372 } 02373 02374 // now just to be safe, add the hexes bridge-adjacent across vertices 02375 // to the hexes we already have 02376 Range tmp_hexes; 02377 MeshTopoUtil mtu( mbImpl ); 02378 for( Range::iterator rit = hexes.begin(); rit != hexes.end(); ++rit ) 02379 { 02380 result = mtu.get_bridge_adjacencies( *rit, 0, 3, tmp_hexes ); 02381 if( MB_SUCCESS != result ) return result; 02382 } 02383 hexes.merge( tmp_hexes ); 02384 02385 return MB_SUCCESS; 02386 } 02387 02388 //! returns true if all vertices are dual to hexes (not faces) 02389 bool DualTool::is_blind( const EntityHandle chord_or_sheet ) 02390 { 02391 // must be an entity set 02392 if( TYPE_FROM_HANDLE( chord_or_sheet ) != MBENTITYSET ) return false; 02393 02394 // get the vertices 02395 Range verts, ents; 02396 ErrorCode result = mbImpl->get_entities_by_handle( chord_or_sheet, ents ); 02397 if( MB_SUCCESS != result || ents.empty() ) return false; 02398 02399 result = mbImpl->get_adjacencies( ents, 0, false, verts, Interface::UNION ); 02400 if( MB_SUCCESS != result || verts.empty() ) return false; 02401 02402 for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit ) 02403 { 02404 // get dual entity for this vertex 02405 EntityHandle dual_ent = get_dual_entity( *rit ); 02406 if( 0 == dual_ent ) continue; 02407 if( TYPE_FROM_HANDLE( dual_ent ) == MBQUAD ) return false; 02408 } 02409 02410 // if none of the vertices' duals were quads, chord_or_sheet must be blind 02411 return true; 02412 } 02413 02414 //! given a 1-cell and a chord, return the neighboring vertices on the 02415 //! chord, in the same order as the 1-cell's vertices 02416 ErrorCode DualTool::get_opposite_verts( const EntityHandle middle_edge, const EntityHandle chord, EntityHandle* verts ) 02417 { 02418 // get the edges on the chord, in order, and move to middle_edge 02419 std::vector< EntityHandle > chord_edges; 02420 const EntityHandle* connect; 02421 int num_connect; 02422 02423 ErrorCode result = mbImpl->get_entities_by_handle( chord, chord_edges );RR; 02424 std::vector< EntityHandle >::iterator vit = std::find( chord_edges.begin(), chord_edges.end(), middle_edge ); 02425 result = mbImpl->get_connectivity( middle_edge, connect, num_connect );RR; 02426 02427 if( 02428 // middle_edge isn't on this chord 02429 vit == chord_edges.end() || 02430 // chord only has 1 edge 02431 chord_edges.size() == 1 || 02432 // middle_edge is at beginning or end and chord isn't blind 02433 ( ( vit == chord_edges.begin() || vit == chord_edges.end() - 1 ) && !is_blind( chord ) ) ) 02434 return MB_FAILURE; 02435 02436 else if( chord_edges.size() == 2 ) 02437 { 02438 // easier if it's a 2-edge blind chord, just get vertices in opposite order 02439 verts[0] = connect[1]; 02440 verts[1] = connect[0]; 02441 return MB_SUCCESS; 02442 } 02443 02444 // get vertices with the prev edge & subtract vertices of 1-cell 02445 if( vit == chord_edges.begin() ) 02446 vit = chord_edges.end() - 1; 02447 else 02448 --vit; 02449 Range dum_connect, middle_connect; 02450 result = mbImpl->get_connectivity( &middle_edge, 1, middle_connect );RR; 02451 result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR; 02452 dum_connect = subtract( dum_connect, middle_connect ); 02453 if( dum_connect.size() != 1 ) 02454 { 02455 std::cerr << "Trouble traversing chord." << std::endl; 02456 return MB_FAILURE; 02457 } 02458 02459 // put in verts[0] 02460 verts[0] = *dum_connect.begin(); 02461 02462 // same with prev edge 02463 ++vit; 02464 if( vit == chord_edges.end() ) vit = chord_edges.begin(); 02465 ++vit; 02466 dum_connect.clear(); 02467 result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR; 02468 dum_connect = subtract( dum_connect, middle_connect ); 02469 if( dum_connect.size() != 1 ) 02470 { 02471 std::cerr << "Trouble traversing chord." << std::endl; 02472 return MB_FAILURE; 02473 } 02474 02475 // put in verts[1] 02476 verts[1] = *dum_connect.begin(); 02477 02478 // if verts[0] and 1st vertex of 1cell don't have common edge, switch verts 02479 MeshTopoUtil mtu( mbImpl ); 02480 if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) ) 02481 { 02482 EntityHandle dum_h = verts[0]; 02483 verts[0] = verts[1]; 02484 verts[1] = dum_h; 02485 } 02486 02487 if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) ) 02488 { 02489 std::cerr << "Trouble traversing chord." << std::endl; 02490 return MB_FAILURE; 02491 } 02492 02493 return MB_SUCCESS; 02494 } 02495 02496 ErrorCode DualTool::get_dual_entities( const EntityHandle dual_ent, 02497 Range* dcells, 02498 Range* dedges, 02499 Range* dverts, 02500 Range* dverts_loop, 02501 Range* dedges_loop ) 02502 { 02503 ErrorCode result = MB_SUCCESS; 02504 02505 if( NULL != dcells ) 02506 { 02507 result = mbImpl->get_entities_by_type( dual_ent, MBPOLYGON, *dcells ); 02508 if( MB_SUCCESS != result ) return result; 02509 } 02510 02511 if( NULL != dedges ) 02512 { 02513 if( NULL != dcells ) 02514 result = mbImpl->get_adjacencies( *dcells, 1, false, *dedges, Interface::UNION ); 02515 else 02516 result = mbImpl->get_entities_by_type( dual_ent, MBEDGE, *dedges ); 02517 02518 if( MB_SUCCESS != result ) return result; 02519 } 02520 02521 if( NULL != dverts ) 02522 { 02523 if( NULL != dcells ) 02524 result = mbImpl->get_adjacencies( *dcells, 0, false, *dverts, Interface::UNION ); 02525 else if( NULL != dedges ) 02526 result = mbImpl->get_adjacencies( *dedges, 0, false, *dverts, Interface::UNION ); 02527 else 02528 { 02529 Range all_ents; 02530 result = mbImpl->get_entities_by_handle( dual_ent, all_ents );RR; 02531 result = mbImpl->get_adjacencies( all_ents, 0, false, *dverts, Interface::UNION ); 02532 } 02533 02534 if( MB_SUCCESS != result ) return result; 02535 } 02536 02537 if( NULL != dverts_loop && NULL != dverts ) 02538 { 02539 static std::vector< EntityHandle > dual_ents; 02540 dual_ents.resize( dverts->size() ); 02541 result = mbImpl->tag_get_data( dualEntity_tag(), *dverts, &dual_ents[0] ); 02542 if( MB_SUCCESS != result ) return result; 02543 Range::iterator rit; 02544 unsigned int i; 02545 for( rit = dverts->begin(), i = 0; rit != dverts->end(); ++rit, i++ ) 02546 if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBQUAD ) dverts_loop->insert( *rit ); 02547 } 02548 02549 if( NULL != dedges_loop && NULL != dedges ) 02550 { 02551 static std::vector< EntityHandle > dual_ents; 02552 dual_ents.resize( dedges->size() ); 02553 result = mbImpl->tag_get_data( dualEntity_tag(), *dedges, &dual_ents[0] ); 02554 if( MB_SUCCESS != result ) return result; 02555 Range::iterator rit; 02556 unsigned int i; 02557 for( rit = dedges->begin(), i = 0; rit != dedges->end(); ++rit, i++ ) 02558 if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBEDGE ) dedges_loop->insert( *rit ); 02559 } 02560 02561 return result; 02562 } 02563 02564 ErrorCode DualTool::list_entities( const EntityHandle* entities, const int num_entities ) const 02565 { 02566 Range temp_range; 02567 ErrorCode result; 02568 if( NULL == entities && 0 == num_entities ) 02569 return mbImpl->list_entities( entities, num_entities ); 02570 02571 else if( NULL == entities && 0 < num_entities ) 02572 { 02573 02574 // list all entities of all types 02575 std::cout << std::endl; 02576 for( EntityType this_type = MBVERTEX; this_type < MBMAXTYPE; this_type++ ) 02577 { 02578 result = mbImpl->get_entities_by_type( 0, this_type, temp_range ); 02579 if( MB_SUCCESS != result ) return result; 02580 } 02581 } 02582 02583 else 02584 { 02585 std::copy( entities, entities + num_entities, range_inserter( temp_range ) ); 02586 } 02587 02588 return list_entities( temp_range ); 02589 } 02590 02591 ErrorCode DualTool::list_entities( const Range& entities ) const 02592 { 02593 // now print each entity, listing the dual information first then calling Interface to do 02594 // the rest 02595 ErrorCode result = MB_SUCCESS, tmp_result; 02596 for( Range::const_iterator iter = entities.begin(); iter != entities.end(); ++iter ) 02597 { 02598 EntityType this_type = TYPE_FROM_HANDLE( *iter ); 02599 std::cout << CN::EntityTypeName( this_type ) << " " << ID_FROM_HANDLE( *iter ) << ":" << std::endl; 02600 02601 EntityHandle dual_ent = get_dual_entity( *iter ); 02602 if( 0 != dual_ent ) 02603 { 02604 std::cout << "Dual to " << CN::EntityTypeName( mbImpl->type_from_handle( dual_ent ) ) << " " 02605 << mbImpl->id_from_handle( dual_ent ) << std::endl; 02606 } 02607 02608 if( TYPE_FROM_HANDLE( *iter ) == MBENTITYSET ) 02609 { 02610 EntityHandle chord = 0, sheet = 0; 02611 int id; 02612 result = mbImpl->tag_get_data( dualCurve_tag(), &( *iter ), 1, &chord ); 02613 if( MB_SUCCESS != result ) return result; 02614 result = mbImpl->tag_get_data( dualSurface_tag(), &( *iter ), 1, &sheet ); 02615 if( MB_SUCCESS != result ) return result; 02616 result = mbImpl->tag_get_data( globalId_tag(), &( *iter ), 1, &id ); 02617 if( MB_SUCCESS != result ) return result; 02618 02619 if( 0 != chord ) std::cout << "(Dual chord " << id << ")" << std::endl; 02620 if( 0 != sheet ) std::cout << "(Dual sheet " << id << ")" << std::endl; 02621 } 02622 02623 tmp_result = mbImpl->list_entity( *iter ); 02624 if( MB_SUCCESS != tmp_result ) result = tmp_result; 02625 } 02626 02627 return result; 02628 } 02629 02630 ErrorCode DualTool::face_shrink( EntityHandle odedge ) 02631 { 02632 // some preliminary checking 02633 if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE; 02634 02635 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 02636 02637 std::cout << "FS("; 02638 print_cell( odedge ); 02639 std::cout << ")" << std::endl; 02640 02641 EntityHandle quads[4], hexes[2]; 02642 std::vector< EntityHandle > connects[4], side_quads[2]; 02643 02644 // get the quads along the chord through the 2 hexes, and the vertices 02645 // for those quads 02646 ErrorCode result = fs_get_quads( odedge, quads, hexes, connects ); 02647 if( MB_SUCCESS != result ) return result; 02648 02649 // flip/rotate connect arrays so they align & are same sense 02650 result = fs_check_quad_sense( hexes[0], quads[0], connects ); 02651 if( MB_SUCCESS != result ) return result; 02652 02653 // get the quad loops along the "side" surfaces 02654 result = fs_get_quad_loops( hexes, connects, side_quads ); 02655 if( MB_SUCCESS != result ) return result; 02656 02657 // ok, done with setup; now delete dual entities affected by this operation, 02658 // which is all the entities adjacent to vertices of dual edge 02659 Range adj_verts, adj_edges, dual_ents, cells1or2; 02660 MeshTopoUtil mtu( mbImpl ); 02661 result = mtu.get_bridge_adjacencies( odedge, 0, 1, adj_edges ); 02662 if( MB_SUCCESS != result ) return result; 02663 result = mbImpl->get_adjacencies( adj_edges, 0, false, adj_verts, Interface::UNION ); 02664 if( MB_SUCCESS != result ) return result; 02665 for( int i = 1; i <= 3; i++ ) 02666 { 02667 result = mbImpl->get_adjacencies( adj_verts, i, false, dual_ents, Interface::UNION ); 02668 if( MB_SUCCESS != result ) return result; 02669 } 02670 02671 // before deleting dual, grab the 1- and 2-cells 02672 for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit ) 02673 { 02674 int dim = mbImpl->dimension_from_handle( *rit ); 02675 if( 1 == dim || 2 == dim ) cells1or2.insert( *rit ); 02676 } 02677 Range dual_hps; 02678 for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit ) 02679 dual_hps.insert( get_dual_hyperplane( *rit ) ); 02680 02681 dual_ents.insert( odedge ); 02682 result = delete_dual_entities( dual_ents ); 02683 if( MB_SUCCESS != result ) return result; 02684 02685 // after deleting cells, check for empty chords & sheets, and delete those too 02686 for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit ) 02687 { 02688 Range tmp_ents; 02689 result = mbImpl->get_entities_by_handle( *rit, tmp_ents ); 02690 if( MB_SUCCESS != result ) return result; 02691 if( tmp_ents.empty() ) 02692 { 02693 result = mbImpl->delete_entities( &( *rit ), 1 ); 02694 if( MB_SUCCESS != result ) return result; 02695 } 02696 } 02697 02698 // remove any explicit adjacencies between side quads and hexes; don't check 02699 // for error, since there might not be adjacencies 02700 for( int i = 0; i < 4; i++ ) 02701 { 02702 for( int j = 0; j < 2; j++ ) 02703 { 02704 result = mbImpl->remove_adjacencies( side_quads[j][i], &hexes[j], 1 ); 02705 } 02706 } 02707 02708 // make inner ring of vertices 02709 // get centroid of quad2 02710 double q2coords[12], avg[3] = { 0.0, 0.0, 0.0 }; 02711 result = mbImpl->get_coords( &connects[1][0], 4, q2coords ); 02712 if( MB_SUCCESS != result ) return result; 02713 for( int i = 0; i < 4; i++ ) 02714 { 02715 avg[0] += q2coords[3 * i]; 02716 avg[1] += q2coords[3 * i + 1]; 02717 avg[2] += q2coords[3 * i + 2]; 02718 } 02719 avg[0] *= .25; 02720 avg[1] *= .25; 02721 avg[2] *= .25; 02722 // position new vertices 02723 connects[3].resize( 4 ); 02724 for( int i = 0; i < 4; i++ ) 02725 { 02726 q2coords[3 * i] = avg[0] + .25 * ( q2coords[3 * i] - avg[0] ); 02727 q2coords[3 * i + 1] = avg[1] + .25 * ( q2coords[3 * i + 1] - avg[1] ); 02728 q2coords[3 * i + 2] = avg[2] + .25 * ( q2coords[3 * i + 2] - avg[2] ); 02729 result = mbImpl->create_vertex( &q2coords[3 * i], connects[3][i] ); 02730 if( MB_SUCCESS != result ) return result; 02731 } 02732 02733 // ok, now have the 4 connectivity arrays for 4 quads; construct hexes 02734 EntityHandle hconnect[8], new_hexes[4]; 02735 int new_hex_ids[4]; 02736 02737 for( int i = 0; i < 4; i++ ) 02738 { 02739 int i1 = i, i2 = ( i + 1 ) % 4; 02740 hconnect[0] = connects[0][i1]; 02741 hconnect[1] = connects[0][i2]; 02742 hconnect[2] = connects[3][i2]; 02743 hconnect[3] = connects[3][i1]; 02744 02745 hconnect[4] = connects[1][i1]; 02746 hconnect[5] = connects[1][i2]; 02747 hconnect[6] = connects[2][i2]; 02748 hconnect[7] = connects[2][i1]; 02749 02750 result = mbImpl->create_element( MBHEX, hconnect, 8, new_hexes[i] ); 02751 if( MB_SUCCESS != result ) return result; 02752 02753 // test for equiv entities from the side quads, and make explicit adjacencies 02754 // if there are any 02755 for( int j = 0; j < 2; j++ ) 02756 { 02757 if( mtu.equivalent_entities( side_quads[j][i] ) ) 02758 { 02759 result = mbImpl->add_adjacencies( side_quads[j][i], &new_hexes[i], 1, false ); 02760 if( MB_SUCCESS != result ) return result; 02761 } 02762 } 02763 02764 new_hex_ids[i] = ++maxHexId; 02765 } 02766 02767 // set the global id tag on the new hexes 02768 result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 4, new_hex_ids ); 02769 if( MB_SUCCESS != result ) return result; 02770 02771 // now fixup other two hexes; start by getting hex through quads 0, 1 02772 // make this first hex switch to the other side, to make the dual look like 02773 // a hex push 02774 int tmp_ids[2]; 02775 result = mbImpl->tag_get_data( globalId_tag(), hexes, 2, tmp_ids ); 02776 if( MB_SUCCESS != result ) return result; 02777 02778 result = mbImpl->delete_entities( hexes, 2 ); 02779 if( MB_SUCCESS != result ) return result; 02780 result = mbImpl->delete_entities( &quads[1], 1 ); 02781 if( MB_SUCCESS != result ) return result; 02782 for( int i = 0; i < 4; i++ ) 02783 { 02784 hconnect[i] = connects[3][i]; 02785 hconnect[4 + i] = connects[2][i]; 02786 } 02787 result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[0] ); 02788 if( MB_SUCCESS != result ) return result; 02789 02790 for( int i = 0; i < 4; i++ ) 02791 { 02792 hconnect[i] = connects[0][i]; 02793 hconnect[4 + i] = connects[3][i]; 02794 } 02795 result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[1] ); 02796 if( MB_SUCCESS != result ) return result; 02797 02798 // check for and fix any explicit adjacencies on either end quad 02799 if( mtu.equivalent_entities( quads[0] ) ) mbImpl->add_adjacencies( quads[0], &hexes[1], 1, false ); 02800 if( mtu.equivalent_entities( quads[2] ) ) mbImpl->add_adjacencies( quads[2], &hexes[0], 1, false ); 02801 02802 // re-set the global ids for the hexes to what they were 02803 result = mbImpl->tag_set_data( globalId_tag(), hexes, 2, tmp_ids ); 02804 if( MB_SUCCESS != result ) return result; 02805 02806 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 02807 02808 // now update the dual 02809 Range tmph; 02810 result = mtu.get_bridge_adjacencies( hexes[0], 0, 3, tmph ); 02811 if( MB_SUCCESS != result ) return result; 02812 result = mtu.get_bridge_adjacencies( hexes[1], 0, 3, tmph ); 02813 if( MB_SUCCESS != result ) return result; 02814 tmph.insert( hexes[1] ); 02815 result = construct_hex_dual( tmph ); 02816 if( MB_SUCCESS != result ) return result; 02817 02818 return result; 02819 } 02820 02821 ErrorCode DualTool::fs_get_quad_loops( EntityHandle* hexes, 02822 std::vector< EntityHandle >* connects, 02823 std::vector< EntityHandle >* side_quads ) 02824 { 02825 for( int i = 0; i < 4; i++ ) 02826 { 02827 for( int j = 0; j < 2; j++ ) 02828 { 02829 Range adj_ents, dum_quads; 02830 adj_ents.insert( hexes[j] ); 02831 adj_ents.insert( connects[j][i] ); 02832 adj_ents.insert( connects[j][( i + 1 ) % 4] ); 02833 adj_ents.insert( connects[j + 1][i] ); 02834 adj_ents.insert( connects[j + 1][( i + 1 ) % 4] ); 02835 02836 ErrorCode result = mbImpl->get_adjacencies( adj_ents, 2, false, dum_quads ); 02837 if( MB_SUCCESS != result ) return result; 02838 assert( 1 == dum_quads.size() ); 02839 side_quads[j].push_back( *dum_quads.begin() ); 02840 } 02841 } 02842 02843 return MB_SUCCESS; 02844 } 02845 02846 ErrorCode DualTool::fs_check_quad_sense( EntityHandle hex0, EntityHandle quad0, std::vector< EntityHandle >* connects ) 02847 { 02848 // check sense of 0th quad wrt hex; since sense is out of element, 02849 // switch if quad is NOT reversed wrt hex 02850 int dum1, dum2, sense = 0; 02851 ErrorCode result = mbImpl->side_number( hex0, quad0, dum1, sense, dum2 ); 02852 if( MB_SUCCESS != result ) return result; 02853 assert( 0 != sense ); 02854 if( 1 == sense ) 02855 { 02856 // just switch sense of this one; others will get switched next 02857 EntityHandle dum = connects[0][0]; 02858 connects[0][0] = connects[0][2]; 02859 connects[0][2] = dum; 02860 } 02861 02862 // check sense of 1st, 2nd quads, rotate if necessary to align connect arrays 02863 int index0 = -1, index2 = -1, sense0 = 0, sense2 = 0; 02864 MeshTopoUtil mtu( mbImpl ); 02865 for( int i = 0; i < 4; i++ ) 02866 { 02867 if( 0 != mtu.common_entity( connects[0][0], connects[1][i], 1 ) ) 02868 { 02869 index0 = i; 02870 if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 1 ) % 4], 1 ) ) 02871 sense0 = 1; 02872 else if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 4 - 1 ) % 4], 1 ) ) 02873 sense0 = -1; 02874 break; 02875 } 02876 } 02877 02878 assert( index0 != -1 && sense0 != 0 ); 02879 02880 if( sense0 == -1 ) 02881 { 02882 EntityHandle dumh = connects[1][0]; 02883 connects[1][0] = connects[1][2]; 02884 connects[1][2] = dumh; 02885 if( index0 % 2 == 0 ) index0 = ( index0 + 2 ) % 4; 02886 } 02887 02888 if( index0 != 0 ) 02889 { 02890 std::vector< EntityHandle > tmpc; 02891 for( int i = 0; i < 4; i++ ) 02892 tmpc.push_back( connects[1][( index0 + i ) % 4] ); 02893 connects[1].swap( tmpc ); 02894 } 02895 02896 for( int i = 0; i < 4; i++ ) 02897 { 02898 if( 0 != mtu.common_entity( connects[1][0], connects[2][i], 1 ) ) 02899 { 02900 index2 = i; 02901 if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 1 ) % 4], 1 ) ) 02902 sense2 = 1; 02903 else if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 4 - 1 ) % 4], 1 ) ) 02904 sense2 = -1; 02905 break; 02906 } 02907 } 02908 02909 assert( index2 != -1 && sense2 != 0 ); 02910 02911 if( sense2 == -1 ) 02912 { 02913 EntityHandle dumh = connects[2][0]; 02914 connects[2][0] = connects[2][2]; 02915 connects[2][2] = dumh; 02916 if( index2 % 2 == 0 ) index2 = ( index2 + 2 ) % 4; 02917 } 02918 02919 if( index2 != 0 ) 02920 { 02921 std::vector< EntityHandle > tmpc; 02922 for( int i = 0; i < 4; i++ ) 02923 tmpc.push_back( connects[2][( index2 + i ) % 4] ); 02924 connects[2].swap( tmpc ); 02925 } 02926 02927 return MB_SUCCESS; 02928 } 02929 02930 //! effect reverse face shrink operation 02931 ErrorCode DualTool::rev_face_shrink( EntityHandle odedge ) 02932 { 02933 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 02934 02935 // some preliminary checking 02936 if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE; 02937 02938 std::cout << "-FS("; 02939 print_cell( odedge ); 02940 std::cout << ")" << std::endl; 02941 02942 EntityHandle quads[4], hexes[2]; 02943 std::vector< EntityHandle > connects[4], side_quads[2]; 02944 02945 // get three quads (shared quad & 2 end quads), hexes, and quad 02946 // connects 02947 ErrorCode result = fs_get_quads( odedge, quads, hexes, connects ); 02948 if( MB_SUCCESS != result ) return result; 02949 02950 // adjust sense & rotation so they're aligned, together & wrt first 02951 // hex 02952 result = fs_check_quad_sense( hexes[0], quads[0], connects ); 02953 if( MB_SUCCESS != result ) return result; 02954 02955 result = fsr_get_fourth_quad( connects, side_quads ); 02956 if( MB_SUCCESS != result ) 02957 { 02958 std::cout << "Can't do -FS here, two hexes must be adjacent to ring of 4 hexes." << std::endl; 02959 return result; 02960 } 02961 02962 Range adj_ents, outer_hexes, all_adjs; 02963 02964 // first get the entities connected to interior 4 verts 02965 for( int i = 1; i <= 3; i++ ) 02966 { 02967 result = mbImpl->get_adjacencies( &connects[1][0], 4, i, false, adj_ents, Interface::UNION ); 02968 if( MB_SUCCESS != result ) return result; 02969 } 02970 02971 // next get all entities adjacent to those; these will have their dual 02972 // entities deleted 02973 for( int i = 0; i < 3; i++ ) 02974 { 02975 result = mbImpl->get_adjacencies( adj_ents, i, false, all_adjs, Interface::UNION ); 02976 if( MB_SUCCESS != result ) return result; 02977 } 02978 02979 // get the dual entities and delete them 02980 Range dual_ents, dual_hps; 02981 for( Range::iterator rit = all_adjs.begin(); rit != all_adjs.end(); ++rit ) 02982 { 02983 EntityHandle this_ent = get_dual_entity( *rit ); 02984 dual_ents.insert( this_ent ); 02985 } 02986 02987 // before deleting dual, grab the 1- and 2-cells 02988 for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit ) 02989 { 02990 int dim = mbImpl->dimension_from_handle( *rit ); 02991 if( 1 == dim || 2 == dim ) dual_hps.insert( get_dual_hyperplane( *rit ) ); 02992 } 02993 02994 result = delete_dual_entities( dual_ents ); 02995 if( MB_SUCCESS != result ) return result; 02996 02997 // after deleting cells, check for empty chords & sheets, and delete those too 02998 for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit ) 02999 { 03000 Range tmp_ents; 03001 result = mbImpl->get_entities_by_handle( *rit, tmp_ents ); 03002 if( MB_SUCCESS != result ) return result; 03003 if( tmp_ents.empty() ) 03004 { 03005 result = mbImpl->delete_entities( &( *rit ), 1 ); 03006 if( MB_SUCCESS != result ) return result; 03007 } 03008 } 03009 03010 // before re-connecting two hexes, check for existing quad on 4th quad vertices; 03011 // if there is a quad there, need to add explicit adjs to any adj hexes, since 03012 // by definition there'll be another quad on those vertices 03013 bool need_explicit = false; 03014 Range adj_quads; 03015 result = mbImpl->get_adjacencies( &connects[3][0], 4, 2, false, adj_quads ); 03016 if( MB_MULTIPLE_ENTITIES_FOUND == result || !adj_quads.empty() ) 03017 { 03018 // there's already a quad for these 4 vertices; by definition, 03019 // we'll be creating equivalent entities, so that original quad 03020 // needs explicit adj's to its bounding elements 03021 need_explicit = true; 03022 for( Range::iterator rit = adj_quads.begin(); rit != adj_quads.end(); ++rit ) 03023 { 03024 Range adj_hexes; 03025 result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, adj_hexes );RR; 03026 result = mbImpl->add_adjacencies( *rit, adj_hexes, false );RR; 03027 } 03028 } 03029 03030 // re-connect the two hexes 03031 std::vector< EntityHandle > new_connect; 03032 std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) ); 03033 std::copy( connects[2].begin(), connects[2].end(), std::back_inserter( new_connect ) ); 03034 result = mbImpl->set_connectivity( hexes[0], &new_connect[0], 8 ); 03035 if( MB_SUCCESS != result ) return result; 03036 03037 new_connect.clear(); 03038 std::copy( connects[0].begin(), connects[0].end(), std::back_inserter( new_connect ) ); 03039 std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) ); 03040 result = mbImpl->set_connectivity( hexes[1], &new_connect[0], 8 ); 03041 if( MB_SUCCESS != result ) return result; 03042 03043 // test for equiv entities from the side quads, and make explicit 03044 // adjacencies if there are any 03045 MeshTopoUtil mtu( mbImpl ); 03046 for( int j = 0; j < 2; j++ ) 03047 { 03048 for( int i = 0; i < 4; i++ ) 03049 { 03050 if( mtu.equivalent_entities( side_quads[j][i] ) ) 03051 { 03052 result = mbImpl->add_adjacencies( side_quads[j][i], &hexes[j], 1, false ); 03053 if( MB_SUCCESS != result ) return result; 03054 } 03055 } 03056 } 03057 03058 // remove hexes we want to keep 03059 adj_ents.erase( hexes[0] ); 03060 adj_ents.erase( hexes[1] ); 03061 03062 // delete the other interior entities 03063 result = mbImpl->delete_entities( adj_ents ); 03064 if( MB_SUCCESS != result ) return result; 03065 03066 EntityHandle new_quad; 03067 result = mbImpl->create_element( MBQUAD, &connects[3][0], 4, new_quad );RR; 03068 if( need_explicit ) 03069 { 03070 result = mbImpl->add_adjacencies( new_quad, hexes, 2, false );RR; 03071 } 03072 03073 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies(); 03074 03075 // now update the dual 03076 result = construct_hex_dual( hexes, 2 ); 03077 if( MB_SUCCESS != result ) return result; 03078 03079 return MB_SUCCESS; 03080 } 03081 03082 ErrorCode DualTool::fsr_get_fourth_quad( std::vector< EntityHandle >* connects, 03083 std::vector< EntityHandle >* side_quads ) 03084 { 03085 // given the first three quad connectivities in ordered vectors, get the fourth, 03086 // where the fourth is really the 4 vertices originally shared by the 2 hexes 03087 // before the face shrink on them 03088 03089 // vertex on 4th quad is in quad adj to other 3 verts 03090 for( int i = 0; i < 4; i++ ) 03091 { 03092 Range start_verts, tmp_verts, quads; 03093 for( int j = 0; j < 3; j++ ) 03094 start_verts.insert( connects[j][i] ); 03095 ErrorCode result = mbImpl->get_adjacencies( start_verts, 2, false, quads ); 03096 if( MB_SUCCESS != result ) return result; 03097 assert( quads.size() == 1 ); 03098 result = mbImpl->get_adjacencies( &( *quads.begin() ), 1, 0, false, tmp_verts );RR; 03099 tmp_verts = subtract( tmp_verts, start_verts ); 03100 assert( 1 == tmp_verts.size() ); 03101 connects[3].push_back( *tmp_verts.begin() ); 03102 } 03103 03104 // now get the side quads 03105 for( int i = 0; i < 4; i++ ) 03106 { 03107 Range dum_ents, hexes; 03108 dum_ents.insert( connects[1][i] ); 03109 dum_ents.insert( connects[1][( i + 1 ) % 4] ); 03110 dum_ents.insert( connects[3][i] ); 03111 ErrorCode result = mbImpl->get_adjacencies( dum_ents, 3, false, hexes ); 03112 if( MB_SUCCESS != result ) return result; 03113 assert( 1 == hexes.size() ); 03114 03115 hexes.insert( connects[0][i] ); 03116 hexes.insert( connects[0][( i + 1 ) % 4] ); 03117 hexes.insert( connects[3][i] ); 03118 hexes.insert( connects[3][( i + 1 ) % 4] ); 03119 dum_ents.clear(); 03120 result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents ); 03121 if( MB_SUCCESS != result ) return result; 03122 assert( dum_ents.size() == 1 ); 03123 side_quads[0].push_back( *dum_ents.begin() ); 03124 03125 hexes.erase( connects[0][i] ); 03126 hexes.erase( connects[0][( i + 1 ) % 4] ); 03127 hexes.insert( connects[2][i] ); 03128 hexes.insert( connects[2][( i + 1 ) % 4] ); 03129 dum_ents.clear(); 03130 result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents ); 03131 if( MB_SUCCESS != result ) return result; 03132 side_quads[1].push_back( *dum_ents.begin() ); 03133 } 03134 03135 return MB_SUCCESS; 03136 } 03137 03138 ErrorCode DualTool::fs_get_quads( EntityHandle odedge, 03139 EntityHandle* quads, 03140 EntityHandle* hexes, 03141 std::vector< EntityHandle >* connects ) 03142 { 03143 // need to get the three quads along the chord 03144 EntityHandle chord = get_dual_hyperplane( odedge ); 03145 if( 0 == chord ) return MB_FAILURE; 03146 03147 std::vector< EntityHandle > edges; 03148 ErrorCode result = mbImpl->get_entities_by_handle( chord, edges ); 03149 if( MB_FAILURE == result ) return result; 03150 03151 std::vector< EntityHandle >::iterator vit = std::find( edges.begin(), edges.end(), odedge ); 03152 // shouldn't be first or last edge on chord 03153 if( vit == edges.end() || *edges.begin() == *vit || *edges.rbegin() == *vit ) return MB_FAILURE; 03154 03155 // get quads/connectivity for first 3 quads 03156 quads[0] = get_dual_entity( *( vit - 1 ) ); 03157 quads[1] = get_dual_entity( *vit ); 03158 quads[2] = get_dual_entity( *( vit + 1 ) ); 03159 for( int i = 0; i < 3; i++ ) 03160 { 03161 result = mbImpl->get_connectivity( &quads[i], 1, connects[i], true ); 03162 if( MB_SUCCESS != result ) return result; 03163 } 03164 03165 Range tmph; 03166 result = mbImpl->get_adjacencies( quads, 2, 3, false, tmph ); 03167 if( MB_SUCCESS != result ) return result; 03168 assert( tmph.size() == 1 ); 03169 hexes[0] = *tmph.begin(); 03170 03171 tmph.clear(); 03172 result = mbImpl->get_adjacencies( &quads[1], 2, 3, false, tmph ); 03173 if( MB_SUCCESS != result ) return result; 03174 assert( tmph.size() == 1 ); 03175 hexes[1] = *tmph.begin(); 03176 03177 return MB_SUCCESS; 03178 } 03179 03180 ErrorCode DualTool::delete_whole_dual() 03181 { 03182 // delete dual hyperplanes 03183 Range dual_surfs, dual_curves; 03184 ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs );RR; 03185 result = mbImpl->delete_entities( dual_surfs );RR; 03186 result = this->get_dual_hyperplanes( mbImpl, 1, dual_curves );RR; 03187 result = mbImpl->delete_entities( dual_curves );RR; 03188 03189 // gather up all dual entities 03190 Range dual_ents; 03191 result = mbImpl->get_entities_by_type_and_tag( 0, MBVERTEX, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR; 03192 result = mbImpl->get_entities_by_type_and_tag( 0, MBEDGE, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR; 03193 result = mbImpl->get_entities_by_type_and_tag( 0, MBPOLYGON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR; 03194 result = 03195 mbImpl->get_entities_by_type_and_tag( 0, MBPOLYHEDRON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR; 03196 03197 // delete them, in reverse order of dimension 03198 ErrorCode tmp_result; 03199 for( Range::reverse_iterator rit = dual_ents.rbegin(); rit != dual_ents.rend(); ++rit ) 03200 { 03201 tmp_result = mbImpl->delete_entities( &( *rit ), 1 ); 03202 if( MB_SUCCESS != tmp_result ) result = tmp_result; 03203 } 03204 RR; 03205 03206 // delete dual-related tags 03207 if( 0 != dualSurfaceTag ) 03208 { 03209 tmp_result = mbImpl->tag_delete( dualSurfaceTag ); 03210 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result; 03211 } 03212 if( 0 != dualCurveTag ) 03213 { 03214 tmp_result = mbImpl->tag_delete( dualCurveTag ); 03215 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result; 03216 } 03217 if( 0 != dualEntityTag ) 03218 { 03219 tmp_result = mbImpl->tag_delete( dualEntityTag ); 03220 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result; 03221 } 03222 if( 0 != extraDualEntityTag ) 03223 { 03224 tmp_result = mbImpl->tag_delete( extraDualEntityTag ); 03225 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result; 03226 } 03227 if( 0 != dualGraphicsPointTag ) 03228 { 03229 tmp_result = mbImpl->tag_delete( dualGraphicsPointTag ); 03230 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result; 03231 } 03232 03233 return MB_SUCCESS; 03234 } 03235 03236 ErrorCode DualTool::check_dual_adjs() 03237 { 03238 // check primal-dual correspondence 03239 03240 // get the primal entities 03241 Range pents[4]; 03242 ErrorCode result = mbImpl->get_entities_by_type( 0, MBHEX, pents[3] );RR; 03243 for( int i = 2; i >= 0; i-- ) 03244 { 03245 result = mbImpl->get_adjacencies( pents[3], 2, false, pents[2], Interface::UNION );RR; 03246 } 03247 03248 // for each primal entity of dimension pd 03249 #define PRENT( ent ) CN::EntityTypeName( TYPE_FROM_HANDLE( ent ) ) << " " << ID_FROM_HANDLE( ent ) 03250 ErrorCode overall_result = MB_SUCCESS; 03251 for( int pd = 1; pd <= 3; pd++ ) 03252 { 03253 for( Range::iterator prit = pents[pd].begin(); prit != pents[pd].end(); ++prit ) 03254 { 03255 // get corresponding dual entity of dimension dd = 3-pd 03256 EntityHandle dual_ent = get_dual_entity( *prit ); 03257 if( 0 == dual_ent ) std::cerr << "Problem getting dual entity for " << PRENT( *prit ) << std::endl; 03258 03259 // for each sub dimension sd = 0..pd-1 03260 for( int sd = 0; sd < pd; sd++ ) 03261 { 03262 Range R1, R2, R3; 03263 // R1 = entities bounding primal entity of dim sd 03264 result = mbImpl->get_adjacencies( &( *prit ), 1, sd, false, R1 );RR; 03265 03266 // R2 = entities bounded by dual entity, of dim 3-sd 03267 result = mbImpl->get_adjacencies( &dual_ent, 1, 3 - sd, false, R2 );RR; 03268 03269 if( R1.size() != R2.size() ) 03270 { 03271 std::cerr << PRENT( *prit ) << ": number of adj ents in " 03272 << "primal/dual don't agree for dimension " << sd << "." << std::endl; 03273 overall_result = MB_FAILURE; 03274 } 03275 03276 // for each entity in R1, get its dual and look for it in R2 03277 for( Range::iterator r1it = R1.begin(); r1it != R1.end(); ++r1it ) 03278 { 03279 EntityHandle tmp_dual = get_dual_entity( *r1it ); 03280 if( R2.find( tmp_dual ) == R2.end() ) 03281 { 03282 std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r1it ) << " isn't adjacent in dual." 03283 << std::endl; 03284 overall_result = MB_FAILURE; 03285 } 03286 } 03287 // ditto for R2 03288 for( Range::iterator r2it = R2.begin(); r2it != R2.end(); ++r2it ) 03289 { 03290 EntityHandle tmp_prim = get_dual_entity( *r2it ); 03291 if( R1.find( tmp_prim ) == R1.end() ) 03292 { 03293 std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r2it ) << " isn't adjacent in primal." 03294 << std::endl; 03295 overall_result = MB_FAILURE; 03296 } 03297 } 03298 } 03299 } 03300 } 03301 03302 return overall_result; 03303 } 03304 03305 } // namespace moab