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