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