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