MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 #include <iostream> 00002 #include <map> 00003 00004 #include "moab/FBEngine.hpp" 00005 #include "moab/Interface.hpp" 00006 #include "moab/GeomTopoTool.hpp" 00007 #include "moab/GeomUtil.hpp" 00008 #include "moab/OrientedBoxTreeTool.hpp" 00009 00010 #include <cstdlib> 00011 #include <cstring> 00012 #include <map> 00013 #include <set> 00014 #include <queue> 00015 #include <algorithm> 00016 #include <cassert> 00017 00018 #include "SmoothCurve.hpp" 00019 #include "SmoothFace.hpp" 00020 00021 // this is just to replace MBI with moab interface, which is _mbImpl in this class 00022 #define MBI _mbImpl 00023 #define MBERRORR( rval, STR ) \ 00024 { \ 00025 if( MB_SUCCESS != ( rval ) ) \ 00026 { \ 00027 std::cout << ( STR ) << std::endl; \ 00028 return rval; \ 00029 } \ 00030 } 00031 00032 namespace moab 00033 { 00034 00035 // some tolerances for ray tracing and geometry intersections 00036 // these are involved in ray tracing, at least 00037 00038 double tolerance = 0.01; // TODO: how is this used ???? 00039 double tolerance_segment = 0.01; // for segments intersection, points collapse, area coordinates for triangles 00040 // it should be a relative, not an absolute value 00041 // as this splitting operation can create small edges, it should be relatively high 00042 // or, it should be coordinated with the decimation errors 00043 // we really just want to preserve the integrity of the mesh, we should avoid creating small edges 00044 // or angles 00045 const bool Debug_surf_eval = false; 00046 bool debug_splits = false; 00047 00048 // will compute intersection between a segment and slice of a plane 00049 // output is the intersection point 00050 bool intersect_segment_and_plane_slice( CartVect& from, CartVect& to, CartVect& p1, CartVect& p2, CartVect&, 00051 CartVect& normPlane, CartVect& intx_point, double& parPos ) 00052 { 00053 // 00054 // plane eq is normPlane % r + d = 0, or normPlane % r - normPlane%p1 = 0 00055 double dd = -normPlane % p1; 00056 double valFrom = normPlane % from + dd; 00057 double valTo = normPlane % to + dd; 00058 00059 if( fabs( valFrom ) < tolerance_segment ) 00060 { 00061 intx_point = from; 00062 parPos = 0.; 00063 double proj1 = ( intx_point - p1 ) % ( p2 - p1 ); 00064 double proj2 = ( intx_point - p2 ) % ( p1 - p2 ); 00065 if( proj1 <= -tolerance_segment || proj2 <= -tolerance_segment ) return false; 00066 if( debug_splits ) std::cout << "intx : " << intx_point << "\n"; 00067 return true; 00068 } 00069 if( fabs( valTo ) < tolerance_segment ) 00070 { 00071 intx_point = to; 00072 parPos = 1; 00073 double proj1 = ( intx_point - p1 ) % ( p2 - p1 ); 00074 double proj2 = ( intx_point - p2 ) % ( p1 - p2 ); 00075 if( proj1 <= -tolerance_segment || proj2 <= -tolerance_segment ) return false; 00076 if( debug_splits ) std::cout << "intx : " << intx_point << "\n"; 00077 return true; 00078 } 00079 if( valFrom * valTo > 0 ) return false; // no intersection, although it could be very close 00080 // else, it could intersect the plane; check for the slice too. 00081 parPos = valFrom / ( valFrom - valTo ); // this is 0 for valFrom 0, 1 for valTo 0 00082 intx_point = from + ( to - from ) * parPos; 00083 // now check if the intx_point is indeed between p1 and p2 in the slice. 00084 double proj1 = ( intx_point - p1 ) % ( p2 - p1 ); 00085 double proj2 = ( intx_point - p2 ) % ( p1 - p2 ); 00086 if( proj1 <= -tolerance_segment || proj2 <= -tolerance_segment ) return false; 00087 00088 if( debug_splits ) std::cout << "intx : " << intx_point << "\n"; 00089 return true; 00090 } 00091 00092 ErrorCode area_coordinates( Interface* mbi, EntityHandle tri, CartVect& pnt, double* area_coord, 00093 EntityHandle& boundary_handle ) 00094 { 00095 00096 int nnodes; 00097 const EntityHandle* conn3; 00098 ErrorCode rval = mbi->get_connectivity( tri, conn3, nnodes ); 00099 MBERRORR( rval, "Failed to get connectivity" ); 00100 assert( 3 == nnodes ); 00101 CartVect P[3]; 00102 rval = mbi->get_coords( conn3, nnodes, (double*)&P[0] ); 00103 MBERRORR( rval, "Failed to get coordinates" ); 00104 00105 CartVect r0( P[0] - pnt ); 00106 CartVect r1( P[1] - pnt ); 00107 CartVect r2( P[2] - pnt ); 00108 if( debug_splits ) 00109 { 00110 std::cout << " nodes:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n"; 00111 std::cout << " distances: " << r0.length() << " " << r1.length() << " " << r2.length() << "\n"; 00112 } 00113 if( r0.length() < tolerance_segment ) 00114 { 00115 area_coord[0] = 1.; 00116 area_coord[1] = 0.; 00117 area_coord[2] = 0.; 00118 boundary_handle = conn3[0]; 00119 return MB_SUCCESS; 00120 } 00121 if( r1.length() < tolerance_segment ) 00122 { 00123 area_coord[0] = 0.; 00124 area_coord[1] = 1.; 00125 area_coord[2] = 0.; 00126 boundary_handle = conn3[1]; 00127 return MB_SUCCESS; 00128 } 00129 if( r2.length() < tolerance_segment ) 00130 { 00131 area_coord[0] = 0.; 00132 area_coord[1] = 0.; 00133 area_coord[2] = 1.; 00134 boundary_handle = conn3[2]; 00135 return MB_SUCCESS; 00136 } 00137 00138 CartVect v1( P[1] - P[0] ); 00139 CartVect v2( P[2] - P[0] ); 00140 00141 double areaDouble = ( v1 * v2 ).length(); // the same for CartVect 00142 if( areaDouble < tolerance_segment * tolerance_segment ) { MBERRORR( MB_FAILURE, "area of triangle too small" ); } 00143 area_coord[0] = ( r1 * r2 ).length() / areaDouble; 00144 area_coord[1] = ( r2 * r0 ).length() / areaDouble; 00145 area_coord[2] = ( r0 * r1 ).length() / areaDouble; 00146 00147 if( fabs( area_coord[0] + area_coord[1] + area_coord[2] - 1 ) > tolerance_segment ) 00148 { 00149 MBERRORR( MB_FAILURE, "point outside triangle" ); 00150 } 00151 // the tolerance is used here for area coordinates (0 to 1), and in other 00152 // parts it is used as an absolute distance; pretty inconsistent. 00153 bool side0 = ( area_coord[0] < tolerance_segment ); 00154 bool side1 = ( area_coord[1] < tolerance_segment ); 00155 bool side2 = ( area_coord[2] < tolerance_segment ); 00156 if( !side0 && !side1 && !side2 ) return MB_SUCCESS; // interior point 00157 // now, find out what boundary is in question 00158 // first, get all edges, in order 00159 std::vector< EntityHandle > edges; 00160 EntityHandle nn2[2]; 00161 for( int i = 0; i < 3; i++ ) 00162 { 00163 nn2[0] = conn3[( i + 1 ) % 3]; 00164 nn2[1] = conn3[( i + 2 ) % 3]; 00165 std::vector< EntityHandle > adjacent; 00166 rval = mbi->get_adjacencies( nn2, 2, 1, false, adjacent, Interface::INTERSECT ); 00167 MBERRORR( rval, "Failed to get edges" ); 00168 if( adjacent.size() != 1 ) MBERRORR( MB_FAILURE, "Failed to get adjacent edges" ); 00169 // should be only one edge here 00170 edges.push_back( adjacent[0] ); 00171 } 00172 00173 if( side0 ) boundary_handle = edges[0]; 00174 if( side1 ) boundary_handle = edges[1]; 00175 if( side2 ) boundary_handle = edges[2]; 00176 00177 return MB_SUCCESS; 00178 } 00179 00180 FBEngine::FBEngine( Interface* impl, GeomTopoTool* topoTool, const bool smooth ) 00181 : _mbImpl( impl ), _my_geomTopoTool( topoTool ), _t_created( false ), _smooth( smooth ), _initialized( false ), 00182 _smthFace( NULL ), _smthCurve( NULL ) 00183 { 00184 if( !_my_geomTopoTool ) 00185 { 00186 _my_geomTopoTool = new GeomTopoTool( _mbImpl ); 00187 _t_created = true; 00188 } 00189 // should this be part of the constructor or not? 00190 // Init(); 00191 } 00192 FBEngine::~FBEngine() 00193 { 00194 clean(); 00195 _smooth = false; 00196 } 00197 00198 void FBEngine::clean() 00199 { 00200 if( _smooth ) 00201 { 00202 _faces.clear(); 00203 _edges.clear(); 00204 int size1 = _my_gsets[1].size(); 00205 int i = 0; 00206 for( i = 0; i < size1; i++ ) 00207 delete _smthCurve[i]; 00208 delete[] _smthCurve; 00209 _smthCurve = NULL; 00210 size1 = _my_gsets[2].size(); 00211 for( i = 0; i < size1; i++ ) 00212 delete _smthFace[i]; 00213 delete[] _smthFace; 00214 _smthFace = NULL; 00215 //_smooth = false; 00216 } 00217 00218 for( int j = 0; j < 5; j++ ) 00219 _my_gsets[j].clear(); 00220 if( _t_created ) delete _my_geomTopoTool; 00221 _my_geomTopoTool = NULL; 00222 _t_created = false; 00223 } 00224 00225 ErrorCode FBEngine::Init() 00226 { 00227 if( !_initialized ) 00228 { 00229 if( !_my_geomTopoTool ) return MB_FAILURE; 00230 00231 ErrorCode rval = _my_geomTopoTool->find_geomsets( _my_gsets ); 00232 assert( rval == MB_SUCCESS ); 00233 if( MB_SUCCESS != rval ) { return rval; } 00234 00235 rval = split_quads(); 00236 assert( rval == MB_SUCCESS ); 00237 00238 rval = _my_geomTopoTool->construct_obb_trees(); 00239 assert( rval == MB_SUCCESS ); 00240 00241 if( _smooth ) rval = initializeSmoothing(); 00242 assert( rval == MB_SUCCESS ); 00243 00244 _initialized = true; 00245 } 00246 return MB_SUCCESS; 00247 } 00248 ErrorCode FBEngine::initializeSmoothing() 00249 { 00250 // 00251 /*ErrorCode rval = Init(); 00252 MBERRORR(rval, "failed initialize");*/ 00253 // first of all, we need to retrieve all the surfaces from the (root) set 00254 // in icesheet_test we use iGeom, but maybe that is a stretch 00255 // get directly the sets with geom dim 2, and from there create the SmoothFace 00256 Tag geom_tag; 00257 ErrorCode rval = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag ); 00258 MBERRORR( rval, "can't get geom tag" ); 00259 00260 int numSurfaces = _my_gsets[2].size(); 00261 // SmoothFace ** smthFace = new SmoothFace *[numSurfaces]; 00262 _smthFace = new SmoothFace*[numSurfaces]; 00263 00264 // there should also be a map from surfaces to evaluators 00265 // std::map<MBEntityHandle, SmoothFace*> mapSurfaces; 00266 00267 int i = 0; 00268 Range::iterator it; 00269 for( it = _my_gsets[2].begin(); it != _my_gsets[2].end(); ++it, i++ ) 00270 { 00271 EntityHandle face = *it; 00272 _smthFace[i] = new SmoothFace( MBI, face, _my_geomTopoTool ); // geom topo tool will be used for searching, 00273 // among other things; also for senses in edge sets... 00274 _faces[face] = _smthFace[i]; 00275 } 00276 00277 int numCurves = _my_gsets[1].size(); // csets.size(); 00278 // SmoothCurve ** smthCurve = new SmoothCurve *[numCurves]; 00279 _smthCurve = new SmoothCurve*[numCurves]; 00280 // there should also be a map from surfaces to evaluators 00281 // std::map<MBEntityHandle, SmoothCurve*> mapCurves; 00282 00283 i = 0; 00284 for( it = _my_gsets[1].begin(); it != _my_gsets[1].end(); ++it, i++ ) 00285 { 00286 EntityHandle curve = *it; 00287 _smthCurve[i] = new SmoothCurve( MBI, curve, _my_geomTopoTool ); 00288 _edges[curve] = _smthCurve[i]; 00289 } 00290 00291 for( i = 0; i < numSurfaces; i++ ) 00292 { 00293 _smthFace[i]->init_gradient(); // this will also retrieve the triangles in each surface 00294 _smthFace[i]->compute_tangents_for_each_edge(); // this one will consider all edges 00295 // internal, so the 00296 // tangents are all in the direction of the edge; a little bit of waste, as we store 00297 // one tangent for each edge node , even though they are equal here... 00298 // no loops are considered 00299 } 00300 00301 // this will be used to mark boundary edges, so for them the control points are computed earlier 00302 unsigned char value = 0; // default value is "not used"=0 for the tag 00303 // unsigned char def_data_bit = 1;// valid by default 00304 // rval = mb->tag_create("valid", 1, MB_TAG_BIT, validTag, &def_data_bit); 00305 Tag markTag; 00306 rval = MBI->tag_get_handle( "MARKER", 1, MB_TYPE_BIT, markTag, MB_TAG_EXCL | MB_TAG_BIT, 00307 &value ); // default value : 0 = not computed yet 00308 // each feature edge will need to have a way to retrieve at every moment the surfaces it belongs 00309 // to from edge sets, using the sense tag, we can get faces, and from each face, using the map, 00310 // we can get the SmoothFace (surface evaluator), that has everything, including the normals!!! 00311 assert( rval == MB_SUCCESS ); 00312 00313 // create the tag also for control points on the edges 00314 double defCtrlPoints[9] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. }; 00315 Tag edgeCtrlTag; 00316 rval = MBI->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, edgeCtrlTag, MB_TAG_DENSE | MB_TAG_CREAT, 00317 &defCtrlPoints ); 00318 assert( rval == MB_SUCCESS ); 00319 00320 Tag facetCtrlTag; 00321 double defControls[18] = { 0. }; 00322 rval = MBI->tag_get_handle( "CONTROLFACE", 18, MB_TYPE_DOUBLE, facetCtrlTag, MB_TAG_CREAT | MB_TAG_DENSE, 00323 &defControls ); 00324 assert( rval == MB_SUCCESS ); 00325 00326 Tag facetEdgeCtrlTag; 00327 double defControls2[27] = { 0. }; // corresponding to 9 control points on edges, in order 00328 // from edge 0, 1, 2 ( 1-2, 2-0, 0-1 ) 00329 rval = MBI->tag_get_handle( "CONTROLEDGEFACE", 27, MB_TYPE_DOUBLE, facetEdgeCtrlTag, MB_TAG_CREAT | MB_TAG_DENSE, 00330 &defControls2 ); 00331 assert( rval == MB_SUCCESS ); 00332 // if the 00333 double min_dot = -1.0; // depends on _angle, but now we ignore it, for the time being 00334 for( i = 0; i < numCurves; i++ ) 00335 { 00336 _smthCurve[i]->compute_tangents_for_each_edge(); // do we need surfaces now? or just the chains? 00337 // the computed edges will be marked accordingly; later one, only internal edges to surfaces 00338 // are left 00339 _smthCurve[i]->compute_control_points_on_boundary_edges( min_dot, _faces, edgeCtrlTag, markTag ); 00340 } 00341 00342 // when done with boundary edges, compute the control points on all edges in the surfaces 00343 00344 for( i = 0; i < numSurfaces; i++ ) 00345 { 00346 // also pass the tags for 00347 _smthFace[i]->compute_control_points_on_edges( min_dot, edgeCtrlTag, markTag ); 00348 } 00349 00350 // now we should be able to compute the control points for the facets 00351 00352 for( i = 0; i < numSurfaces; i++ ) 00353 { 00354 // also pass the tags for edge and facet control points 00355 _smthFace[i]->compute_internal_control_points_on_facets( min_dot, facetCtrlTag, facetEdgeCtrlTag ); 00356 } 00357 // we will need to compute the tangents for all edges in the model 00358 // they will be needed for control points for each edge 00359 // the boundary edges and the feature edges are more complicated 00360 // the boundary edges need to consider their loops, but feature edges need to consider loops and 00361 // the normals on each connected surface 00362 00363 // some control points 00364 if( Debug_surf_eval ) 00365 for( i = 0; i < numSurfaces; i++ ) 00366 _smthFace[i]->DumpModelControlPoints(); 00367 00368 return MB_SUCCESS; 00369 } 00370 00371 // clean up the smooth tags data if created, so the files will be smaller 00372 // if saved 00373 // also, recompute the tags if topology is modified 00374 void FBEngine::delete_smooth_tags() 00375 { 00376 // get all tags from database that are created for smooth data, and 00377 // delete them; it will delete all data associated with them 00378 // first tags from faces, edges: 00379 std::vector< Tag > smoothTags; 00380 int size1 = (int)_my_gsets[2].size(); 00381 00382 for( int i = 0; i < size1; i++ ) 00383 { 00384 // these 2 will append gradient tag and plane tag 00385 _smthFace[i]->append_smooth_tags( smoothTags ); 00386 } 00387 // then , get other tags: 00388 // "TANGENTS", "MARKER", "CONTROLEDGE", "CONTROLFACE", "CONTROLEDGEFACE" 00389 Tag tag_handle; 00390 ErrorCode rval = _mbImpl->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tag_handle ); 00391 if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle ); 00392 00393 rval = _mbImpl->tag_get_handle( "MARKER", 1, MB_TYPE_BIT, tag_handle ); 00394 if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle ); 00395 00396 rval = _mbImpl->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, tag_handle ); 00397 if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle ); 00398 00399 rval = _mbImpl->tag_get_handle( "CONTROLFACE", 18, MB_TYPE_DOUBLE, tag_handle ); 00400 if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle ); 00401 00402 rval = _mbImpl->tag_get_handle( "CONTROLEDGEFACE", 27, MB_TYPE_DOUBLE, tag_handle ); 00403 if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle ); 00404 00405 // a lot of tags, delete them 00406 for( unsigned int k = 0; k < smoothTags.size(); k++ ) 00407 { 00408 // could be a lot of data 00409 _mbImpl->tag_delete( smoothTags[k] ); 00410 } 00411 } 00412 /* 00413 #define COPY_RANGE(r, vec) { \ 00414 EntityHandle *tmp_ptr = reinterpret_cast<EntityHandle*>(vec); \ 00415 std::copy(r.begin(), r.end(), tmp_ptr);} 00416 */ 00417 00418 /*static inline void 00419 ProcessError(const char* desc);*/ 00420 00421 ErrorCode FBEngine::getRootSet( EntityHandle* root_set ) 00422 { 00423 *root_set = _my_geomTopoTool->get_root_model_set(); 00424 return MB_SUCCESS; 00425 } 00426 00427 ErrorCode FBEngine::getNumEntSets( EntityHandle set, int num_hops, int* all_sets ) 00428 { 00429 ErrorCode rval = MBI->num_contained_meshsets( set, all_sets, num_hops + 1 ); 00430 return rval; 00431 } 00432 00433 ErrorCode FBEngine::createEntSet( int isList, EntityHandle* pSet ) 00434 { 00435 ErrorCode rval; 00436 00437 if( isList ) 00438 rval = MBI->create_meshset( MESHSET_ORDERED, *pSet ); 00439 else 00440 rval = MBI->create_meshset( MESHSET_SET, *pSet ); 00441 00442 return rval; 00443 } 00444 00445 ErrorCode FBEngine::getEntities( EntityHandle set_handle, int entity_type, Range& gentities ) 00446 { 00447 int i; 00448 if( 0 > entity_type || 4 < entity_type ) { return MB_FAILURE; } 00449 else if( entity_type < 4 ) 00450 { // 4 means all entities 00451 gentities = _my_geomTopoTool->geoRanges()[entity_type]; // all from root set! 00452 } 00453 else 00454 { 00455 gentities.clear(); 00456 for( i = 0; i < 4; i++ ) 00457 { 00458 gentities.merge( _my_geomTopoTool->geoRanges()[i] ); 00459 } 00460 } 00461 Range sets; 00462 // see now if they are in the set passed as input or not 00463 ErrorCode rval = MBI->get_entities_by_type( set_handle, MBENTITYSET, sets ); 00464 MBERRORR( rval, "can't get sets in the initial set" ); 00465 gentities = intersect( gentities, sets ); 00466 00467 return MB_SUCCESS; 00468 } 00469 00470 ErrorCode FBEngine::addEntArrToSet( const Range& entities, EntityHandle set ) 00471 { 00472 return MBI->add_entities( set, entities ); 00473 } 00474 00475 ErrorCode FBEngine::addEntSet( EntityHandle entity_set_to_add, EntityHandle entity_set_handle ) 00476 { 00477 return MBI->add_entities( entity_set_handle, &entity_set_to_add, 1 ); 00478 } 00479 00480 ErrorCode FBEngine::getNumOfType( EntityHandle set, int ent_type, int* pNum ) 00481 { 00482 if( 0 > ent_type || 4 < ent_type ) 00483 { 00484 std::cout << "Invalid type\n"; 00485 return MB_FAILURE; 00486 } 00487 // get sets of geom dimension tag from here, and intersect with the gentities from geo 00488 // ranges 00489 00490 // get the geom dimensions sets in the set (AKA gentities) 00491 Range geom_sets; 00492 Tag geom_tag; 00493 ErrorCode rval = 00494 _mbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag, MB_TAG_SPARSE | MB_TAG_CREAT ); 00495 MBERRORR( rval, "Failed to get geom tag." ); 00496 rval = _mbImpl->get_entities_by_type_and_tag( set, MBENTITYSET, &geom_tag, NULL, 1, geom_sets, Interface::UNION ); 00497 MBERRORR( rval, "Failed to get gentities from set" ); 00498 00499 if( ent_type == 4 ) 00500 { 00501 *pNum = 0; 00502 for( int k = 0; k <= 3; k++ ) 00503 { 00504 Range gEntsOfTypeK = intersect( geom_sets, _my_geomTopoTool->geoRanges()[k] ); 00505 *pNum += (int)gEntsOfTypeK.size(); 00506 } 00507 } 00508 else 00509 { 00510 Range gEntsOfType = intersect( geom_sets, _my_geomTopoTool->geoRanges()[ent_type] ); 00511 *pNum = (int)gEntsOfType.size(); 00512 } 00513 // we do not really check if it is in the set or not; 00514 // _my_gsets[i].find(gent) != _my_gsets[i].end() 00515 return MB_SUCCESS; 00516 } 00517 00518 ErrorCode FBEngine::getEntType( EntityHandle gent, int* type ) 00519 { 00520 for( int i = 0; i < 4; i++ ) 00521 { 00522 if( _my_geomTopoTool->geoRanges()[i].find( gent ) != _my_geomTopoTool->geoRanges()[i].end() ) 00523 { 00524 *type = i; 00525 return MB_SUCCESS; 00526 } 00527 } 00528 *type = -1; // failure 00529 return MB_FAILURE; 00530 } 00531 ErrorCode FBEngine::getEntBoundBox( EntityHandle gent, double* min_x, double* min_y, double* min_z, double* max_x, 00532 double* max_y, double* max_z ) 00533 { 00534 ErrorCode rval; 00535 int type; 00536 rval = getEntType( gent, &type ); 00537 MBERRORR( rval, "Failed to get entity type." ); 00538 00539 if( type == 0 ) 00540 { 00541 rval = getVtxCoord( gent, min_x, min_y, min_z ); 00542 MBERRORR( rval, "Failed to get vertex coordinates." ); 00543 max_x = min_x; 00544 max_y = min_y; 00545 max_z = min_z; 00546 } 00547 else if( type == 1 ) 00548 { 00549 MBERRORR( MB_FAILURE, "iGeom_getEntBoundBox is not supported for Edge entity type." ); 00550 } 00551 else if( type == 2 || type == 3 ) 00552 { 00553 00554 EntityHandle root; 00555 CartVect center, axis[3]; 00556 rval = _my_geomTopoTool->get_root( gent, root ); 00557 MBERRORR( rval, "Failed to get tree root in iGeom_getEntBoundBox." ); 00558 rval = _my_geomTopoTool->obb_tree()->box( root, center.array(), axis[0].array(), axis[1].array(), 00559 axis[2].array() ); 00560 MBERRORR( rval, "Failed to get closest point in iGeom_getEntBoundBox." ); 00561 00562 CartVect absv[3]; 00563 for( int i = 0; i < 3; i++ ) 00564 { 00565 absv[i] = CartVect( fabs( axis[i][0] ), fabs( axis[i][1] ), fabs( axis[i][2] ) ); 00566 } 00567 CartVect min, max; 00568 min = center - absv[0] - absv[1] - absv[2]; 00569 max = center + absv[0] + absv[1] + absv[2]; 00570 *min_x = min[0]; 00571 *min_y = min[1]; 00572 *min_z = min[2]; 00573 *max_x = max[0]; 00574 *max_y = max[1]; 00575 *max_z = max[2]; 00576 } 00577 else 00578 return MB_FAILURE; 00579 00580 return MB_SUCCESS; 00581 } 00582 ErrorCode FBEngine::getEntClosestPt( EntityHandle this_gent, double near_x, double near_y, double near_z, double* on_x, 00583 double* on_y, double* on_z ) 00584 { 00585 ErrorCode rval; 00586 int type; 00587 rval = getEntType( this_gent, &type ); 00588 MBERRORR( rval, "Failed to get entity type." ); 00589 00590 if( type == 0 ) 00591 { 00592 rval = getVtxCoord( this_gent, on_x, on_y, on_z ); 00593 MBERRORR( rval, "Failed to get vertex coordinates." ); 00594 } 00595 else if( _smooth && type == 1 ) 00596 { 00597 *on_x = near_x; 00598 *on_y = near_y; 00599 *on_z = near_z; 00600 SmoothCurve* smthcurve = _edges[this_gent]; 00601 // call the new method from smooth edge 00602 smthcurve->move_to_curve( *on_x, *on_y, *on_z ); 00603 } 00604 else if( type == 2 || type == 3 ) 00605 { 00606 double point[3] = { near_x, near_y, near_z }; 00607 double point_out[3]; 00608 EntityHandle root, facet_out; 00609 if( _smooth && 2 == type ) 00610 { 00611 SmoothFace* smthFace = _faces[this_gent]; 00612 *on_x = near_x; 00613 *on_y = near_y; 00614 *on_z = near_z; 00615 smthFace->move_to_surface( *on_x, *on_y, *on_z ); 00616 } 00617 else 00618 { 00619 rval = _my_geomTopoTool->get_root( this_gent, root ); 00620 MBERRORR( rval, "Failed to get tree root in iGeom_getEntClosestPt." ); 00621 rval = _my_geomTopoTool->obb_tree()->closest_to_location( point, root, point_out, facet_out ); 00622 MBERRORR( rval, "Failed to get closest point in iGeom_getEntClosestPt." ); 00623 00624 *on_x = point_out[0]; 00625 *on_y = point_out[1]; 00626 *on_z = point_out[2]; 00627 } 00628 } 00629 else 00630 return MB_TYPE_OUT_OF_RANGE; 00631 00632 return MB_SUCCESS; 00633 } 00634 00635 ErrorCode FBEngine::getVtxCoord( EntityHandle vertex_handle, double* x0, double* y0, double* z0 ) 00636 { 00637 int type; 00638 ErrorCode rval = getEntType( vertex_handle, &type ); 00639 MBERRORR( rval, "Failed to get entity type in getVtxCoord." ); 00640 00641 if( type != 0 ) { MBERRORR( MB_FAILURE, "Entity is not a vertex type." ); } 00642 00643 Range entities; 00644 rval = MBI->get_entities_by_type( vertex_handle, MBVERTEX, entities ); 00645 MBERRORR( rval, "can't get nodes in vertex set." ); 00646 00647 if( entities.size() != 1 ) { MBERRORR( MB_FAILURE, "Vertex has multiple points." ); } 00648 double coords[3]; 00649 EntityHandle node = entities[0]; 00650 rval = MBI->get_coords( &node, 1, coords ); 00651 MBERRORR( rval, "can't get coordinates." ); 00652 *x0 = coords[0]; 00653 *y0 = coords[1]; 00654 *z0 = coords[2]; 00655 00656 return MB_SUCCESS; 00657 } 00658 00659 ErrorCode FBEngine::gsubtract( EntityHandle entity_set_1, EntityHandle entity_set_2, EntityHandle result_entity_set ) 00660 { 00661 /*result_entity_set = subtract(entity_set_1, entity_set_2);*/ 00662 Range ents1, ents2; 00663 ErrorCode rval = MBI->get_entities_by_type( entity_set_1, MBENTITYSET, ents1 ); 00664 MBERRORR( rval, "can't get entities from set 1." ); 00665 00666 rval = MBI->get_entities_by_type( entity_set_2, MBENTITYSET, ents2 ); 00667 MBERRORR( rval, "can't get entities from set 2." ); 00668 00669 ents1 = subtract( ents1, ents2 ); 00670 rval = MBI->clear_meshset( &result_entity_set, 1 ); 00671 MBERRORR( rval, "can't empty set." ); 00672 00673 rval = MBI->add_entities( result_entity_set, ents1 ); 00674 MBERRORR( rval, "can't add result to set." ); 00675 00676 return rval; 00677 } 00678 00679 ErrorCode FBEngine::getEntNrmlXYZ( EntityHandle entity_handle, double x, double y, double z, double* nrml_i, 00680 double* nrml_j, double* nrml_k ) 00681 { 00682 // just do for surface and volume 00683 int type; 00684 ErrorCode rval = getEntType( entity_handle, &type ); 00685 MBERRORR( rval, "Failed to get entity type in iGeom_getEntNrmlXYZ." ); 00686 00687 if( type != 2 && type != 3 ) 00688 { 00689 MBERRORR( MB_FAILURE, "Entities passed into gentityNormal must be face or volume." ); 00690 } 00691 00692 if( _smooth && 2 == type ) 00693 { 00694 SmoothFace* smthFace = _faces[entity_handle]; 00695 //*on_x = near_x; *on_y = near_y; *on_z = near_z; 00696 smthFace->normal_at( x, y, z, *nrml_i, *nrml_j, *nrml_k ); 00697 } 00698 else 00699 { 00700 // get closest location and facet 00701 double point[3] = { x, y, z }; 00702 double point_out[3]; 00703 EntityHandle root, facet_out; 00704 _my_geomTopoTool->get_root( entity_handle, root ); 00705 rval = _my_geomTopoTool->obb_tree()->closest_to_location( point, root, point_out, facet_out ); 00706 MBERRORR( rval, "Failed to get closest location in iGeom_getEntNrmlXYZ." ); 00707 00708 // get facet normal 00709 const EntityHandle* conn; 00710 int len; 00711 CartVect coords[3], normal; 00712 rval = MBI->get_connectivity( facet_out, conn, len ); 00713 MBERRORR( rval, "Failed to get triangle connectivity in iGeom_getEntNrmlXYZ." ); 00714 if( len != 3 ) MBERRORR( MB_FAILURE, " not a triangle, error " ); 00715 00716 rval = MBI->get_coords( conn, len, coords[0].array() ); 00717 MBERRORR( rval, "Failed to get triangle coordinates in iGeom_getEntNrmlXYZ." ); 00718 00719 coords[1] -= coords[0]; 00720 coords[2] -= coords[0]; 00721 normal = coords[1] * coords[2]; 00722 normal.normalize(); 00723 *nrml_i = normal[0]; 00724 *nrml_j = normal[1]; 00725 *nrml_k = normal[2]; 00726 } 00727 return MB_SUCCESS; 00728 } 00729 00730 ErrorCode FBEngine::getPntRayIntsct( double x, double y, double z, double dir_x, double dir_y, double dir_z, 00731 std::vector< EntityHandle >& intersect_entity_handles, 00732 /* int storage_order,*/ 00733 std::vector< double >& intersect_coords, std::vector< double >& param_coords ) 00734 { 00735 // this is pretty cool 00736 // we will return only surfaces (gentities ) 00737 // 00738 ErrorCode rval; 00739 00740 unsigned int numfaces = _my_gsets[2].size(); 00741 // do ray fire 00742 const double point[] = { x, y, z }; 00743 const double dir[] = { dir_x, dir_y, dir_z }; 00744 CartVect P( point ); 00745 CartVect V( dir ); 00746 00747 // std::vector<double> distances; 00748 std::vector< EntityHandle > facets; 00749 // std::vector<EntityHandle> sets; 00750 unsigned int i; 00751 for( i = 0; i < numfaces; i++ ) 00752 { 00753 EntityHandle face = _my_gsets[2][i]; 00754 EntityHandle rootForFace; 00755 rval = _my_geomTopoTool->get_root( face, rootForFace ); 00756 MBERRORR( rval, "Failed to get root of face." ); 00757 std::vector< double > distances_out; 00758 std::vector< EntityHandle > sets_out; 00759 std::vector< EntityHandle > facets_out; 00760 rval = _my_geomTopoTool->obb_tree()->ray_intersect_sets( distances_out, sets_out, facets_out, rootForFace, 00761 tolerance, point, dir ); 00762 unsigned int j; 00763 for( j = 0; j < distances_out.size(); j++ ) 00764 param_coords.push_back( distances_out[j] ); 00765 for( j = 0; j < sets_out.size(); j++ ) 00766 intersect_entity_handles.push_back( sets_out[j] ); 00767 for( j = 0; j < facets_out.size(); j++ ) 00768 facets.push_back( facets_out[j] ); 00769 00770 MBERRORR( rval, "Failed to get ray intersections." ); 00771 } 00772 // facets.size == distances.size()!! 00773 for( i = 0; i < param_coords.size(); i++ ) 00774 { 00775 CartVect intx = P + param_coords[i] * V; 00776 for( int j = 0; j < 3; j++ ) 00777 intersect_coords.push_back( intx[j] ); 00778 } 00779 if( _smooth ) 00780 { 00781 // correct the intersection point and the distance for smooth surfaces 00782 for( i = 0; i < intersect_entity_handles.size(); i++ ) 00783 { 00784 // EntityHandle geoSet = MBH_cast(sets[i]); 00785 SmoothFace* sFace = _faces[intersect_entity_handles[i]]; 00786 // correct coordinates and distance from point 00787 /*moab::ErrorCode ray_intersection_correct(moab::EntityHandle facet, // (IN) the facet 00788 where the patch is defined moab::CartVect &pt, // (IN) shoot from moab::CartVect &ray, 00789 // (IN) ray direction moab::CartVect &eval_pt, // (INOUT) The intersection point double 00790 & distance, // (IN OUT) the new distance bool &outside);*/ 00791 CartVect pos( &( intersect_coords[3 * i] ) ); 00792 double dist = param_coords[i]; 00793 bool outside = false; 00794 rval = sFace->ray_intersection_correct( facets[i], P, V, pos, dist, outside ); 00795 MBERRORR( rval, "Failed to get better point on ray." ); 00796 param_coords[i] = dist; 00797 00798 for( int j = 0; j < 3; j++ ) 00799 intersect_coords[3 * i + j] = pos[j]; 00800 } 00801 } 00802 return MB_SUCCESS; 00803 } 00804 00805 ErrorCode FBEngine::getAdjacentEntities( const EntityHandle from, const int to_dim, Range& adjs ) 00806 { 00807 int this_dim = -1; 00808 for( int i = 0; i < 4; i++ ) 00809 { 00810 if( _my_geomTopoTool->geoRanges()[i].find( from ) != _my_geomTopoTool->geoRanges()[i].end() ) 00811 { 00812 this_dim = i; 00813 break; 00814 } 00815 } 00816 00817 // check target dimension 00818 if( -1 == this_dim ) 00819 { 00820 // ProcessError(iBase_FAILURE, "Entity not a geometry entity."); 00821 return MB_FAILURE; 00822 } 00823 else if( 0 > to_dim || 3 < to_dim ) 00824 { 00825 // ProcessError(iBase_FAILURE, "To dimension must be between 0 and 3."); 00826 return MB_FAILURE; 00827 } 00828 else if( to_dim == this_dim ) 00829 { 00830 // ProcessError(iBase_FAILURE, 00831 // "To dimension must be different from entity dimension."); 00832 return MB_FAILURE; 00833 } 00834 00835 ErrorCode rval = MB_SUCCESS; 00836 adjs.clear(); 00837 if( to_dim > this_dim ) 00838 { 00839 int diffDim = to_dim - this_dim; 00840 rval = MBI->get_parent_meshsets( from, adjs, diffDim ); 00841 if( MB_SUCCESS != rval ) return rval; 00842 if( diffDim > 1 ) 00843 { 00844 // subtract the parents that come with diffDim-1 hops 00845 Range extra; 00846 rval = MBI->get_parent_meshsets( from, extra, diffDim - 1 ); 00847 if( MB_SUCCESS != rval ) return rval; 00848 adjs = subtract( adjs, extra ); 00849 } 00850 } 00851 else 00852 { 00853 int diffDim = this_dim - to_dim; 00854 rval = MBI->get_child_meshsets( from, adjs, diffDim ); 00855 if( MB_SUCCESS != rval ) return rval; 00856 if( diffDim > 1 ) 00857 { 00858 // subtract the children that come with diffDim-1 hops 00859 Range extra; 00860 rval = MBI->get_child_meshsets( from, extra, diffDim - 1 ); 00861 if( MB_SUCCESS != rval ) return rval; 00862 adjs = subtract( adjs, extra ); 00863 } 00864 } 00865 00866 return rval; 00867 } 00868 00869 // so far, this one is 00870 // used only for __MKModelEntityGeo tag 00871 00872 ErrorCode FBEngine::createTag( const char* tag_name, int tag_size, int tag_type, Tag& tag_handle_out ) 00873 { 00874 // this is copied from iMesh_MOAB.cpp; different name to not have trouble 00875 // with it 00876 // also, we do not want to depend on iMesh.h... 00877 // iMesh is more complicated, because of the options passed 00878 00879 DataType mb_data_type_table2[] = { MB_TYPE_OPAQUE, MB_TYPE_INTEGER, MB_TYPE_DOUBLE, MB_TYPE_HANDLE, 00880 MB_TYPE_HANDLE }; 00881 moab::TagType storage = MB_TAG_SPARSE; 00882 ErrorCode result; 00883 00884 result = 00885 MBI->tag_get_handle( tag_name, tag_size, mb_data_type_table2[tag_type], tag_handle_out, storage | MB_TAG_EXCL ); 00886 00887 if( MB_SUCCESS != result ) 00888 { 00889 std::string msg( "iMesh_createTag: " ); 00890 if( MB_ALREADY_ALLOCATED == result ) 00891 { 00892 msg += "Tag already exists with name: \""; 00893 msg += tag_name; 00894 std::cout << msg << "\n"; 00895 } 00896 else 00897 { 00898 std::cout << "Failed to create tag with name: " << tag_name << "\n"; 00899 return MB_FAILURE; 00900 } 00901 } 00902 00903 // end copy 00904 return MB_SUCCESS; 00905 } 00906 00907 ErrorCode FBEngine::getArrData( const EntityHandle* entity_handles, int entity_handles_size, Tag tag_handle, 00908 void* tag_values_out ) 00909 { 00910 // responsibility of the user to have tag_values_out properly allocated 00911 // only some types of Tags are possible (double, int, etc) 00912 return MBI->tag_get_data( tag_handle, entity_handles, entity_handles_size, tag_values_out ); 00913 } 00914 00915 ErrorCode FBEngine::setArrData( const EntityHandle* entity_handles, int entity_handles_size, Tag tag_handle, 00916 const void* tag_values ) 00917 { 00918 // responsibility of the user to have tag_values_out properly allocated 00919 // only some types of Tags are possible (double, int, etc) 00920 return MBI->tag_set_data( tag_handle, entity_handles, entity_handles_size, tag_values ); 00921 } 00922 00923 ErrorCode FBEngine::getEntAdj( EntityHandle handle, int type_requested, Range& adjEnts ) 00924 { 00925 return getAdjacentEntities( handle, type_requested, adjEnts ); 00926 } 00927 00928 ErrorCode FBEngine::getEgFcSense( EntityHandle mbedge, EntityHandle mbface, int& sense_out ) 00929 { 00930 00931 // this one is important, for establishing the orientation of the edges in faces 00932 // use senses 00933 std::vector< EntityHandle > faces; 00934 std::vector< int > senses; // 0 is forward and 1 is backward 00935 ErrorCode rval = _my_geomTopoTool->get_senses( mbedge, faces, senses ); 00936 if( MB_SUCCESS != rval ) return rval; 00937 00938 for( unsigned int i = 0; i < faces.size(); i++ ) 00939 { 00940 if( faces[i] == mbface ) 00941 { 00942 sense_out = senses[i]; 00943 return MB_SUCCESS; 00944 } 00945 } 00946 return MB_FAILURE; 00947 } 00948 // we assume the measures array was allocated correctly 00949 ErrorCode FBEngine::measure( const EntityHandle* moab_entities, int entities_size, double* measures ) 00950 { 00951 ErrorCode rval; 00952 for( int i = 0; i < entities_size; i++ ) 00953 { 00954 measures[i] = 0.; 00955 00956 int type; 00957 EntityHandle gset = moab_entities[i]; 00958 rval = getEntType( gset, &type ); 00959 if( MB_SUCCESS != rval ) return rval; 00960 if( type == 1 ) 00961 { // edge: get all edges part of the edge set 00962 Range entities; 00963 rval = MBI->get_entities_by_type( gset, MBEDGE, entities ); 00964 if( MB_SUCCESS != rval ) return rval; 00965 00966 for( Range::iterator it = entities.begin(); it != entities.end(); ++it ) 00967 { 00968 EntityHandle edge = *it; 00969 CartVect vv[2]; 00970 const EntityHandle* conn2 = NULL; 00971 int num_nodes; 00972 rval = MBI->get_connectivity( edge, conn2, num_nodes ); 00973 if( MB_SUCCESS != rval || num_nodes != 2 ) return MB_FAILURE; 00974 rval = MBI->get_coords( conn2, 2, (double*)&( vv[0][0] ) ); 00975 if( MB_SUCCESS != rval ) return rval; 00976 00977 vv[0] = vv[1] - vv[0]; 00978 measures[i] += vv[0].length(); 00979 } 00980 } 00981 if( type == 2 ) 00982 { // surface 00983 // get triangles in surface; TODO: quads! 00984 Range entities; 00985 rval = MBI->get_entities_by_type( gset, MBTRI, entities ); 00986 if( MB_SUCCESS != rval ) return rval; 00987 00988 for( Range::iterator it = entities.begin(); it != entities.end(); ++it ) 00989 { 00990 EntityHandle tri = *it; 00991 CartVect vv[3]; 00992 const EntityHandle* conn3 = NULL; 00993 int num_nodes; 00994 rval = MBI->get_connectivity( tri, conn3, num_nodes ); 00995 if( MB_SUCCESS != rval || num_nodes != 3 ) return MB_FAILURE; 00996 rval = MBI->get_coords( conn3, 3, (double*)&( vv[0][0] ) ); 00997 if( MB_SUCCESS != rval ) return rval; 00998 00999 vv[1] = vv[1] - vv[0]; 01000 vv[2] = vv[2] - vv[0]; 01001 vv[0] = vv[1] * vv[2]; 01002 measures[i] += vv[0].length() / 2; // area of triangle 01003 } 01004 } 01005 } 01006 return MB_SUCCESS; 01007 } 01008 01009 ErrorCode FBEngine::getEntNrmlSense( EntityHandle /*face*/, EntityHandle /*region*/, int& /*sense*/ ) 01010 { 01011 return MB_NOT_IMPLEMENTED; // not implemented 01012 } 01013 01014 ErrorCode FBEngine::getEgEvalXYZ( EntityHandle /*edge*/, double /*x*/, double /*y*/, double /*z*/, double& /*on_x*/, 01015 double& /*on_y*/, double& /*on_z*/, double& /*tngt_i*/, double& /*tngt_j*/, 01016 double& /*tngt_k*/, double& /*cvtr_i*/, double& /*cvtr_j*/, double& /*cvtr_k*/ ) 01017 { 01018 return MB_NOT_IMPLEMENTED; // not implemented 01019 } 01020 ErrorCode FBEngine::getFcEvalXYZ( EntityHandle /*face*/, double /*x*/, double /*y*/, double /*z*/, double& /*on_x*/, 01021 double& /*on_y*/, double& /*on_z*/, double& /*nrml_i*/, double& /*nrml_j*/, 01022 double& /*nrml_k*/, double& /*cvtr1_i*/, double& /*cvtr1_j*/, double& /*cvtr1_k*/, 01023 double& /*cvtr2_i*/, double& /*cvtr2_j*/, double& /*cvtr2_k*/ ) 01024 { 01025 return MB_NOT_IMPLEMENTED; // not implemented 01026 } 01027 01028 ErrorCode FBEngine::getEgVtxSense( EntityHandle edge, EntityHandle vtx1, EntityHandle vtx2, int& sense ) 01029 { 01030 // need to decide first or second vertex 01031 // important for moab 01032 int type; 01033 01034 EntityHandle v1, v2; 01035 ErrorCode rval = getEntType( vtx1, &type ); 01036 if( MB_SUCCESS != rval || type != 0 ) return MB_FAILURE; 01037 // edge: get one vertex as part of the vertex set 01038 Range entities; 01039 rval = MBI->get_entities_by_type( vtx1, MBVERTEX, entities ); 01040 if( MB_SUCCESS != rval ) return rval; 01041 if( entities.size() < 1 ) return MB_FAILURE; 01042 v1 = entities[0]; // the first vertex 01043 entities.clear(); 01044 rval = getEntType( vtx2, &type ); 01045 if( MB_SUCCESS != rval || type != 0 ) return MB_FAILURE; 01046 rval = MBI->get_entities_by_type( vtx2, MBVERTEX, entities ); 01047 if( MB_SUCCESS != rval ) return rval; 01048 if( entities.size() < 1 ) return MB_FAILURE; 01049 v2 = entities[0]; // the first vertex 01050 entities.clear(); 01051 // now get the edges, and get the first node and the last node in sequence of edges 01052 // the order is important... 01053 // these are ordered sets !! 01054 std::vector< EntityHandle > ents; 01055 rval = MBI->get_entities_by_type( edge, MBEDGE, ents ); 01056 if( MB_SUCCESS != rval ) return rval; 01057 if( ents.size() < 1 ) return MB_FAILURE; 01058 01059 const EntityHandle* conn = NULL; 01060 int len; 01061 EntityHandle startNode, endNode; 01062 rval = MBI->get_connectivity( ents[0], conn, len ); 01063 if( MB_SUCCESS != rval ) return rval; 01064 startNode = conn[0]; 01065 rval = MBI->get_connectivity( ents[ents.size() - 1], conn, len ); 01066 if( MB_SUCCESS != rval ) return rval; 01067 01068 endNode = conn[1]; 01069 sense = 1; // 01070 if( ( startNode == endNode ) && ( v1 == startNode ) ) 01071 { 01072 sense = 0; // periodic 01073 } 01074 if( ( startNode == v1 ) && ( endNode == v2 ) ) 01075 { 01076 sense = 1; // forward 01077 } 01078 if( ( startNode == v2 ) && ( endNode == v1 ) ) 01079 { 01080 sense = -1; // reverse 01081 } 01082 return MB_SUCCESS; 01083 } 01084 01085 ErrorCode FBEngine::getEntURange( EntityHandle edge, double& u_min, double& u_max ) 01086 { 01087 SmoothCurve* smoothCurve = _edges[edge]; // this is a map 01088 // now, call smoothCurve methods 01089 smoothCurve->get_param_range( u_min, u_max ); 01090 return MB_SUCCESS; 01091 } 01092 01093 ErrorCode FBEngine::getEntUtoXYZ( EntityHandle edge, double u, double& x, double& y, double& z ) 01094 { 01095 SmoothCurve* smoothCurve = _edges[edge]; // this is a map 01096 // now, call smoothCurve methods 01097 smoothCurve->position_from_u( u, x, y, z ); 01098 return MB_SUCCESS; 01099 } 01100 01101 ErrorCode FBEngine::getEntTgntU( EntityHandle edge, double u, double& i, double& j, double& k ) 01102 { 01103 SmoothCurve* smoothCurve = _edges[edge]; // this is a map 01104 // now, call smoothCurve methods 01105 double tg[3]; 01106 double x, y, z; 01107 smoothCurve->position_from_u( u, x, y, z, tg ); 01108 i = tg[0]; 01109 j = tg[1]; 01110 k = tg[2]; 01111 return MB_SUCCESS; 01112 } 01113 ErrorCode FBEngine::isEntAdj( EntityHandle entity1, EntityHandle entity2, bool& adjacent_out ) 01114 { 01115 int type1, type2; 01116 ErrorCode rval = getEntType( entity1, &type1 ); 01117 if( MB_SUCCESS != rval ) return rval; 01118 rval = getEntType( entity2, &type2 ); 01119 if( MB_SUCCESS != rval ) return rval; 01120 01121 Range adjs; 01122 if( type1 < type2 ) 01123 { 01124 rval = MBI->get_parent_meshsets( entity1, adjs, type2 - type1 ); 01125 if( MB_SUCCESS != rval ) return rval; // MBERRORR("Failed to get parent meshsets in iGeom_isEntAdj."); 01126 } 01127 else 01128 { 01129 // note: if they ave the same type, they will not be adjacent, in our definition 01130 rval = MBI->get_child_meshsets( entity1, adjs, type1 - type2 ); 01131 if( MB_SUCCESS != rval ) return rval; // MBERRORR("Failed to get child meshsets in iGeom_isEntAdj."); 01132 } 01133 01134 // adjacent_out = adjs.find(entity2) != _my_gsets[type2].end(); 01135 // hmmm, possible bug here; is this called? 01136 adjacent_out = adjs.find( entity2 ) != adjs.end(); 01137 01138 return MB_SUCCESS; 01139 } 01140 01141 ErrorCode FBEngine::split_surface_with_direction( EntityHandle face, std::vector< double >& xyz, double* direction, 01142 int closed, double min_dot, EntityHandle& oNewFace ) 01143 { 01144 01145 // first of all, find all intersection points (piercing in the face along the direction) 01146 // assume it is robust; what if it is not sufficiently robust? 01147 // if the polyline is open, find the intersection with the boundary edges, of the 01148 // polyline extruded at ends 01149 01150 ErrorCode rval; 01151 01152 // then find the position 01153 int numIniPoints = (int)xyz.size() / 3; 01154 if( ( closed && numIniPoints < 3 ) || ( !closed && numIniPoints < 2 ) ) 01155 MBERRORR( MB_FAILURE, "not enough polyline points " ); 01156 EntityHandle rootForFace; 01157 01158 rval = _my_geomTopoTool->get_root( face, rootForFace ); 01159 MBERRORR( rval, "Failed to get root of face." ); 01160 01161 const double dir[] = { direction[0], direction[1], direction[2] }; 01162 std::vector< EntityHandle > nodes; // get the nodes closest to the ray traces of interest 01163 01164 // these are nodes on the boundary of original face; 01165 // if the cutting line is not closed, the starting - ending vertices of the 01166 // polygonal line must come from this list 01167 01168 std::vector< CartVect > b_pos; 01169 std::vector< EntityHandle > boundary_nodes; 01170 std::vector< EntityHandle > splittingNodes; 01171 Range boundary_mesh_edges; 01172 if( !closed ) 01173 { 01174 rval = boundary_nodes_on_face( face, boundary_nodes ); 01175 MBERRORR( rval, "Failed to get boundary nodes." ); 01176 b_pos.resize( boundary_nodes.size() ); 01177 rval = _mbImpl->get_coords( &( boundary_nodes[0] ), boundary_nodes.size(), (double*)( &b_pos[0][0] ) ); 01178 MBERRORR( rval, "Failed to get coordinates for boundary nodes." ); 01179 rval = boundary_mesh_edges_on_face( face, boundary_mesh_edges ); 01180 MBERRORR( rval, "Failed to get mesh boundary edges for face." ); 01181 } 01182 // 01183 int i = 0; 01184 CartVect dirct( direction ); 01185 dirct.normalize(); // maybe an overkill? 01186 for( ; i < numIniPoints; i++ ) 01187 { 01188 01189 const double point[] = { xyz[3 * i], xyz[3 * i + 1], xyz[3 * i + 2] }; // or even point( &(xyz[3*i]) ); // 01190 CartVect p1( point ); 01191 if( !closed && ( ( 0 == i ) || ( numIniPoints - 1 == i ) ) ) 01192 { 01193 01194 // find the intersection point between a plane and boundary mesh edges 01195 // this will be the closest point on the boundary of face 01196 /// the segment is the first or last segment in the polyline 01197 int i1 = i + 1; 01198 if( i == numIniPoints - 1 ) i1 = i - 1; // previous point if the last 01199 // the direction is from point to point1 01200 const double point1[] = { xyz[3 * i1], xyz[3 * i1 + 1], xyz[3 * i1 + 2] }; 01201 CartVect p2( point1 ); 01202 CartVect normPlane = ( p2 - p1 ) * dirct; 01203 normPlane.normalize(); 01204 //(roughly, from p1 to p2, perpendicular to dirct, in the "xy" plane 01205 // if the intx point is "outside" p1 - p2, skip if the intx point is closer to p2 01206 CartVect perpDir = dirct * normPlane; 01207 Range::iterator ite = boundary_mesh_edges.begin(); 01208 // do a linear search for the best intersection point position (on a boundary edge) 01209 if( debug_splits ) 01210 { 01211 std::cout << " p1:" << p1 << "\n"; 01212 std::cout << " p2:" << p2 << "\n"; 01213 std::cout << " perpDir:" << perpDir << "\n"; 01214 std::cout << " boundary edges size:" << boundary_mesh_edges.size() << "\n"; 01215 } 01216 for( ; ite != boundary_mesh_edges.end(); ++ite ) 01217 { 01218 EntityHandle candidateEdge = *ite; 01219 const EntityHandle* conn2; 01220 int nno; 01221 rval = _mbImpl->get_connectivity( candidateEdge, conn2, nno ); 01222 MBERRORR( rval, "Failed to get conn for boundary edge" ); 01223 CartVect pts[2]; 01224 rval = _mbImpl->get_coords( conn2, 2, &( pts[0][0] ) ); 01225 MBERRORR( rval, "Failed to get coords of nodes for boundary edge" ); 01226 CartVect intx_point; 01227 double parPos; 01228 bool intersect = 01229 intersect_segment_and_plane_slice( pts[0], pts[1], p1, p2, dirct, normPlane, intx_point, parPos ); 01230 if( debug_splits ) 01231 { 01232 std::cout << " Edge:" << _mbImpl->id_from_handle( candidateEdge ) << "\n"; 01233 std::cout << " Node 1:" << _mbImpl->id_from_handle( conn2[0] ) << pts[0] << "\n"; 01234 std::cout << " Node 2:" << _mbImpl->id_from_handle( conn2[1] ) << pts[1] << "\n"; 01235 std::cout << " Intersect bool:" << intersect << "\n"; 01236 } 01237 if( intersect ) 01238 { 01239 double proj1 = ( intx_point - p1 ) % perpDir; 01240 double proj2 = ( intx_point - p2 ) % perpDir; 01241 if( ( fabs( proj1 ) > fabs( proj2 ) ) // this means it is closer to p2 than p1 01242 ) 01243 continue; // basically, this means the intersection point is with a 01244 // boundary edge on the other side, closer to p2 than p1, so we 01245 // skip it 01246 if( parPos == 0 ) 01247 { 01248 // close to vertex 1, nothing to do 01249 nodes.push_back( conn2[0] ); 01250 splittingNodes.push_back( conn2[0] ); 01251 } 01252 else if( parPos == 1. ) 01253 { 01254 // close to vertex 2, nothing to do 01255 nodes.push_back( conn2[1] ); 01256 splittingNodes.push_back( conn2[1] ); 01257 } 01258 else 01259 { 01260 // break the edge, create a new node at intersection point (will be smoothed 01261 // out) 01262 EntityHandle newVertex; 01263 rval = _mbImpl->create_vertex( &( intx_point[0] ), newVertex ); 01264 MBERRORR( rval, "can't create vertex" ); 01265 nodes.push_back( newVertex ); 01266 split_internal_edge( candidateEdge, newVertex ); 01267 splittingNodes.push_back( newVertex ); 01268 _brokenEdges[newVertex] = candidateEdge; 01269 _piercedEdges.insert( candidateEdge ); 01270 } 01271 break; // break from the loop over boundary edges, we are interested in the 01272 // first 01273 // split (hopefully, the only split) 01274 } 01275 } 01276 if( ite == boundary_mesh_edges.end() ) 01277 MBERRORR( MB_FAILURE, "Failed to find boundary intersection edge. Bail out" ); 01278 } 01279 else 01280 { 01281 std::vector< double > distances_out; 01282 std::vector< EntityHandle > sets_out; 01283 std::vector< EntityHandle > facets_out; 01284 rval = _my_geomTopoTool->obb_tree()->ray_intersect_sets( distances_out, sets_out, facets_out, rootForFace, 01285 tolerance, point, dir ); 01286 MBERRORR( rval, "Failed to get ray intersections." ); 01287 if( distances_out.size() < 1 ) 01288 MBERRORR( MB_FAILURE, "Failed to get one intersection point, bad direction." ); 01289 01290 if( distances_out.size() > 1 ) 01291 { 01292 std::cout << " too many intersection points. Only the first one considered\n"; 01293 } 01294 std::vector< EntityHandle >::iterator pFace = std::find( sets_out.begin(), sets_out.end(), face ); 01295 01296 if( pFace == sets_out.end() ) MBERRORR( MB_FAILURE, "Failed to intersect given face, bad direction." ); 01297 unsigned int index = pFace - sets_out.begin(); 01298 // get the closest node of the triangle, and modify locally the triangle(s), so the 01299 // intersection point is at a new vertex, if needed 01300 CartVect P( point ); 01301 CartVect Dir( dir ); 01302 CartVect newPoint = P + distances_out[index] * Dir; 01303 // get the triangle coordinates 01304 01305 double area_coord[3]; 01306 EntityHandle boundary_handle = 0; // if 0, not on a boundary 01307 rval = area_coordinates( _mbImpl, facets_out[index], newPoint, area_coord, boundary_handle ); 01308 MBERRORR( rval, "Failed to get area coordinates" ); 01309 01310 if( debug_splits ) 01311 { 01312 std::cout << " int point:" << newPoint << " area coord " << area_coord[0] << " " << area_coord[1] << " " 01313 << area_coord[2] << "\n"; 01314 std::cout << " triangle: " << _mbImpl->id_from_handle( facets_out[index] ) 01315 << " boundary:" << boundary_handle << "\n"; 01316 } 01317 EntityType type; 01318 if( boundary_handle ) type = _mbImpl->type_from_handle( boundary_handle ); 01319 if( boundary_handle && ( type == MBVERTEX ) ) 01320 { 01321 // nothing to do, we are happy 01322 nodes.push_back( boundary_handle ); 01323 } 01324 else 01325 { 01326 // for an edge, we will split 2 triangles 01327 // for interior point, we will create 3 triangles out of one 01328 // create a new vertex 01329 EntityHandle newVertex; 01330 rval = _mbImpl->create_vertex( &( newPoint[0] ), newVertex ); 01331 if( boundary_handle ) // this is edge 01332 { 01333 split_internal_edge( boundary_handle, newVertex ); 01334 _piercedEdges.insert( boundary_handle ); // to be removed at the end 01335 } 01336 else 01337 divide_triangle( facets_out[index], newVertex ); 01338 01339 nodes.push_back( newVertex ); 01340 } 01341 } 01342 } 01343 // now, we have to find more intersection points, either interior to triangles, or on edges, or 01344 // on vertices use the same tolerance as before starting from 2 points on 2 triangles, and 01345 // having the direction, get more intersection points between the plane formed by direction and 01346 // those 2 points, and edges from triangulation (the triangles involved will be part of the same 01347 // gentity , original face ( moab set) 01348 // int closed = 1;// closed = 0 if the polyline is not closed 01349 01350 CartVect Dir( direction ); 01351 std::vector< EntityHandle > chainedEdges; 01352 01353 for( i = 0; i < numIniPoints - 1 + closed; i++ ) 01354 { 01355 int nextIndex = ( i + 1 ) % numIniPoints; 01356 std::vector< EntityHandle > trianglesAlong; 01357 std::vector< CartVect > points; 01358 // otherwise to edges or even nodes 01359 std::vector< EntityHandle > entities; 01360 // start with initial points, intersect along the direction, find the facets 01361 rval = compute_intersection_points( face, nodes[i], nodes[nextIndex], Dir, points, entities, trianglesAlong ); 01362 MBERRORR( rval, "can't get intersection points along a line" ); 01363 std::vector< EntityHandle > nodesAlongPolyline; 01364 // refactor code; move some edge creation for each 2 intersection points 01365 nodesAlongPolyline.push_back( entities[0] ); // it is for sure a node 01366 int num_points = (int)points.size(); // it should be num_triangles + 1 01367 for( int j = 0; j < num_points - 1; j++ ) 01368 { 01369 EntityHandle tri = trianglesAlong[j]; // this is happening in trianglesAlong i 01370 EntityHandle e1 = entities[j]; 01371 EntityHandle e2 = entities[j + 1]; 01372 EntityType et1 = _mbImpl->type_from_handle( e1 ); 01373 // EntityHandle vertex1 = nodesAlongPolyline[i];// irrespective of the entity type i, 01374 // we already have the vertex there 01375 EntityType et2 = _mbImpl->type_from_handle( e2 ); 01376 if( et2 == MBVERTEX ) { nodesAlongPolyline.push_back( e2 ); } 01377 else // if (et2==MBEDGE) 01378 { 01379 CartVect coord_vert = points[j + 1]; 01380 EntityHandle newVertex; 01381 rval = _mbImpl->create_vertex( (double*)&coord_vert, newVertex ); 01382 MBERRORR( rval, "can't create vertex" ); 01383 nodesAlongPolyline.push_back( newVertex ); 01384 } 01385 // if vertices, do not split anything, just get the edge for polyline 01386 if( et2 == MBVERTEX && et1 == MBVERTEX ) 01387 { 01388 // nothing to do, just continue; 01389 continue; // continue the for loop 01390 } 01391 01392 if( debug_splits ) 01393 { 01394 std::cout << "tri: type: " << _mbImpl->type_from_handle( tri ) 01395 << " id:" << _mbImpl->id_from_handle( tri ) << "\n e1:" << e1 01396 << " id:" << _mbImpl->id_from_handle( e1 ) << " e2:" << e2 01397 << " id:" << _mbImpl->id_from_handle( e2 ) << "\n"; 01398 } 01399 // here, at least one is an edge 01400 rval = BreakTriangle2( tri, e1, e2, nodesAlongPolyline[j], nodesAlongPolyline[j + 1] ); 01401 MBERRORR( rval, "can't break triangle 2" ); 01402 if( et2 == MBEDGE ) _piercedEdges.insert( e2 ); 01403 _piercedTriangles.insert( tri ); 01404 } 01405 // nodesAlongPolyline will define a new geometric edge 01406 if( debug_splits ) 01407 { 01408 std::cout << "nodesAlongPolyline: " << nodesAlongPolyline.size() << "\n"; 01409 std::cout << "points: " << num_points << "\n"; 01410 } 01411 // if needed, create edges along polyline, or revert the existing ones, to 01412 // put them in a new edge set 01413 EntityHandle new_geo_edge; 01414 rval = create_new_gedge( nodesAlongPolyline, new_geo_edge ); 01415 MBERRORR( rval, "can't create a new edge" ); 01416 chainedEdges.push_back( new_geo_edge ); 01417 // end copy 01418 } 01419 // the segment between point_i and point_i+1 is in trianglesAlong_i 01420 // points_i is on entities_i 01421 // all these edges are oriented correctly 01422 rval = split_surface( face, chainedEdges, splittingNodes, oNewFace ); 01423 MBERRORR( rval, "can't split surface" ); 01424 // 01425 rval = chain_edges( min_dot ); // acos(0.8)~= 36 degrees 01426 MBERRORR( rval, "can't chain edges" ); 01427 return MB_SUCCESS; 01428 } 01429 /** 01430 * this method splits along the polyline defined by points and entities 01431 * the polyline will be defined with 01432 * // the entities are now only nodes and edges, no triangles!!! 01433 * the first and last ones are also nodes for sure 01434 */ 01435 ErrorCode FBEngine::split_surface( EntityHandle face, std::vector< EntityHandle >& chainedEdges, 01436 std::vector< EntityHandle >& splittingNodes, EntityHandle& newFace ) 01437 { 01438 // use now the chained edges to create a new face (loop or clean split) 01439 // use a fill to determine the new sets, up to the polyline 01440 // points on the polyline will be moved to the closest point location, with some constraints 01441 // then the sets will be reset, geometry recomputed. new vertices, new edges, etc. 01442 01443 Range iniTris; 01444 ErrorCode rval; 01445 rval = _mbImpl->get_entities_by_type( face, MBTRI, iniTris ); 01446 MBERRORR( rval, "can't get initial triangles" ); 01447 01448 // start from a triangle that is not in the triangles to delete 01449 // flood fill 01450 01451 bool closed = splittingNodes.size() == 0; 01452 if( !closed ) 01453 { 01454 // 01455 if( splittingNodes.size() != 2 ) MBERRORR( MB_FAILURE, "need to have exactly 2 nodes for splitting" ); 01456 // we will have to split the boundary edges 01457 // first, find the actual boundary, and try to split with the 2 end points (nodes) 01458 // get the adjacent edges, and see which one has the end nodes 01459 rval = split_boundary( face, splittingNodes[0] ); 01460 MBERRORR( rval, "can't split with first node" ); 01461 rval = split_boundary( face, splittingNodes[1] ); 01462 MBERRORR( rval, "can't split with second node)" ); 01463 } 01464 // we will separate triangles to delete, unaffected, new_triangles, 01465 // nodesAlongPolyline, 01466 Range first, second; 01467 rval = separate( face, chainedEdges, first, second ); 01468 01469 // now, we are done with the computations; 01470 // we need to put the new nodes on the smooth surface 01471 if( this->_smooth ) 01472 { 01473 rval = smooth_new_intx_points( face, chainedEdges ); 01474 MBERRORR( rval, "can't smooth new points" ); 01475 } 01476 01477 // create the new set 01478 rval = _mbImpl->create_meshset( MESHSET_SET, newFace ); 01479 MBERRORR( rval, "can't create a new face" ); 01480 01481 _my_geomTopoTool->add_geo_set( newFace, 2 ); 01482 01483 // the new face will have the first set (positive sense triangles, to the left) 01484 rval = _mbImpl->add_entities( newFace, first ); 01485 MBERRORR( rval, "can't add first range triangles to new face" ); 01486 01487 for( unsigned int j = 0; j < chainedEdges.size(); j++ ) 01488 { 01489 EntityHandle new_geo_edge = chainedEdges[j]; 01490 // both faces will have the edge now 01491 rval = _mbImpl->add_parent_child( face, new_geo_edge ); 01492 MBERRORR( rval, "can't add parent child relations for new edge" ); 01493 01494 rval = _mbImpl->add_parent_child( newFace, new_geo_edge ); 01495 MBERRORR( rval, "can't add parent child relations for new edge" ); 01496 // add senses 01497 // sense newFace is 1, old face is -1 01498 rval = _my_geomTopoTool->set_sense( new_geo_edge, newFace, 1 ); 01499 MBERRORR( rval, "can't set sense for new edge" ); 01500 01501 rval = _my_geomTopoTool->set_sense( new_geo_edge, face, -1 ); 01502 MBERRORR( rval, "can't set sense for new edge in original face" ); 01503 } 01504 01505 rval = set_neumann_tags( face, newFace ); 01506 MBERRORR( rval, "can't set NEUMANN set tags" ); 01507 01508 // now, we should remove from the original set all tris, and put the "second" range 01509 rval = _mbImpl->remove_entities( face, iniTris ); 01510 MBERRORR( rval, "can't remove original tris from initial face set" ); 01511 01512 rval = _mbImpl->add_entities( face, second ); 01513 MBERRORR( rval, "can't add second range to the original set" ); 01514 01515 if( !closed ) 01516 { 01517 rval = redistribute_boundary_edges_to_faces( face, newFace, chainedEdges ); 01518 MBERRORR( rval, "fail to reset the proper boundary faces" ); 01519 } 01520 01521 /*if (_smooth) 01522 delete_smooth_tags();// they need to be recomputed, anyway 01523 // this will remove the extra smooth faces and edges 01524 clean();*/ 01525 // also, these nodes need to be moved to the smooth surface, sometimes before deleting the old 01526 // triangles 01527 // remove the triangles from the set, then delete triangles (also some edges need to be 01528 // deleted!) 01529 rval = _mbImpl->delete_entities( _piercedTriangles ); 01530 01531 MBERRORR( rval, "can't delete triangles" ); 01532 _piercedTriangles.clear(); 01533 // delete edges that are broke up in 2 01534 rval = _mbImpl->delete_entities( _piercedEdges ); 01535 MBERRORR( rval, "can't delete edges" ); 01536 _piercedEdges.clear(); 01537 01538 if( debug_splits ) 01539 { 01540 _mbImpl->write_file( "newFace.vtk", "vtk", 0, &newFace, 1 ); 01541 _mbImpl->write_file( "leftoverFace.vtk", "vtk", 0, &face, 1 ); 01542 } 01543 return MB_SUCCESS; 01544 } 01545 01546 ErrorCode FBEngine::smooth_new_intx_points( EntityHandle face, std::vector< EntityHandle >& chainedEdges ) 01547 { 01548 01549 // do not move nodes from the original face 01550 // first get all triangles, and then all nodes from those triangles 01551 01552 Range tris; 01553 ErrorCode rval = _mbImpl->get_entities_by_type( face, MBTRI, tris ); 01554 MBERRORR( rval, "can't get triangles" ); 01555 01556 Range ini_nodes; 01557 rval = _mbImpl->get_connectivity( tris, ini_nodes ); 01558 MBERRORR( rval, "can't get connectivities" ); 01559 01560 SmoothFace* smthFace = _faces[face]; 01561 01562 // get all nodes from chained edges 01563 Range mesh_edges; 01564 for( unsigned int j = 0; j < chainedEdges.size(); j++ ) 01565 { 01566 // keep adding to the range of mesh edges 01567 rval = _mbImpl->get_entities_by_dimension( chainedEdges[j], 1, mesh_edges ); 01568 MBERRORR( rval, "can't get mesh edges" ); 01569 } 01570 // nodes along polyline 01571 Range nodes_on_polyline; 01572 rval = _mbImpl->get_connectivity( mesh_edges, nodes_on_polyline, true ); // corners only 01573 MBERRORR( rval, "can't get nodes on the polyline" ); 01574 01575 Range new_intx_nodes = subtract( nodes_on_polyline, ini_nodes ); 01576 01577 std::vector< double > ini_coords; 01578 int num_points = (int)new_intx_nodes.size(); 01579 ini_coords.resize( 3 * num_points ); 01580 rval = _mbImpl->get_coords( new_intx_nodes, &( ini_coords[0] ) ); 01581 MBERRORR( rval, "can't get coordinates" ); 01582 01583 int i = 0; 01584 for( Range::iterator it = new_intx_nodes.begin(); it != new_intx_nodes.end(); ++it ) 01585 { 01586 /*EntityHandle node = *it;*/ 01587 int i3 = 3 * i; 01588 smthFace->move_to_surface( ini_coords[i3], ini_coords[i3 + 1], ini_coords[i3 + 2] ); 01589 // reset the coordinates of this node 01590 ++i; 01591 } 01592 rval = _mbImpl->set_coords( new_intx_nodes, &( ini_coords[0] ) ); 01593 MBERRORR( rval, "can't set smoothed coordinates for the new nodes" ); 01594 01595 return MB_SUCCESS; 01596 } 01597 // we will use the fact that the splitting edge is oriented right now 01598 // to the left will be new face, to the right, old face 01599 // (to the left, positively oriented triangles) 01600 ErrorCode FBEngine::separate( EntityHandle face, std::vector< EntityHandle >& chainedEdges, Range& first, 01601 Range& second ) 01602 { 01603 // Range unaffectedTriangles = subtract(iniTriangles, _piercedTriangles); 01604 // insert in each 01605 // start with a new triangle, and flood to get the first range; what is left is the 01606 // second range 01607 // flood fill is considering edges adjacent to seed triangles; if there is 01608 // an edge in the new_geo_edge, it is skipped; triangles in the 01609 // triangles to delete are not added 01610 // first, create all edges of the new triangles 01611 01612 // 01613 // new face will have the new edge oriented positively 01614 // get mesh edges from geo edge (splitting gedge); 01615 01616 Range mesh_edges; 01617 ErrorCode rval; 01618 // mesh_edges 01619 for( unsigned int j = 0; j < chainedEdges.size(); j++ ) 01620 { 01621 // this will keep adding edges to the mesh_edges range 01622 // when we get out, the mesh_edges will be in this range, but not ordered 01623 rval = _mbImpl->get_entities_by_type( chainedEdges[j], MBEDGE, mesh_edges ); 01624 MBERRORR( rval, "can't get new polyline edges" ); 01625 if( debug_splits ) 01626 { 01627 std::cout << " At chained edge " << j << " " << _mbImpl->id_from_handle( chainedEdges[j] ) 01628 << " mesh_edges Range size:" << mesh_edges.size() << "\n"; 01629 } 01630 } 01631 01632 // get a positive triangle adjacent to mesh_edge[0] 01633 // add to first triangles to the left, second triangles to the right of the mesh_edges ; 01634 01635 // create a temp tag, and when done, delete it 01636 // default value: 0 01637 // 3 to be deleted, pierced 01638 // 1 first set 01639 // 2 second set 01640 // create the tag also for control points on the edges 01641 int defVal = 0; 01642 Tag separateTag; 01643 rval = MBI->tag_get_handle( "SEPARATE_TAG", 1, MB_TYPE_INTEGER, separateTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal ); 01644 MBERRORR( rval, "can't create temp tag for separation" ); 01645 // the deleted triangles will get a value 3, from start 01646 int delVal = 3; 01647 for( Range::iterator it1 = this->_piercedTriangles.begin(); it1 != _piercedTriangles.end(); ++it1 ) 01648 { 01649 EntityHandle trToDelete = *it1; 01650 rval = _mbImpl->tag_set_data( separateTag, &trToDelete, 1, &delVal ); 01651 MBERRORR( rval, "can't set delete tag value" ); 01652 } 01653 01654 // find a triangle that will be in the first range, positively oriented about the splitting edge 01655 EntityHandle seed1 = 0; 01656 for( Range::iterator it = mesh_edges.begin(); it != mesh_edges.end() && !seed1; ++it ) 01657 { 01658 EntityHandle meshEdge = *it; 01659 Range adj_tri; 01660 rval = _mbImpl->get_adjacencies( &meshEdge, 1, 2, false, adj_tri ); 01661 MBERRORR( rval, "can't get adj_tris to mesh edge" ); 01662 01663 for( Range::iterator it2 = adj_tri.begin(); it2 != adj_tri.end(); ++it2 ) 01664 { 01665 EntityHandle tr = *it2; 01666 if( _piercedTriangles.find( tr ) != _piercedTriangles.end() ) 01667 continue; // do not attach pierced triangles, they are not good 01668 int num1, sense, offset; 01669 rval = _mbImpl->side_number( tr, meshEdge, num1, sense, offset ); 01670 MBERRORR( rval, "edge not adjacent" ); 01671 if( sense == 1 ) 01672 { 01673 // firstSet.insert(tr); 01674 if( !seed1 ) 01675 { 01676 seed1 = tr; 01677 break; 01678 } 01679 } 01680 } 01681 } 01682 01683 // flood fill first set, the rest will be in second set 01684 // the edges from new_geo_edge will not be crossed 01685 01686 // get edges of face (adjacencies) 01687 // also get the old boundary edges, from face; they will be edges to not cross, too 01688 Range bound_edges; 01689 rval = getAdjacentEntities( face, 1, bound_edges ); 01690 MBERRORR( rval, "can't get boundary edges" ); 01691 01692 // add to the do not cross edges range, all edges from initial boundary 01693 Range initialBoundaryEdges; 01694 for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it ) 01695 { 01696 EntityHandle bound_edge = *it; 01697 rval = _mbImpl->get_entities_by_dimension( bound_edge, 1, initialBoundaryEdges ); 01698 } 01699 01700 Range doNotCrossEdges = unite( initialBoundaryEdges, mesh_edges ); // add the splitting edges ! 01701 01702 // use a second method, with tags 01703 // 01704 std::queue< EntityHandle > queue1; 01705 queue1.push( seed1 ); 01706 std::vector< EntityHandle > arr1; 01707 while( !queue1.empty() ) 01708 { 01709 // start copy 01710 EntityHandle currentTriangle = queue1.front(); 01711 queue1.pop(); 01712 arr1.push_back( currentTriangle ); 01713 // add new triangles that share an edge 01714 Range currentEdges; 01715 rval = _mbImpl->get_adjacencies( ¤tTriangle, 1, 1, true, currentEdges, Interface::UNION ); 01716 MBERRORR( rval, "can't get adjacencies" ); 01717 for( Range::iterator it = currentEdges.begin(); it != currentEdges.end(); ++it ) 01718 { 01719 EntityHandle frontEdge = *it; 01720 if( doNotCrossEdges.find( frontEdge ) == doNotCrossEdges.end() ) 01721 { 01722 // this is an edge that can be crossed 01723 Range adj_tri; 01724 rval = _mbImpl->get_adjacencies( &frontEdge, 1, 2, false, adj_tri, Interface::UNION ); 01725 MBERRORR( rval, "can't get adj_tris" ); 01726 // if the triangle is not in first range, add it to the queue 01727 for( Range::iterator it2 = adj_tri.begin(); it2 != adj_tri.end(); ++it2 ) 01728 { 01729 EntityHandle tri2 = *it2; 01730 int val = 0; 01731 rval = _mbImpl->tag_get_data( separateTag, &tri2, 1, &val ); 01732 MBERRORR( rval, "can't get tag value" ); 01733 if( val ) continue; 01734 // else, set it to 1 01735 val = 1; 01736 rval = _mbImpl->tag_set_data( separateTag, &tri2, 1, &val ); 01737 MBERRORR( rval, "can't get tag value" ); 01738 01739 queue1.push( tri2 ); 01740 } 01741 } // end edge do not cross 01742 } // end while 01743 } 01744 01745 std::sort( arr1.begin(), arr1.end() ); 01746 // Range first1; 01747 std::copy( arr1.rbegin(), arr1.rend(), range_inserter( first ) ); 01748 01749 // std::cout<< "\n first1.size() " << first1.size() << " first.size(): " << first.size() << 01750 // "\n"; 01751 if( debug_splits ) 01752 { 01753 EntityHandle tmpSet; 01754 _mbImpl->create_meshset( MESHSET_SET, tmpSet ); 01755 _mbImpl->add_entities( tmpSet, first ); 01756 _mbImpl->write_file( "dbg1.vtk", "vtk", 0, &tmpSet, 1 ); 01757 } 01758 // now, decide the set 2: 01759 // first, get all ini tris 01760 Range initr; 01761 rval = _mbImpl->get_entities_by_type( face, MBTRI, initr ); 01762 MBERRORR( rval, "can't get tris " ); 01763 second = unite( initr, _newTriangles ); 01764 Range second2 = subtract( second, _piercedTriangles ); 01765 second = subtract( second2, first ); 01766 _newTriangles.clear(); 01767 if( debug_splits ) 01768 { 01769 std::cout << "\n second.size() " << second.size() << " first.size(): " << first.size() << "\n"; 01770 // debugging code 01771 EntityHandle tmpSet2; 01772 _mbImpl->create_meshset( MESHSET_SET, tmpSet2 ); 01773 _mbImpl->add_entities( tmpSet2, second ); 01774 _mbImpl->write_file( "dbg2.vtk", "vtk", 0, &tmpSet2, 1 ); 01775 } 01776 /*Range intex = intersect(first, second); 01777 if (!intex.empty() && debug_splits) 01778 { 01779 std::cout << "error, the sets should be disjoint\n"; 01780 for (Range::iterator it1=intex.begin(); it1!=intex.end(); ++it1) 01781 { 01782 std::cout<<_mbImpl->id_from_handle(*it1) << "\n"; 01783 } 01784 }*/ 01785 rval = _mbImpl->tag_delete( separateTag ); 01786 MBERRORR( rval, "can't delete tag " ); 01787 return MB_SUCCESS; 01788 } 01789 // if there is an edge between 2 nodes, then check it's orientation, and revert it if needed 01790 ErrorCode FBEngine::create_new_gedge( std::vector< EntityHandle >& nodesAlongPolyline, EntityHandle& new_geo_edge ) 01791 { 01792 01793 ErrorCode rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_geo_edge ); 01794 MBERRORR( rval, "can't create geo edge" ); 01795 01796 // now, get the edges, or create if not existing 01797 std::vector< EntityHandle > mesh_edges; 01798 for( unsigned int i = 0; i < nodesAlongPolyline.size() - 1; i++ ) 01799 { 01800 EntityHandle n1 = nodesAlongPolyline[i], n2 = nodesAlongPolyline[i + 1]; 01801 01802 EntityHandle nn2[2]; 01803 nn2[0] = n1; 01804 nn2[1] = n2; 01805 01806 std::vector< EntityHandle > adjacent; 01807 rval = _mbImpl->get_adjacencies( nn2, 2, 1, false, adjacent, Interface::INTERSECT ); 01808 // see if there is an edge between those 2 already, and if it is oriented as we like 01809 bool new_edge = true; 01810 if( adjacent.size() >= 1 ) 01811 { 01812 // check the orientation 01813 const EntityHandle* conn2 = NULL; 01814 int nnodes = 0; 01815 rval = _mbImpl->get_connectivity( adjacent[0], conn2, nnodes ); 01816 MBERRORR( rval, "can't get connectivity" ); 01817 if( conn2[0] == nn2[0] && conn2[1] == nn2[1] ) 01818 { 01819 // everything is fine 01820 mesh_edges.push_back( adjacent[0] ); 01821 new_edge = false; // we found one that's good, no need to create a new one 01822 } 01823 else 01824 { 01825 _piercedEdges.insert( adjacent[0] ); // we want to remove this one, it will be not needed 01826 } 01827 } 01828 if( new_edge ) 01829 { 01830 // there is no edge between n1 and n2, create one 01831 EntityHandle mesh_edge; 01832 rval = _mbImpl->create_element( MBEDGE, nn2, 2, mesh_edge ); 01833 MBERRORR( rval, "Failed to create a new edge" ); 01834 mesh_edges.push_back( mesh_edge ); 01835 } 01836 } 01837 01838 // add loops edges to the edge set 01839 rval = _mbImpl->add_entities( new_geo_edge, &mesh_edges[0], 01840 mesh_edges.size() ); // only one edge 01841 MBERRORR( rval, "can't add edges to new_geo_set" ); 01842 // check vertex sets for vertex 1 and vertex 2? 01843 // get all sets of dimension 0 from database, and see if our ends are here or not 01844 01845 Range ends_geo_edge; 01846 ends_geo_edge.insert( nodesAlongPolyline[0] ); 01847 ends_geo_edge.insert( nodesAlongPolyline[nodesAlongPolyline.size() - 1] ); 01848 01849 for( unsigned int k = 0; k < ends_geo_edge.size(); k++ ) 01850 { 01851 EntityHandle node = ends_geo_edge[k]; 01852 EntityHandle nodeSet; 01853 bool found = find_vertex_set_for_node( node, nodeSet ); 01854 01855 if( !found ) 01856 { 01857 // create a node set and add the node 01858 01859 rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet ); 01860 MBERRORR( rval, "Failed to create a new vertex set" ); 01861 01862 rval = _mbImpl->add_entities( nodeSet, &node, 1 ); 01863 MBERRORR( rval, "Failed to add the node to the set" ); 01864 01865 rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 ); // 01866 MBERRORR( rval, "Failed to commit the node set" ); 01867 01868 if( debug_splits ) 01869 { 01870 std::cout << " create a vertex set " << _mbImpl->id_from_handle( nodeSet ) 01871 << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << " for node " << node 01872 << "\n"; 01873 } 01874 } 01875 01876 rval = _mbImpl->add_parent_child( new_geo_edge, nodeSet ); 01877 MBERRORR( rval, "Failed to add parent child relation" ); 01878 } 01879 // finally, put the edge in the range of edges 01880 rval = _my_geomTopoTool->add_geo_set( new_geo_edge, 1 );MB_CHK_ERR( rval ); 01881 01882 return rval; 01883 } 01884 01885 void FBEngine::print_debug_triangle( EntityHandle t ) 01886 { 01887 std::cout << " triangle id:" << _mbImpl->id_from_handle( t ) << "\n"; 01888 const EntityHandle* conn3 = NULL; 01889 int nnodes = 0; 01890 _mbImpl->get_connectivity( t, conn3, nnodes ); 01891 // get coords 01892 CartVect P[3]; 01893 _mbImpl->get_coords( conn3, 3, (double*)&P[0] ); 01894 std::cout << " nodes:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n"; 01895 CartVect PP[3]; 01896 PP[0] = P[1] - P[0]; 01897 PP[1] = P[2] - P[1]; 01898 PP[2] = P[0] - P[2]; 01899 01900 std::cout << " pos:" << P[0] << " " << P[1] << " " << P[2] << "\n"; 01901 std::cout << " x,y diffs " << PP[0][0] << " " << PP[0][1] << ", " << PP[1][0] << " " << PP[1][1] << ", " 01902 << PP[2][0] << " " << PP[2][1] << "\n"; 01903 return; 01904 } 01905 // actual breaking of triangles 01906 // case 1: n2 interior to triangle 01907 ErrorCode FBEngine::BreakTriangle( EntityHandle, EntityHandle, EntityHandle, EntityHandle, EntityHandle, EntityHandle ) 01908 { 01909 std::cout << "FBEngine::BreakTriangle not implemented yet\n"; 01910 return MB_FAILURE; 01911 } 01912 // case 2, n1 and n2 on boundary 01913 ErrorCode FBEngine::BreakTriangle2( EntityHandle tri, EntityHandle e1, EntityHandle e2, EntityHandle n1, 01914 EntityHandle n2 ) // nodesAlongPolyline are on entities! 01915 { 01916 // we have the nodes, we just need to reconnect to form new triangles 01917 ErrorCode rval; 01918 const EntityHandle* conn3 = NULL; 01919 int nnodes = 0; 01920 rval = _mbImpl->get_connectivity( tri, conn3, nnodes ); 01921 MBERRORR( rval, "Failed to get connectivity" ); 01922 assert( 3 == nnodes ); 01923 01924 EntityType et1 = _mbImpl->type_from_handle( e1 ); 01925 EntityType et2 = _mbImpl->type_from_handle( e2 ); 01926 01927 if( MBVERTEX == et1 ) 01928 { 01929 // find the vertex in conn3, and form 2 other triangles 01930 int index = -1; 01931 for( index = 0; index < 3; index++ ) 01932 { 01933 if( conn3[index] == e1 ) // also n1 01934 break; 01935 } 01936 if( index == 3 ) return MB_FAILURE; 01937 // 2 triangles: n1, index+1, n2, and n1, n2, index+2 01938 EntityHandle conn[6] = { n1, conn3[( index + 1 ) % 3], n2, n1, n2, conn3[( index + 2 ) % 3] }; 01939 EntityHandle newTriangle; 01940 rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle ); 01941 MBERRORR( rval, "Failed to create a new triangle" ); 01942 _newTriangles.insert( newTriangle ); 01943 if( debug_splits ) print_debug_triangle( newTriangle ); 01944 rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle ); // the second triangle 01945 MBERRORR( rval, "Failed to create a new triangle" ); 01946 _newTriangles.insert( newTriangle ); 01947 if( debug_splits ) print_debug_triangle( newTriangle ); 01948 } 01949 else if( MBVERTEX == et2 ) 01950 { 01951 int index = -1; 01952 for( index = 0; index < 3; index++ ) 01953 { 01954 if( conn3[index] == e2 ) // also n2 01955 break; 01956 } 01957 if( index == 3 ) return MB_FAILURE; 01958 // 2 triangles: n1, index+1, n2, and n1, n2, index+2 01959 EntityHandle conn[6] = { n2, conn3[( index + 1 ) % 3], n1, n2, n1, conn3[( index + 2 ) % 3] }; 01960 EntityHandle newTriangle; 01961 rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle ); 01962 MBERRORR( rval, "Failed to create a new triangle" ); 01963 _newTriangles.insert( newTriangle ); 01964 if( debug_splits ) print_debug_triangle( newTriangle ); 01965 rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle ); // the second triangle 01966 MBERRORR( rval, "Failed to create a new triangle" ); 01967 _newTriangles.insert( newTriangle ); 01968 if( debug_splits ) print_debug_triangle( newTriangle ); 01969 } 01970 else 01971 { 01972 // both are edges adjacent to triangle tri 01973 // there are several configurations possible for n1, n2, between conn3 nodes. 01974 int num1, num2, sense, offset; 01975 rval = _mbImpl->side_number( tri, e1, num1, sense, offset ); 01976 MBERRORR( rval, "edge not adjacent" ); 01977 01978 rval = _mbImpl->side_number( tri, e2, num2, sense, offset ); 01979 MBERRORR( rval, "edge not adjacent" ); 01980 01981 const EntityHandle* conn12; // connectivity for edge 1 01982 const EntityHandle* conn22; // connectivity for edge 2 01983 // int nnodes; 01984 rval = _mbImpl->get_connectivity( e1, conn12, nnodes ); 01985 MBERRORR( rval, "Failed to get connectivity of edge 1" ); 01986 assert( 2 == nnodes ); 01987 rval = _mbImpl->get_connectivity( e2, conn22, nnodes ); 01988 MBERRORR( rval, "Failed to get connectivity of edge 2" ); 01989 assert( 2 == nnodes ); 01990 // now, having the side number, work through 01991 if( debug_splits ) 01992 { 01993 std::cout << "tri conn3:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n"; 01994 std::cout << " edge1: conn12:" << conn12[0] << " " << conn12[1] << " side: " << num1 << "\n"; 01995 std::cout << " edge2: conn22:" << conn22[0] << " " << conn22[1] << " side: " << num2 << "\n"; 01996 } 01997 int unaffectedSide = 3 - num1 - num2; 01998 int i3 = ( unaffectedSide + 2 ) % 3; // to 0 is 2, to 1 is 0, to 2 is 1 01999 // triangles will be formed with triVertexIndex , n1, n2 (in what order?) 02000 EntityHandle v1, v2; // to hold the 2 nodes on edges 02001 if( num1 == i3 ) 02002 { 02003 v1 = n1; 02004 v2 = n2; 02005 } 02006 else // if (num2==i3) 02007 { 02008 v1 = n2; 02009 v2 = n1; 02010 } 02011 // three triangles are formed 02012 int i1 = ( i3 + 1 ) % 3; 02013 int i2 = ( i3 + 2 ) % 3; 02014 // we could break the surface differently 02015 EntityHandle conn[9] = { conn3[i3], v1, v2, v1, conn3[i1], conn3[i2], v2, v1, conn3[i2] }; 02016 EntityHandle newTriangle; 02017 if( debug_splits ) std::cout << "Split 2 edges :\n"; 02018 rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle ); 02019 MBERRORR( rval, "Failed to create a new triangle" ); 02020 _newTriangles.insert( newTriangle ); 02021 if( debug_splits ) print_debug_triangle( newTriangle ); 02022 rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle ); // the second triangle 02023 MBERRORR( rval, "Failed to create a new triangle" ); 02024 _newTriangles.insert( newTriangle ); 02025 if( debug_splits ) print_debug_triangle( newTriangle ); 02026 rval = _mbImpl->create_element( MBTRI, conn + 6, 3, newTriangle ); // the second triangle 02027 MBERRORR( rval, "Failed to create a new triangle" ); 02028 _newTriangles.insert( newTriangle ); 02029 if( debug_splits ) print_debug_triangle( newTriangle ); 02030 } 02031 02032 return MB_SUCCESS; 02033 } 02034 02035 // build the list of intx points and entities from database involved 02036 // vertices, edges, triangles 02037 // it could be just a list of vertices (easiest case to handle after) 02038 02039 ErrorCode FBEngine::compute_intersection_points( EntityHandle&, EntityHandle from, EntityHandle to, CartVect& Dir, 02040 std::vector< CartVect >& points, std::vector< EntityHandle >& entities, 02041 std::vector< EntityHandle >& triangles ) 02042 { 02043 // keep a stack of triangles to process, and do not add those already processed 02044 // either mark them, or maybe keep them in a local set? 02045 // from and to are now nodes, start from them 02046 CartVect p1, p2; // the position of from and to 02047 ErrorCode rval = _mbImpl->get_coords( &from, 1, (double*)&p1 ); 02048 MBERRORR( rval, "failed to get 'from' coordinates" ); 02049 rval = _mbImpl->get_coords( &to, 1, (double*)&p2 ); 02050 MBERRORR( rval, "failed to get 'from' coordinates" ); 02051 02052 CartVect vect( p2 - p1 ); 02053 double dist2 = vect.length(); 02054 if( dist2 < tolerance_segment ) 02055 { 02056 // we are done, return 02057 return MB_SUCCESS; 02058 } 02059 CartVect normPlane = Dir * vect; 02060 normPlane.normalize(); 02061 std::set< EntityHandle > visitedTriangles; 02062 CartVect currentPoint = p1; 02063 // push the current point if it is empty 02064 if( points.size() == 0 ) 02065 { 02066 points.push_back( p1 ); 02067 entities.push_back( from ); // this is a node now 02068 } 02069 02070 // these will be used in many loops 02071 CartVect intx = p1; // somewhere to start 02072 double param = -1.; 02073 02074 // first intersection 02075 EntityHandle currentBoundary = from; // it is a node, in the beginning 02076 02077 vect = p2 - currentPoint; 02078 while( vect.length() > 0. ) 02079 { 02080 // advance towards "to" node, from boundary handle 02081 EntityType etype = _mbImpl->type_from_handle( currentBoundary ); 02082 // if vertex, look for other triangles connected which intersect our plane (defined by p1, 02083 // p2, dir) 02084 std::vector< EntityHandle > adj_tri; 02085 rval = _mbImpl->get_adjacencies( ¤tBoundary, 1, 2, false, adj_tri ); 02086 unsigned int j = 0; 02087 EntityHandle tri; 02088 for( ; j < adj_tri.size(); j++ ) 02089 { 02090 tri = adj_tri[j]; 02091 if( visitedTriangles.find( tri ) != visitedTriangles.end() ) 02092 continue; // get another triangle, this one was already visited 02093 // check if it is one of the triangles that was pierced already 02094 if( _piercedTriangles.find( tri ) != _piercedTriangles.end() ) continue; 02095 // if vertex, look for opposite edge 02096 // if edge, look for 2 opposite edges 02097 // get vertices 02098 int nnodes; 02099 const EntityHandle* conn3; 02100 rval = _mbImpl->get_connectivity( tri, conn3, nnodes ); 02101 MBERRORR( rval, "Failed to get connectivity" ); 02102 // if one of the nodes is to, stop right there 02103 { 02104 if( conn3[0] == to || conn3[1] == to || conn3[2] == to ) 02105 { 02106 visitedTriangles.insert( tri ); 02107 triangles.push_back( tri ); 02108 currentPoint = p2; 02109 points.push_back( p2 ); 02110 entities.push_back( to ); // we are done 02111 break; // this is break from for loop, we still need to get out of while 02112 // we will get out, because vect will become 0, (p2-p2) 02113 } 02114 } 02115 EntityHandle nn2[2]; 02116 if( MBVERTEX == etype ) 02117 { 02118 nn2[0] = conn3[0]; 02119 nn2[1] = conn3[1]; 02120 if( nn2[0] == currentBoundary ) nn2[0] = conn3[2]; 02121 if( nn2[1] == currentBoundary ) nn2[1] = conn3[2]; 02122 // get coordinates 02123 CartVect Pt[2]; 02124 02125 rval = _mbImpl->get_coords( nn2, 2, (double*)&Pt[0] ); 02126 MBERRORR( rval, "Failed to get coordinates" ); 02127 // see the intersection 02128 if( intersect_segment_and_plane_slice( Pt[0], Pt[1], currentPoint, p2, Dir, normPlane, intx, param ) ) 02129 { 02130 // we should stop for loop, and decide if it is edge or vertex 02131 if( param == 0. ) 02132 currentBoundary = nn2[0]; 02133 else 02134 { 02135 if( param == 1. ) 02136 currentBoundary = nn2[1]; 02137 else // param between 0 and 1, so edge 02138 { 02139 // find the edge between vertices 02140 std::vector< EntityHandle > edges1; 02141 // change the create flag to true, because that edge must exist in 02142 // current triangle if we want to advance; nn2 are 2 nodes in current 02143 // triangle!! 02144 rval = _mbImpl->get_adjacencies( nn2, 2, 1, true, edges1, Interface::INTERSECT ); 02145 MBERRORR( rval, "Failed to get edges" ); 02146 if( edges1.size() != 1 ) MBERRORR( MB_FAILURE, "Failed to get adjacent edges to 2 nodes" ); 02147 currentBoundary = edges1[0]; 02148 } 02149 } 02150 visitedTriangles.insert( tri ); 02151 currentPoint = intx; 02152 points.push_back( intx ); 02153 entities.push_back( currentBoundary ); 02154 triangles.push_back( tri ); 02155 if( debug_splits ) 02156 std::cout << "vtx new tri : " << _mbImpl->id_from_handle( tri ) 02157 << " type bdy:" << _mbImpl->type_from_handle( currentBoundary ) << "\n"; 02158 break; // out of for loop over triangles 02159 } 02160 } 02161 else // this is MBEDGE, we have the other triangle to try out 02162 { 02163 // first find the nodes from existing boundary 02164 int nnodes2; 02165 const EntityHandle* conn2; 02166 rval = _mbImpl->get_connectivity( currentBoundary, conn2, nnodes2 ); 02167 MBERRORR( rval, "Failed to get connectivity" ); 02168 int thirdIndex = -1; 02169 for( int tj = 0; tj < 3; tj++ ) 02170 { 02171 if( ( conn3[tj] != conn2[0] ) && ( conn3[tj] != conn2[1] ) ) 02172 { 02173 thirdIndex = tj; 02174 break; 02175 } 02176 } 02177 if( -1 == thirdIndex ) MBERRORR( MB_FAILURE, " can't get third node" ); 02178 CartVect Pt[3]; 02179 rval = _mbImpl->get_coords( conn3, 3, (double*)&Pt[0] ); 02180 MBERRORR( rval, "Failed to get coords" ); 02181 int indexFirst = ( thirdIndex + 1 ) % 3; 02182 int indexSecond = ( thirdIndex + 2 ) % 3; 02183 int index[2] = { indexFirst, indexSecond }; 02184 for( int k = 0; k < 2; k++ ) 02185 { 02186 nn2[0] = conn3[index[k]], nn2[1] = conn3[thirdIndex]; 02187 if( intersect_segment_and_plane_slice( Pt[index[k]], Pt[thirdIndex], currentPoint, p2, Dir, 02188 normPlane, intx, param ) ) 02189 { 02190 // we should stop for loop, and decide if it is edge or vertex 02191 if( param == 0. ) 02192 currentBoundary = conn3[index[k]]; // it is not really possible 02193 else 02194 { 02195 if( param == 1. ) 02196 currentBoundary = conn3[thirdIndex]; 02197 else // param between 0 and 1, so edge is fine 02198 { 02199 // find the edge between vertices 02200 std::vector< EntityHandle > edges1; 02201 // change the create flag to true, because that edge must exist in 02202 // current triangle if we want to advance; nn2 are 2 nodes in 02203 // current triangle!! 02204 rval = _mbImpl->get_adjacencies( nn2, 2, 1, true, edges1, Interface::INTERSECT ); 02205 MBERRORR( rval, "Failed to get edges" ); 02206 if( edges1.size() != 1 ) 02207 MBERRORR( MB_FAILURE, "Failed to get adjacent edges to 2 nodes" ); 02208 currentBoundary = edges1[0]; 02209 } 02210 } 02211 visitedTriangles.insert( tri ); 02212 currentPoint = intx; 02213 points.push_back( intx ); 02214 entities.push_back( currentBoundary ); 02215 triangles.push_back( tri ); 02216 if( debug_splits ) 02217 { 02218 std::cout << "edge new tri : " << _mbImpl->id_from_handle( tri ) 02219 << " type bdy: " << _mbImpl->type_from_handle( currentBoundary ) 02220 << " id: " << _mbImpl->id_from_handle( currentBoundary ) << "\n"; 02221 _mbImpl->list_entity( currentBoundary ); 02222 } 02223 break; // out of for loop over triangles 02224 } 02225 // we should not reach here 02226 } 02227 } 02228 } 02229 /*if (j==adj_tri.size()) 02230 { 02231 MBERRORR(MB_FAILURE, "did not advance"); 02232 }*/ 02233 vect = p2 - currentPoint; 02234 } 02235 02236 if( debug_splits ) 02237 std::cout << "nb entities: " << entities.size() << " triangles:" << triangles.size() 02238 << " points.size(): " << points.size() << "\n"; 02239 02240 return MB_SUCCESS; 02241 } 02242 02243 ErrorCode FBEngine::split_edge_at_point( EntityHandle edge, CartVect& point, EntityHandle& new_edge ) 02244 { 02245 // return MB_NOT_IMPLEMENTED; 02246 // first, we need to find the closest point on the smooth edge, then 02247 // split the smooth edge, then call the split_edge_at_mesh_node 02248 // or maybe just find the closest node?? 02249 // first of all, we need to find a point on the smooth edge, close by 02250 // then split the mesh edge (if needed) 02251 // this would be quite a drama, as the new edge has to be inserted in 02252 // order for proper geo edge definition 02253 02254 // first of all, find a closest point 02255 // usually called for 02256 if( debug_splits ) 02257 { 02258 std::cout << "Split edge " << _mbImpl->id_from_handle( edge ) << " at point:" << point << "\n"; 02259 } 02260 int dim = _my_geomTopoTool->dimension( edge ); 02261 if( dim != 1 ) return MB_FAILURE; 02262 if( !_smooth ) return MB_FAILURE; // call it only for smooth option... 02263 // maybe we should do something for linear too 02264 02265 SmoothCurve* curve = this->_edges[edge]; 02266 EntityHandle closeNode; 02267 int edgeIndex; 02268 double u = curve->u_from_position( point[0], point[1], point[2], closeNode, edgeIndex ); 02269 if( 0 == closeNode ) 02270 { 02271 // we really need to split an existing edge 02272 // do not do that yet 02273 // evaluate from u: 02274 /*double pos[3]; 02275 curve->position_from_u(u, pos[0], pos[1], pos[2] );*/ 02276 // create a new node here, and split one edge 02277 // change also connectivity, create new triangles on both sides ... 02278 std::cout << "not found a close node, u is: " << u << " edge index: " << edgeIndex << "\n"; 02279 return MB_FAILURE; // not implemented yet 02280 } 02281 return split_edge_at_mesh_node( edge, closeNode, new_edge ); 02282 } 02283 02284 ErrorCode FBEngine::split_edge_at_mesh_node( EntityHandle edge, EntityHandle node, EntityHandle& new_edge ) 02285 { 02286 // the node should be in the list of nodes 02287 02288 int dim = _my_geomTopoTool->dimension( edge ); 02289 if( dim != 1 ) return MB_FAILURE; // not an edge 02290 02291 if( debug_splits ) 02292 { 02293 std::cout << "Split edge " << _mbImpl->id_from_handle( edge ) 02294 << " with global id: " << _my_geomTopoTool->global_id( edge ) 02295 << " at node:" << _mbImpl->id_from_handle( node ) << "\n"; 02296 } 02297 02298 // now get the edges 02299 // the order is important... 02300 // these are ordered sets !! 02301 std::vector< EntityHandle > ents; 02302 ErrorCode rval = _mbImpl->get_entities_by_type( edge, MBEDGE, ents ); 02303 if( MB_SUCCESS != rval ) return rval; 02304 if( ents.size() < 1 ) return MB_FAILURE; // no edges 02305 02306 const EntityHandle* conn = NULL; 02307 int len; 02308 // find the edge connected to the splitting node 02309 int num_mesh_edges = (int)ents.size(); 02310 int index_edge; 02311 EntityHandle firstNode = 0; 02312 for( index_edge = 0; index_edge < num_mesh_edges; index_edge++ ) 02313 { 02314 rval = MBI->get_connectivity( ents[index_edge], conn, len ); 02315 if( MB_SUCCESS != rval ) return rval; 02316 if( index_edge == 0 ) firstNode = conn[0]; // will be used to decide vertex sets adjacencies 02317 if( conn[0] == node ) 02318 { 02319 if( index_edge == 0 ) 02320 { 02321 new_edge = 0; // no need to split, it is the first node 02322 return MB_SUCCESS; // no split 02323 } 02324 else 02325 return MB_FAILURE; // we should have found the node already , wrong 02326 // connectivity 02327 } 02328 else if( conn[1] == node ) 02329 { 02330 // we found the index of interest 02331 break; 02332 } 02333 } 02334 if( index_edge == num_mesh_edges - 1 ) 02335 { 02336 new_edge = 0; // no need to split, it is the last node 02337 return MB_SUCCESS; // no split 02338 } 02339 02340 // here, we have 0 ... index_edge edges in the first set, 02341 // create a vertex set and add the node to it 02342 02343 if( debug_splits ) 02344 { 02345 std::cout << "Split edge with " << num_mesh_edges << " mesh edges, at index (0 based) " << index_edge << "\n"; 02346 } 02347 02348 // at this moment, the node set should have been already created in new_geo_edge 02349 EntityHandle nodeSet; // the node set that has the node (find it!) 02350 bool found = find_vertex_set_for_node( node, nodeSet ); 02351 02352 if( !found ) 02353 { 02354 // create a node set and add the node 02355 02356 // must be an error, but create one nevertheless 02357 rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet ); 02358 MBERRORR( rval, "Failed to create a new vertex set" ); 02359 02360 rval = _mbImpl->add_entities( nodeSet, &node, 1 ); 02361 MBERRORR( rval, "Failed to add the node to the set" ); 02362 02363 rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 ); // 02364 MBERRORR( rval, "Failed to commit the node set" ); 02365 02366 if( debug_splits ) 02367 { 02368 std::cout << " create a vertex set (this should have been created before!)" 02369 << _mbImpl->id_from_handle( nodeSet ) 02370 << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << "\n"; 02371 } 02372 } 02373 02374 // we need to remove the remaining mesh edges from first set, and add it 02375 // to the second set, in order 02376 02377 rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_edge ); 02378 MBERRORR( rval, "can't create geo edge" ); 02379 02380 int remaining = num_mesh_edges - 1 - index_edge; 02381 02382 // add edges to the edge set 02383 rval = _mbImpl->add_entities( new_edge, &ents[index_edge + 1], remaining ); 02384 MBERRORR( rval, "can't add edges to the new edge" ); 02385 02386 // also, remove the second node set from old edge 02387 // remove the edges index_edge+1 and up 02388 02389 rval = _mbImpl->remove_entities( edge, &ents[index_edge + 1], remaining ); 02390 MBERRORR( rval, "can't remove edges from the old edge" ); 02391 02392 // need to find the adjacent vertex sets 02393 Range vertexRange; 02394 rval = getAdjacentEntities( edge, 0, vertexRange ); 02395 02396 EntityHandle secondSet; 02397 if( vertexRange.size() == 1 ) 02398 { 02399 // initially a periodic edge, OK to add the new set to both edges, and the 02400 // second set 02401 secondSet = vertexRange[0]; 02402 } 02403 else 02404 { 02405 if( vertexRange.size() > 2 ) return MB_FAILURE; // something must be wrong with too many vertices 02406 // find first node 02407 int k; 02408 for( k = 0; k < 2; k++ ) 02409 { 02410 Range verts; 02411 rval = _mbImpl->get_entities_by_type( vertexRange[k], MBVERTEX, verts ); 02412 02413 MBERRORR( rval, "can't get vertices from vertex set" ); 02414 if( verts.size() != 1 ) MBERRORR( MB_FAILURE, " node set not defined well" ); 02415 if( firstNode == verts[0] ) 02416 { 02417 secondSet = vertexRange[1 - k]; // the other set; it is 1 or 0 02418 break; 02419 } 02420 } 02421 if( k >= 2 ) 02422 { 02423 // it means we didn't find the right node set 02424 MBERRORR( MB_FAILURE, " can't find the right vertex" ); 02425 } 02426 // remove the second set from the connectivity with the 02427 // edge (parent-child relation) 02428 // remove_parent_child( EntityHandle parent, 02429 // EntityHandle child ) 02430 rval = _mbImpl->remove_parent_child( edge, secondSet ); 02431 MBERRORR( rval, " can't remove the second vertex from edge" ); 02432 } 02433 // at this point, we just need to add to both edges the new node set (vertex) 02434 rval = _mbImpl->add_parent_child( edge, nodeSet ); 02435 MBERRORR( rval, " can't add new vertex to old edge" ); 02436 02437 rval = _mbImpl->add_parent_child( new_edge, nodeSet ); 02438 MBERRORR( rval, " can't add new vertex to new edge" ); 02439 02440 // one thing that I forgot: add the secondSet as a child to new edge!!! 02441 // (so far, the new edge has only one end vertex!) 02442 rval = _mbImpl->add_parent_child( new_edge, secondSet ); 02443 MBERRORR( rval, " can't add second vertex to new edge" ); 02444 02445 // now, add the edge and vertex to geo tool 02446 02447 rval = _my_geomTopoTool->add_geo_set( new_edge, 1 ); 02448 MBERRORR( rval, " can't add new edge" ); 02449 02450 // next, get the adjacent faces to initial edge, and add them as parents 02451 // to the new edge 02452 02453 // need to find the adjacent face sets 02454 Range faceRange; 02455 rval = getAdjacentEntities( edge, 2, faceRange ); 02456 02457 // these faces will be adjacent to the new edge too! 02458 // need to set the senses of new edge within faces 02459 02460 for( Range::iterator it = faceRange.begin(); it != faceRange.end(); ++it ) 02461 { 02462 EntityHandle face = *it; 02463 rval = _mbImpl->add_parent_child( face, new_edge ); 02464 MBERRORR( rval, " can't add new edge - face parent relation" ); 02465 int sense; 02466 rval = _my_geomTopoTool->get_sense( edge, face, sense ); 02467 MBERRORR( rval, " can't get initial sense of edge in the adjacent face" ); 02468 // keep the same sense for the new edge within the faces 02469 rval = _my_geomTopoTool->set_sense( new_edge, face, sense ); 02470 MBERRORR( rval, " can't set sense of new edge in the adjacent face" ); 02471 } 02472 02473 return MB_SUCCESS; 02474 } 02475 02476 ErrorCode FBEngine::split_bedge_at_new_mesh_node( EntityHandle edge, EntityHandle node, EntityHandle brokenEdge, 02477 EntityHandle& new_edge ) 02478 { 02479 // the node should be in the list of nodes 02480 02481 int dim = _my_geomTopoTool->dimension( edge ); 02482 if( dim != 1 ) return MB_FAILURE; // not an edge 02483 02484 if( debug_splits ) 02485 { 02486 std::cout << "Split edge " << _mbImpl->id_from_handle( edge ) 02487 << " with global id: " << _my_geomTopoTool->global_id( edge ) 02488 << " at new node:" << _mbImpl->id_from_handle( node ) << "breaking mesh edge" 02489 << _mbImpl->id_from_handle( brokenEdge ) << "\n"; 02490 } 02491 02492 // now get the edges 02493 // the order is important... 02494 // these are ordered sets !! 02495 std::vector< EntityHandle > ents; 02496 ErrorCode rval = _mbImpl->get_entities_by_type( edge, MBEDGE, ents ); 02497 if( MB_SUCCESS != rval ) return rval; 02498 if( ents.size() < 1 ) return MB_FAILURE; // no edges 02499 02500 const EntityHandle* conn = NULL; 02501 int len; 02502 // the edge connected to the splitting node is brokenEdge 02503 // find the small edges it is broken into, which are connected to 02504 // new vertex (node) and its ends; also, if necessary, reorient these small edges 02505 // for proper orientation 02506 rval = _mbImpl->get_connectivity( brokenEdge, conn, len ); 02507 MBERRORR( rval, "Failed to get connectivity of broken edge" ); 02508 02509 // we already created the new edges, make sure they are oriented fine; if not, revert them 02510 EntityHandle conn02[] = { conn[0], node }; // first node and new node 02511 // default, intersect 02512 std::vector< EntityHandle > adj_edges; 02513 rval = _mbImpl->get_adjacencies( conn02, 2, 1, false, adj_edges ); 02514 if( adj_edges.size() < 1 || rval != MB_SUCCESS ) MBERRORR( MB_FAILURE, " Can't find small split edge" ); 02515 02516 // get this edge connectivity; 02517 EntityHandle firstEdge = adj_edges[0]; 02518 const EntityHandle* connActual = NULL; 02519 rval = _mbImpl->get_connectivity( firstEdge, connActual, len ); 02520 MBERRORR( rval, "Failed to get connectivity of first split edge" ); 02521 // if it is the same as conn02, we are happy 02522 if( conn02[0] != connActual[0] ) 02523 { 02524 // reset connectivity of edge 02525 rval = _mbImpl->set_connectivity( firstEdge, conn02, 2 ); 02526 MBERRORR( rval, "Failed to reset connectivity of first split edge" ); 02527 } 02528 // now treat the second edge 02529 adj_edges.clear(); // 02530 EntityHandle conn21[] = { node, conn[1] }; // new node and second node 02531 rval = _mbImpl->get_adjacencies( conn21, 2, 1, false, adj_edges ); 02532 if( adj_edges.size() < 1 || rval != MB_SUCCESS ) MBERRORR( MB_FAILURE, " Can't find second small split edge" ); 02533 02534 // get this edge connectivity; 02535 EntityHandle secondEdge = adj_edges[0]; 02536 rval = _mbImpl->get_connectivity( firstEdge, connActual, len ); 02537 MBERRORR( rval, "Failed to get connectivity of first split edge" ); 02538 // if it is the same as conn21, we are happy 02539 if( conn21[0] != connActual[0] ) 02540 { 02541 // reset connectivity of edge 02542 rval = _mbImpl->set_connectivity( secondEdge, conn21, 2 ); 02543 MBERRORR( rval, "Failed to reset connectivity of first split edge" ); 02544 } 02545 02546 int num_mesh_edges = (int)ents.size(); 02547 int index_edge; // this is the index of the edge that will be removed (because it is split) 02548 // the rest of edges will be put in order in the (remaining) edge and new edge 02549 02550 for( index_edge = 0; index_edge < num_mesh_edges; index_edge++ ) 02551 if( brokenEdge == ents[index_edge] ) break; 02552 if( index_edge >= num_mesh_edges ) MBERRORR( MB_FAILURE, "can't find the broken edge" ); 02553 02554 // so the edges up to index_edge and firstEdge, will form the "old" edge 02555 // the edges secondEdge and from index_edge+1 to end will form the new_edge 02556 02557 // here, we have 0 ... index_edge edges in the first set, 02558 // create a vertex set and add the node to it 02559 02560 if( debug_splits ) 02561 { 02562 std::cout << "Split edge with " << num_mesh_edges << " mesh edges, at index (0 based) " << index_edge << "\n"; 02563 } 02564 02565 // at this moment, the node set should have been already created in new_geo_edge 02566 EntityHandle nodeSet; // the node set that has the node (find it!) 02567 bool found = find_vertex_set_for_node( node, nodeSet ); 02568 02569 if( !found ) 02570 { 02571 // create a node set and add the node 02572 02573 // must be an error, but create one nevertheless 02574 rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet ); 02575 MBERRORR( rval, "Failed to create a new vertex set" ); 02576 02577 rval = _mbImpl->add_entities( nodeSet, &node, 1 ); 02578 MBERRORR( rval, "Failed to add the node to the set" ); 02579 02580 rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 ); // 02581 MBERRORR( rval, "Failed to commit the node set" ); 02582 02583 if( debug_splits ) 02584 { 02585 std::cout << " create a vertex set (this should have been created before!)" 02586 << _mbImpl->id_from_handle( nodeSet ) 02587 << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << "\n"; 02588 } 02589 } 02590 02591 // we need to remove the remaining mesh edges from first set, and add it 02592 // to the second set, in order 02593 02594 rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_edge ); 02595 MBERRORR( rval, "can't create geo edge" ); 02596 02597 int remaining = num_mesh_edges - 1 - index_edge; 02598 02599 // add edges to the new edge set 02600 rval = _mbImpl->add_entities( new_edge, &secondEdge, 1 ); // add the second split edge to new edge 02601 MBERRORR( rval, "can't add second split edge to the new edge" ); 02602 // then add the rest 02603 rval = _mbImpl->add_entities( new_edge, &ents[index_edge + 1], remaining ); 02604 MBERRORR( rval, "can't add edges to the new edge" ); 02605 02606 // also, remove the second node set from old edge 02607 // remove the edges index_edge and up 02608 02609 rval = _mbImpl->remove_entities( edge, &ents[index_edge], remaining + 1 ); // include the 02610 MBERRORR( rval, "can't remove edges from the old edge" ); 02611 // add the firstEdge too 02612 rval = _mbImpl->add_entities( edge, &firstEdge, 1 ); // add the second split edge to new edge 02613 MBERRORR( rval, "can't add first split edge to the old edge" ); 02614 02615 // need to find the adjacent vertex sets 02616 Range vertexRange; 02617 rval = getAdjacentEntities( edge, 0, vertexRange ); 02618 02619 EntityHandle secondSet; 02620 if( vertexRange.size() == 1 ) 02621 { 02622 // initially a periodic edge, OK to add the new set to both edges, and the 02623 // second set 02624 secondSet = vertexRange[0]; 02625 } 02626 else 02627 { 02628 if( vertexRange.size() > 2 ) return MB_FAILURE; // something must be wrong with too many vertices 02629 // find first node 02630 EntityHandle firstNode; 02631 02632 rval = MBI->get_connectivity( ents[0], conn, len ); 02633 if( MB_SUCCESS != rval ) return rval; 02634 firstNode = conn[0]; // this is the first node of the initial edge 02635 // we will use it to identify the vertex set associated with it 02636 int k; 02637 for( k = 0; k < 2; k++ ) 02638 { 02639 Range verts; 02640 rval = _mbImpl->get_entities_by_type( vertexRange[k], MBVERTEX, verts ); 02641 02642 MBERRORR( rval, "can't get vertices from vertex set" ); 02643 if( verts.size() != 1 ) MBERRORR( MB_FAILURE, " node set not defined well" ); 02644 if( firstNode == verts[0] ) 02645 { 02646 secondSet = vertexRange[1 - k]; // the other set; it is 1 or 0 02647 break; 02648 } 02649 } 02650 if( k >= 2 ) 02651 { 02652 // it means we didn't find the right node set 02653 MBERRORR( MB_FAILURE, " can't find the right vertex" ); 02654 } 02655 // remove the second set from the connectivity with the 02656 // edge (parent-child relation) 02657 // remove_parent_child( EntityHandle parent, 02658 // EntityHandle child ) 02659 rval = _mbImpl->remove_parent_child( edge, secondSet ); 02660 MBERRORR( rval, " can't remove the second vertex from edge" ); 02661 } 02662 // at this point, we just need to add to both edges the new node set (vertex) 02663 rval = _mbImpl->add_parent_child( edge, nodeSet ); 02664 MBERRORR( rval, " can't add new vertex to old edge" ); 02665 02666 rval = _mbImpl->add_parent_child( new_edge, nodeSet ); 02667 MBERRORR( rval, " can't add new vertex to new edge" ); 02668 02669 // one thing that I forgot: add the secondSet as a child to new edge!!! 02670 // (so far, the new edge has only one end vertex!) 02671 rval = _mbImpl->add_parent_child( new_edge, secondSet ); 02672 MBERRORR( rval, " can't add second vertex to new edge" ); 02673 02674 // now, add the edge and vertex to geo tool 02675 02676 rval = _my_geomTopoTool->add_geo_set( new_edge, 1 ); 02677 MBERRORR( rval, " can't add new edge" ); 02678 02679 // next, get the adjacent faces to initial edge, and add them as parents 02680 // to the new edge 02681 02682 // need to find the adjacent face sets 02683 Range faceRange; 02684 rval = getAdjacentEntities( edge, 2, faceRange ); 02685 02686 // these faces will be adjacent to the new edge too! 02687 // need to set the senses of new edge within faces 02688 02689 for( Range::iterator it = faceRange.begin(); it != faceRange.end(); ++it ) 02690 { 02691 EntityHandle face = *it; 02692 rval = _mbImpl->add_parent_child( face, new_edge ); 02693 MBERRORR( rval, " can't add new edge - face parent relation" ); 02694 int sense; 02695 rval = _my_geomTopoTool->get_sense( edge, face, sense ); 02696 MBERRORR( rval, " can't get initial sense of edge in the adjacent face" ); 02697 // keep the same sense for the new edge within the faces 02698 rval = _my_geomTopoTool->set_sense( new_edge, face, sense ); 02699 MBERRORR( rval, " can't set sense of new edge in the adjacent face" ); 02700 } 02701 02702 return MB_SUCCESS; 02703 } 02704 02705 ErrorCode FBEngine::split_boundary( EntityHandle face, EntityHandle atNode ) 02706 { 02707 // find the boundary edges, and split the one that we find it is a part of 02708 if( debug_splits ) 02709 { 02710 std::cout << "Split face " << _mbImpl->id_from_handle( face ) 02711 << " at node:" << _mbImpl->id_from_handle( atNode ) << "\n"; 02712 } 02713 Range bound_edges; 02714 ErrorCode rval = getAdjacentEntities( face, 1, bound_edges ); 02715 MBERRORR( rval, " can't get boundary edges" ); 02716 bool brokEdge = _brokenEdges.find( atNode ) != _brokenEdges.end(); 02717 02718 for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it ) 02719 { 02720 EntityHandle b_edge = *it; 02721 // get all edges in range 02722 Range mesh_edges; 02723 rval = _mbImpl->get_entities_by_dimension( b_edge, 1, mesh_edges ); 02724 MBERRORR( rval, " can't get mesh edges" ); 02725 if( brokEdge ) 02726 { 02727 EntityHandle brokenEdge = _brokenEdges[atNode]; 02728 if( mesh_edges.find( brokenEdge ) != mesh_edges.end() ) 02729 { 02730 EntityHandle new_edge; 02731 rval = split_bedge_at_new_mesh_node( b_edge, atNode, brokenEdge, new_edge ); 02732 return rval; 02733 } 02734 } 02735 else 02736 { 02737 Range nodes; 02738 rval = _mbImpl->get_connectivity( mesh_edges, nodes ); 02739 MBERRORR( rval, " can't get nodes from mesh edges" ); 02740 02741 if( nodes.find( atNode ) != nodes.end() ) 02742 { 02743 // we found our boundary edge candidate 02744 EntityHandle new_edge; 02745 rval = split_edge_at_mesh_node( b_edge, atNode, new_edge ); 02746 return rval; 02747 } 02748 } 02749 } 02750 // if the node was not found in any "current" boundary, it broke an existing 02751 // boundary edge 02752 MBERRORR( MB_FAILURE, " we did not find an appropriate boundary edge" ); 02753 ; // 02754 } 02755 02756 bool FBEngine::find_vertex_set_for_node( EntityHandle iNode, EntityHandle& oVertexSet ) 02757 { 02758 bool found = false; 02759 Range vertex_sets; 02760 02761 const int zero = 0; 02762 const void* const zero_val[] = { &zero }; 02763 Tag geom_tag; 02764 ErrorCode rval = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag ); 02765 if( MB_SUCCESS != rval ) return false; 02766 rval = _mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom_tag, zero_val, 1, vertex_sets ); 02767 if( MB_SUCCESS != rval ) return false; 02768 // local _gsets, as we might have not updated the local lists 02769 // see if ends of geo edge generated is in a geo set 0 or not 02770 02771 for( Range::iterator vsit = vertex_sets.begin(); vsit != vertex_sets.end(); ++vsit ) 02772 { 02773 EntityHandle vset = *vsit; 02774 // is the node part of this set? 02775 if( _mbImpl->contains_entities( vset, &iNode, 1 ) ) 02776 { 02777 02778 found = true; 02779 oVertexSet = vset; 02780 break; 02781 } 02782 } 02783 return found; 02784 } 02785 ErrorCode FBEngine::redistribute_boundary_edges_to_faces( EntityHandle face, EntityHandle newFace, 02786 std::vector< EntityHandle >& chainedEdges ) 02787 { 02788 02789 // so far, original boundary edges are all parent/child relations for face 02790 // we should get them all, and see if they are truly adjacent to face or newFace 02791 // decide also on the orientation/sense with respect to the triangles 02792 Range r1; // range in old face 02793 Range r2; // range of tris in new face 02794 ErrorCode rval = _mbImpl->get_entities_by_dimension( face, 2, r1 ); 02795 MBERRORR( rval, " can't get triangles from old face" ); 02796 rval = _mbImpl->get_entities_by_dimension( newFace, 2, r2 ); 02797 MBERRORR( rval, " can't get triangles from new face" ); 02798 // right now, all boundary edges are children of face 02799 // we need to get them all, and verify if they indeed are adjacent to face 02800 Range children; 02801 rval = _mbImpl->get_child_meshsets( face, children ); // only direct children are of interest 02802 MBERRORR( rval, " can't get children sets from face" ); 02803 02804 for( Range::iterator it = children.begin(); it != children.end(); ++it ) 02805 { 02806 EntityHandle edge = *it; 02807 if( std::find( chainedEdges.begin(), chainedEdges.end(), edge ) != chainedEdges.end() ) 02808 continue; // we already set this one fine 02809 // get one mesh edge from the edge; we have to get all of them!! 02810 if( _my_geomTopoTool->dimension( edge ) != 1 ) continue; // not an edge 02811 Range mesh_edges; 02812 rval = _mbImpl->get_entities_by_handle( edge, mesh_edges ); 02813 MBERRORR( rval, " can't get mesh edges from edge" ); 02814 if( mesh_edges.empty() ) MBERRORR( MB_FAILURE, " no mesh edges" ); 02815 EntityHandle mesh_edge = mesh_edges[0]; // first one is enough 02816 // get its triangles; see which one is in r1 or r2; it should not be in both 02817 Range triangles; 02818 rval = _mbImpl->get_adjacencies( &mesh_edge, 1, 2, false, triangles ); 02819 MBERRORR( rval, " can't get adjacent triangles" ); 02820 Range intx1 = intersect( triangles, r1 ); 02821 Range intx2 = intersect( triangles, r2 ); 02822 if( !intx1.empty() && !intx2.empty() ) MBERRORR( MB_FAILURE, " at least one should be empty" ); 02823 02824 if( intx2.empty() ) 02825 { 02826 // it means it is in the range r1; the sense should be fine, no need to reset 02827 // the sense should have been fine, also 02828 continue; 02829 } 02830 // so it must be a triangle in r2; 02831 EntityHandle triangle = intx2[0]; // one triangle only 02832 // remove the edge from boundary of face, and add it to the boundary of newFace 02833 // remove_parent_child( EntityHandle parent, EntityHandle child ) 02834 rval = _mbImpl->remove_parent_child( face, edge ); 02835 MBERRORR( rval, " can't remove parent child relation for edge" ); 02836 // add to the new face 02837 rval = _mbImpl->add_parent_child( newFace, edge ); 02838 MBERRORR( rval, " can't add parent child relation for edge" ); 02839 02840 // set some sense, based on the sense of the mesh_edge in triangle 02841 int num1, sense, offset; 02842 rval = _mbImpl->side_number( triangle, mesh_edge, num1, sense, offset ); 02843 MBERRORR( rval, "mesh edge not adjacent to triangle" ); 02844 02845 rval = this->_my_geomTopoTool->set_sense( edge, newFace, sense ); 02846 MBERRORR( rval, "can't set proper sense of edge in face" ); 02847 } 02848 02849 return MB_SUCCESS; 02850 } 02851 02852 ErrorCode FBEngine::set_neumann_tags( EntityHandle face, EntityHandle newFace ) 02853 { 02854 // these are for debugging purposes only 02855 // check the initial tag, then 02856 Tag ntag; 02857 ErrorCode rval = _mbImpl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, ntag ); 02858 MBERRORR( rval, "can't get tag handle" ); 02859 // check the value for face 02860 int nval; 02861 rval = _mbImpl->tag_get_data( ntag, &face, 1, &nval ); 02862 if( MB_SUCCESS == rval ) { nval++; } 02863 else 02864 { 02865 nval = 1; 02866 rval = _mbImpl->tag_set_data( ntag, &face, 1, &nval ); 02867 MBERRORR( rval, "can't set tag" ); 02868 nval = 2; 02869 } 02870 rval = _mbImpl->tag_set_data( ntag, &newFace, 1, &nval ); 02871 MBERRORR( rval, "can't set tag" ); 02872 02873 return MB_SUCCESS; 02874 } 02875 02876 // split the quads if needed; it will create a new gtt, which will 02877 // contain triangles instead of quads 02878 ErrorCode FBEngine::split_quads() 02879 { 02880 // first see if we have any quads in the 2d faces 02881 // _my_gsets[2] is a range of surfaces (moab sets of dimension 2) 02882 int num_quads = 0; 02883 for( Range::iterator it = _my_gsets[2].begin(); it != _my_gsets[2].end(); ++it ) 02884 { 02885 EntityHandle surface = *it; 02886 int num_q = 0; 02887 _mbImpl->get_number_entities_by_type( surface, MBQUAD, num_q ); 02888 num_quads += num_q; 02889 } 02890 02891 if( num_quads == 0 ) return MB_SUCCESS; // nothing to do 02892 02893 GeomTopoTool* new_gtt = NULL; 02894 ErrorCode rval = _my_geomTopoTool->duplicate_model( new_gtt ); 02895 MBERRORR( rval, "can't duplicate model" ); 02896 if( this->_t_created ) delete _my_geomTopoTool; 02897 02898 _t_created = true; // this one is for sure created here, even if the original gtt was not 02899 02900 // if we were using smart pointers, we would decrease the reference to the _my_geomTopoTool, at 02901 // least 02902 _my_geomTopoTool = new_gtt; 02903 02904 // replace the _my_gsets with the new ones, from the new set 02905 _my_geomTopoTool->find_geomsets( _my_gsets ); 02906 02907 // we have duplicated now the model, and the quads are part of the new _my_gsets[2] 02908 // we will split them now, by creating triangles along the smallest diagonal 02909 // maybe we will come up with a better criteria, but for the time being, it should be fine. 02910 // what we will do: we will remove all quads from the surface sets, and split them 02911 02912 for( Range::iterator it2 = _my_gsets[2].begin(); it2 != _my_gsets[2].end(); ++it2 ) 02913 { 02914 EntityHandle surface = *it2; 02915 Range quads; 02916 rval = _mbImpl->get_entities_by_type( surface, MBQUAD, quads ); 02917 MBERRORR( rval, "can't get quads from the surface set" ); 02918 rval = _mbImpl->remove_entities( surface, quads ); 02919 MBERRORR( rval, "can't remove quads from the surface set" ); // they are not deleted, just 02920 // removed from the set 02921 for( Range::iterator it = quads.begin(); it != quads.end(); ++it ) 02922 { 02923 EntityHandle quad = *it; 02924 int nnodes; 02925 const EntityHandle* conn; 02926 rval = _mbImpl->get_connectivity( quad, conn, nnodes ); 02927 MBERRORR( rval, "can't get quad connectivity" ); 02928 // get all vertices position, to see which one is the shortest diagonal 02929 CartVect pos[4]; 02930 rval = _mbImpl->get_coords( conn, 4, (double*)&pos[0] ); 02931 MBERRORR( rval, "can't get node coordinates" ); 02932 bool diag1 = ( ( pos[2] - pos[0] ).length_squared() < ( pos[3] - pos[1] ).length_squared() ); 02933 EntityHandle newTris[2]; 02934 EntityHandle tri1[3] = { conn[0], conn[1], conn[2] }; 02935 EntityHandle tri2[3] = { conn[0], conn[2], conn[3] }; 02936 if( !diag1 ) 02937 { 02938 tri1[2] = conn[3]; 02939 tri2[0] = conn[1]; 02940 } 02941 rval = _mbImpl->create_element( MBTRI, tri1, 3, newTris[0] ); 02942 MBERRORR( rval, "can't create triangle 1" ); 02943 rval = _mbImpl->create_element( MBTRI, tri2, 3, newTris[1] ); 02944 MBERRORR( rval, "can't create triangle 2" ); 02945 rval = _mbImpl->add_entities( surface, newTris, 2 ); 02946 MBERRORR( rval, "can't add triangles to the set" ); 02947 } 02948 // 02949 } 02950 return MB_SUCCESS; 02951 } 02952 ErrorCode FBEngine::boundary_mesh_edges_on_face( EntityHandle face, Range& boundary_mesh_edges ) 02953 { 02954 // this list is used only for finding the intersecting mesh edge for starting the 02955 // polygonal cutting line at boundary (if !closed) 02956 Range bound_edges; 02957 ErrorCode rval = getAdjacentEntities( face, 1, bound_edges ); 02958 MBERRORR( rval, " can't get boundary edges" ); 02959 for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it ) 02960 { 02961 EntityHandle b_edge = *it; 02962 // get all edges in range 02963 // Range mesh_edges; 02964 rval = _mbImpl->get_entities_by_dimension( b_edge, 1, boundary_mesh_edges ); 02965 MBERRORR( rval, " can't get mesh edges" ); 02966 } 02967 return MB_SUCCESS; 02968 } 02969 ErrorCode FBEngine::boundary_nodes_on_face( EntityHandle face, std::vector< EntityHandle >& boundary_nodes ) 02970 { 02971 // even if we repeat some nodes, it is OK 02972 // this list is used only for finding the closest boundary node for starting the 02973 // polygonal cutting line at boundary (if !closed) 02974 Range bound_edges; 02975 ErrorCode rval = getAdjacentEntities( face, 1, bound_edges ); 02976 MBERRORR( rval, " can't get boundary edges" ); 02977 Range b_nodes; 02978 for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it ) 02979 { 02980 EntityHandle b_edge = *it; 02981 // get all edges in range 02982 Range mesh_edges; 02983 rval = _mbImpl->get_entities_by_dimension( b_edge, 1, mesh_edges ); 02984 MBERRORR( rval, " can't get mesh edges" ); 02985 rval = _mbImpl->get_connectivity( mesh_edges, b_nodes ); 02986 MBERRORR( rval, " can't get nodes from mesh edges" ); 02987 } 02988 // create now a vector based on Range of nodes 02989 boundary_nodes.resize( b_nodes.size() ); 02990 std::copy( b_nodes.begin(), b_nodes.end(), boundary_nodes.begin() ); 02991 return MB_SUCCESS; 02992 } 02993 // used for splitting an edge 02994 ErrorCode FBEngine::split_internal_edge( EntityHandle& edge, EntityHandle& newVertex ) 02995 { 02996 // split the edge, and form 4 new triangles and 2 edges 02997 // get 2 triangles adjacent to edge 02998 Range adj_tri; 02999 ErrorCode rval = _mbImpl->get_adjacencies( &edge, 1, 2, false, adj_tri ); 03000 MBERRORR( rval, "can't get adj_tris" ); 03001 adj_tri = subtract( adj_tri, _piercedTriangles ); 03002 if( adj_tri.size() >= 3 ) { std::cout << "WARNING: non manifold geometry. Are you sure?"; } 03003 for( Range::iterator it = adj_tri.begin(); it != adj_tri.end(); ++it ) 03004 { 03005 EntityHandle tri = *it; 03006 _piercedTriangles.insert( tri ); 03007 const EntityHandle* conn3; 03008 int nnodes; 03009 rval = _mbImpl->get_connectivity( tri, conn3, nnodes ); 03010 MBERRORR( rval, "can't get nodes" ); 03011 int num1, sense, offset; 03012 rval = _mbImpl->side_number( tri, edge, num1, sense, offset ); 03013 MBERRORR( rval, "can't get side number" ); 03014 // after we know the side number, we can split in 2 triangles 03015 // node i is opposite to edge i 03016 int num2 = ( num1 + 1 ) % 3; 03017 int num3 = ( num2 + 1 ) % 3; 03018 // the edge from num1 to num2 is split into 2 edges 03019 EntityHandle t1[] = { conn3[num2], conn3[num3], newVertex }; 03020 EntityHandle t2[] = { conn3[num1], newVertex, conn3[num3] }; 03021 EntityHandle newTriangle, newTriangle2; 03022 rval = _mbImpl->create_element( MBTRI, t1, 3, newTriangle ); 03023 MBERRORR( rval, "can't create triangle" ); 03024 _newTriangles.insert( newTriangle ); 03025 rval = _mbImpl->create_element( MBTRI, t2, 3, newTriangle2 ); 03026 MBERRORR( rval, "can't create triangle" ); 03027 _newTriangles.insert( newTriangle2 ); 03028 // create edges with this, indirectly 03029 std::vector< EntityHandle > edges0; 03030 rval = _mbImpl->get_adjacencies( &newTriangle, 1, 1, true, edges0 ); 03031 MBERRORR( rval, "can't get new edges" ); 03032 edges0.clear(); 03033 rval = _mbImpl->get_adjacencies( &newTriangle2, 1, 1, true, edges0 ); 03034 MBERRORR( rval, "can't get new edges" ); 03035 if( debug_splits ) 03036 { 03037 std::cout << "2 (out of 4) triangles formed:\n"; 03038 _mbImpl->list_entity( newTriangle ); 03039 print_debug_triangle( newTriangle ); 03040 _mbImpl->list_entity( newTriangle2 ); 03041 print_debug_triangle( newTriangle2 ); 03042 } 03043 } 03044 return MB_SUCCESS; 03045 } 03046 // triangle split 03047 ErrorCode FBEngine::divide_triangle( EntityHandle triangle, EntityHandle& newVertex ) 03048 { 03049 // 03050 _piercedTriangles.insert( triangle ); 03051 int nnodes = 0; 03052 const EntityHandle* conn3 = NULL; 03053 ErrorCode rval = _mbImpl->get_connectivity( triangle, conn3, nnodes ); 03054 MBERRORR( rval, "can't get nodes" ); 03055 EntityHandle t1[] = { conn3[0], conn3[1], newVertex }; 03056 EntityHandle t2[] = { conn3[1], conn3[2], newVertex }; 03057 EntityHandle t3[] = { conn3[2], conn3[0], newVertex }; 03058 EntityHandle newTriangle, newTriangle2, newTriangle3; 03059 rval = _mbImpl->create_element( MBTRI, t1, 3, newTriangle ); 03060 MBERRORR( rval, "can't create triangle" ); 03061 _newTriangles.insert( newTriangle ); 03062 rval = _mbImpl->create_element( MBTRI, t2, 3, newTriangle3 ); 03063 MBERRORR( rval, "can't create triangle" ); 03064 _newTriangles.insert( newTriangle3 ); 03065 rval = _mbImpl->create_element( MBTRI, t3, 3, newTriangle2 ); 03066 MBERRORR( rval, "can't create triangle" ); 03067 _newTriangles.insert( newTriangle2 ); 03068 03069 // create all edges 03070 std::vector< EntityHandle > edges0; 03071 rval = _mbImpl->get_adjacencies( &newTriangle, 1, 1, true, edges0 ); 03072 MBERRORR( rval, "can't get new edges" ); 03073 edges0.clear(); 03074 rval = _mbImpl->get_adjacencies( &newTriangle2, 1, 1, true, edges0 ); 03075 MBERRORR( rval, "can't get new edges" ); 03076 if( debug_splits ) 03077 { 03078 std::cout << "3 triangles formed:\n"; 03079 _mbImpl->list_entity( newTriangle ); 03080 print_debug_triangle( newTriangle ); 03081 _mbImpl->list_entity( newTriangle3 ); 03082 print_debug_triangle( newTriangle3 ); 03083 _mbImpl->list_entity( newTriangle2 ); 03084 print_debug_triangle( newTriangle2 ); 03085 std::cout << "original nodes in tri:\n"; 03086 _mbImpl->list_entity( conn3[0] ); 03087 _mbImpl->list_entity( conn3[1] ); 03088 _mbImpl->list_entity( conn3[2] ); 03089 } 03090 return MB_SUCCESS; 03091 } 03092 03093 ErrorCode FBEngine::create_volume_with_direction( EntityHandle newFace1, EntityHandle newFace2, double* direction, 03094 EntityHandle& volume ) 03095 { 03096 03097 // MESHSET 03098 // ErrorCode rval = _mbImpl->create_meshset(MESHSET_ORDERED, new_geo_edge); 03099 ErrorCode rval = _mbImpl->create_meshset( MESHSET_SET, volume ); 03100 MBERRORR( rval, "can't create volume" ); 03101 03102 int volumeMatId = 1; // just give a mat id, for debugging, mostly 03103 Tag matTag; 03104 rval = _mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, matTag ); 03105 MBERRORR( rval, "can't get material tag" ); 03106 03107 rval = _mbImpl->tag_set_data( matTag, &volume, 1, &volumeMatId ); 03108 MBERRORR( rval, "can't set material tag value on volume" ); 03109 03110 // get the edges of those 2 faces, and get the vertices of those edges 03111 // in order, they should be created in the same order (?); check for that, anyway 03112 rval = _mbImpl->add_parent_child( volume, newFace1 ); 03113 MBERRORR( rval, "can't add first face to volume" ); 03114 03115 rval = _mbImpl->add_parent_child( volume, newFace2 ); 03116 MBERRORR( rval, "can't add second face to volume" ); 03117 03118 // first is bottom, so it is negatively oriented 03119 rval = _my_geomTopoTool->add_geo_set( volume, 3 ); 03120 MBERRORR( rval, "can't add volume to the gtt" ); 03121 03122 // set senses 03123 // bottom face is negatively oriented, its normal is toward interior of the volume 03124 rval = _my_geomTopoTool->set_sense( newFace1, volume, -1 ); 03125 MBERRORR( rval, "can't set bottom face sense to the volume" ); 03126 03127 // the top face is positively oriented 03128 rval = _my_geomTopoTool->set_sense( newFace2, volume, 1 ); 03129 MBERRORR( rval, "can't set top face sense to the volume" ); 03130 03131 // the children should be in the same direction 03132 // get the side edges of each face, and form lateral faces, along direction 03133 std::vector< EntityHandle > edges1; 03134 std::vector< EntityHandle > edges2; 03135 03136 rval = _mbImpl->get_child_meshsets( newFace1, edges1 ); // no hops 03137 MBERRORR( rval, "can't get children edges or first face, bottom" ); 03138 03139 rval = _mbImpl->get_child_meshsets( newFace2, edges2 ); // no hops 03140 MBERRORR( rval, "can't get children edges for second face, top" ); 03141 03142 if( edges1.size() != edges2.size() ) MBERRORR( MB_FAILURE, "wrong correspondence " ); 03143 03144 for( unsigned int i = 0; i < edges1.size(); ++i ) 03145 { 03146 EntityHandle newLatFace; 03147 rval = weave_lateral_face_from_edges( edges1[i], edges2[i], direction, newLatFace ); 03148 MBERRORR( rval, "can't weave lateral face" ); 03149 rval = _mbImpl->add_parent_child( volume, newLatFace ); 03150 MBERRORR( rval, "can't add lateral face to volume" ); 03151 03152 // set sense as positive 03153 rval = _my_geomTopoTool->set_sense( newLatFace, volume, 1 ); 03154 MBERRORR( rval, "can't set lateral face sense to the volume" ); 03155 } 03156 03157 rval = set_default_neumann_tags(); 03158 MBERRORR( rval, "can't set new neumann tags" ); 03159 03160 return MB_SUCCESS; 03161 } 03162 03163 ErrorCode FBEngine::get_nodes_from_edge( EntityHandle gedge, std::vector< EntityHandle >& nodes ) 03164 { 03165 std::vector< EntityHandle > ents; 03166 ErrorCode rval = _mbImpl->get_entities_by_type( gedge, MBEDGE, ents ); 03167 if( MB_SUCCESS != rval ) return rval; 03168 if( ents.size() < 1 ) return MB_FAILURE; 03169 03170 nodes.resize( ents.size() + 1 ); 03171 const EntityHandle* conn = NULL; 03172 int len; 03173 for( unsigned int i = 0; i < ents.size(); ++i ) 03174 { 03175 rval = _mbImpl->get_connectivity( ents[i], conn, len ); 03176 MBERRORR( rval, "can't get edge connectivity" ); 03177 nodes[i] = conn[0]; 03178 } 03179 // the last one is conn[1] 03180 nodes[ents.size()] = conn[1]; 03181 return MB_SUCCESS; 03182 } 03183 ErrorCode FBEngine::weave_lateral_face_from_edges( EntityHandle bEdge, EntityHandle tEdge, double* direction, 03184 EntityHandle& newLatFace ) 03185 { 03186 // in weird cases might need to create new vertices in the interior; 03187 // in most cases, it is OK 03188 03189 ErrorCode rval = _mbImpl->create_meshset( MESHSET_SET, newLatFace ); 03190 MBERRORR( rval, "can't create new lateral face" ); 03191 03192 EntityHandle v[4]; // vertex sets 03193 // bot edge will be v1->v2 03194 // top edge will be v3->v4 03195 // we need to create edges from v1 to v3 and from v2 to v4 03196 std::vector< EntityHandle > adj; 03197 rval = _mbImpl->get_child_meshsets( bEdge, adj ); 03198 MBERRORR( rval, "can't get children nodes" ); 03199 bool periodic = false; 03200 if( adj.size() == 1 ) 03201 { 03202 v[0] = v[1] = adj[0]; 03203 periodic = true; 03204 } 03205 else 03206 { 03207 v[0] = adj[0]; 03208 v[1] = adj[1]; 03209 } 03210 int senseB; 03211 rval = getEgVtxSense( bEdge, v[0], v[1], senseB ); 03212 MBERRORR( rval, "can't get bottom edge sense" ); 03213 if( -1 == senseB ) 03214 { 03215 v[1] = adj[0]; // so , bEdge will be oriented from v[0] to v[1], and will start at 03216 // nodes1, coords1.. 03217 v[0] = adj[1]; 03218 } 03219 adj.clear(); 03220 rval = _mbImpl->get_child_meshsets( tEdge, adj ); 03221 MBERRORR( rval, "can't get children nodes" ); 03222 if( adj.size() == 1 ) 03223 { 03224 v[2] = v[3] = adj[0]; 03225 if( !periodic ) MBERRORR( MB_FAILURE, "top edge is periodic, but bottom edge is not" ); 03226 } 03227 else 03228 { 03229 v[2] = adj[0]; 03230 v[3] = adj[1]; 03231 if( periodic ) MBERRORR( MB_FAILURE, "top edge is not periodic, but bottom edge is" ); 03232 } 03233 03234 // now, using nodes on bottom edge and top edge, create triangles, oriented outwards the 03235 // volume (sense positive on bottom edge) 03236 std::vector< EntityHandle > nodes1; 03237 rval = get_nodes_from_edge( bEdge, nodes1 ); 03238 MBERRORR( rval, "can't get nodes from bott edge" ); 03239 03240 std::vector< EntityHandle > nodes2; 03241 rval = get_nodes_from_edge( tEdge, nodes2 ); 03242 MBERRORR( rval, "can't get nodes from top edge" ); 03243 03244 std::vector< CartVect > coords1, coords2; 03245 coords1.resize( nodes1.size() ); 03246 coords2.resize( nodes2.size() ); 03247 03248 int N1 = (int)nodes1.size(); 03249 int N2 = (int)nodes2.size(); 03250 03251 rval = _mbImpl->get_coords( &( nodes1[0] ), nodes1.size(), (double*)&( coords1[0] ) ); 03252 MBERRORR( rval, "can't get coords of nodes from bott edge" ); 03253 03254 rval = _mbImpl->get_coords( &( nodes2[0] ), nodes2.size(), (double*)&( coords2[0] ) ); 03255 MBERRORR( rval, "can't get coords of nodes from top edge" ); 03256 CartVect up( direction ); 03257 03258 // see if the start and end coordinates are matching, if not, reverse edge 2 nodes and 03259 // coordinates 03260 CartVect v1 = ( coords1[0] - coords2[0] ) * up; 03261 CartVect v2 = ( coords1[0] - coords2[N2 - 1] ) * up; 03262 if( v1.length_squared() > v2.length_squared() ) 03263 { 03264 // we need to reverse coords2 and node2, as nodes are not above each other 03265 // the edges need to be found between v0 and v3, v1 and v2! 03266 for( unsigned int k = 0; k < nodes2.size() / 2; k++ ) 03267 { 03268 EntityHandle tmp = nodes2[k]; 03269 nodes2[k] = nodes2[N2 - 1 - k]; 03270 nodes2[N2 - 1 - k] = tmp; 03271 CartVect tv = coords2[k]; 03272 coords2[k] = coords2[N2 - 1 - k]; 03273 coords2[N2 - 1 - k] = tv; 03274 } 03275 } 03276 // make sure v[2] has nodes2[0], if not, reverse v[2] and v[3] 03277 if( !_mbImpl->contains_entities( v[2], &( nodes2[0] ), 1 ) ) 03278 { 03279 // reverse v[2] and v[3], so v[2] will be above v[0] 03280 EntityHandle tmp = v[2]; 03281 v[2] = v[3]; 03282 v[3] = tmp; 03283 } 03284 // now we know that v[0]--v[3] will be vertex sets in the order we want 03285 EntityHandle nd[4] = { nodes1[0], nodes1[N1 - 1], nodes2[0], nodes2[N2 - 1] }; 03286 03287 adj.clear(); 03288 EntityHandle e1, e2; 03289 // find edge 1 between v[0] and v[2], and e2 between v[1] and v[3] 03290 rval = _mbImpl->get_parent_meshsets( v[0], adj ); 03291 MBERRORR( rval, "can't get edges connected to vertex set 1" ); 03292 bool found = false; 03293 for( unsigned int j = 0; j < adj.size(); j++ ) 03294 { 03295 EntityHandle ed = adj[j]; 03296 Range vertices; 03297 rval = _mbImpl->get_child_meshsets( ed, vertices ); 03298 if( vertices.find( v[2] ) != vertices.end() ) 03299 { 03300 found = true; 03301 e1 = ed; 03302 break; 03303 } 03304 } 03305 if( !found ) 03306 { 03307 // create an edge from v[0] to v[2] 03308 rval = _mbImpl->create_meshset( MESHSET_SET, e1 ); 03309 MBERRORR( rval, "can't create edge 1" ); 03310 03311 rval = _mbImpl->add_parent_child( e1, v[0] ); 03312 MBERRORR( rval, "can't add parent - child relation" ); 03313 03314 rval = _mbImpl->add_parent_child( e1, v[2] ); 03315 MBERRORR( rval, "can't add parent - child relation" ); 03316 03317 EntityHandle nn2[2] = { nd[0], nd[2] }; 03318 EntityHandle medge; 03319 rval = _mbImpl->create_element( MBEDGE, nn2, 2, medge ); 03320 MBERRORR( rval, "can't create mesh edge" ); 03321 03322 rval = _mbImpl->add_entities( e1, &medge, 1 ); 03323 MBERRORR( rval, "can't add mesh edge to geo edge" ); 03324 03325 rval = this->_my_geomTopoTool->add_geo_set( e1, 1 ); 03326 MBERRORR( rval, "can't add edge to gtt" ); 03327 } 03328 03329 // find the edge from v2 to v4 (if not, create one) 03330 rval = _mbImpl->get_parent_meshsets( v[1], adj ); 03331 MBERRORR( rval, "can't get edges connected to vertex set 2" ); 03332 found = false; 03333 for( unsigned int i = 0; i < adj.size(); i++ ) 03334 { 03335 EntityHandle ed = adj[i]; 03336 Range vertices; 03337 rval = _mbImpl->get_child_meshsets( ed, vertices ); 03338 if( vertices.find( v[3] ) != vertices.end() ) 03339 { 03340 found = true; 03341 e2 = ed; 03342 break; 03343 } 03344 } 03345 if( !found ) 03346 { 03347 // create an edge from v2 to v4 03348 rval = _mbImpl->create_meshset( MESHSET_SET, e2 ); 03349 MBERRORR( rval, "can't create edge 1" ); 03350 03351 rval = _mbImpl->add_parent_child( e2, v[1] ); 03352 MBERRORR( rval, "can't add parent - child relation" ); 03353 03354 rval = _mbImpl->add_parent_child( e2, v[3] ); 03355 MBERRORR( rval, "can't add parent - child relation" ); 03356 03357 EntityHandle nn2[2] = { nd[1], nd[3] }; 03358 EntityHandle medge; 03359 rval = _mbImpl->create_element( MBEDGE, nn2, 2, medge ); 03360 MBERRORR( rval, "can't create mesh edge" ); 03361 03362 rval = _mbImpl->add_entities( e2, &medge, 1 ); 03363 MBERRORR( rval, "can't add mesh edge to geo edge" ); 03364 03365 rval = _my_geomTopoTool->add_geo_set( e2, 1 ); 03366 MBERRORR( rval, "can't add edge to gtt" ); 03367 } 03368 03369 // now we have the four edges, add them to the face, as children 03370 03371 // add children to face 03372 rval = _mbImpl->add_parent_child( newLatFace, bEdge ); 03373 MBERRORR( rval, "can't add parent - child relation" ); 03374 03375 rval = _mbImpl->add_parent_child( newLatFace, tEdge ); 03376 MBERRORR( rval, "can't add parent - child relation" ); 03377 03378 rval = _mbImpl->add_parent_child( newLatFace, e1 ); 03379 MBERRORR( rval, "can't add parent - child relation" ); 03380 03381 rval = _mbImpl->add_parent_child( newLatFace, e2 ); 03382 MBERRORR( rval, "can't add parent - child relation" ); 03383 03384 rval = _my_geomTopoTool->add_geo_set( newLatFace, 2 ); 03385 MBERRORR( rval, "can't add face to gtt" ); 03386 // add senses 03387 // 03388 rval = _my_geomTopoTool->set_sense( bEdge, newLatFace, 1 ); 03389 MBERRORR( rval, "can't set bottom edge sense to the lateral face" ); 03390 03391 int Tsense; 03392 rval = getEgVtxSense( tEdge, v[3], v[2], Tsense ); 03393 MBERRORR( rval, "can't get top edge sense" ); 03394 // we need to see what sense has topEdge in face 03395 rval = _my_geomTopoTool->set_sense( tEdge, newLatFace, Tsense ); 03396 MBERRORR( rval, "can't set top edge sense to the lateral face" ); 03397 03398 rval = _my_geomTopoTool->set_sense( e1, newLatFace, -1 ); 03399 MBERRORR( rval, "can't set first vert edge sense" ); 03400 03401 rval = _my_geomTopoTool->set_sense( e2, newLatFace, 1 ); 03402 MBERRORR( rval, "can't set second edge sense to the lateral face" ); 03403 // first, create edges along direction, for the 03404 03405 int indexB = 0, indexT = 0; // indices of the current nodes in the weaving process 03406 // weaving is either up or down; the triangles are oriented positively either way 03407 // up is nodes1[indexB], nodes2[indexT+1], nodes2[indexT] 03408 // down is nodes1[indexB], nodes1[indexB+1], nodes2[indexT] 03409 // the decision to weave up or down is based on which one is closer to the direction normal 03410 /* 03411 * 03412 * --------*------*-----------* ^ 03413 * / . \ . . ------> dir1 | up 03414 * /. \ . . | 03415 * -----*------------*----* 03416 * 03417 */ 03418 // we have to change the logic to account for cases when the curve in xy plane is not straight 03419 // (for example, when merging curves with a min_dot < -0.5, which means that the edges 03420 // could be even closed (periodic), with one vertex 03421 // the top and bottom curves should follow the same path in the "xy" plane (better to say 03422 // the plane perpendicular to the direction of weaving) 03423 // in this logic, the vector dir1 varies along the curve !!! 03424 CartVect dir1 = coords1[1] - coords1[0]; // we should have at least 2 nodes, N1>=2!! 03425 03426 CartVect planeNormal = dir1 * up; 03427 dir1 = up * planeNormal; 03428 dir1.normalize(); 03429 // this direction will be updated constantly with the direction of last edge added 03430 bool weaveDown = true; 03431 03432 CartVect startP = coords1[0]; // used for origin of comparisons 03433 while( 1 ) 03434 { 03435 if( ( indexB == N1 - 1 ) && ( indexT == N2 - 1 ) ) break; // we cannot advance anymore 03436 if( indexB == N1 - 1 ) { weaveDown = false; } 03437 else if( indexT == N2 - 1 ) 03438 { 03439 weaveDown = true; 03440 } 03441 else 03442 { 03443 // none are at the end, we have to decide which way to go, based on which index + 1 is 03444 // closer 03445 double proj1 = ( coords1[indexB + 1] - startP ) % dir1; 03446 double proj2 = ( coords2[indexT + 1] - startP ) % dir1; 03447 if( proj1 < proj2 ) 03448 weaveDown = true; 03449 else 03450 weaveDown = false; 03451 } 03452 EntityHandle nTri[3] = { nodes1[indexB], 0, nodes2[indexT] }; 03453 if( weaveDown ) 03454 { 03455 nTri[1] = nodes1[indexB + 1]; 03456 nTri[2] = nodes2[indexT]; 03457 indexB++; 03458 } 03459 else 03460 { 03461 nTri[1] = nodes2[indexT + 1]; 03462 indexT++; 03463 } 03464 EntityHandle triangle; 03465 rval = _mbImpl->create_element( MBTRI, nTri, 3, triangle ); 03466 MBERRORR( rval, "can't create triangle" ); 03467 03468 rval = _mbImpl->add_entities( newLatFace, &triangle, 1 ); 03469 MBERRORR( rval, "can't add triangle to face set" ); 03470 if( weaveDown ) 03471 { 03472 // increase was from nodes1[indexB-1] to nodes1[indexb] 03473 dir1 = coords1[indexB] - coords1[indexB - 1]; // we should have at least 2 nodes, N1>=2!! 03474 } 03475 else 03476 { 03477 dir1 = coords2[indexT] - coords2[indexT - 1]; 03478 } 03479 dir1 = up * ( dir1 * up ); 03480 dir1.normalize(); 03481 } 03482 // we do not check yet if the triangles are inverted 03483 // if yes, we should correct them. HOW? 03484 // we probably need a better meshing strategy, one that can overcome really bad meshes. 03485 // again, these faces are not what we should use for geometry, maybe we should use the 03486 // extruded quads, identified AFTER hexa are created. 03487 // right now, I see only a cosmetic use of these lateral faces 03488 // the whole idea of volume maybe is overrated 03489 // volume should be just quads extruded from bottom ? 03490 // 03491 return MB_SUCCESS; 03492 } 03493 // this will be used during creation of the volume, to set unique 03494 // tags on all surfaces 03495 // it is changing original tags, so do not use it if you want to preserve 03496 // initial neumann tags 03497 ErrorCode FBEngine::set_default_neumann_tags() 03498 { 03499 // these are for debugging purposes only 03500 // check the initial tag, then 03501 Tag ntag; 03502 ErrorCode rval = _mbImpl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, ntag ); 03503 MBERRORR( rval, "can't get tag handle" ); 03504 // get all surfaces in the model now 03505 Range sets[5]; 03506 rval = _my_geomTopoTool->find_geomsets( sets ); 03507 MBERRORR( rval, "can't get geo sets" ); 03508 int nfaces = (int)sets[2].size(); 03509 int* vals = new int[nfaces]; 03510 for( int i = 0; i < nfaces; i++ ) 03511 { 03512 vals[i] = i + 1; 03513 } 03514 rval = _mbImpl->tag_set_data( ntag, sets[2], (void*)vals ); 03515 MBERRORR( rval, "can't set tag values for neumann sets" ); 03516 03517 delete[] vals; 03518 03519 return MB_SUCCESS; 03520 } 03521 // a reverse operation for splitting an gedge at a mesh node 03522 ErrorCode FBEngine::chain_edges( double min_dot ) 03523 { 03524 Range sets[5]; 03525 ErrorCode rval; 03526 while( 1 ) // break out only if no edges are chained 03527 { 03528 rval = _my_geomTopoTool->find_geomsets( sets ); 03529 MBERRORR( rval, "can't get geo sets" ); 03530 // these gentities are "always" current, while the ones in this-> _my_gsets[0:4] are 03531 // the "originals" before FBEngine modifications 03532 int nedges = (int)sets[1].size(); 03533 // as long as we have chainable edges, continue; 03534 bool chain = false; 03535 for( int i = 0; i < nedges; i++ ) 03536 { 03537 EntityHandle edge = sets[1][i]; 03538 EntityHandle next_edge; 03539 bool chainable = false; 03540 rval = chain_able_edge( edge, min_dot, next_edge, chainable ); 03541 MBERRORR( rval, "can't determine chain-ability" ); 03542 if( chainable ) 03543 { 03544 rval = chain_two_edges( edge, next_edge ); 03545 MBERRORR( rval, "can't chain 2 edges" ); 03546 chain = true; 03547 break; // interrupt for loop 03548 } 03549 } 03550 if( !chain ) 03551 { 03552 break; // break out of while loop 03553 } 03554 } 03555 return MB_SUCCESS; 03556 } 03557 03558 // determine if from the end of edge we can extend with another edge; return also the 03559 // extension edge (next_edge) 03560 ErrorCode FBEngine::chain_able_edge( EntityHandle edge, double min_dot, EntityHandle& next_edge, bool& chainable ) 03561 { 03562 // get the end, then get all parents of end 03563 // see if some are the starts of 03564 chainable = false; 03565 EntityHandle v1, v2; 03566 ErrorCode rval = get_vert_edges( edge, v1, v2 ); 03567 MBERRORR( rval, "can't get vertices" ); 03568 if( v1 == v2 ) return MB_SUCCESS; // it is periodic, can't chain it with another edge! 03569 03570 // v2 is a set, get its parents, which should be edges 03571 Range edges; 03572 rval = _mbImpl->get_parent_meshsets( v2, edges ); 03573 MBERRORR( rval, "can't get parents of vertex set" ); 03574 // get parents of current edge (faces) 03575 Range faces; 03576 rval = _mbImpl->get_parent_meshsets( edge, faces ); 03577 MBERRORR( rval, "can't get parents of edge set" ); 03578 // get also the last edge "tangent" at the vertex 03579 std::vector< EntityHandle > mesh_edges; 03580 rval = _mbImpl->get_entities_by_type( edge, MBEDGE, mesh_edges ); 03581 MBERRORR( rval, "can't get mesh edges from edge set" ); 03582 EntityHandle lastMeshEdge = mesh_edges[mesh_edges.size() - 1]; 03583 const EntityHandle* conn2 = NULL; 03584 int len = 0; 03585 rval = _mbImpl->get_connectivity( lastMeshEdge, conn2, len ); 03586 MBERRORR( rval, "can't connectivity of last mesh edge" ); 03587 // get the coordinates of last edge 03588 if( len != 2 ) MBERRORR( MB_FAILURE, "bad number of vertices" ); 03589 CartVect P[2]; 03590 rval = _mbImpl->get_coords( conn2, len, (double*)&P[0] ); 03591 MBERRORR( rval, "Failed to get coordinates" ); 03592 03593 CartVect vec1( P[1] - P[0] ); 03594 vec1.normalize(); 03595 for( Range::iterator edgeIter = edges.begin(); edgeIter != edges.end(); ++edgeIter ) 03596 { 03597 EntityHandle otherEdge = *edgeIter; 03598 if( edge == otherEdge ) continue; 03599 // get faces parents of this edge 03600 Range faces2; 03601 rval = _mbImpl->get_parent_meshsets( otherEdge, faces2 ); 03602 MBERRORR( rval, "can't get parents of other edge set" ); 03603 if( faces != faces2 ) continue; 03604 // now, if the first mesh edge is within given angle, we can go on 03605 std::vector< EntityHandle > mesh_edges2; 03606 rval = _mbImpl->get_entities_by_type( otherEdge, MBEDGE, mesh_edges2 ); 03607 MBERRORR( rval, "can't get mesh edges from other edge set" ); 03608 EntityHandle firstMeshEdge = mesh_edges2[0]; 03609 const EntityHandle* conn22; 03610 int len2; 03611 rval = _mbImpl->get_connectivity( firstMeshEdge, conn22, len2 ); 03612 MBERRORR( rval, "can't connectivity of first mesh edge" ); 03613 if( len2 != 2 ) MBERRORR( MB_FAILURE, "bad number of vertices" ); 03614 if( conn2[1] != conn22[0] ) continue; // the mesh edges are not one after the other 03615 // get the coordinates of first edge 03616 03617 // CartVect P2[2]; 03618 rval = _mbImpl->get_coords( conn22, len, (double*)&P[0] ); 03619 CartVect vec2( P[1] - P[0] ); 03620 vec2.normalize(); 03621 if( vec1 % vec2 < min_dot ) continue; 03622 // we found our edge, victory! we can get out 03623 next_edge = otherEdge; 03624 chainable = true; 03625 return MB_SUCCESS; 03626 } 03627 03628 return MB_SUCCESS; // in general, hard to come by chain-able edges 03629 } 03630 ErrorCode FBEngine::chain_two_edges( EntityHandle edge, EntityHandle next_edge ) 03631 { 03632 // the biggest thing is to see the sense tags; or maybe not... 03633 // they should be correct :) 03634 // get the vertex sets 03635 EntityHandle v11, v12, v21, v22; 03636 ErrorCode rval = get_vert_edges( edge, v11, v12 ); 03637 MBERRORR( rval, "can't get vert sets" ); 03638 rval = get_vert_edges( next_edge, v21, v22 ); 03639 MBERRORR( rval, "can't get vert sets" ); 03640 assert( v12 == v21 ); 03641 std::vector< EntityHandle > mesh_edges; 03642 rval = MBI->get_entities_by_type( next_edge, MBEDGE, mesh_edges ); 03643 MBERRORR( rval, "can't get mesh edges" ); 03644 03645 rval = _mbImpl->add_entities( edge, &mesh_edges[0], (int)mesh_edges.size() ); 03646 MBERRORR( rval, "can't add new mesh edges" ); 03647 // remove the child - parent relation for second vertex of first edge 03648 rval = _mbImpl->remove_parent_child( edge, v12 ); 03649 MBERRORR( rval, "can't remove parent - child relation between first edge and middle vertex" ); 03650 03651 if( v22 != v11 ) // the edge would become periodic, do not add again the relationship 03652 { 03653 rval = _mbImpl->add_parent_child( edge, v22 ); 03654 MBERRORR( rval, "can't add second vertex to edge " ); 03655 } 03656 // we can now safely eliminate next_edge 03657 rval = _mbImpl->remove_parent_child( next_edge, v21 ); 03658 MBERRORR( rval, "can't remove child - parent relation " ); 03659 03660 rval = _mbImpl->remove_parent_child( next_edge, v22 ); 03661 MBERRORR( rval, "can't remove child - parent relation " ); 03662 03663 // remove the next_edge relation to the faces 03664 Range faces; 03665 rval = _mbImpl->get_parent_meshsets( next_edge, faces ); 03666 MBERRORR( rval, "can't get parent faces " ); 03667 03668 for( Range::iterator it = faces.begin(); it != faces.end(); ++it ) 03669 { 03670 EntityHandle ff = *it; 03671 rval = _mbImpl->remove_parent_child( ff, next_edge ); 03672 MBERRORR( rval, "can't remove parent-edge rel " ); 03673 } 03674 03675 rval = _mbImpl->delete_entities( &next_edge, 1 ); 03676 MBERRORR( rval, "can't remove edge set " ); 03677 03678 // delete the vertex set that is idle now (v12 = v21) 03679 rval = _mbImpl->delete_entities( &v12, 1 ); 03680 MBERRORR( rval, "can't remove edge set " ); 03681 return MB_SUCCESS; 03682 } 03683 ErrorCode FBEngine::get_vert_edges( EntityHandle edge, EntityHandle& v1, EntityHandle& v2 ) 03684 { 03685 // need to decide first or second vertex 03686 // important for moab 03687 03688 Range children; 03689 // EntityHandle v1, v2; 03690 ErrorCode rval = _mbImpl->get_child_meshsets( edge, children ); 03691 MBERRORR( rval, "can't get child meshsets" ); 03692 if( children.size() == 1 ) 03693 { 03694 // this is periodic edge, get out early 03695 v1 = children[0]; 03696 v2 = v1; 03697 return MB_SUCCESS; 03698 } 03699 else if( children.size() > 2 ) 03700 MBERRORR( MB_FAILURE, "too many vertices in one edge" ); 03701 // edge: get one vertex as part of the vertex set 03702 Range entities; 03703 rval = MBI->get_entities_by_type( children[0], MBVERTEX, entities ); 03704 MBERRORR( rval, "can't get entities from vertex set" ); 03705 if( entities.size() < 1 ) MBERRORR( MB_FAILURE, "no mesh nodes in vertex set" ); 03706 EntityHandle node0 = entities[0]; // the first vertex 03707 entities.clear(); 03708 03709 // now get the edges, and get the first node and the last node in sequence of edges 03710 // the order is important... 03711 // these are ordered sets !! 03712 std::vector< EntityHandle > ents; 03713 rval = MBI->get_entities_by_type( edge, MBEDGE, ents ); 03714 MBERRORR( rval, "can't get mesh edges" ); 03715 if( ents.size() < 1 ) MBERRORR( MB_FAILURE, "no mesh edges in edge set" ); 03716 03717 const EntityHandle* conn = NULL; 03718 int len; 03719 rval = MBI->get_connectivity( ents[0], conn, len ); 03720 MBERRORR( rval, "can't connectivity of first mesh edge" ); 03721 03722 if( conn[0] == node0 ) 03723 { 03724 v1 = children[0]; 03725 v2 = children[1]; 03726 } 03727 else // the other way around, although we should check (we are paranoid) 03728 { 03729 v2 = children[0]; 03730 v1 = children[1]; 03731 } 03732 03733 return MB_SUCCESS; 03734 } 03735 } // namespace moab