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