Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
moab::FBEngine Class Reference

#include <FBEngine.hpp>

+ Collaboration diagram for moab::FBEngine:

Public Member Functions

 FBEngine (Interface *impl, GeomTopoTool *geomTopoTool=NULL, const bool smooth=false)
 ~FBEngine ()
ErrorCode Init ()
ErrorCode getRootSet (EntityHandle *root_set)
ErrorCode getNumEntSets (EntityHandle set, int num_hops, int *all_sets)
ErrorCode createEntSet (int isList, EntityHandle *pSet)
ErrorCode addEntSet (EntityHandle entity_set_to_add, EntityHandle entity_set_handle)
ErrorCode getEntities (EntityHandle root_set, int ent_type, Range &gentities)
ErrorCode addEntArrToSet (const Range &entities, EntityHandle set)
ErrorCode getNumOfType (EntityHandle set, int ent_type, int *pNum)
ErrorCode getEntType (EntityHandle gent, int *type)
ErrorCode getEntBoundBox (EntityHandle this_gent, double *x0, double *y0, double *z0, double *x1, double *y1, double *z1)
ErrorCode getEntClosestPt (EntityHandle this_gent, double x, double y, double z, double *x1, double *y1, double *y3)
ErrorCode getVtxCoord (EntityHandle this_gent, double *x0, double *y0, double *z0)
ErrorCode gsubtract (EntityHandle entity_set_1, EntityHandle entity_set_2, EntityHandle result_entity_set)
ErrorCode getEntNrmlXYZ (EntityHandle entity_handle, double x, double y, double z, double *nrml_i, double *nrml_j, double *nrml_k)
ErrorCode getPntRayIntsct (double x, double y, double z, double dir_x, double dir_y, double dir_z, std::vector< EntityHandle > &intersect_entity_handles, std::vector< double > &intersect_coords, std::vector< double > &param_coords)
ErrorCode createTag (const char *tag_name, int tag_num_type_values, int tag_type, Tag &tag_handle_out)
Interfacemoab_instance ()
ErrorCode getArrData (const moab::EntityHandle *entity_handles, int entity_handles_size, Tag tag_handle, void *tag_values_out)
ErrorCode setArrData (const EntityHandle *entity_handles, int entity_handles_size, Tag tag_handle, const void *tag_values)
ErrorCode getEntAdj (EntityHandle handle, int type_requested, Range &adjEnts)
ErrorCode getEgFcSense (EntityHandle mbedge, EntityHandle mbface, int &sense)
ErrorCode measure (const EntityHandle *moab_entities, int entities_size, double *measures)
ErrorCode getEntNrmlSense (EntityHandle face, EntityHandle region, int &sense)
ErrorCode getEgEvalXYZ (EntityHandle edge, double x, double y, double z, double &on_x, double &on_y, double &on_z, double &tngt_i, double &tngt_j, double &tngt_k, double &cvtr_i, double &cvtr_j, double &cvtr_k)
ErrorCode getFcEvalXYZ (EntityHandle face, double x, double y, double z, double &on_x, double &on_y, double &on_z, double &nrml_i, double &nrml_j, double &nrml_k, double &cvtr1_i, double &cvtr1_j, double &cvtr1_k, double &cvtr2_i, double &cvtr2_j, double &cvtr2_k)
ErrorCode getEgVtxSense (EntityHandle edge, EntityHandle vtx1, EntityHandle vtx2, int &sense)
ErrorCode getEntURange (EntityHandle edge, double &u_min, double &u_max)
ErrorCode getEntUtoXYZ (EntityHandle edge, double u, double &x, double &y, double &z)
ErrorCode getEntTgntU (EntityHandle edge, double u, double &i, double &j, double &k)
ErrorCode isEntAdj (EntityHandle entity1, EntityHandle entity2, bool &adjacent_out)
ErrorCode split_surface_with_direction (EntityHandle face, std::vector< double > &xyz, double *direction, int closed, double min_dot, EntityHandle &oNewFace)
ErrorCode split_surface (EntityHandle face, std::vector< EntityHandle > &chainedEdges, std::vector< EntityHandle > &splittingNodes, EntityHandle &newFace)
ErrorCode split_edge_at_point (EntityHandle edge, CartVect &point, EntityHandle &new_edge)
ErrorCode split_edge_at_mesh_node (EntityHandle edge, EntityHandle node, EntityHandle &new_edge)
ErrorCode split_bedge_at_new_mesh_node (EntityHandle b_edge, EntityHandle atNode, EntityHandle brokenEdge, EntityHandle &new_edge)
void clean ()
void delete_smooth_tags ()
GeomTopoToolget_gtt ()
ErrorCode create_volume_with_direction (EntityHandle newFace1, EntityHandle newFace2, double *direction, EntityHandle &volume)
ErrorCode get_nodes_from_edge (EntityHandle gedge, std::vector< EntityHandle > &nodes)
ErrorCode weave_lateral_face_from_edges (EntityHandle bEdge, EntityHandle tEdge, double *direction, EntityHandle &newLatFace)
ErrorCode chain_edges (double min_dot)
ErrorCode chain_two_edges (EntityHandle edge, EntityHandle next_edge)
ErrorCode get_vert_edges (EntityHandle edge, EntityHandle &v1, EntityHandle &v2)
void set_smooth ()

Private Member Functions

ErrorCode initializeSmoothing ()
ErrorCode getAdjacentEntities (const EntityHandle from, const int to_dim, Range &adj_ents)
ErrorCode compute_intersection_points (EntityHandle &face, EntityHandle from, EntityHandle to, CartVect &Dir, std::vector< CartVect > &points, std::vector< EntityHandle > &entities, std::vector< EntityHandle > &triangles)
ErrorCode BreakTriangle (EntityHandle tri, EntityHandle e1, EntityHandle e3, EntityHandle n1, EntityHandle n2, EntityHandle n3)
ErrorCode BreakTriangle2 (EntityHandle tri, EntityHandle e1, EntityHandle e2, EntityHandle n1, EntityHandle n2)
void print_debug_triangle (EntityHandle triangle)
ErrorCode create_new_gedge (std::vector< EntityHandle > &nodesAlongPolyline, EntityHandle &new_geo_edge)
ErrorCode separate (EntityHandle face, std::vector< EntityHandle > &chainedEdges, Range &first, Range &second)
ErrorCode smooth_new_intx_points (EntityHandle face, std::vector< EntityHandle > &chainedEdges)
ErrorCode split_boundary (EntityHandle face, EntityHandle atNode)
bool find_vertex_set_for_node (EntityHandle iNode, EntityHandle &oVertexSet)
ErrorCode redistribute_boundary_edges_to_faces (EntityHandle face, EntityHandle newFace, std::vector< EntityHandle > &chainedEdges)
ErrorCode set_neumann_tags (EntityHandle face, EntityHandle newFace)
ErrorCode split_quads ()
ErrorCode boundary_nodes_on_face (EntityHandle face, std::vector< EntityHandle > &boundary_nodes)
ErrorCode boundary_mesh_edges_on_face (EntityHandle face, Range &boundary_mesh_edges)
ErrorCode split_internal_edge (EntityHandle &edge, EntityHandle &newVertex)
ErrorCode divide_triangle (EntityHandle triangle, EntityHandle &newVertex)
ErrorCode set_default_neumann_tags ()
ErrorCode chain_able_edge (EntityHandle edge, double min_dot, EntityHandle &next_edge, bool &chainable)

Private Attributes

Interface_mbImpl
GeomTopoTool_my_geomTopoTool
bool _t_created
bool _smooth
bool _initialized
Range _my_gsets [5]
std::map< EntityHandle,
SmoothFace * > 
_faces
std::map< EntityHandle,
SmoothCurve * > 
_edges
SmoothFace ** _smthFace
SmoothCurve ** _smthCurve
Range _piercedTriangles
Range _newTriangles
Range _piercedEdges
std::map< EntityHandle,
EntityHandle
_brokenEdges

Detailed Description

Definition at line 24 of file FBEngine.hpp.


Constructor & Destructor Documentation

moab::FBEngine::FBEngine ( Interface impl,
GeomTopoTool geomTopoTool = NULL,
const bool  smooth = false 
)

Definition at line 192 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, and _t_created.

    : _mbImpl( impl ), _my_geomTopoTool( topoTool ), _t_created( false ), _smooth( smooth ), _initialized( false ),
      _smthFace( NULL ), _smthCurve( NULL )
{
    if( !_my_geomTopoTool )
    {
        _my_geomTopoTool = new GeomTopoTool( _mbImpl );
        _t_created       = true;
    }
    // should this be part of the constructor or not?
    // Init();
}

Definition at line 204 of file FBEngine.cpp.

References _smooth, and clean().

{
    clean();
    _smooth = false;
}

Member Function Documentation

Definition at line 488 of file FBEngine.cpp.

References MBI.

{
    return MBI->add_entities( set, entities );
}
ErrorCode moab::FBEngine::addEntSet ( EntityHandle  entity_set_to_add,
EntityHandle  entity_set_handle 
)

Definition at line 493 of file FBEngine.cpp.

References MBI.

{
    return MBI->add_entities( entity_set_handle, &entity_set_to_add, 1 );
}
ErrorCode moab::FBEngine::boundary_mesh_edges_on_face ( EntityHandle  face,
Range boundary_mesh_edges 
) [private]

Definition at line 3047 of file FBEngine.cpp.

References _mbImpl, moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_entities_by_dimension(), getAdjacentEntities(), MB_SUCCESS, and MBERRORR.

Referenced by split_surface_with_direction().

{
    // this list is used only for finding the intersecting mesh edge for starting the
    // polygonal cutting line at boundary (if !closed)
    Range bound_edges;
    ErrorCode rval = getAdjacentEntities( face, 1, bound_edges );
    MBERRORR( rval, " can't get boundary edges" );
    for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
    {
        EntityHandle b_edge = *it;
        // get all edges in range
        // Range mesh_edges;
        rval = _mbImpl->get_entities_by_dimension( b_edge, 1, boundary_mesh_edges );
        MBERRORR( rval, " can't get mesh edges" );
    }
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::boundary_nodes_on_face ( EntityHandle  face,
std::vector< EntityHandle > &  boundary_nodes 
) [private]

Definition at line 3064 of file FBEngine.cpp.

References _mbImpl, moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_entities_by_dimension(), getAdjacentEntities(), MB_SUCCESS, MBERRORR, and moab::Range::size().

Referenced by split_surface_with_direction().

{
    // even if we repeat some nodes, it is OK
    // this list is used only for finding the closest boundary node for starting the
    // polygonal cutting line at boundary (if !closed)
    Range bound_edges;
    ErrorCode rval = getAdjacentEntities( face, 1, bound_edges );
    MBERRORR( rval, " can't get boundary edges" );
    Range b_nodes;
    for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
    {
        EntityHandle b_edge = *it;
        // get all edges in range
        Range mesh_edges;
        rval = _mbImpl->get_entities_by_dimension( b_edge, 1, mesh_edges );
        MBERRORR( rval, " can't get mesh edges" );
        rval = _mbImpl->get_connectivity( mesh_edges, b_nodes );
        MBERRORR( rval, " can't get nodes from mesh edges" );
    }
    // create now a vector based on Range of nodes
    boundary_nodes.resize( b_nodes.size() );
    std::copy( b_nodes.begin(), b_nodes.end(), boundary_nodes.begin() );
    return MB_SUCCESS;
}

Definition at line 1989 of file FBEngine.cpp.

{
    std::cout << "FBEngine::BreakTriangle not implemented yet\n";
    return MB_FAILURE;
}

Definition at line 1995 of file FBEngine.cpp.

References _mbImpl, _newTriangles, moab::Interface::create_element(), moab::debug_splits, ErrorCode, moab::Interface::get_connectivity(), moab::Range::insert(), MB_SUCCESS, MBERRORR, MBTRI, MBVERTEX, print_debug_triangle(), moab::Interface::side_number(), and moab::Interface::type_from_handle().

Referenced by split_surface_with_direction().

{
    // we have the nodes, we just need to reconnect to form new triangles
    ErrorCode rval;
    const EntityHandle* conn3 = NULL;
    int nnodes                = 0;
    rval                      = _mbImpl->get_connectivity( tri, conn3, nnodes );
    MBERRORR( rval, "Failed to get connectivity" );
    assert( 3 == nnodes );

    EntityType et1 = _mbImpl->type_from_handle( e1 );
    EntityType et2 = _mbImpl->type_from_handle( e2 );

    if( MBVERTEX == et1 )
    {
        // find the vertex in conn3, and form 2 other triangles
        int index = -1;
        for( index = 0; index < 3; index++ )
        {
            if( conn3[index] == e1 )  // also n1
                break;
        }
        if( index == 3 ) return MB_FAILURE;
        // 2 triangles: n1, index+1, n2, and n1, n2, index+2
        EntityHandle conn[6] = { n1, conn3[( index + 1 ) % 3], n2, n1, n2, conn3[( index + 2 ) % 3] };
        EntityHandle newTriangle;
        rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle );
        MBERRORR( rval, "Failed to create a new triangle" );
        _newTriangles.insert( newTriangle );
        if( debug_splits ) print_debug_triangle( newTriangle );
        rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle );  // the second triangle
        MBERRORR( rval, "Failed to create a new triangle" );
        _newTriangles.insert( newTriangle );
        if( debug_splits ) print_debug_triangle( newTriangle );
    }
    else if( MBVERTEX == et2 )
    {
        int index = -1;
        for( index = 0; index < 3; index++ )
        {
            if( conn3[index] == e2 )  // also n2
                break;
        }
        if( index == 3 ) return MB_FAILURE;
        // 2 triangles: n1, index+1, n2, and n1, n2, index+2
        EntityHandle conn[6] = { n2, conn3[( index + 1 ) % 3], n1, n2, n1, conn3[( index + 2 ) % 3] };
        EntityHandle newTriangle;
        rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle );
        MBERRORR( rval, "Failed to create a new triangle" );
        _newTriangles.insert( newTriangle );
        if( debug_splits ) print_debug_triangle( newTriangle );
        rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle );  // the second triangle
        MBERRORR( rval, "Failed to create a new triangle" );
        _newTriangles.insert( newTriangle );
        if( debug_splits ) print_debug_triangle( newTriangle );
    }
    else
    {
        // both are edges adjacent to triangle tri
        // there are several configurations possible for n1, n2, between conn3 nodes.
        int num1, num2, sense, offset;
        rval = _mbImpl->side_number( tri, e1, num1, sense, offset );
        MBERRORR( rval, "edge not adjacent" );

        rval = _mbImpl->side_number( tri, e2, num2, sense, offset );
        MBERRORR( rval, "edge not adjacent" );

        const EntityHandle* conn12;  // connectivity for edge 1
        const EntityHandle* conn22;  // connectivity for edge 2
        // int nnodes;
        rval = _mbImpl->get_connectivity( e1, conn12, nnodes );
        MBERRORR( rval, "Failed to get connectivity of edge 1" );
        assert( 2 == nnodes );
        rval = _mbImpl->get_connectivity( e2, conn22, nnodes );
        MBERRORR( rval, "Failed to get connectivity of edge 2" );
        assert( 2 == nnodes );
        // now, having the side number, work through
        if( debug_splits )
        {
            std::cout << "tri conn3:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n";
            std::cout << " edge1: conn12:" << conn12[0] << " " << conn12[1] << "  side: " << num1 << "\n";
            std::cout << " edge2: conn22:" << conn22[0] << " " << conn22[1] << "  side: " << num2 << "\n";
        }
        int unaffectedSide = 3 - num1 - num2;
        int i3             = ( unaffectedSide + 2 ) % 3;  // to 0 is 2, to 1 is 0, to 2 is 1
        // triangles will be formed with triVertexIndex , n1, n2 (in what order?)
        EntityHandle v1, v2;  // to hold the 2 nodes on edges
        if( num1 == i3 )
        {
            v1 = n1;
            v2 = n2;
        }
        else  // if (num2==i3)
        {
            v1 = n2;
            v2 = n1;
        }
        // three triangles are formed
        int i1 = ( i3 + 1 ) % 3;
        int i2 = ( i3 + 2 ) % 3;
        // we could break the surface differently
        EntityHandle conn[9] = { conn3[i3], v1, v2, v1, conn3[i1], conn3[i2], v2, v1, conn3[i2] };
        EntityHandle newTriangle;
        if( debug_splits ) std::cout << "Split 2 edges :\n";
        rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle );
        MBERRORR( rval, "Failed to create a new triangle" );
        _newTriangles.insert( newTriangle );
        if( debug_splits ) print_debug_triangle( newTriangle );
        rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle );  // the second triangle
        MBERRORR( rval, "Failed to create a new triangle" );
        _newTriangles.insert( newTriangle );
        if( debug_splits ) print_debug_triangle( newTriangle );
        rval = _mbImpl->create_element( MBTRI, conn + 6, 3, newTriangle );  // the second triangle
        MBERRORR( rval, "Failed to create a new triangle" );
        _newTriangles.insert( newTriangle );
        if( debug_splits ) print_debug_triangle( newTriangle );
    }

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::chain_able_edge ( EntityHandle  edge,
double  min_dot,
EntityHandle next_edge,
bool &  chainable 
) [private]

Definition at line 3665 of file FBEngine.cpp.

References _mbImpl, moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), moab::Interface::get_parent_meshsets(), get_vert_edges(), MB_SUCCESS, MBEDGE, MBERRORR, and moab::CartVect::normalize().

Referenced by chain_edges().

{
    // get the end, then get all parents of end
    // see if some are the starts of
    chainable = false;
    EntityHandle v1, v2;
    ErrorCode rval = get_vert_edges( edge, v1, v2 );
    MBERRORR( rval, "can't get vertices" );
    if( v1 == v2 ) return MB_SUCCESS;  // it is periodic, can't chain it with another edge!

    // v2 is a set, get its parents, which should be edges
    Range edges;
    rval = _mbImpl->get_parent_meshsets( v2, edges );
    MBERRORR( rval, "can't get parents of vertex set" );
    // get parents of current edge (faces)
    Range faces;
    rval = _mbImpl->get_parent_meshsets( edge, faces );
    MBERRORR( rval, "can't get parents of edge set" );
    // get also the last edge "tangent" at the vertex
    std::vector< EntityHandle > mesh_edges;
    rval = _mbImpl->get_entities_by_type( edge, MBEDGE, mesh_edges );
    MBERRORR( rval, "can't get mesh edges from edge set" );
    EntityHandle lastMeshEdge = mesh_edges[mesh_edges.size() - 1];
    const EntityHandle* conn2 = NULL;
    int len                   = 0;
    rval                      = _mbImpl->get_connectivity( lastMeshEdge, conn2, len );
    MBERRORR( rval, "can't connectivity of last mesh edge" );
    // get the coordinates of last edge
    if( len != 2 ) MBERRORR( MB_FAILURE, "bad number of vertices" );
    CartVect P[2];
    rval = _mbImpl->get_coords( conn2, len, (double*)&P[0] );
    MBERRORR( rval, "Failed to get coordinates" );

    CartVect vec1( P[1] - P[0] );
    vec1.normalize();
    for( Range::iterator edgeIter = edges.begin(); edgeIter != edges.end(); ++edgeIter )
    {
        EntityHandle otherEdge = *edgeIter;
        if( edge == otherEdge ) continue;
        // get faces parents of this edge
        Range faces2;
        rval = _mbImpl->get_parent_meshsets( otherEdge, faces2 );
        MBERRORR( rval, "can't get parents of other edge set" );
        if( faces != faces2 ) continue;
        // now, if the first mesh edge is within given angle, we can go on
        std::vector< EntityHandle > mesh_edges2;
        rval = _mbImpl->get_entities_by_type( otherEdge, MBEDGE, mesh_edges2 );
        MBERRORR( rval, "can't get mesh edges from other edge set" );
        EntityHandle firstMeshEdge = mesh_edges2[0];
        const EntityHandle* conn22;
        int len2;
        rval = _mbImpl->get_connectivity( firstMeshEdge, conn22, len2 );
        MBERRORR( rval, "can't connectivity of first mesh edge" );
        if( len2 != 2 ) MBERRORR( MB_FAILURE, "bad number of vertices" );
        if( conn2[1] != conn22[0] ) continue;  // the mesh edges are not one after the other
        // get the coordinates of first edge

        // CartVect P2[2];
        rval = _mbImpl->get_coords( conn22, len, (double*)&P[0] );
        CartVect vec2( P[1] - P[0] );
        vec2.normalize();
        if( vec1 % vec2 < min_dot ) continue;
        // we found our edge, victory! we can get out
        next_edge = otherEdge;
        chainable = true;
        return MB_SUCCESS;
    }

    return MB_SUCCESS;  // in general, hard to come by chain-able edges
}

Definition at line 3627 of file FBEngine.cpp.

References _my_geomTopoTool, chain_able_edge(), chain_two_edges(), ErrorCode, moab::GeomTopoTool::find_geomsets(), MB_SUCCESS, MBERRORR, and size.

Referenced by split_surface_with_direction().

{
    Range sets[5];
    ErrorCode rval;
    while( 1 )  // break out only if no edges are chained
    {
        rval = _my_geomTopoTool->find_geomsets( sets );
        MBERRORR( rval, "can't get geo sets" );
        // these gentities are "always" current, while the ones in this-> _my_gsets[0:4] are
        // the "originals" before FBEngine modifications
        int nedges = (int)sets[1].size();
        // as long as we have chainable edges, continue;
        bool chain = false;
        for( int i = 0; i < nedges; i++ )
        {
            EntityHandle edge = sets[1][i];
            EntityHandle next_edge;
            bool chainable = false;
            rval           = chain_able_edge( edge, min_dot, next_edge, chainable );
            MBERRORR( rval, "can't determine chain-ability" );
            if( chainable )
            {
                rval = chain_two_edges( edge, next_edge );
                MBERRORR( rval, "can't chain 2 edges" );
                chain = true;
                break;  // interrupt for loop
            }
        }
        if( !chain )
        {
            break;  // break out of while loop
        }
    }
    return MB_SUCCESS;
}

Definition at line 3735 of file FBEngine.cpp.

References _mbImpl, moab::Interface::add_entities(), moab::Interface::add_parent_child(), moab::Range::begin(), moab::Interface::delete_entities(), moab::Range::end(), ErrorCode, moab::Interface::get_parent_meshsets(), get_vert_edges(), MB_SUCCESS, MBEDGE, MBERRORR, MBI, and moab::Interface::remove_parent_child().

Referenced by chain_edges().

{
    // the biggest thing is to see the sense tags; or maybe not...
    // they should be correct :)
    // get the vertex sets
    EntityHandle v11, v12, v21, v22;
    ErrorCode rval = get_vert_edges( edge, v11, v12 );
    MBERRORR( rval, "can't get vert sets" );
    rval = get_vert_edges( next_edge, v21, v22 );
    MBERRORR( rval, "can't get vert sets" );
    assert( v12 == v21 );
    std::vector< EntityHandle > mesh_edges;
    rval = MBI->get_entities_by_type( next_edge, MBEDGE, mesh_edges );
    MBERRORR( rval, "can't get mesh edges" );

    rval = _mbImpl->add_entities( edge, &mesh_edges[0], (int)mesh_edges.size() );
    MBERRORR( rval, "can't add new mesh edges" );
    // remove the child - parent relation for second vertex of first edge
    rval = _mbImpl->remove_parent_child( edge, v12 );
    MBERRORR( rval, "can't remove parent - child relation between first edge and middle vertex" );

    if( v22 != v11 )  // the edge would become periodic, do not add again the relationship
    {
        rval = _mbImpl->add_parent_child( edge, v22 );
        MBERRORR( rval, "can't add second vertex to edge " );
    }
    // we can now safely eliminate next_edge
    rval = _mbImpl->remove_parent_child( next_edge, v21 );
    MBERRORR( rval, "can't remove child - parent relation " );

    rval = _mbImpl->remove_parent_child( next_edge, v22 );
    MBERRORR( rval, "can't remove child - parent relation " );

    // remove the next_edge relation to the faces
    Range faces;
    rval = _mbImpl->get_parent_meshsets( next_edge, faces );
    MBERRORR( rval, "can't get parent faces " );

    for( Range::iterator it = faces.begin(); it != faces.end(); ++it )
    {
        EntityHandle ff = *it;
        rval            = _mbImpl->remove_parent_child( ff, next_edge );
        MBERRORR( rval, "can't remove parent-edge rel " );
    }

    rval = _mbImpl->delete_entities( &next_edge, 1 );
    MBERRORR( rval, "can't remove edge set " );

    // delete the vertex set that is idle now (v12 = v21)
    rval = _mbImpl->delete_entities( &v12, 1 );
    MBERRORR( rval, "can't remove edge set " );
    return MB_SUCCESS;
}

Definition at line 210 of file FBEngine.cpp.

References _edges, _faces, _my_geomTopoTool, _my_gsets, _smooth, _smthCurve, _smthFace, _t_created, and moab::Range::size().

Referenced by ~FBEngine().

{
    if( _smooth )
    {
        _faces.clear();
        _edges.clear();
        int size1 = _my_gsets[1].size();
        int i     = 0;
        for( i = 0; i < size1; i++ )
            delete _smthCurve[i];
        delete[] _smthCurve;
        _smthCurve = NULL;
        size1      = _my_gsets[2].size();
        for( i = 0; i < size1; i++ )
            delete _smthFace[i];
        delete[] _smthFace;
        _smthFace = NULL;
        //_smooth = false;
    }

    for( int j = 0; j < 5; j++ )
        _my_gsets[j].clear();
    if( _t_created ) delete _my_geomTopoTool;
    _my_geomTopoTool = NULL;
    _t_created       = false;
}
ErrorCode moab::FBEngine::compute_intersection_points ( EntityHandle face,
EntityHandle  from,
EntityHandle  to,
CartVect Dir,
std::vector< CartVect > &  points,
std::vector< EntityHandle > &  entities,
std::vector< EntityHandle > &  triangles 
) [private]

Definition at line 2124 of file FBEngine.cpp.

References _mbImpl, _piercedTriangles, moab::debug_splits, moab::Range::end(), ErrorCode, moab::Range::find(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::id_from_handle(), moab::Interface::INTERSECT, moab::intersect_segment_and_plane_slice(), moab::CartVect::length(), moab::Interface::list_entity(), MB_SUCCESS, MBERRORR, MBVERTEX, moab::CartVect::normalize(), moab::tolerance_segment, and moab::Interface::type_from_handle().

Referenced by split_surface_with_direction().

{
    // keep a stack of triangles to process, and do not add those already processed
    // either mark them, or maybe keep them in a local set?
    // from and to are now nodes, start from them
    CartVect p1, p2;  // the position of from and to
    ErrorCode rval = _mbImpl->get_coords( &from, 1, (double*)&p1 );
    MBERRORR( rval, "failed to get 'from' coordinates" );
    rval = _mbImpl->get_coords( &to, 1, (double*)&p2 );
    MBERRORR( rval, "failed to get 'from' coordinates" );

    CartVect vect( p2 - p1 );
    double dist2 = vect.length();
    if( dist2 < tolerance_segment )
    {
        // we are done, return
        return MB_SUCCESS;
    }
    CartVect normPlane = Dir * vect;
    normPlane.normalize();
    std::set< EntityHandle > visitedTriangles;
    CartVect currentPoint = p1;
    // push the current point if it is empty
    if( points.size() == 0 )
    {
        points.push_back( p1 );
        entities.push_back( from );  // this is a node now
    }

    // these will be used in many loops
    CartVect intx = p1;  // somewhere to start
    double param  = -1.;

    // first intersection
    EntityHandle currentBoundary = from;  // it is a node, in the beginning

    vect = p2 - currentPoint;
    while( vect.length() > 0. )
    {
        // advance towards "to" node, from boundary handle
        EntityType etype = _mbImpl->type_from_handle( currentBoundary );
        // if vertex, look for other triangles connected which intersect our plane (defined by p1,
        // p2, dir)
        std::vector< EntityHandle > adj_tri;
        rval           = _mbImpl->get_adjacencies( &currentBoundary, 1, 2, false, adj_tri );
        unsigned int j = 0;
        EntityHandle tri;
        for( ; j < adj_tri.size(); j++ )
        {
            tri = adj_tri[j];
            if( visitedTriangles.find( tri ) != visitedTriangles.end() )
                continue;  // get another triangle, this one was already visited
            // check if it is one of the triangles that was pierced already
            if( _piercedTriangles.find( tri ) != _piercedTriangles.end() ) continue;
            // if vertex, look for opposite edge
            // if edge, look for 2 opposite edges
            // get vertices
            int nnodes;
            const EntityHandle* conn3;
            rval = _mbImpl->get_connectivity( tri, conn3, nnodes );
            MBERRORR( rval, "Failed to get connectivity" );
            // if one of the nodes is to, stop right there
            {
                if( conn3[0] == to || conn3[1] == to || conn3[2] == to )
                {
                    visitedTriangles.insert( tri );
                    triangles.push_back( tri );
                    currentPoint = p2;
                    points.push_back( p2 );
                    entities.push_back( to );  // we are done
                    break;                     // this is break from for loop, we still need to get out of while
                    // we will get out, because vect will become 0, (p2-p2)
                }
            }
            EntityHandle nn2[2];
            if( MBVERTEX == etype )
            {
                nn2[0] = conn3[0];
                nn2[1] = conn3[1];
                if( nn2[0] == currentBoundary ) nn2[0] = conn3[2];
                if( nn2[1] == currentBoundary ) nn2[1] = conn3[2];
                // get coordinates
                CartVect Pt[2];

                rval = _mbImpl->get_coords( nn2, 2, (double*)&Pt[0] );
                MBERRORR( rval, "Failed to get coordinates" );
                // see the intersection
                if( intersect_segment_and_plane_slice( Pt[0], Pt[1], currentPoint, p2, Dir, normPlane, intx, param ) )
                {
                    // we should stop for loop, and decide if it is edge or vertex
                    if( param == 0. )
                        currentBoundary = nn2[0];
                    else
                    {
                        if( param == 1. )
                            currentBoundary = nn2[1];
                        else  // param between 0 and 1, so edge
                        {
                            // find the edge between vertices
                            std::vector< EntityHandle > edges1;
                            // change the create flag to true, because that edge must exist in
                            // current triangle if we want to advance; nn2 are 2 nodes in current
                            // triangle!!
                            rval = _mbImpl->get_adjacencies( nn2, 2, 1, true, edges1, Interface::INTERSECT );
                            MBERRORR( rval, "Failed to get edges" );
                            if( edges1.size() != 1 ) MBERRORR( MB_FAILURE, "Failed to get adjacent edges to 2 nodes" );
                            currentBoundary = edges1[0];
                        }
                    }
                    visitedTriangles.insert( tri );
                    currentPoint = intx;
                    points.push_back( intx );
                    entities.push_back( currentBoundary );
                    triangles.push_back( tri );
                    if( debug_splits )
                        std::cout << "vtx new tri : " << _mbImpl->id_from_handle( tri )
                                  << " type bdy:" << _mbImpl->type_from_handle( currentBoundary ) << "\n";
                    break;  // out of for loop over triangles
                }
            }
            else  // this is MBEDGE, we have the other triangle to try out
            {
                // first find the nodes from existing boundary
                int nnodes2;
                const EntityHandle* conn2;
                rval = _mbImpl->get_connectivity( currentBoundary, conn2, nnodes2 );
                MBERRORR( rval, "Failed to get connectivity" );
                int thirdIndex = -1;
                for( int tj = 0; tj < 3; tj++ )
                {
                    if( ( conn3[tj] != conn2[0] ) && ( conn3[tj] != conn2[1] ) )
                    {
                        thirdIndex = tj;
                        break;
                    }
                }
                if( -1 == thirdIndex ) MBERRORR( MB_FAILURE, " can't get third node" );
                CartVect Pt[3];
                rval = _mbImpl->get_coords( conn3, 3, (double*)&Pt[0] );
                MBERRORR( rval, "Failed to get coords" );
                int indexFirst  = ( thirdIndex + 1 ) % 3;
                int indexSecond = ( thirdIndex + 2 ) % 3;
                int index[2]    = { indexFirst, indexSecond };
                for( int k = 0; k < 2; k++ )
                {
                    nn2[0] = conn3[index[k]], nn2[1] = conn3[thirdIndex];
                    if( intersect_segment_and_plane_slice( Pt[index[k]], Pt[thirdIndex], currentPoint, p2, Dir,
                                                           normPlane, intx, param ) )
                    {
                        // we should stop for loop, and decide if it is edge or vertex
                        if( param == 0. )
                            currentBoundary = conn3[index[k]];  // it is not really possible
                        else
                        {
                            if( param == 1. )
                                currentBoundary = conn3[thirdIndex];
                            else  // param between 0 and 1, so edge is fine
                            {
                                // find the edge between vertices
                                std::vector< EntityHandle > edges1;
                                // change the create flag to true, because that edge must exist in
                                // current triangle if we want to advance; nn2 are 2 nodes in
                                // current triangle!!
                                rval = _mbImpl->get_adjacencies( nn2, 2, 1, true, edges1, Interface::INTERSECT );
                                MBERRORR( rval, "Failed to get edges" );
                                if( edges1.size() != 1 )
                                    MBERRORR( MB_FAILURE, "Failed to get adjacent edges to 2 nodes" );
                                currentBoundary = edges1[0];
                            }
                        }
                        visitedTriangles.insert( tri );
                        currentPoint = intx;
                        points.push_back( intx );
                        entities.push_back( currentBoundary );
                        triangles.push_back( tri );
                        if( debug_splits )
                        {
                            std::cout << "edge new tri : " << _mbImpl->id_from_handle( tri )
                                      << "  type bdy: " << _mbImpl->type_from_handle( currentBoundary )
                                      << " id: " << _mbImpl->id_from_handle( currentBoundary ) << "\n";
                            _mbImpl->list_entity( currentBoundary );
                        }
                        break;  // out of for loop over triangles
                    }
                    // we should not reach here
                }
            }
        }
        /*if (j==adj_tri.size())
         {
         MBERRORR(MB_FAILURE, "did not advance");
         }*/
        vect = p2 - currentPoint;
    }

    if( debug_splits )
        std::cout << "nb entities: " << entities.size() << " triangles:" << triangles.size()
                  << " points.size(): " << points.size() << "\n";

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::create_new_gedge ( std::vector< EntityHandle > &  nodesAlongPolyline,
EntityHandle new_geo_edge 
) [private]

Definition at line 1872 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, _piercedEdges, moab::Interface::add_entities(), moab::GeomTopoTool::add_geo_set(), moab::Interface::add_parent_child(), moab::Interface::create_element(), moab::Interface::create_meshset(), moab::debug_splits, ErrorCode, find_vertex_set_for_node(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::GeomTopoTool::global_id(), moab::Interface::id_from_handle(), moab::Range::insert(), moab::Interface::INTERSECT, MB_CHK_ERR, MBEDGE, MBERRORR, MESHSET_SET, and moab::Range::size().

Referenced by split_surface_with_direction().

{

    ErrorCode rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_geo_edge );
    MBERRORR( rval, "can't create geo edge" );

    // now, get the edges, or create if not existing
    std::vector< EntityHandle > mesh_edges;
    for( unsigned int i = 0; i < nodesAlongPolyline.size() - 1; i++ )
    {
        EntityHandle n1 = nodesAlongPolyline[i], n2 = nodesAlongPolyline[i + 1];

        EntityHandle nn2[2];
        nn2[0] = n1;
        nn2[1] = n2;

        std::vector< EntityHandle > adjacent;
        rval = _mbImpl->get_adjacencies( nn2, 2, 1, false, adjacent, Interface::INTERSECT );
        // see if there is an edge between those 2 already, and if it is oriented as we like
        bool new_edge = true;
        if( adjacent.size() >= 1 )
        {
            // check the orientation
            const EntityHandle* conn2 = NULL;
            int nnodes                = 0;
            rval                      = _mbImpl->get_connectivity( adjacent[0], conn2, nnodes );
            MBERRORR( rval, "can't get connectivity" );
            if( conn2[0] == nn2[0] && conn2[1] == nn2[1] )
            {
                // everything is fine
                mesh_edges.push_back( adjacent[0] );
                new_edge = false;  // we found one that's good, no need to create a new one
            }
            else
            {
                _piercedEdges.insert( adjacent[0] );  // we want to remove this one, it will be not needed
            }
        }
        if( new_edge )
        {
            // there is no edge between n1 and n2, create one
            EntityHandle mesh_edge;
            rval = _mbImpl->create_element( MBEDGE, nn2, 2, mesh_edge );
            MBERRORR( rval, "Failed to create a new edge" );
            mesh_edges.push_back( mesh_edge );
        }
    }

    // add loops edges to the edge set
    rval = _mbImpl->add_entities( new_geo_edge, &mesh_edges[0],
                                  mesh_edges.size() );  // only one edge
    MBERRORR( rval, "can't add edges to new_geo_set" );
    // check vertex sets for vertex 1 and vertex 2?
    // get all sets of dimension 0 from database, and see if our ends are here or not

    Range ends_geo_edge;
    ends_geo_edge.insert( nodesAlongPolyline[0] );
    ends_geo_edge.insert( nodesAlongPolyline[nodesAlongPolyline.size() - 1] );

    for( unsigned int k = 0; k < ends_geo_edge.size(); k++ )
    {
        EntityHandle node = ends_geo_edge[k];
        EntityHandle nodeSet;
        bool found = find_vertex_set_for_node( node, nodeSet );

        if( !found )
        {
            // create a node set and add the node

            rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet );
            MBERRORR( rval, "Failed to create a new vertex set" );

            rval = _mbImpl->add_entities( nodeSet, &node, 1 );
            MBERRORR( rval, "Failed to add the node to the set" );

            rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 );  //
            MBERRORR( rval, "Failed to commit the node set" );

            if( debug_splits )
            {
                std::cout << " create a vertex set " << _mbImpl->id_from_handle( nodeSet )
                          << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << " for node " << node
                          << "\n";
            }
        }

        rval = _mbImpl->add_parent_child( new_geo_edge, nodeSet );
        MBERRORR( rval, "Failed to add parent child relation" );
    }
    // finally, put the edge in the range of edges
    rval = _my_geomTopoTool->add_geo_set( new_geo_edge, 1 );MB_CHK_ERR( rval );

    return rval;
}
ErrorCode moab::FBEngine::create_volume_with_direction ( EntityHandle  newFace1,
EntityHandle  newFace2,
double *  direction,
EntityHandle volume 
)

Definition at line 3191 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, moab::GeomTopoTool::add_geo_set(), moab::Interface::add_parent_child(), moab::Interface::create_meshset(), ErrorCode, moab::Interface::get_child_meshsets(), MATERIAL_SET_TAG_NAME, MB_SUCCESS, MB_TYPE_INTEGER, MBERRORR, MESHSET_SET, set_default_neumann_tags(), moab::GeomTopoTool::set_sense(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), and weave_lateral_face_from_edges().

{

    // MESHSET
    // ErrorCode rval = _mbImpl->create_meshset(MESHSET_ORDERED, new_geo_edge);
    ErrorCode rval = _mbImpl->create_meshset( MESHSET_SET, volume );
    MBERRORR( rval, "can't create volume" );

    int volumeMatId = 1;  // just give a mat id, for debugging, mostly
    Tag matTag;
    rval = _mbImpl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, matTag );
    MBERRORR( rval, "can't get material tag" );

    rval = _mbImpl->tag_set_data( matTag, &volume, 1, &volumeMatId );
    MBERRORR( rval, "can't set material tag value on volume" );

    // get the edges of those 2 faces, and get the vertices of those edges
    // in order, they should be created in the same order (?); check for that, anyway
    rval = _mbImpl->add_parent_child( volume, newFace1 );
    MBERRORR( rval, "can't add first face to volume" );

    rval = _mbImpl->add_parent_child( volume, newFace2 );
    MBERRORR( rval, "can't add second face to volume" );

    // first is bottom, so it is negatively oriented
    rval = _my_geomTopoTool->add_geo_set( volume, 3 );
    MBERRORR( rval, "can't add volume to the gtt" );

    // set senses
    // bottom face is negatively oriented, its normal is toward interior of the volume
    rval = _my_geomTopoTool->set_sense( newFace1, volume, -1 );
    MBERRORR( rval, "can't set bottom face sense to the volume" );

    // the top face is positively oriented
    rval = _my_geomTopoTool->set_sense( newFace2, volume, 1 );
    MBERRORR( rval, "can't set top face sense to the volume" );

    // the children should be in the same direction
    //   get the side edges of each face, and form lateral faces, along direction
    std::vector< EntityHandle > edges1;
    std::vector< EntityHandle > edges2;

    rval = _mbImpl->get_child_meshsets( newFace1, edges1 );  // no hops
    MBERRORR( rval, "can't get children edges or first face, bottom" );

    rval = _mbImpl->get_child_meshsets( newFace2, edges2 );  // no hops
    MBERRORR( rval, "can't get children edges for second face, top" );

    if( edges1.size() != edges2.size() ) MBERRORR( MB_FAILURE, "wrong correspondence " );

    for( unsigned int i = 0; i < edges1.size(); ++i )
    {
        EntityHandle newLatFace;
        rval = weave_lateral_face_from_edges( edges1[i], edges2[i], direction, newLatFace );
        MBERRORR( rval, "can't weave lateral face" );
        rval = _mbImpl->add_parent_child( volume, newLatFace );
        MBERRORR( rval, "can't add lateral face to volume" );

        // set sense as positive
        rval = _my_geomTopoTool->set_sense( newLatFace, volume, 1 );
        MBERRORR( rval, "can't set lateral face sense to the volume" );
    }

    rval = set_default_neumann_tags();
    MBERRORR( rval, "can't set new neumann tags" );

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::createEntSet ( int  isList,
EntityHandle pSet 
)

Definition at line 448 of file FBEngine.cpp.

References ErrorCode, MBI, and MESHSET_SET.

{
    ErrorCode rval;

    if( isList )
        rval = MBI->create_meshset( MESHSET_ORDERED, *pSet );
    else
        rval = MBI->create_meshset( MESHSET_SET, *pSet );

    return rval;
}
ErrorCode moab::FBEngine::createTag ( const char *  tag_name,
int  tag_num_type_values,
int  tag_type,
Tag tag_handle_out 
)

Definition at line 917 of file FBEngine.cpp.

References ErrorCode, MB_ALREADY_ALLOCATED, MB_SUCCESS, MB_TAG_EXCL, MB_TAG_SPARSE, MB_TYPE_DOUBLE, MB_TYPE_HANDLE, MB_TYPE_INTEGER, MB_TYPE_OPAQUE, MBI, and TagType.

{
    // this is copied from iMesh_MOAB.cpp; different name to not have trouble
    // with it
    // also, we do not want to depend on iMesh.h...
    // iMesh is more complicated, because of the options passed

    DataType mb_data_type_table2[] = { MB_TYPE_OPAQUE, MB_TYPE_INTEGER, MB_TYPE_DOUBLE, MB_TYPE_HANDLE,
                                       MB_TYPE_HANDLE };
    moab::TagType storage          = MB_TAG_SPARSE;
    ErrorCode result;

    result =
        MBI->tag_get_handle( tag_name, tag_size, mb_data_type_table2[tag_type], tag_handle_out, storage | MB_TAG_EXCL );

    if( MB_SUCCESS != result )
    {
        std::string msg( "iMesh_createTag: " );
        if( MB_ALREADY_ALLOCATED == result )
        {
            msg += "Tag already exists with name: \"";
            msg += tag_name;
            std::cout << msg << "\n";
        }
        else
        {
            std::cout << "Failed to create tag with name: " << tag_name << "\n";
            return MB_FAILURE;
        }
    }

    // end copy
    return MB_SUCCESS;
}

Definition at line 389 of file FBEngine.cpp.

References _mbImpl, _my_gsets, _smthFace, moab::SmoothFace::append_smooth_tags(), ErrorCode, MB_TAG_NOT_FOUND, MB_TYPE_BIT, MB_TYPE_DOUBLE, size, moab::Interface::tag_delete(), and moab::Interface::tag_get_handle().

{
    // get all tags from database that are created for smooth data, and
    // delete them; it will delete all data associated with them
    // first tags from faces, edges:
    std::vector< Tag > smoothTags;
    int size1 = (int)_my_gsets[2].size();

    for( int i = 0; i < size1; i++ )
    {
        // these 2 will append gradient tag and plane tag
        _smthFace[i]->append_smooth_tags( smoothTags );
    }
    // then , get other tags:
    // "TANGENTS", "MARKER", "CONTROLEDGE", "CONTROLFACE", "CONTROLEDGEFACE"
    Tag tag_handle;
    ErrorCode rval = _mbImpl->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tag_handle );
    if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );

    rval = _mbImpl->tag_get_handle( "MARKER", 1, MB_TYPE_BIT, tag_handle );
    if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );

    rval = _mbImpl->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, tag_handle );
    if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );

    rval = _mbImpl->tag_get_handle( "CONTROLFACE", 18, MB_TYPE_DOUBLE, tag_handle );
    if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );

    rval = _mbImpl->tag_get_handle( "CONTROLEDGEFACE", 27, MB_TYPE_DOUBLE, tag_handle );
    if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );

    // a lot of tags, delete them
    for( unsigned int k = 0; k < smoothTags.size(); k++ )
    {
        // could be a lot of data
        _mbImpl->tag_delete( smoothTags[k] );
    }
}
ErrorCode moab::FBEngine::divide_triangle ( EntityHandle  triangle,
EntityHandle newVertex 
) [private]

Definition at line 3145 of file FBEngine.cpp.

References _mbImpl, _newTriangles, _piercedTriangles, moab::Interface::create_element(), moab::debug_splits, ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Range::insert(), moab::Interface::list_entity(), MB_SUCCESS, MBERRORR, MBTRI, and print_debug_triangle().

Referenced by split_surface_with_direction().

{
    //
    _piercedTriangles.insert( triangle );
    int nnodes                = 0;
    const EntityHandle* conn3 = NULL;
    ErrorCode rval            = _mbImpl->get_connectivity( triangle, conn3, nnodes );
    MBERRORR( rval, "can't get nodes" );
    EntityHandle t1[] = { conn3[0], conn3[1], newVertex };
    EntityHandle t2[] = { conn3[1], conn3[2], newVertex };
    EntityHandle t3[] = { conn3[2], conn3[0], newVertex };
    EntityHandle newTriangle, newTriangle2, newTriangle3;
    rval = _mbImpl->create_element( MBTRI, t1, 3, newTriangle );
    MBERRORR( rval, "can't create triangle" );
    _newTriangles.insert( newTriangle );
    rval = _mbImpl->create_element( MBTRI, t2, 3, newTriangle3 );
    MBERRORR( rval, "can't create triangle" );
    _newTriangles.insert( newTriangle3 );
    rval = _mbImpl->create_element( MBTRI, t3, 3, newTriangle2 );
    MBERRORR( rval, "can't create triangle" );
    _newTriangles.insert( newTriangle2 );

    // create all edges
    std::vector< EntityHandle > edges0;
    rval = _mbImpl->get_adjacencies( &newTriangle, 1, 1, true, edges0 );
    MBERRORR( rval, "can't get new edges" );
    edges0.clear();
    rval = _mbImpl->get_adjacencies( &newTriangle2, 1, 1, true, edges0 );
    MBERRORR( rval, "can't get new edges" );
    if( debug_splits )
    {
        std::cout << "3 triangles formed:\n";
        _mbImpl->list_entity( newTriangle );
        print_debug_triangle( newTriangle );
        _mbImpl->list_entity( newTriangle3 );
        print_debug_triangle( newTriangle3 );
        _mbImpl->list_entity( newTriangle2 );
        print_debug_triangle( newTriangle2 );
        std::cout << "original nodes in tri:\n";
        _mbImpl->list_entity( conn3[0] );
        _mbImpl->list_entity( conn3[1] );
        _mbImpl->list_entity( conn3[2] );
    }
    return MB_SUCCESS;
}
bool moab::FBEngine::find_vertex_set_for_node ( EntityHandle  iNode,
EntityHandle oVertexSet 
) [private]

Definition at line 2847 of file FBEngine.cpp.

References _mbImpl, moab::Range::begin(), moab::Interface::contains_entities(), moab::Range::end(), ErrorCode, GEOM_DIMENSION_TAG_NAME, moab::Interface::get_entities_by_type_and_tag(), MB_SUCCESS, MB_TYPE_INTEGER, MBENTITYSET, and MBI.

Referenced by create_new_gedge(), split_bedge_at_new_mesh_node(), and split_edge_at_mesh_node().

{
    bool found = false;
    Range vertex_sets;

    const int zero               = 0;
    const void* const zero_val[] = { &zero };
    Tag geom_tag;
    ErrorCode rval = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag );
    if( MB_SUCCESS != rval ) return false;
    rval = _mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom_tag, zero_val, 1, vertex_sets );
    if( MB_SUCCESS != rval ) return false;
    // local _gsets, as we might have not updated the local lists
    // see if ends of geo edge generated is in a geo set 0 or not

    for( Range::iterator vsit = vertex_sets.begin(); vsit != vertex_sets.end(); ++vsit )
    {
        EntityHandle vset = *vsit;
        // is the node part of this set?
        if( _mbImpl->contains_entities( vset, &iNode, 1 ) )
        {

            found      = true;
            oVertexSet = vset;
            break;
        }
    }
    return found;
}

Definition at line 183 of file FBEngine.hpp.

References _my_geomTopoTool.

Referenced by FBiGeom_dtor2().

    {
        return this->_my_geomTopoTool;
    }

Definition at line 3263 of file FBEngine.cpp.

References _mbImpl, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_entities_by_type(), MB_SUCCESS, MBEDGE, and MBERRORR.

Referenced by weave_lateral_face_from_edges().

{
    std::vector< EntityHandle > ents;
    ErrorCode rval = _mbImpl->get_entities_by_type( gedge, MBEDGE, ents );
    if( MB_SUCCESS != rval ) return rval;
    if( ents.size() < 1 ) return MB_FAILURE;

    nodes.resize( ents.size() + 1 );
    const EntityHandle* conn = NULL;
    int len;
    for( unsigned int i = 0; i < ents.size(); ++i )
    {
        rval = _mbImpl->get_connectivity( ents[i], conn, len );
        MBERRORR( rval, "can't get edge connectivity" );
        nodes[i] = conn[0];
    }
    // the last one is conn[1]
    nodes[ents.size()] = conn[1];
    return MB_SUCCESS;
}

Definition at line 3788 of file FBEngine.cpp.

References _mbImpl, children, moab::Range::clear(), entities, ErrorCode, moab::Interface::get_child_meshsets(), MB_SUCCESS, MBEDGE, MBERRORR, MBI, MBVERTEX, and moab::Range::size().

Referenced by chain_able_edge(), and chain_two_edges().

{
    // need to decide first or second vertex
    // important for moab

    Range children;
    // EntityHandle v1, v2;
    ErrorCode rval = _mbImpl->get_child_meshsets( edge, children );
    MBERRORR( rval, "can't get child meshsets" );
    if( children.size() == 1 )
    {
        // this is periodic edge, get out early
        v1 = children[0];
        v2 = v1;
        return MB_SUCCESS;
    }
    else if( children.size() > 2 )
        MBERRORR( MB_FAILURE, "too many vertices in one edge" );
    // edge: get one vertex as part of the vertex set
    Range entities;
    rval = MBI->get_entities_by_type( children[0], MBVERTEX, entities );
    MBERRORR( rval, "can't get entities from vertex set" );
    if( entities.size() < 1 ) MBERRORR( MB_FAILURE, "no mesh nodes in vertex set" );
    EntityHandle node0 = entities[0];  // the first vertex
    entities.clear();

    // now get the edges, and get the first node and the last node in sequence of edges
    // the order is important...
    // these are ordered sets !!
    std::vector< EntityHandle > ents;
    rval = MBI->get_entities_by_type( edge, MBEDGE, ents );
    MBERRORR( rval, "can't get mesh edges" );
    if( ents.size() < 1 ) MBERRORR( MB_FAILURE, "no mesh edges in edge set" );

    const EntityHandle* conn = NULL;
    int len;
    rval = MBI->get_connectivity( ents[0], conn, len );
    MBERRORR( rval, "can't connectivity of first mesh edge" );

    if( conn[0] == node0 )
    {
        v1 = children[0];
        v2 = children[1];
    }
    else  // the other way around, although we should check (we are paranoid)
    {
        v2 = children[0];
        v1 = children[1];
    }

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getAdjacentEntities ( const EntityHandle  from,
const int  to_dim,
Range adj_ents 
) [private]

Definition at line 850 of file FBEngine.cpp.

References _my_geomTopoTool, moab::Range::clear(), moab::Range::end(), ErrorCode, moab::Range::find(), moab::GeomTopoTool::geoRanges(), MB_SUCCESS, MBI, and moab::subtract().

Referenced by boundary_mesh_edges_on_face(), boundary_nodes_on_face(), getEntAdj(), separate(), split_bedge_at_new_mesh_node(), split_boundary(), and split_edge_at_mesh_node().

{
    int this_dim = -1;
    for( int i = 0; i < 4; i++ )
    {
        if( _my_geomTopoTool->geoRanges()[i].find( from ) != _my_geomTopoTool->geoRanges()[i].end() )
        {
            this_dim = i;
            break;
        }
    }

    // check target dimension
    if( -1 == this_dim )
    {
        // ProcessError(iBase_FAILURE, "Entity not a geometry entity.");
        return MB_FAILURE;
    }
    else if( 0 > to_dim || 3 < to_dim )
    {
        // ProcessError(iBase_FAILURE, "To dimension must be between 0 and 3.");
        return MB_FAILURE;
    }
    else if( to_dim == this_dim )
    {
        // ProcessError(iBase_FAILURE,
        //      "To dimension must be different from entity dimension.");
        return MB_FAILURE;
    }

    ErrorCode rval = MB_SUCCESS;
    adjs.clear();
    if( to_dim > this_dim )
    {
        int diffDim = to_dim - this_dim;
        rval        = MBI->get_parent_meshsets( from, adjs, diffDim );
        if( MB_SUCCESS != rval ) return rval;
        if( diffDim > 1 )
        {
            // subtract the parents that come with diffDim-1 hops
            Range extra;
            rval = MBI->get_parent_meshsets( from, extra, diffDim - 1 );
            if( MB_SUCCESS != rval ) return rval;
            adjs = subtract( adjs, extra );
        }
    }
    else
    {
        int diffDim = this_dim - to_dim;
        rval        = MBI->get_child_meshsets( from, adjs, diffDim );
        if( MB_SUCCESS != rval ) return rval;
        if( diffDim > 1 )
        {
            // subtract the children that come with diffDim-1 hops
            Range extra;
            rval = MBI->get_child_meshsets( from, extra, diffDim - 1 );
            if( MB_SUCCESS != rval ) return rval;
            adjs = subtract( adjs, extra );
        }
    }

    return rval;
}
ErrorCode moab::FBEngine::getArrData ( const moab::EntityHandle entity_handles,
int  entity_handles_size,
Tag  tag_handle,
void *  tag_values_out 
)

Definition at line 952 of file FBEngine.cpp.

References MBI.

{
    // responsibility of the user to have tag_values_out properly allocated
    // only some types of Tags are possible (double, int, etc)
    return MBI->tag_get_data( tag_handle, entity_handles, entity_handles_size, tag_values_out );
}
ErrorCode moab::FBEngine::getEgEvalXYZ ( EntityHandle  edge,
double  x,
double  y,
double  z,
double &  on_x,
double &  on_y,
double &  on_z,
double &  tngt_i,
double &  tngt_j,
double &  tngt_k,
double &  cvtr_i,
double &  cvtr_j,
double &  cvtr_k 
)

Definition at line 1063 of file FBEngine.cpp.

References MB_NOT_IMPLEMENTED.

Referenced by FBiGeom_getEgEvalXYZ().

{
    return MB_NOT_IMPLEMENTED;  // not implemented
}
ErrorCode moab::FBEngine::getEgFcSense ( EntityHandle  mbedge,
EntityHandle  mbface,
int &  sense 
)

Definition at line 977 of file FBEngine.cpp.

References _my_geomTopoTool, ErrorCode, moab::GeomTopoTool::get_senses(), and MB_SUCCESS.

Referenced by FBiGeom_getEgFcSense(), and FBiGeom_getEntNrmlSense().

{

    // this one is important, for establishing the orientation of the edges in faces
    // use senses
    std::vector< EntityHandle > faces;
    std::vector< int > senses;  // 0 is forward and 1 is backward
    ErrorCode rval = _my_geomTopoTool->get_senses( mbedge, faces, senses );
    if( MB_SUCCESS != rval ) return rval;

    for( unsigned int i = 0; i < faces.size(); i++ )
    {
        if( faces[i] == mbface )
        {
            sense_out = senses[i];
            return MB_SUCCESS;
        }
    }
    return MB_FAILURE;
}
ErrorCode moab::FBEngine::getEgVtxSense ( EntityHandle  edge,
EntityHandle  vtx1,
EntityHandle  vtx2,
int &  sense 
)

Definition at line 1099 of file FBEngine.cpp.

References moab::Range::clear(), entities, ErrorCode, getEntType(), MB_SUCCESS, MBEDGE, MBI, MBVERTEX, and moab::Range::size().

Referenced by FBiGeom_getEgVtxSense(), and weave_lateral_face_from_edges().

{
    // need to decide first or second vertex
    // important for moab
    int type;

    EntityHandle v1, v2;
    ErrorCode rval = getEntType( vtx1, &type );
    if( MB_SUCCESS != rval || type != 0 ) return MB_FAILURE;
    // edge: get one vertex as part of the vertex set
    Range entities;
    rval = MBI->get_entities_by_type( vtx1, MBVERTEX, entities );
    if( MB_SUCCESS != rval ) return rval;
    if( entities.size() < 1 ) return MB_FAILURE;
    v1 = entities[0];  // the first vertex
    entities.clear();
    rval = getEntType( vtx2, &type );
    if( MB_SUCCESS != rval || type != 0 ) return MB_FAILURE;
    rval = MBI->get_entities_by_type( vtx2, MBVERTEX, entities );
    if( MB_SUCCESS != rval ) return rval;
    if( entities.size() < 1 ) return MB_FAILURE;
    v2 = entities[0];  // the first vertex
    entities.clear();
    // now get the edges, and get the first node and the last node in sequence of edges
    // the order is important...
    // these are ordered sets !!
    std::vector< EntityHandle > ents;
    rval = MBI->get_entities_by_type( edge, MBEDGE, ents );
    if( MB_SUCCESS != rval ) return rval;
    if( ents.size() < 1 ) return MB_FAILURE;

    const EntityHandle* conn = NULL;
    int len;
    EntityHandle startNode, endNode;
    rval = MBI->get_connectivity( ents[0], conn, len );
    if( MB_SUCCESS != rval ) return rval;
    startNode = conn[0];
    rval      = MBI->get_connectivity( ents[ents.size() - 1], conn, len );
    if( MB_SUCCESS != rval ) return rval;

    endNode = conn[1];
    sense   = 1;  //
    if( ( startNode == endNode ) && ( v1 == startNode ) )
    {
        sense = 0;  // periodic
    }
    if( ( startNode == v1 ) && ( endNode == v2 ) )
    {
        sense = 1;  // forward
    }
    if( ( startNode == v2 ) && ( endNode == v1 ) )
    {
        sense = -1;  // reverse
    }
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getEntAdj ( EntityHandle  handle,
int  type_requested,
Range adjEnts 
)

Definition at line 972 of file FBEngine.cpp.

References getAdjacentEntities().

Referenced by FBiGeom_getArrAdj(), FBiGeom_getEnt2ndAdj(), and FBiGeom_getEntAdj().

{
    return getAdjacentEntities( handle, type_requested, adjEnts );
}
ErrorCode moab::FBEngine::getEntBoundBox ( EntityHandle  this_gent,
double *  x0,
double *  y0,
double *  z0,
double *  x1,
double *  y1,
double *  z1 
)

Definition at line 549 of file FBEngine.cpp.

References _my_geomTopoTool, moab::CartVect::array(), moab::OrientedBoxTreeTool::box(), center(), ErrorCode, moab::GeomTopoTool::get_root(), getEntType(), getVtxCoord(), MB_SUCCESS, MBERRORR, and moab::GeomTopoTool::obb_tree().

{
    ErrorCode rval;
    int type;
    rval = getEntType( gent, &type );
    MBERRORR( rval, "Failed to get entity type." );

    if( type == 0 )
    {
        rval = getVtxCoord( gent, min_x, min_y, min_z );
        MBERRORR( rval, "Failed to get vertex coordinates." );
        max_x = min_x;
        max_y = min_y;
        max_z = min_z;
    }
    else if( type == 1 )
    {
        MBERRORR( MB_FAILURE, "iGeom_getEntBoundBox is not supported for Edge entity type." );
    }
    else if( type == 2 || type == 3 )
    {

        EntityHandle root;
        CartVect center, axis[3];
        rval = _my_geomTopoTool->get_root( gent, root );
        MBERRORR( rval, "Failed to get tree root in iGeom_getEntBoundBox." );
        rval = _my_geomTopoTool->obb_tree()->box( root, center.array(), axis[0].array(), axis[1].array(),
                                                  axis[2].array() );
        MBERRORR( rval, "Failed to get closest point in iGeom_getEntBoundBox." );

        CartVect absv[3];
        for( int i = 0; i < 3; i++ )
        {
            absv[i] = CartVect( fabs( axis[i][0] ), fabs( axis[i][1] ), fabs( axis[i][2] ) );
        }
        CartVect min, max;
        min    = center - absv[0] - absv[1] - absv[2];
        max    = center + absv[0] + absv[1] + absv[2];
        *min_x = min[0];
        *min_y = min[1];
        *min_z = min[2];
        *max_x = max[0];
        *max_y = max[1];
        *max_z = max[2];
    }
    else
        return MB_FAILURE;

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getEntClosestPt ( EntityHandle  this_gent,
double  x,
double  y,
double  z,
double *  x1,
double *  y1,
double *  y3 
)

Definition at line 605 of file FBEngine.cpp.

References _edges, _faces, _my_geomTopoTool, _smooth, moab::OrientedBoxTreeTool::closest_to_location(), ErrorCode, moab::GeomTopoTool::get_root(), getEntType(), getVtxCoord(), MB_SUCCESS, MB_TYPE_OUT_OF_RANGE, MBERRORR, moab::SmoothCurve::move_to_curve(), moab::SmoothFace::move_to_surface(), and moab::GeomTopoTool::obb_tree().

Referenced by FBiGeom_getEntClosestPt(), and FBiGeom_getFcEvalXYZ().

{
    ErrorCode rval;
    int type;
    rval = getEntType( this_gent, &type );
    MBERRORR( rval, "Failed to get entity type." );

    if( type == 0 )
    {
        rval = getVtxCoord( this_gent, on_x, on_y, on_z );
        MBERRORR( rval, "Failed to get vertex coordinates." );
    }
    else if( _smooth && type == 1 )
    {
        *on_x                  = near_x;
        *on_y                  = near_y;
        *on_z                  = near_z;
        SmoothCurve* smthcurve = _edges[this_gent];
        // call the new method from smooth edge
        smthcurve->move_to_curve( *on_x, *on_y, *on_z );
    }
    else if( type == 2 || type == 3 )
    {
        double point[3] = { near_x, near_y, near_z };
        double point_out[3];
        EntityHandle root, facet_out;
        if( _smooth && 2 == type )
        {
            SmoothFace* smthFace = _faces[this_gent];
            *on_x                = near_x;
            *on_y                = near_y;
            *on_z                = near_z;
            smthFace->move_to_surface( *on_x, *on_y, *on_z );
        }
        else
        {
            rval = _my_geomTopoTool->get_root( this_gent, root );
            MBERRORR( rval, "Failed to get tree root in iGeom_getEntClosestPt." );
            rval = _my_geomTopoTool->obb_tree()->closest_to_location( point, root, point_out, facet_out );
            MBERRORR( rval, "Failed to get closest point in iGeom_getEntClosestPt." );

            *on_x = point_out[0];
            *on_y = point_out[1];
            *on_z = point_out[2];
        }
    }
    else
        return MB_TYPE_OUT_OF_RANGE;

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getEntities ( EntityHandle  root_set,
int  ent_type,
Range gentities 
)

Definition at line 460 of file FBEngine.cpp.

References _my_geomTopoTool, moab::Range::clear(), ErrorCode, moab::GeomTopoTool::geoRanges(), moab::intersect(), MB_SUCCESS, MBENTITYSET, MBERRORR, MBI, and moab::Range::merge().

Referenced by FBiGeom_getEntities().

{
    int i;
    if( 0 > entity_type || 4 < entity_type )
    {
        return MB_FAILURE;
    }
    else if( entity_type < 4 )
    {                                                            // 4 means all entities
        gentities = _my_geomTopoTool->geoRanges()[entity_type];  // all from root set!
    }
    else
    {
        gentities.clear();
        for( i = 0; i < 4; i++ )
        {
            gentities.merge( _my_geomTopoTool->geoRanges()[i] );
        }
    }
    Range sets;
    // see now if they are in the set passed as input or not
    ErrorCode rval = MBI->get_entities_by_type( set_handle, MBENTITYSET, sets );
    MBERRORR( rval, "can't get sets in the initial set" );
    gentities = intersect( gentities, sets );

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getEntNrmlSense ( EntityHandle  face,
EntityHandle  region,
int &  sense 
)

Definition at line 1058 of file FBEngine.cpp.

References MB_NOT_IMPLEMENTED.

{
    return MB_NOT_IMPLEMENTED;  // not implemented
}
ErrorCode moab::FBEngine::getEntNrmlXYZ ( EntityHandle  entity_handle,
double  x,
double  y,
double  z,
double *  nrml_i,
double *  nrml_j,
double *  nrml_k 
)

Definition at line 713 of file FBEngine.cpp.

References _faces, _my_geomTopoTool, _smooth, moab::OrientedBoxTreeTool::closest_to_location(), ErrorCode, moab::GeomTopoTool::get_root(), getEntType(), MB_SUCCESS, MBERRORR, MBI, moab::SmoothFace::normal_at(), moab::CartVect::normalize(), and moab::GeomTopoTool::obb_tree().

Referenced by FBiGeom_getEntNrmlXYZ(), and FBiGeom_getFcEvalXYZ().

{
    // just do for surface and volume
    int type;
    ErrorCode rval = getEntType( entity_handle, &type );
    MBERRORR( rval, "Failed to get entity type in iGeom_getEntNrmlXYZ." );

    if( type != 2 && type != 3 )
    {
        MBERRORR( MB_FAILURE, "Entities passed into gentityNormal must be face or volume." );
    }

    if( _smooth && 2 == type )
    {
        SmoothFace* smthFace = _faces[entity_handle];
        //*on_x = near_x; *on_y = near_y; *on_z = near_z;
        smthFace->normal_at( x, y, z, *nrml_i, *nrml_j, *nrml_k );
    }
    else
    {
        // get closest location and facet
        double point[3] = { x, y, z };
        double point_out[3];
        EntityHandle root, facet_out;
        _my_geomTopoTool->get_root( entity_handle, root );
        rval = _my_geomTopoTool->obb_tree()->closest_to_location( point, root, point_out, facet_out );
        MBERRORR( rval, "Failed to get closest location in iGeom_getEntNrmlXYZ." );

        // get facet normal
        const EntityHandle* conn;
        int len;
        CartVect coords[3], normal;
        rval = MBI->get_connectivity( facet_out, conn, len );
        MBERRORR( rval, "Failed to get triangle connectivity in iGeom_getEntNrmlXYZ." );
        if( len != 3 ) MBERRORR( MB_FAILURE, " not a triangle, error " );

        rval = MBI->get_coords( conn, len, coords[0].array() );
        MBERRORR( rval, "Failed to get triangle coordinates in iGeom_getEntNrmlXYZ." );

        coords[1] -= coords[0];
        coords[2] -= coords[0];
        normal = coords[1] * coords[2];
        normal.normalize();
        *nrml_i = normal[0];
        *nrml_j = normal[1];
        *nrml_k = normal[2];
    }
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getEntTgntU ( EntityHandle  edge,
double  u,
double &  i,
double &  j,
double &  k 
)

Definition at line 1172 of file FBEngine.cpp.

References _edges, MB_SUCCESS, and moab::SmoothCurve::position_from_u().

Referenced by FBiGeom_getEntTgntU().

{
    SmoothCurve* smoothCurve = _edges[edge];  // this is a map
    // now, call smoothCurve methods
    double tg[3];
    double x, y, z;
    smoothCurve->position_from_u( u, x, y, z, tg );
    i = tg[0];
    j = tg[1];
    k = tg[2];
    return MB_SUCCESS;
}

Definition at line 536 of file FBEngine.cpp.

References _my_geomTopoTool, moab::Range::end(), moab::Range::find(), moab::GeomTopoTool::geoRanges(), and MB_SUCCESS.

Referenced by FBiGeom_getEntType(), getEgVtxSense(), getEntBoundBox(), getEntClosestPt(), getEntNrmlXYZ(), getVtxCoord(), isEntAdj(), and measure().

{
    for( int i = 0; i < 4; i++ )
    {
        if( _my_geomTopoTool->geoRanges()[i].find( gent ) != _my_geomTopoTool->geoRanges()[i].end() )
        {
            *type = i;
            return MB_SUCCESS;
        }
    }
    *type = -1;  // failure
    return MB_FAILURE;
}
ErrorCode moab::FBEngine::getEntURange ( EntityHandle  edge,
double &  u_min,
double &  u_max 
)

Definition at line 1156 of file FBEngine.cpp.

References _edges, moab::SmoothCurve::get_param_range(), and MB_SUCCESS.

Referenced by FBiGeom_getEntURange().

{
    SmoothCurve* smoothCurve = _edges[edge];  // this is a map
    // now, call smoothCurve methods
    smoothCurve->get_param_range( u_min, u_max );
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getEntUtoXYZ ( EntityHandle  edge,
double  u,
double &  x,
double &  y,
double &  z 
)

Definition at line 1164 of file FBEngine.cpp.

References _edges, MB_SUCCESS, and moab::SmoothCurve::position_from_u().

Referenced by FBiGeom_getEntUtoXYZ().

{
    SmoothCurve* smoothCurve = _edges[edge];  // this is a map
    // now, call smoothCurve methods
    smoothCurve->position_from_u( u, x, y, z );
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getFcEvalXYZ ( EntityHandle  face,
double  x,
double  y,
double  z,
double &  on_x,
double &  on_y,
double &  on_z,
double &  nrml_i,
double &  nrml_j,
double &  nrml_k,
double &  cvtr1_i,
double &  cvtr1_j,
double &  cvtr1_k,
double &  cvtr2_i,
double &  cvtr2_j,
double &  cvtr2_k 
)

Definition at line 1079 of file FBEngine.cpp.

References MB_NOT_IMPLEMENTED.

{
    return MB_NOT_IMPLEMENTED;  // not implemented
}
ErrorCode moab::FBEngine::getNumEntSets ( EntityHandle  set,
int  num_hops,
int *  all_sets 
)

Definition at line 442 of file FBEngine.cpp.

References ErrorCode, and MBI.

{
    ErrorCode rval = MBI->num_contained_meshsets( set, all_sets, num_hops + 1 );
    return rval;
}
ErrorCode moab::FBEngine::getNumOfType ( EntityHandle  set,
int  ent_type,
int *  pNum 
)

Definition at line 498 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, ErrorCode, GEOM_DIMENSION_TAG_NAME, moab::GeomTopoTool::geoRanges(), moab::Interface::get_entities_by_type_and_tag(), moab::intersect(), MB_SUCCESS, MB_TAG_CREAT, MB_TAG_SPARSE, MB_TYPE_INTEGER, MBENTITYSET, MBERRORR, moab::Range::size(), moab::Interface::tag_get_handle(), and moab::Interface::UNION.

Referenced by FBiGeom_getNumOfType().

{
    if( 0 > ent_type || 4 < ent_type )
    {
        std::cout << "Invalid type\n";
        return MB_FAILURE;
    }
    // get sets of geom dimension tag from here, and intersect with the gentities from geo
    // ranges

    // get the geom dimensions sets in the set (AKA gentities)
    Range geom_sets;
    Tag geom_tag;
    ErrorCode rval =
        _mbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
    MBERRORR( rval, "Failed to get geom tag." );
    rval = _mbImpl->get_entities_by_type_and_tag( set, MBENTITYSET, &geom_tag, NULL, 1, geom_sets, Interface::UNION );
    MBERRORR( rval, "Failed to get gentities from set" );

    if( ent_type == 4 )
    {
        *pNum = 0;
        for( int k = 0; k <= 3; k++ )
        {
            Range gEntsOfTypeK = intersect( geom_sets, _my_geomTopoTool->geoRanges()[k] );
            *pNum += (int)gEntsOfTypeK.size();
        }
    }
    else
    {
        Range gEntsOfType = intersect( geom_sets, _my_geomTopoTool->geoRanges()[ent_type] );
        *pNum             = (int)gEntsOfType.size();
    }
    // we do not really check if it is in the set or not;
    // _my_gsets[i].find(gent) != _my_gsets[i].end()
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getPntRayIntsct ( double  x,
double  y,
double  z,
double  dir_x,
double  dir_y,
double  dir_z,
std::vector< EntityHandle > &  intersect_entity_handles,
std::vector< double > &  intersect_coords,
std::vector< double > &  param_coords 
)

Definition at line 769 of file FBEngine.cpp.

References _faces, _my_geomTopoTool, _my_gsets, _smooth, ErrorCode, moab::GeomTopoTool::get_root(), MB_SUCCESS, MBERRORR, moab::GeomTopoTool::obb_tree(), moab::OrientedBoxTreeTool::ray_intersect_sets(), moab::SmoothFace::ray_intersection_correct(), moab::Range::size(), and moab::tolerance.

Referenced by FBiGeom_getPntRayIntsct().

{
    // this is pretty cool
    // we will return only surfaces (gentities )
    //
    ErrorCode rval;

    unsigned int numfaces = _my_gsets[2].size();
    // do ray fire
    const double point[] = { x, y, z };
    const double dir[]   = { dir_x, dir_y, dir_z };
    CartVect P( point );
    CartVect V( dir );

    // std::vector<double> distances;
    std::vector< EntityHandle > facets;
    // std::vector<EntityHandle> sets;
    unsigned int i;
    for( i = 0; i < numfaces; i++ )
    {
        EntityHandle face = _my_gsets[2][i];
        EntityHandle rootForFace;
        rval = _my_geomTopoTool->get_root( face, rootForFace );
        MBERRORR( rval, "Failed to get root of face." );
        std::vector< double > distances_out;
        std::vector< EntityHandle > sets_out;
        std::vector< EntityHandle > facets_out;
        rval = _my_geomTopoTool->obb_tree()->ray_intersect_sets( distances_out, sets_out, facets_out, rootForFace,
                                                                 tolerance, point, dir );
        unsigned int j;
        for( j = 0; j < distances_out.size(); j++ )
            param_coords.push_back( distances_out[j] );
        for( j = 0; j < sets_out.size(); j++ )
            intersect_entity_handles.push_back( sets_out[j] );
        for( j = 0; j < facets_out.size(); j++ )
            facets.push_back( facets_out[j] );

        MBERRORR( rval, "Failed to get ray intersections." );
    }
    // facets.size == distances.size()!!
    for( i = 0; i < param_coords.size(); i++ )
    {
        CartVect intx = P + param_coords[i] * V;
        for( int j = 0; j < 3; j++ )
            intersect_coords.push_back( intx[j] );
    }
    if( _smooth )
    {
        // correct the intersection point and the distance for smooth surfaces
        for( i = 0; i < intersect_entity_handles.size(); i++ )
        {
            // EntityHandle geoSet = MBH_cast(sets[i]);
            SmoothFace* sFace = _faces[intersect_entity_handles[i]];
            // correct coordinates and distance from point
            /*moab::ErrorCode ray_intersection_correct(moab::EntityHandle facet, // (IN) the facet
             where the patch is defined moab::CartVect &pt, // (IN) shoot from moab::CartVect &ray,
             // (IN) ray direction moab::CartVect &eval_pt, // (INOUT) The intersection point double
             & distance, // (IN OUT) the new distance bool &outside);*/
            CartVect pos( &( intersect_coords[3 * i] ) );
            double dist  = param_coords[i];
            bool outside = false;
            rval         = sFace->ray_intersection_correct( facets[i], P, V, pos, dist, outside );
            MBERRORR( rval, "Failed to get better point on ray." );
            param_coords[i] = dist;

            for( int j = 0; j < 3; j++ )
                intersect_coords[3 * i + j] = pos[j];
        }
    }
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::getVtxCoord ( EntityHandle  this_gent,
double *  x0,
double *  y0,
double *  z0 
)

Definition at line 663 of file FBEngine.cpp.

References entities, ErrorCode, getEntType(), MB_SUCCESS, MBERRORR, MBI, MBVERTEX, and moab::Range::size().

Referenced by FBiGeom_getVtxCoord(), getEntBoundBox(), and getEntClosestPt().

{
    int type;
    ErrorCode rval = getEntType( vertex_handle, &type );
    MBERRORR( rval, "Failed to get entity type in getVtxCoord." );

    if( type != 0 )
    {
        MBERRORR( MB_FAILURE, "Entity is not a vertex type." );
    }

    Range entities;
    rval = MBI->get_entities_by_type( vertex_handle, MBVERTEX, entities );
    MBERRORR( rval, "can't get nodes in vertex set." );

    if( entities.size() != 1 )
    {
        MBERRORR( MB_FAILURE, "Vertex has multiple points." );
    }
    double coords[3];
    EntityHandle node = entities[0];
    rval              = MBI->get_coords( &node, 1, coords );
    MBERRORR( rval, "can't get coordinates." );
    *x0 = coords[0];
    *y0 = coords[1];
    *z0 = coords[2];

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::gsubtract ( EntityHandle  entity_set_1,
EntityHandle  entity_set_2,
EntityHandle  result_entity_set 
)

Definition at line 693 of file FBEngine.cpp.

References ErrorCode, MBENTITYSET, MBERRORR, MBI, and moab::subtract().

{
    /*result_entity_set = subtract(entity_set_1, entity_set_2);*/
    Range ents1, ents2;
    ErrorCode rval = MBI->get_entities_by_type( entity_set_1, MBENTITYSET, ents1 );
    MBERRORR( rval, "can't get entities from set 1." );

    rval = MBI->get_entities_by_type( entity_set_2, MBENTITYSET, ents2 );
    MBERRORR( rval, "can't get entities from set 2." );

    ents1 = subtract( ents1, ents2 );
    rval  = MBI->clear_meshset( &result_entity_set, 1 );
    MBERRORR( rval, "can't empty set." );

    rval = MBI->add_entities( result_entity_set, ents1 );
    MBERRORR( rval, "can't add result to set." );

    return rval;
}

Definition at line 237 of file FBEngine.cpp.

References _initialized, _my_geomTopoTool, _my_gsets, _smooth, moab::GeomTopoTool::construct_obb_trees(), ErrorCode, moab::GeomTopoTool::find_geomsets(), initializeSmoothing(), MB_SUCCESS, and split_quads().

Referenced by FBiGeom_load(), and FBiGeom_newGeomFromMesh().

{
    if( !_initialized )
    {
        if( !_my_geomTopoTool ) return MB_FAILURE;

        ErrorCode rval = _my_geomTopoTool->find_geomsets( _my_gsets );
        assert( rval == MB_SUCCESS );
        if( MB_SUCCESS != rval )
        {
            return rval;
        }

        rval = split_quads();
        assert( rval == MB_SUCCESS );

        rval = _my_geomTopoTool->construct_obb_trees();
        assert( rval == MB_SUCCESS );

        if( _smooth ) rval = initializeSmoothing();
        assert( rval == MB_SUCCESS );

        _initialized = true;
    }
    return MB_SUCCESS;
}

Definition at line 263 of file FBEngine.cpp.

References _edges, _faces, _my_geomTopoTool, _my_gsets, _smthCurve, _smthFace, moab::SmoothCurve::compute_control_points_on_boundary_edges(), moab::SmoothFace::compute_control_points_on_edges(), moab::SmoothFace::compute_internal_control_points_on_facets(), moab::SmoothFace::compute_tangents_for_each_edge(), moab::SmoothCurve::compute_tangents_for_each_edge(), moab::Debug_surf_eval, moab::Range::end(), ErrorCode, GEOM_DIMENSION_TAG_NAME, moab::SmoothFace::init_gradient(), MB_SUCCESS, MB_TAG_BIT, MB_TAG_CREAT, MB_TAG_DENSE, MB_TAG_EXCL, MB_TYPE_BIT, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MBERRORR, MBI, and moab::Range::size().

Referenced by Init().

{
    //
    /*ErrorCode rval = Init();
     MBERRORR(rval, "failed initialize");*/
    // first of all, we need to retrieve all the surfaces from the (root) set
    // in icesheet_test we use iGeom, but maybe that is a stretch
    // get directly the sets with geom dim 2, and from there create the SmoothFace
    Tag geom_tag;
    ErrorCode rval = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag );
    MBERRORR( rval, "can't get geom tag" );

    int numSurfaces = _my_gsets[2].size();
    // SmoothFace ** smthFace = new SmoothFace *[numSurfaces];
    _smthFace = new SmoothFace*[numSurfaces];

    // there should also be a map from surfaces to evaluators
    // std::map<MBEntityHandle, SmoothFace*> mapSurfaces;

    int i = 0;
    Range::iterator it;
    for( it = _my_gsets[2].begin(); it != _my_gsets[2].end(); ++it, i++ )
    {
        EntityHandle face = *it;
        _smthFace[i] = new SmoothFace( MBI, face, _my_geomTopoTool );  // geom topo tool will be used for searching,
        // among other things; also for senses in edge sets...
        _faces[face] = _smthFace[i];
    }

    int numCurves = _my_gsets[1].size();  // csets.size();
    // SmoothCurve ** smthCurve = new SmoothCurve *[numCurves];
    _smthCurve = new SmoothCurve*[numCurves];
    // there should also be a map from surfaces to evaluators
    // std::map<MBEntityHandle, SmoothCurve*> mapCurves;

    i = 0;
    for( it = _my_gsets[1].begin(); it != _my_gsets[1].end(); ++it, i++ )
    {
        EntityHandle curve = *it;
        _smthCurve[i]      = new SmoothCurve( MBI, curve, _my_geomTopoTool );
        _edges[curve]      = _smthCurve[i];
    }

    for( i = 0; i < numSurfaces; i++ )
    {
        _smthFace[i]->init_gradient();                   // this will also retrieve the triangles in each surface
        _smthFace[i]->compute_tangents_for_each_edge();  // this one will consider all edges
                                                         // internal, so the
        // tangents are all in the direction of the edge; a little bit of waste, as we store
        // one tangent for each edge node , even though they are equal here...
        // no loops are considered
    }

    // this will be used to mark boundary edges, so for them the control points are computed earlier
    unsigned char value = 0;  // default value is "not used"=0 for the tag
    // unsigned char def_data_bit = 1;// valid by default
    // rval = mb->tag_create("valid", 1, MB_TAG_BIT, validTag, &def_data_bit);
    Tag markTag;
    rval = MBI->tag_get_handle( "MARKER", 1, MB_TYPE_BIT, markTag, MB_TAG_EXCL | MB_TAG_BIT,
                                &value );  // default value : 0 = not computed yet
    // each feature edge will need to have a way to retrieve at every moment the surfaces it belongs
    // to from edge sets, using the sense tag, we can get faces, and from each face, using the map,
    // we can get the SmoothFace (surface evaluator), that has everything, including the normals!!!
    assert( rval == MB_SUCCESS );

    // create the tag also for control points on the edges
    double defCtrlPoints[9] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. };
    Tag edgeCtrlTag;
    rval = MBI->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, edgeCtrlTag, MB_TAG_DENSE | MB_TAG_CREAT,
                                &defCtrlPoints );
    assert( rval == MB_SUCCESS );

    Tag facetCtrlTag;
    double defControls[18] = { 0. };
    rval = MBI->tag_get_handle( "CONTROLFACE", 18, MB_TYPE_DOUBLE, facetCtrlTag, MB_TAG_CREAT | MB_TAG_DENSE,
                                &defControls );
    assert( rval == MB_SUCCESS );

    Tag facetEdgeCtrlTag;
    double defControls2[27] = { 0. };  // corresponding to 9 control points on edges, in order
                                       // from edge 0, 1, 2 ( 1-2, 2-0, 0-1 )
    rval = MBI->tag_get_handle( "CONTROLEDGEFACE", 27, MB_TYPE_DOUBLE, facetEdgeCtrlTag, MB_TAG_CREAT | MB_TAG_DENSE,
                                &defControls2 );
    assert( rval == MB_SUCCESS );
    // if the
    double min_dot = -1.0;  // depends on _angle, but now we ignore it, for the time being
    for( i = 0; i < numCurves; i++ )
    {
        _smthCurve[i]->compute_tangents_for_each_edge();  // do we need surfaces now? or just the chains?
        // the computed edges will be marked accordingly; later one, only internal edges to surfaces
        // are left
        _smthCurve[i]->compute_control_points_on_boundary_edges( min_dot, _faces, edgeCtrlTag, markTag );
    }

    // when done with boundary edges, compute the control points on all edges in the surfaces

    for( i = 0; i < numSurfaces; i++ )
    {
        // also pass the tags for
        _smthFace[i]->compute_control_points_on_edges( min_dot, edgeCtrlTag, markTag );
    }

    // now we should be able to compute the control points for the facets

    for( i = 0; i < numSurfaces; i++ )
    {
        // also pass the tags for edge and facet control points
        _smthFace[i]->compute_internal_control_points_on_facets( min_dot, facetCtrlTag, facetEdgeCtrlTag );
    }
    // we will need to compute the tangents for all edges in the model
    // they will be needed for control points for each edge
    // the boundary edges and the feature edges are more complicated
    // the boundary edges need to consider their loops, but feature edges need to consider loops and
    // the normals on each connected surface

    // some control points
    if( Debug_surf_eval )
        for( i = 0; i < numSurfaces; i++ )
            _smthFace[i]->DumpModelControlPoints();

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::isEntAdj ( EntityHandle  entity1,
EntityHandle  entity2,
bool &  adjacent_out 
)

Definition at line 1184 of file FBEngine.cpp.

References moab::Range::end(), ErrorCode, moab::Range::find(), getEntType(), MB_SUCCESS, and MBI.

Referenced by FBiGeom_isEntAdj().

{
    int type1, type2;
    ErrorCode rval = getEntType( entity1, &type1 );
    if( MB_SUCCESS != rval ) return rval;
    rval = getEntType( entity2, &type2 );
    if( MB_SUCCESS != rval ) return rval;

    Range adjs;
    if( type1 < type2 )
    {
        rval = MBI->get_parent_meshsets( entity1, adjs, type2 - type1 );
        if( MB_SUCCESS != rval ) return rval;  // MBERRORR("Failed to get parent meshsets in iGeom_isEntAdj.");
    }
    else
    {
        // note: if they ave the same type, they will not be adjacent, in our definition
        rval = MBI->get_child_meshsets( entity1, adjs, type1 - type2 );
        if( MB_SUCCESS != rval ) return rval;  // MBERRORR("Failed to get child meshsets in iGeom_isEntAdj.");
    }

    // adjacent_out = adjs.find(entity2) != _my_gsets[type2].end();
    // hmmm, possible bug here; is this called?
    adjacent_out = adjs.find( entity2 ) != adjs.end();

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::measure ( const EntityHandle moab_entities,
int  entities_size,
double *  measures 
)

Definition at line 998 of file FBEngine.cpp.

References moab::Range::begin(), moab::Range::end(), entities, ErrorCode, getEntType(), moab::CartVect::length(), MB_SUCCESS, MBEDGE, MBI, and MBTRI.

Referenced by FBiGeom_measure().

{
    ErrorCode rval;
    for( int i = 0; i < entities_size; i++ )
    {
        measures[i] = 0.;

        int type;
        EntityHandle gset = moab_entities[i];
        rval              = getEntType( gset, &type );
        if( MB_SUCCESS != rval ) return rval;
        if( type == 1 )
        {  // edge: get all edges part of the edge set
            Range entities;
            rval = MBI->get_entities_by_type( gset, MBEDGE, entities );
            if( MB_SUCCESS != rval ) return rval;

            for( Range::iterator it = entities.begin(); it != entities.end(); ++it )
            {
                EntityHandle edge = *it;
                CartVect vv[2];
                const EntityHandle* conn2 = NULL;
                int num_nodes;
                rval = MBI->get_connectivity( edge, conn2, num_nodes );
                if( MB_SUCCESS != rval || num_nodes != 2 ) return MB_FAILURE;
                rval = MBI->get_coords( conn2, 2, (double*)&( vv[0][0] ) );
                if( MB_SUCCESS != rval ) return rval;

                vv[0] = vv[1] - vv[0];
                measures[i] += vv[0].length();
            }
        }
        if( type == 2 )
        {  // surface
            // get triangles in surface; TODO: quads!
            Range entities;
            rval = MBI->get_entities_by_type( gset, MBTRI, entities );
            if( MB_SUCCESS != rval ) return rval;

            for( Range::iterator it = entities.begin(); it != entities.end(); ++it )
            {
                EntityHandle tri = *it;
                CartVect vv[3];
                const EntityHandle* conn3 = NULL;
                int num_nodes;
                rval = MBI->get_connectivity( tri, conn3, num_nodes );
                if( MB_SUCCESS != rval || num_nodes != 3 ) return MB_FAILURE;
                rval = MBI->get_coords( conn3, 3, (double*)&( vv[0][0] ) );
                if( MB_SUCCESS != rval ) return rval;

                vv[1] = vv[1] - vv[0];
                vv[2] = vv[2] - vv[0];
                vv[0] = vv[1] * vv[2];
                measures[i] += vv[0].length() / 2;  // area of triangle
            }
        }
    }
    return MB_SUCCESS;
}

Definition at line 91 of file FBEngine.hpp.

References _mbImpl.

    {
        return _mbImpl;
    }
void moab::FBEngine::print_debug_triangle ( EntityHandle  triangle) [private]

Definition at line 1967 of file FBEngine.cpp.

References _mbImpl, moab::Interface::get_connectivity(), moab::Interface::get_coords(), and moab::Interface::id_from_handle().

Referenced by BreakTriangle2(), divide_triangle(), and split_internal_edge().

{
    std::cout << " triangle id:" << _mbImpl->id_from_handle( t ) << "\n";
    const EntityHandle* conn3 = NULL;
    int nnodes                = 0;
    _mbImpl->get_connectivity( t, conn3, nnodes );
    // get coords
    CartVect P[3];
    _mbImpl->get_coords( conn3, 3, (double*)&P[0] );
    std::cout << "  nodes:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n";
    CartVect PP[3];
    PP[0] = P[1] - P[0];
    PP[1] = P[2] - P[1];
    PP[2] = P[0] - P[2];

    std::cout << "  pos:" << P[0] << " " << P[1] << " " << P[2] << "\n";
    std::cout << "   x,y diffs " << PP[0][0] << " " << PP[0][1] << ",  " << PP[1][0] << " " << PP[1][1] << ",  "
              << PP[2][0] << " " << PP[2][1] << "\n";
    return;
}
ErrorCode moab::FBEngine::redistribute_boundary_edges_to_faces ( EntityHandle  face,
EntityHandle  newFace,
std::vector< EntityHandle > &  chainedEdges 
) [private]

Definition at line 2876 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, moab::Interface::add_parent_child(), moab::Range::begin(), children, moab::GeomTopoTool::dimension(), moab::Range::empty(), moab::Range::end(), ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_child_meshsets(), moab::Interface::get_entities_by_dimension(), moab::Interface::get_entities_by_handle(), moab::intersect(), MB_SUCCESS, MBERRORR, moab::Interface::remove_parent_child(), moab::GeomTopoTool::set_sense(), and moab::Interface::side_number().

Referenced by split_surface().

{

    // so far, original boundary edges are all parent/child relations for face
    // we should get them all, and see if they are truly adjacent to face or newFace
    // decide also on the orientation/sense with respect to the triangles
    Range r1;  // range in old face
    Range r2;  // range of tris in new face
    ErrorCode rval = _mbImpl->get_entities_by_dimension( face, 2, r1 );
    MBERRORR( rval, " can't get triangles from old face" );
    rval = _mbImpl->get_entities_by_dimension( newFace, 2, r2 );
    MBERRORR( rval, " can't get triangles from new face" );
    // right now, all boundary edges are children of face
    // we need to get them all, and verify if they indeed are adjacent to face
    Range children;
    rval = _mbImpl->get_child_meshsets( face, children );  // only direct children are of interest
    MBERRORR( rval, " can't get children sets from face" );

    for( Range::iterator it = children.begin(); it != children.end(); ++it )
    {
        EntityHandle edge = *it;
        if( std::find( chainedEdges.begin(), chainedEdges.end(), edge ) != chainedEdges.end() )
            continue;  // we already set this one fine
        // get one mesh edge from the edge; we have to get all of them!!
        if( _my_geomTopoTool->dimension( edge ) != 1 ) continue;  // not an edge
        Range mesh_edges;
        rval = _mbImpl->get_entities_by_handle( edge, mesh_edges );
        MBERRORR( rval, " can't get mesh edges from edge" );
        if( mesh_edges.empty() ) MBERRORR( MB_FAILURE, " no mesh edges" );
        EntityHandle mesh_edge = mesh_edges[0];  // first one is enough
        // get its triangles; see which one is in r1 or r2; it should not be in both
        Range triangles;
        rval = _mbImpl->get_adjacencies( &mesh_edge, 1, 2, false, triangles );
        MBERRORR( rval, " can't get adjacent triangles" );
        Range intx1 = intersect( triangles, r1 );
        Range intx2 = intersect( triangles, r2 );
        if( !intx1.empty() && !intx2.empty() ) MBERRORR( MB_FAILURE, " at least one should be empty" );

        if( intx2.empty() )
        {
            // it means it is in the range r1; the sense should be fine, no need to reset
            // the sense should have been fine, also
            continue;
        }
        // so it must be a triangle in r2;
        EntityHandle triangle = intx2[0];  // one triangle only
        // remove the edge from boundary of face, and add it to the boundary of newFace
        // remove_parent_child( EntityHandle parent,  EntityHandle child )
        rval = _mbImpl->remove_parent_child( face, edge );
        MBERRORR( rval, " can't remove parent child relation for edge" );
        // add to the new face
        rval = _mbImpl->add_parent_child( newFace, edge );
        MBERRORR( rval, " can't add parent child relation for edge" );

        // set some sense, based on the sense of the mesh_edge in triangle
        int num1, sense, offset;
        rval = _mbImpl->side_number( triangle, mesh_edge, num1, sense, offset );
        MBERRORR( rval, "mesh edge not adjacent to triangle" );

        rval = this->_my_geomTopoTool->set_sense( edge, newFace, sense );
        MBERRORR( rval, "can't set proper sense of edge in face" );
    }

    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::separate ( EntityHandle  face,
std::vector< EntityHandle > &  chainedEdges,
Range first,
Range second 
) [private]

Definition at line 1680 of file FBEngine.cpp.

References _mbImpl, _newTriangles, _piercedTriangles, moab::Interface::add_entities(), moab::Range::begin(), moab::Range::clear(), moab::Interface::create_meshset(), moab::debug_splits, moab::Range::end(), ErrorCode, moab::Range::find(), moab::Interface::get_adjacencies(), moab::Interface::get_entities_by_dimension(), moab::Interface::get_entities_by_type(), getAdjacentEntities(), moab::Interface::id_from_handle(), MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_INTEGER, MBEDGE, MBERRORR, MBI, MBTRI, MESHSET_SET, moab::Interface::side_number(), moab::Range::size(), moab::subtract(), moab::Interface::tag_delete(), moab::Interface::tag_get_data(), moab::Interface::tag_set_data(), moab::Interface::UNION, moab::unite(), and moab::Interface::write_file().

Referenced by split_surface().

{
    // Range unaffectedTriangles = subtract(iniTriangles, _piercedTriangles);
    // insert in each
    // start with a new triangle, and flood to get the first range; what is left is the
    // second range
    // flood fill is considering edges adjacent to seed triangles; if there is
    //  an edge in the new_geo_edge, it is skipped; triangles in the
    // triangles to delete are not added
    // first, create all edges of the new triangles

    //
    // new face will have the new edge oriented positively
    // get mesh edges from geo edge (splitting gedge);

    Range mesh_edges;
    ErrorCode rval;
    // mesh_edges
    for( unsigned int j = 0; j < chainedEdges.size(); j++ )
    {
        // this will keep adding edges to the mesh_edges range
        // when we get out, the mesh_edges will be in this range, but not ordered
        rval = _mbImpl->get_entities_by_type( chainedEdges[j], MBEDGE, mesh_edges );
        MBERRORR( rval, "can't get new polyline edges" );
        if( debug_splits )
        {
            std::cout << " At chained edge " << j << " " << _mbImpl->id_from_handle( chainedEdges[j] )
                      << " mesh_edges Range size:" << mesh_edges.size() << "\n";
        }
    }

    // get a positive triangle adjacent to mesh_edge[0]
    // add to first triangles to the left, second triangles to the right of the mesh_edges ;

    // create a temp tag, and when done, delete it
    // default value: 0
    // 3 to be deleted, pierced
    // 1 first set
    // 2 second set
    // create the tag also for control points on the edges
    int defVal = 0;
    Tag separateTag;
    rval = MBI->tag_get_handle( "SEPARATE_TAG", 1, MB_TYPE_INTEGER, separateTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );
    MBERRORR( rval, "can't create temp tag for separation" );
    // the deleted triangles will get a value 3, from start
    int delVal = 3;
    for( Range::iterator it1 = this->_piercedTriangles.begin(); it1 != _piercedTriangles.end(); ++it1 )
    {
        EntityHandle trToDelete = *it1;
        rval                    = _mbImpl->tag_set_data( separateTag, &trToDelete, 1, &delVal );
        MBERRORR( rval, "can't set delete tag value" );
    }

    // find a triangle that will be in the first range, positively oriented about the splitting edge
    EntityHandle seed1 = 0;
    for( Range::iterator it = mesh_edges.begin(); it != mesh_edges.end() && !seed1; ++it )
    {
        EntityHandle meshEdge = *it;
        Range adj_tri;
        rval = _mbImpl->get_adjacencies( &meshEdge, 1, 2, false, adj_tri );
        MBERRORR( rval, "can't get adj_tris to mesh edge" );

        for( Range::iterator it2 = adj_tri.begin(); it2 != adj_tri.end(); ++it2 )
        {
            EntityHandle tr = *it2;
            if( _piercedTriangles.find( tr ) != _piercedTriangles.end() )
                continue;  // do not attach pierced triangles, they are not good
            int num1, sense, offset;
            rval = _mbImpl->side_number( tr, meshEdge, num1, sense, offset );
            MBERRORR( rval, "edge not adjacent" );
            if( sense == 1 )
            {
                // firstSet.insert(tr);
                if( !seed1 )
                {
                    seed1 = tr;
                    break;
                }
            }
        }
    }

    // flood fill first set, the rest will be in second set
    // the edges from new_geo_edge will not be crossed

    // get edges of face (adjacencies)
    // also get the old boundary edges, from face; they will be edges to not cross, too
    Range bound_edges;
    rval = getAdjacentEntities( face, 1, bound_edges );
    MBERRORR( rval, "can't get boundary edges" );

    // add to the do not cross edges range, all edges from initial boundary
    Range initialBoundaryEdges;
    for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
    {
        EntityHandle bound_edge = *it;
        rval                    = _mbImpl->get_entities_by_dimension( bound_edge, 1, initialBoundaryEdges );
    }

    Range doNotCrossEdges = unite( initialBoundaryEdges, mesh_edges );  // add the splitting edges !

    // use a second method, with tags
    //
    std::queue< EntityHandle > queue1;
    queue1.push( seed1 );
    std::vector< EntityHandle > arr1;
    while( !queue1.empty() )
    {
        // start copy
        EntityHandle currentTriangle = queue1.front();
        queue1.pop();
        arr1.push_back( currentTriangle );
        // add new triangles that share an edge
        Range currentEdges;
        rval = _mbImpl->get_adjacencies( &currentTriangle, 1, 1, true, currentEdges, Interface::UNION );
        MBERRORR( rval, "can't get adjacencies" );
        for( Range::iterator it = currentEdges.begin(); it != currentEdges.end(); ++it )
        {
            EntityHandle frontEdge = *it;
            if( doNotCrossEdges.find( frontEdge ) == doNotCrossEdges.end() )
            {
                // this is an edge that can be crossed
                Range adj_tri;
                rval = _mbImpl->get_adjacencies( &frontEdge, 1, 2, false, adj_tri, Interface::UNION );
                MBERRORR( rval, "can't get adj_tris" );
                // if the triangle is not in first range, add it to the queue
                for( Range::iterator it2 = adj_tri.begin(); it2 != adj_tri.end(); ++it2 )
                {
                    EntityHandle tri2 = *it2;
                    int val           = 0;
                    rval              = _mbImpl->tag_get_data( separateTag, &tri2, 1, &val );
                    MBERRORR( rval, "can't get tag value" );
                    if( val ) continue;
                    // else, set it to 1
                    val  = 1;
                    rval = _mbImpl->tag_set_data( separateTag, &tri2, 1, &val );
                    MBERRORR( rval, "can't get tag value" );

                    queue1.push( tri2 );
                }
            }  // end edge do not cross
        }      // end while
    }

    std::sort( arr1.begin(), arr1.end() );
    // Range first1;
    std::copy( arr1.rbegin(), arr1.rend(), range_inserter( first ) );

    // std::cout<< "\n first1.size() " << first1.size() << " first.size(): " << first.size() <<
    // "\n";
    if( debug_splits )
    {
        EntityHandle tmpSet;
        _mbImpl->create_meshset( MESHSET_SET, tmpSet );
        _mbImpl->add_entities( tmpSet, first );
        _mbImpl->write_file( "dbg1.vtk", "vtk", 0, &tmpSet, 1 );
    }
    // now, decide the set 2:
    // first, get all ini tris
    Range initr;
    rval = _mbImpl->get_entities_by_type( face, MBTRI, initr );
    MBERRORR( rval, "can't get tris " );
    second        = unite( initr, _newTriangles );
    Range second2 = subtract( second, _piercedTriangles );
    second        = subtract( second2, first );
    _newTriangles.clear();
    if( debug_splits )
    {
        std::cout << "\n second.size() " << second.size() << " first.size(): " << first.size() << "\n";
        // debugging code
        EntityHandle tmpSet2;
        _mbImpl->create_meshset( MESHSET_SET, tmpSet2 );
        _mbImpl->add_entities( tmpSet2, second );
        _mbImpl->write_file( "dbg2.vtk", "vtk", 0, &tmpSet2, 1 );
    }
    /*Range intex = intersect(first, second);
    if (!intex.empty() && debug_splits)
    {
      std::cout << "error, the sets should be disjoint\n";
      for (Range::iterator it1=intex.begin(); it1!=intex.end(); ++it1)
      {
        std::cout<<_mbImpl->id_from_handle(*it1) << "\n";
      }
    }*/
    rval = _mbImpl->tag_delete( separateTag );
    MBERRORR( rval, "can't delete tag " );
    return MB_SUCCESS;
}

Definition at line 3602 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, ErrorCode, moab::GeomTopoTool::find_geomsets(), MB_SUCCESS, MB_TYPE_INTEGER, MBERRORR, NEUMANN_SET_TAG_NAME, size, moab::Interface::tag_get_handle(), and moab::Interface::tag_set_data().

Referenced by create_volume_with_direction().

{
    // these are for debugging purposes only
    // check the initial tag, then
    Tag ntag;
    ErrorCode rval = _mbImpl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, ntag );
    MBERRORR( rval, "can't get tag handle" );
    // get all surfaces in the model now
    Range sets[5];
    rval = _my_geomTopoTool->find_geomsets( sets );
    MBERRORR( rval, "can't get geo sets" );
    int nfaces = (int)sets[2].size();
    int* vals  = new int[nfaces];
    for( int i = 0; i < nfaces; i++ )
    {
        vals[i] = i + 1;
    }
    rval = _mbImpl->tag_set_data( ntag, sets[2], (void*)vals );
    MBERRORR( rval, "can't set tag values for neumann sets" );

    delete[] vals;

    return MB_SUCCESS;
}

Definition at line 2944 of file FBEngine.cpp.

References _mbImpl, ErrorCode, MB_SUCCESS, MB_TYPE_INTEGER, MBERRORR, NEUMANN_SET_TAG_NAME, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), and moab::Interface::tag_set_data().

Referenced by split_surface().

{
    // these are for debugging purposes only
    // check the initial tag, then
    Tag ntag;
    ErrorCode rval = _mbImpl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, ntag );
    MBERRORR( rval, "can't get tag handle" );
    // check the value for face
    int nval;
    rval = _mbImpl->tag_get_data( ntag, &face, 1, &nval );
    if( MB_SUCCESS == rval )
    {
        nval++;
    }
    else
    {
        nval = 1;
        rval = _mbImpl->tag_set_data( ntag, &face, 1, &nval );
        MBERRORR( rval, "can't set tag" );
        nval = 2;
    }
    rval = _mbImpl->tag_set_data( ntag, &newFace, 1, &nval );
    MBERRORR( rval, "can't set tag" );

    return MB_SUCCESS;
}
void moab::FBEngine::set_smooth ( ) [inline]

Definition at line 218 of file FBEngine.hpp.

References _smooth.

Referenced by FBiGeom_load().

    {
        _smooth = true;
    }
ErrorCode moab::FBEngine::setArrData ( const EntityHandle entity_handles,
int  entity_handles_size,
Tag  tag_handle,
const void *  tag_values 
)

Definition at line 962 of file FBEngine.cpp.

References MBI.

{
    // responsibility of the user to have tag_values_out properly allocated
    // only some types of Tags are possible (double, int, etc)
    return MBI->tag_set_data( tag_handle, entity_handles, entity_handles_size, tag_values );
}
ErrorCode moab::FBEngine::smooth_new_intx_points ( EntityHandle  face,
std::vector< EntityHandle > &  chainedEdges 
) [private]

Definition at line 1626 of file FBEngine.cpp.

References _faces, _mbImpl, moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::get_entities_by_type(), MB_SUCCESS, MBERRORR, MBTRI, moab::SmoothFace::move_to_surface(), moab::Interface::set_coords(), moab::Range::size(), and moab::subtract().

Referenced by split_surface().

{

    // do not move nodes from the original face
    // first get all triangles, and then all nodes from those triangles

    Range tris;
    ErrorCode rval = _mbImpl->get_entities_by_type( face, MBTRI, tris );
    MBERRORR( rval, "can't get triangles" );

    Range ini_nodes;
    rval = _mbImpl->get_connectivity( tris, ini_nodes );
    MBERRORR( rval, "can't get connectivities" );

    SmoothFace* smthFace = _faces[face];

    // get all nodes from chained edges
    Range mesh_edges;
    for( unsigned int j = 0; j < chainedEdges.size(); j++ )
    {
        // keep adding to the range of mesh edges
        rval = _mbImpl->get_entities_by_dimension( chainedEdges[j], 1, mesh_edges );
        MBERRORR( rval, "can't get mesh edges" );
    }
    // nodes along polyline
    Range nodes_on_polyline;
    rval = _mbImpl->get_connectivity( mesh_edges, nodes_on_polyline, true );  // corners only
    MBERRORR( rval, "can't get nodes on the polyline" );

    Range new_intx_nodes = subtract( nodes_on_polyline, ini_nodes );

    std::vector< double > ini_coords;
    int num_points = (int)new_intx_nodes.size();
    ini_coords.resize( 3 * num_points );
    rval = _mbImpl->get_coords( new_intx_nodes, &( ini_coords[0] ) );
    MBERRORR( rval, "can't get coordinates" );

    int i = 0;
    for( Range::iterator it = new_intx_nodes.begin(); it != new_intx_nodes.end(); ++it )
    {
        /*EntityHandle node = *it;*/
        int i3 = 3 * i;
        smthFace->move_to_surface( ini_coords[i3], ini_coords[i3 + 1], ini_coords[i3 + 2] );
        // reset the coordinates of this node
        ++i;
    }
    rval = _mbImpl->set_coords( new_intx_nodes, &( ini_coords[0] ) );
    MBERRORR( rval, "can't set smoothed coordinates for the new nodes" );

    return MB_SUCCESS;
}

Definition at line 2565 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, moab::Interface::add_entities(), moab::GeomTopoTool::add_geo_set(), moab::Interface::add_parent_child(), moab::Range::begin(), moab::Interface::create_meshset(), moab::debug_splits, dim, moab::GeomTopoTool::dimension(), moab::Range::end(), ErrorCode, find_vertex_set_for_node(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_entities_by_type(), moab::GeomTopoTool::get_sense(), getAdjacentEntities(), moab::GeomTopoTool::global_id(), moab::Interface::id_from_handle(), MB_SUCCESS, MBEDGE, MBERRORR, MBI, MBVERTEX, MESHSET_SET, moab::Interface::remove_entities(), moab::Interface::remove_parent_child(), moab::Interface::set_connectivity(), moab::GeomTopoTool::set_sense(), and moab::Range::size().

Referenced by split_boundary().

{
    // the node should be in the list of nodes

    int dim = _my_geomTopoTool->dimension( edge );
    if( dim != 1 ) return MB_FAILURE;  // not an edge

    if( debug_splits )
    {
        std::cout << "Split edge " << _mbImpl->id_from_handle( edge )
                  << " with global id: " << _my_geomTopoTool->global_id( edge )
                  << " at new node:" << _mbImpl->id_from_handle( node ) << "breaking mesh edge"
                  << _mbImpl->id_from_handle( brokenEdge ) << "\n";
    }

    // now get the edges
    // the order is important...
    // these are ordered sets !!
    std::vector< EntityHandle > ents;
    ErrorCode rval = _mbImpl->get_entities_by_type( edge, MBEDGE, ents );
    if( MB_SUCCESS != rval ) return rval;
    if( ents.size() < 1 ) return MB_FAILURE;  // no edges

    const EntityHandle* conn = NULL;
    int len;
    // the edge connected to the splitting node is brokenEdge
    // find the small edges it is broken into, which are connected to
    // new vertex (node) and its ends; also, if necessary, reorient these small edges
    // for proper orientation
    rval = _mbImpl->get_connectivity( brokenEdge, conn, len );
    MBERRORR( rval, "Failed to get connectivity of broken edge" );

    // we already created the new edges, make sure they are oriented fine; if not, revert them
    EntityHandle conn02[] = { conn[0], node };  // first node and new node
    // default, intersect
    std::vector< EntityHandle > adj_edges;
    rval = _mbImpl->get_adjacencies( conn02, 2, 1, false, adj_edges );
    if( adj_edges.size() < 1 || rval != MB_SUCCESS ) MBERRORR( MB_FAILURE, " Can't find small split edge" );

    // get this edge connectivity;
    EntityHandle firstEdge         = adj_edges[0];
    const EntityHandle* connActual = NULL;
    rval                           = _mbImpl->get_connectivity( firstEdge, connActual, len );
    MBERRORR( rval, "Failed to get connectivity of first split edge" );
    // if it is the same as conn02, we are happy
    if( conn02[0] != connActual[0] )
    {
        // reset connectivity of edge
        rval = _mbImpl->set_connectivity( firstEdge, conn02, 2 );
        MBERRORR( rval, "Failed to reset connectivity of first split edge" );
    }
    //  now treat the second edge
    adj_edges.clear();                          //
    EntityHandle conn21[] = { node, conn[1] };  //  new node and second node
    rval                  = _mbImpl->get_adjacencies( conn21, 2, 1, false, adj_edges );
    if( adj_edges.size() < 1 || rval != MB_SUCCESS ) MBERRORR( MB_FAILURE, " Can't find second small split edge" );

    // get this edge connectivity;
    EntityHandle secondEdge = adj_edges[0];
    rval                    = _mbImpl->get_connectivity( firstEdge, connActual, len );
    MBERRORR( rval, "Failed to get connectivity of first split edge" );
    // if it is the same as conn21, we are happy
    if( conn21[0] != connActual[0] )
    {
        // reset connectivity of edge
        rval = _mbImpl->set_connectivity( secondEdge, conn21, 2 );
        MBERRORR( rval, "Failed to reset connectivity of first split edge" );
    }

    int num_mesh_edges = (int)ents.size();
    int index_edge;  // this is the index of the edge that will be removed (because it is split)
    // the rest of edges will be put in order in the (remaining) edge and new edge

    for( index_edge = 0; index_edge < num_mesh_edges; index_edge++ )
        if( brokenEdge == ents[index_edge] ) break;
    if( index_edge >= num_mesh_edges ) MBERRORR( MB_FAILURE, "can't find the broken edge" );

    // so the edges up to index_edge and firstEdge, will form the "old" edge
    // the edges secondEdge and from index_edge+1 to end will form the new_edge

    // here, we have 0 ... index_edge edges in the first set,
    // create a vertex set and add the node to it

    if( debug_splits )
    {
        std::cout << "Split edge with " << num_mesh_edges << " mesh edges, at index (0 based) " << index_edge << "\n";
    }

    // at this moment, the node set should have been already created in new_geo_edge
    EntityHandle nodeSet;  // the node set that has the node (find it!)
    bool found = find_vertex_set_for_node( node, nodeSet );

    if( !found )
    {
        // create a node set and add the node

        // must be an error, but create one nevertheless
        rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet );
        MBERRORR( rval, "Failed to create a new vertex set" );

        rval = _mbImpl->add_entities( nodeSet, &node, 1 );
        MBERRORR( rval, "Failed to add the node to the set" );

        rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 );  //
        MBERRORR( rval, "Failed to commit the node set" );

        if( debug_splits )
        {
            std::cout << " create a vertex set (this should have been created before!)"
                      << _mbImpl->id_from_handle( nodeSet )
                      << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << "\n";
        }
    }

    // we need to remove the remaining mesh edges from first set, and add it
    // to the second set, in order

    rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_edge );
    MBERRORR( rval, "can't create geo edge" );

    int remaining = num_mesh_edges - 1 - index_edge;

    // add edges to the new edge set
    rval = _mbImpl->add_entities( new_edge, &secondEdge, 1 );  // add the second split edge to new edge
    MBERRORR( rval, "can't add second split edge to the new edge" );
    // then add the rest
    rval = _mbImpl->add_entities( new_edge, &ents[index_edge + 1], remaining );
    MBERRORR( rval, "can't add edges to the new edge" );

    // also, remove the second node set from old edge
    // remove the edges index_edge and up

    rval = _mbImpl->remove_entities( edge, &ents[index_edge], remaining + 1 );  // include the
    MBERRORR( rval, "can't remove edges from the old edge" );
    // add the firstEdge too
    rval = _mbImpl->add_entities( edge, &firstEdge, 1 );  // add the second split edge to new edge
    MBERRORR( rval, "can't add first split edge to the old edge" );

    // need to find the adjacent vertex sets
    Range vertexRange;
    rval = getAdjacentEntities( edge, 0, vertexRange );

    EntityHandle secondSet;
    if( vertexRange.size() == 1 )
    {
        // initially a periodic edge, OK to add the new set to both edges, and the
        // second set
        secondSet = vertexRange[0];
    }
    else
    {
        if( vertexRange.size() > 2 ) return MB_FAILURE;  // something must be wrong with too many vertices
        // find first node
        EntityHandle firstNode;

        rval = MBI->get_connectivity( ents[0], conn, len );
        if( MB_SUCCESS != rval ) return rval;
        firstNode = conn[0];  // this is the first node of the initial edge
                              // we will use it to identify the vertex set associated with it
        int k;
        for( k = 0; k < 2; k++ )
        {
            Range verts;
            rval = _mbImpl->get_entities_by_type( vertexRange[k], MBVERTEX, verts );

            MBERRORR( rval, "can't get vertices from vertex set" );
            if( verts.size() != 1 ) MBERRORR( MB_FAILURE, " node set not defined well" );
            if( firstNode == verts[0] )
            {
                secondSet = vertexRange[1 - k];  // the other set; it is 1 or 0
                break;
            }
        }
        if( k >= 2 )
        {
            // it means we didn't find the right node set
            MBERRORR( MB_FAILURE, " can't find the right vertex" );
        }
        // remove the second set from the connectivity with the
        //  edge (parent-child relation)
        // remove_parent_child( EntityHandle parent,
        //                                          EntityHandle child )
        rval = _mbImpl->remove_parent_child( edge, secondSet );
        MBERRORR( rval, " can't remove the second vertex from edge" );
    }
    // at this point, we just need to add to both edges the new node set (vertex)
    rval = _mbImpl->add_parent_child( edge, nodeSet );
    MBERRORR( rval, " can't add new vertex to old edge" );

    rval = _mbImpl->add_parent_child( new_edge, nodeSet );
    MBERRORR( rval, " can't add new vertex to new edge" );

    // one thing that I forgot: add the secondSet as a child to new edge!!!
    // (so far, the new edge has only one end vertex!)
    rval = _mbImpl->add_parent_child( new_edge, secondSet );
    MBERRORR( rval, " can't add second vertex to new edge" );

    // now, add the edge and vertex to geo tool

    rval = _my_geomTopoTool->add_geo_set( new_edge, 1 );
    MBERRORR( rval, " can't add new edge" );

    // next, get the adjacent faces to initial edge, and add them as parents
    // to the new edge

    // need to find the adjacent face sets
    Range faceRange;
    rval = getAdjacentEntities( edge, 2, faceRange );

    // these faces will be adjacent to the new edge too!
    // need to set the senses of new edge within faces

    for( Range::iterator it = faceRange.begin(); it != faceRange.end(); ++it )
    {
        EntityHandle face = *it;
        rval              = _mbImpl->add_parent_child( face, new_edge );
        MBERRORR( rval, " can't add new edge - face parent relation" );
        int sense;
        rval = _my_geomTopoTool->get_sense( edge, face, sense );
        MBERRORR( rval, " can't get initial sense of edge in the adjacent face" );
        // keep the same sense for the new edge within the faces
        rval = _my_geomTopoTool->set_sense( new_edge, face, sense );
        MBERRORR( rval, " can't set sense of new edge in the adjacent face" );
    }

    return MB_SUCCESS;
}

Definition at line 2796 of file FBEngine.cpp.

References _brokenEdges, _mbImpl, moab::Range::begin(), moab::debug_splits, moab::Range::end(), ErrorCode, moab::Range::find(), moab::Interface::get_connectivity(), moab::Interface::get_entities_by_dimension(), getAdjacentEntities(), moab::Interface::id_from_handle(), MBERRORR, split_bedge_at_new_mesh_node(), and split_edge_at_mesh_node().

Referenced by split_surface().

{
    // find the boundary edges, and split the one that we find it is a part of
    if( debug_splits )
    {
        std::cout << "Split face " << _mbImpl->id_from_handle( face )
                  << " at node:" << _mbImpl->id_from_handle( atNode ) << "\n";
    }
    Range bound_edges;
    ErrorCode rval = getAdjacentEntities( face, 1, bound_edges );
    MBERRORR( rval, " can't get boundary edges" );
    bool brokEdge = _brokenEdges.find( atNode ) != _brokenEdges.end();

    for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
    {
        EntityHandle b_edge = *it;
        // get all edges in range
        Range mesh_edges;
        rval = _mbImpl->get_entities_by_dimension( b_edge, 1, mesh_edges );
        MBERRORR( rval, " can't get mesh edges" );
        if( brokEdge )
        {
            EntityHandle brokenEdge = _brokenEdges[atNode];
            if( mesh_edges.find( brokenEdge ) != mesh_edges.end() )
            {
                EntityHandle new_edge;
                rval = split_bedge_at_new_mesh_node( b_edge, atNode, brokenEdge, new_edge );
                return rval;
            }
        }
        else
        {
            Range nodes;
            rval = _mbImpl->get_connectivity( mesh_edges, nodes );
            MBERRORR( rval, " can't get nodes from mesh edges" );

            if( nodes.find( atNode ) != nodes.end() )
            {
                // we found our boundary edge candidate
                EntityHandle new_edge;
                rval = split_edge_at_mesh_node( b_edge, atNode, new_edge );
                return rval;
            }
        }
    }
    // if the node was not found in any "current" boundary, it broke an existing
    // boundary edge
    MBERRORR( MB_FAILURE, " we did not find an appropriate boundary edge" );
    ;  //
}

Definition at line 2373 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, moab::Interface::add_entities(), moab::GeomTopoTool::add_geo_set(), moab::Interface::add_parent_child(), moab::Range::begin(), moab::Interface::create_meshset(), moab::debug_splits, dim, moab::GeomTopoTool::dimension(), moab::Range::end(), ErrorCode, find_vertex_set_for_node(), moab::Interface::get_entities_by_type(), moab::GeomTopoTool::get_sense(), getAdjacentEntities(), moab::GeomTopoTool::global_id(), moab::Interface::id_from_handle(), MB_SUCCESS, MBEDGE, MBERRORR, MBI, MBVERTEX, MESHSET_SET, moab::Interface::remove_entities(), moab::Interface::remove_parent_child(), moab::GeomTopoTool::set_sense(), and moab::Range::size().

Referenced by split_boundary(), and split_edge_at_point().

{
    // the node should be in the list of nodes

    int dim = _my_geomTopoTool->dimension( edge );
    if( dim != 1 ) return MB_FAILURE;  // not an edge

    if( debug_splits )
    {
        std::cout << "Split edge " << _mbImpl->id_from_handle( edge )
                  << " with global id: " << _my_geomTopoTool->global_id( edge )
                  << " at node:" << _mbImpl->id_from_handle( node ) << "\n";
    }

    // now get the edges
    // the order is important...
    // these are ordered sets !!
    std::vector< EntityHandle > ents;
    ErrorCode rval = _mbImpl->get_entities_by_type( edge, MBEDGE, ents );
    if( MB_SUCCESS != rval ) return rval;
    if( ents.size() < 1 ) return MB_FAILURE;  // no edges

    const EntityHandle* conn = NULL;
    int len;
    // find the edge connected to the splitting node
    int num_mesh_edges = (int)ents.size();
    int index_edge;
    EntityHandle firstNode = 0;
    for( index_edge = 0; index_edge < num_mesh_edges; index_edge++ )
    {
        rval = MBI->get_connectivity( ents[index_edge], conn, len );
        if( MB_SUCCESS != rval ) return rval;
        if( index_edge == 0 ) firstNode = conn[0];  // will be used to decide vertex sets adjacencies
        if( conn[0] == node )
        {
            if( index_edge == 0 )
            {
                new_edge = 0;       // no need to split, it is the first node
                return MB_SUCCESS;  // no split
            }
            else
                return MB_FAILURE;  // we should have found the node already , wrong
                                    // connectivity
        }
        else if( conn[1] == node )
        {
            // we found the index of interest
            break;
        }
    }
    if( index_edge == num_mesh_edges - 1 )
    {
        new_edge = 0;       // no need to split, it is the last node
        return MB_SUCCESS;  // no split
    }

    // here, we have 0 ... index_edge edges in the first set,
    // create a vertex set and add the node to it

    if( debug_splits )
    {
        std::cout << "Split edge with " << num_mesh_edges << " mesh edges, at index (0 based) " << index_edge << "\n";
    }

    // at this moment, the node set should have been already created in new_geo_edge
    EntityHandle nodeSet;  // the node set that has the node (find it!)
    bool found = find_vertex_set_for_node( node, nodeSet );

    if( !found )
    {
        // create a node set and add the node

        // must be an error, but create one nevertheless
        rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet );
        MBERRORR( rval, "Failed to create a new vertex set" );

        rval = _mbImpl->add_entities( nodeSet, &node, 1 );
        MBERRORR( rval, "Failed to add the node to the set" );

        rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 );  //
        MBERRORR( rval, "Failed to commit the node set" );

        if( debug_splits )
        {
            std::cout << " create a vertex set (this should have been created before!)"
                      << _mbImpl->id_from_handle( nodeSet )
                      << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << "\n";
        }
    }

    // we need to remove the remaining mesh edges from first set, and add it
    // to the second set, in order

    rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_edge );
    MBERRORR( rval, "can't create geo edge" );

    int remaining = num_mesh_edges - 1 - index_edge;

    // add edges to the edge set
    rval = _mbImpl->add_entities( new_edge, &ents[index_edge + 1], remaining );
    MBERRORR( rval, "can't add edges to the new edge" );

    // also, remove the second node set from old edge
    // remove the edges index_edge+1 and up

    rval = _mbImpl->remove_entities( edge, &ents[index_edge + 1], remaining );
    MBERRORR( rval, "can't remove edges from the old edge" );

    // need to find the adjacent vertex sets
    Range vertexRange;
    rval = getAdjacentEntities( edge, 0, vertexRange );

    EntityHandle secondSet;
    if( vertexRange.size() == 1 )
    {
        // initially a periodic edge, OK to add the new set to both edges, and the
        // second set
        secondSet = vertexRange[0];
    }
    else
    {
        if( vertexRange.size() > 2 ) return MB_FAILURE;  // something must be wrong with too many vertices
        // find first node
        int k;
        for( k = 0; k < 2; k++ )
        {
            Range verts;
            rval = _mbImpl->get_entities_by_type( vertexRange[k], MBVERTEX, verts );

            MBERRORR( rval, "can't get vertices from vertex set" );
            if( verts.size() != 1 ) MBERRORR( MB_FAILURE, " node set not defined well" );
            if( firstNode == verts[0] )
            {
                secondSet = vertexRange[1 - k];  // the other set; it is 1 or 0
                break;
            }
        }
        if( k >= 2 )
        {
            // it means we didn't find the right node set
            MBERRORR( MB_FAILURE, " can't find the right vertex" );
        }
        // remove the second set from the connectivity with the
        //  edge (parent-child relation)
        // remove_parent_child( EntityHandle parent,
        //                                          EntityHandle child )
        rval = _mbImpl->remove_parent_child( edge, secondSet );
        MBERRORR( rval, " can't remove the second vertex from edge" );
    }
    // at this point, we just need to add to both edges the new node set (vertex)
    rval = _mbImpl->add_parent_child( edge, nodeSet );
    MBERRORR( rval, " can't add new vertex to old edge" );

    rval = _mbImpl->add_parent_child( new_edge, nodeSet );
    MBERRORR( rval, " can't add new vertex to new edge" );

    // one thing that I forgot: add the secondSet as a child to new edge!!!
    // (so far, the new edge has only one end vertex!)
    rval = _mbImpl->add_parent_child( new_edge, secondSet );
    MBERRORR( rval, " can't add second vertex to new edge" );

    // now, add the edge and vertex to geo tool

    rval = _my_geomTopoTool->add_geo_set( new_edge, 1 );
    MBERRORR( rval, " can't add new edge" );

    // next, get the adjacent faces to initial edge, and add them as parents
    // to the new edge

    // need to find the adjacent face sets
    Range faceRange;
    rval = getAdjacentEntities( edge, 2, faceRange );

    // these faces will be adjacent to the new edge too!
    // need to set the senses of new edge within faces

    for( Range::iterator it = faceRange.begin(); it != faceRange.end(); ++it )
    {
        EntityHandle face = *it;
        rval              = _mbImpl->add_parent_child( face, new_edge );
        MBERRORR( rval, " can't add new edge - face parent relation" );
        int sense;
        rval = _my_geomTopoTool->get_sense( edge, face, sense );
        MBERRORR( rval, " can't get initial sense of edge in the adjacent face" );
        // keep the same sense for the new edge within the faces
        rval = _my_geomTopoTool->set_sense( new_edge, face, sense );
        MBERRORR( rval, " can't set sense of new edge in the adjacent face" );
    }

    return MB_SUCCESS;
}

Definition at line 2332 of file FBEngine.cpp.

References _edges, _mbImpl, _my_geomTopoTool, _smooth, moab::debug_splits, dim, moab::GeomTopoTool::dimension(), moab::Interface::id_from_handle(), split_edge_at_mesh_node(), and moab::SmoothCurve::u_from_position().

{
    // return MB_NOT_IMPLEMENTED;
    // first, we need to find the closest point on the smooth edge, then
    // split the smooth edge, then call the split_edge_at_mesh_node
    // or maybe just find the closest node??
    // first of all, we need to find a point on the smooth edge, close by
    // then split the mesh edge (if needed)
    // this would be quite a drama, as the new edge has to be inserted in
    // order for proper geo edge definition

    // first of all, find a closest point
    // usually called for
    if( debug_splits )
    {
        std::cout << "Split edge " << _mbImpl->id_from_handle( edge ) << " at point:" << point << "\n";
    }
    int dim = _my_geomTopoTool->dimension( edge );
    if( dim != 1 ) return MB_FAILURE;
    if( !_smooth ) return MB_FAILURE;  // call it only for smooth option...
    // maybe we should do something for linear too

    SmoothCurve* curve = this->_edges[edge];
    EntityHandle closeNode;
    int edgeIndex;
    double u = curve->u_from_position( point[0], point[1], point[2], closeNode, edgeIndex );
    if( 0 == closeNode )
    {
        // we really need to split an existing edge
        // do not do that yet
        // evaluate from u:
        /*double pos[3];
        curve->position_from_u(u,  pos[0], pos[1], pos[2] );*/
        // create a new node here, and split one edge
        // change also connectivity, create new triangles on both sides ...
        std::cout << "not found a close node,  u is: " << u << " edge index: " << edgeIndex << "\n";
        return MB_FAILURE;  // not implemented yet
    }
    return split_edge_at_mesh_node( edge, closeNode, new_edge );
}

Definition at line 3089 of file FBEngine.cpp.

References _mbImpl, _newTriangles, _piercedTriangles, moab::Range::begin(), moab::Interface::create_element(), moab::debug_splits, moab::Range::end(), ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Range::insert(), moab::Interface::list_entity(), MB_SUCCESS, MBERRORR, MBTRI, print_debug_triangle(), moab::Interface::side_number(), moab::Range::size(), and moab::subtract().

Referenced by split_surface_with_direction().

{
    // split the edge, and form 4 new triangles and 2 edges
    // get 2 triangles adjacent to edge
    Range adj_tri;
    ErrorCode rval = _mbImpl->get_adjacencies( &edge, 1, 2, false, adj_tri );
    MBERRORR( rval, "can't get adj_tris" );
    adj_tri = subtract( adj_tri, _piercedTriangles );
    if( adj_tri.size() >= 3 )
    {
        std::cout << "WARNING: non manifold geometry. Are you sure?";
    }
    for( Range::iterator it = adj_tri.begin(); it != adj_tri.end(); ++it )
    {
        EntityHandle tri = *it;
        _piercedTriangles.insert( tri );
        const EntityHandle* conn3;
        int nnodes;
        rval = _mbImpl->get_connectivity( tri, conn3, nnodes );
        MBERRORR( rval, "can't get nodes" );
        int num1, sense, offset;
        rval = _mbImpl->side_number( tri, edge, num1, sense, offset );
        MBERRORR( rval, "can't get side number" );
        // after we know the side number, we can split in 2 triangles
        // node i is opposite to edge i
        int num2 = ( num1 + 1 ) % 3;
        int num3 = ( num2 + 1 ) % 3;
        // the edge from num1 to num2 is split into 2 edges
        EntityHandle t1[] = { conn3[num2], conn3[num3], newVertex };
        EntityHandle t2[] = { conn3[num1], newVertex, conn3[num3] };
        EntityHandle newTriangle, newTriangle2;
        rval = _mbImpl->create_element( MBTRI, t1, 3, newTriangle );
        MBERRORR( rval, "can't create triangle" );
        _newTriangles.insert( newTriangle );
        rval = _mbImpl->create_element( MBTRI, t2, 3, newTriangle2 );
        MBERRORR( rval, "can't create triangle" );
        _newTriangles.insert( newTriangle2 );
        // create edges with this, indirectly
        std::vector< EntityHandle > edges0;
        rval = _mbImpl->get_adjacencies( &newTriangle, 1, 1, true, edges0 );
        MBERRORR( rval, "can't get new edges" );
        edges0.clear();
        rval = _mbImpl->get_adjacencies( &newTriangle2, 1, 1, true, edges0 );
        MBERRORR( rval, "can't get new edges" );
        if( debug_splits )
        {
            std::cout << "2 (out of 4) triangles formed:\n";
            _mbImpl->list_entity( newTriangle );
            print_debug_triangle( newTriangle );
            _mbImpl->list_entity( newTriangle2 );
            print_debug_triangle( newTriangle2 );
        }
    }
    return MB_SUCCESS;
}

Definition at line 2973 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, _my_gsets, _t_created, moab::Interface::add_entities(), moab::Range::begin(), moab::Interface::create_element(), moab::GeomTopoTool::duplicate_model(), moab::Range::end(), ErrorCode, moab::GeomTopoTool::find_geomsets(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), moab::Interface::get_number_entities_by_type(), length_squared(), MB_SUCCESS, MBERRORR, MBQUAD, MBTRI, and moab::Interface::remove_entities().

Referenced by Init().

{
    // first see if we have any quads in the 2d faces
    //  _my_gsets[2] is a range of surfaces (moab sets of dimension 2)
    int num_quads = 0;
    for( Range::iterator it = _my_gsets[2].begin(); it != _my_gsets[2].end(); ++it )
    {
        EntityHandle surface = *it;
        int num_q            = 0;
        _mbImpl->get_number_entities_by_type( surface, MBQUAD, num_q );
        num_quads += num_q;
    }

    if( num_quads == 0 ) return MB_SUCCESS;  // nothing to do

    GeomTopoTool* new_gtt = NULL;
    ErrorCode rval        = _my_geomTopoTool->duplicate_model( new_gtt );
    MBERRORR( rval, "can't duplicate model" );
    if( this->_t_created ) delete _my_geomTopoTool;

    _t_created = true;  // this one is for sure created here, even if the original gtt was not

    // if we were using smart pointers, we would decrease the reference to the _my_geomTopoTool, at
    // least
    _my_geomTopoTool = new_gtt;

    // replace the _my_gsets with the new ones, from the new set
    _my_geomTopoTool->find_geomsets( _my_gsets );

    // we have duplicated now the model, and the quads are part of the new _my_gsets[2]
    // we will split them now, by creating triangles along the smallest diagonal
    // maybe we will come up with a better criteria, but for the time being, it should be fine.
    // what we will do: we will remove all quads from the surface sets, and split them

    for( Range::iterator it2 = _my_gsets[2].begin(); it2 != _my_gsets[2].end(); ++it2 )
    {
        EntityHandle surface = *it2;
        Range quads;
        rval = _mbImpl->get_entities_by_type( surface, MBQUAD, quads );
        MBERRORR( rval, "can't get quads from the surface set" );
        rval = _mbImpl->remove_entities( surface, quads );
        MBERRORR( rval, "can't remove quads from the surface set" );  // they are not deleted, just
                                                                      // removed from the set
        for( Range::iterator it = quads.begin(); it != quads.end(); ++it )
        {
            EntityHandle quad = *it;
            int nnodes;
            const EntityHandle* conn;
            rval = _mbImpl->get_connectivity( quad, conn, nnodes );
            MBERRORR( rval, "can't get quad connectivity" );
            // get all vertices position, to see which one is the shortest diagonal
            CartVect pos[4];
            rval = _mbImpl->get_coords( conn, 4, (double*)&pos[0] );
            MBERRORR( rval, "can't get node coordinates" );
            bool diag1 = ( ( pos[2] - pos[0] ).length_squared() < ( pos[3] - pos[1] ).length_squared() );
            EntityHandle newTris[2];
            EntityHandle tri1[3] = { conn[0], conn[1], conn[2] };
            EntityHandle tri2[3] = { conn[0], conn[2], conn[3] };
            if( !diag1 )
            {
                tri1[2] = conn[3];
                tri2[0] = conn[1];
            }
            rval = _mbImpl->create_element( MBTRI, tri1, 3, newTris[0] );
            MBERRORR( rval, "can't create triangle 1" );
            rval = _mbImpl->create_element( MBTRI, tri2, 3, newTris[1] );
            MBERRORR( rval, "can't create triangle 2" );
            rval = _mbImpl->add_entities( surface, newTris, 2 );
            MBERRORR( rval, "can't add triangles to the set" );
        }
        //
    }
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::split_surface ( EntityHandle  face,
std::vector< EntityHandle > &  chainedEdges,
std::vector< EntityHandle > &  splittingNodes,
EntityHandle newFace 
)

this method splits along the polyline defined by points and entities the polyline will be defined with // the entities are now only nodes and edges, no triangles!!! the first and last ones are also nodes for sure

Definition at line 1513 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, _piercedEdges, _piercedTriangles, _smooth, moab::Interface::add_entities(), moab::GeomTopoTool::add_geo_set(), moab::Interface::add_parent_child(), moab::Range::clear(), moab::Interface::create_meshset(), moab::debug_splits, moab::Interface::delete_entities(), ErrorCode, moab::GeomUtil::first(), moab::Interface::get_entities_by_type(), MB_SUCCESS, MBERRORR, MBTRI, MESHSET_SET, redistribute_boundary_edges_to_faces(), moab::Interface::remove_entities(), separate(), set_neumann_tags(), moab::GeomTopoTool::set_sense(), smooth_new_intx_points(), split_boundary(), and moab::Interface::write_file().

Referenced by split_surface_with_direction().

{
    // use now the chained edges to create a new face (loop or clean split)
    // use a fill to determine the new sets, up to the polyline
    // points on the polyline will be moved to the closest point location, with some constraints
    // then the sets will be reset, geometry recomputed. new vertices, new edges, etc.

    Range iniTris;
    ErrorCode rval;
    rval = _mbImpl->get_entities_by_type( face, MBTRI, iniTris );
    MBERRORR( rval, "can't get initial triangles" );

    // start from a triangle that is not in the triangles to delete
    // flood fill

    bool closed = splittingNodes.size() == 0;
    if( !closed )
    {
        //
        if( splittingNodes.size() != 2 ) MBERRORR( MB_FAILURE, "need to have exactly 2 nodes for splitting" );
        // we will have to split the boundary edges
        // first, find the actual boundary, and try to split with the 2 end points (nodes)
        // get the adjacent edges, and see which one has the end nodes
        rval = split_boundary( face, splittingNodes[0] );
        MBERRORR( rval, "can't split with first node" );
        rval = split_boundary( face, splittingNodes[1] );
        MBERRORR( rval, "can't split with second node)" );
    }
    // we will separate triangles to delete, unaffected, new_triangles,
    //  nodesAlongPolyline,
    Range first, second;
    rval = separate( face, chainedEdges, first, second );

    // now, we are done with the computations;
    // we need to put the new nodes on the smooth surface
    if( this->_smooth )
    {
        rval = smooth_new_intx_points( face, chainedEdges );
        MBERRORR( rval, "can't smooth new points" );
    }

    // create the new set
    rval = _mbImpl->create_meshset( MESHSET_SET, newFace );
    MBERRORR( rval, "can't create a new face" );

    _my_geomTopoTool->add_geo_set( newFace, 2 );

    // the new face will have the first set (positive sense triangles, to the left)
    rval = _mbImpl->add_entities( newFace, first );
    MBERRORR( rval, "can't add first range triangles to new face" );

    for( unsigned int j = 0; j < chainedEdges.size(); j++ )
    {
        EntityHandle new_geo_edge = chainedEdges[j];
        // both faces will have the edge now
        rval = _mbImpl->add_parent_child( face, new_geo_edge );
        MBERRORR( rval, "can't add parent child relations for new edge" );

        rval = _mbImpl->add_parent_child( newFace, new_geo_edge );
        MBERRORR( rval, "can't add parent child relations for new edge" );
        // add senses
        // sense newFace is 1, old face is -1
        rval = _my_geomTopoTool->set_sense( new_geo_edge, newFace, 1 );
        MBERRORR( rval, "can't set sense for new edge" );

        rval = _my_geomTopoTool->set_sense( new_geo_edge, face, -1 );
        MBERRORR( rval, "can't set sense for new edge in original face" );
    }

    rval = set_neumann_tags( face, newFace );
    MBERRORR( rval, "can't set NEUMANN set tags" );

    // now, we should remove from the original set all tris, and put the "second" range
    rval = _mbImpl->remove_entities( face, iniTris );
    MBERRORR( rval, "can't remove original tris from initial face set" );

    rval = _mbImpl->add_entities( face, second );
    MBERRORR( rval, "can't add second range to the original set" );

    if( !closed )
    {
        rval = redistribute_boundary_edges_to_faces( face, newFace, chainedEdges );
        MBERRORR( rval, "fail to reset the proper boundary faces" );
    }

    /*if (_smooth)
      delete_smooth_tags();// they need to be recomputed, anyway
    // this will remove the extra smooth faces and edges
    clean();*/
    // also, these nodes need to be moved to the smooth surface, sometimes before deleting the old
    // triangles
    // remove the triangles from the set, then delete triangles (also some edges need to be
    // deleted!)
    rval = _mbImpl->delete_entities( _piercedTriangles );

    MBERRORR( rval, "can't delete triangles" );
    _piercedTriangles.clear();
    // delete edges that are broke up in 2
    rval = _mbImpl->delete_entities( _piercedEdges );
    MBERRORR( rval, "can't delete edges" );
    _piercedEdges.clear();

    if( debug_splits )
    {
        _mbImpl->write_file( "newFace.vtk", "vtk", 0, &newFace, 1 );
        _mbImpl->write_file( "leftoverFace.vtk", "vtk", 0, &face, 1 );
    }
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::split_surface_with_direction ( EntityHandle  face,
std::vector< double > &  xyz,
double *  direction,
int  closed,
double  min_dot,
EntityHandle oNewFace 
)

the segment is the first or last segment in the polyline

Definition at line 1212 of file FBEngine.cpp.

References _brokenEdges, _mbImpl, _my_geomTopoTool, _piercedEdges, _piercedTriangles, moab::area_coordinates(), moab::Range::begin(), boundary_mesh_edges_on_face(), boundary_nodes_on_face(), BreakTriangle2(), chain_edges(), compute_intersection_points(), create_new_gedge(), moab::Interface::create_vertex(), moab::debug_splits, divide_triangle(), moab::Range::end(), entities, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::GeomTopoTool::get_root(), moab::Interface::id_from_handle(), moab::Range::insert(), moab::intersect(), moab::intersect_segment_and_plane_slice(), MB_SUCCESS, MBEDGE, MBERRORR, MBVERTEX, moab::CartVect::normalize(), moab::GeomTopoTool::obb_tree(), moab::OrientedBoxTreeTool::ray_intersect_sets(), moab::Range::size(), split_internal_edge(), split_surface(), moab::tolerance, and moab::Interface::type_from_handle().

{

    // first of all, find all intersection points (piercing in the face along the direction)
    // assume it is robust; what if it is not sufficiently robust?
    // if the polyline is open, find the intersection with the boundary edges, of the
    // polyline extruded at ends

    ErrorCode rval;

    // then find the position
    int numIniPoints = (int)xyz.size() / 3;
    if( ( closed && numIniPoints < 3 ) || ( !closed && numIniPoints < 2 ) )
        MBERRORR( MB_FAILURE, "not enough polyline points " );
    EntityHandle rootForFace;

    rval = _my_geomTopoTool->get_root( face, rootForFace );
    MBERRORR( rval, "Failed to get root of face." );

    const double dir[] = { direction[0], direction[1], direction[2] };
    std::vector< EntityHandle > nodes;  // get the nodes closest to the ray traces of interest

    // these are nodes on the boundary of original face;
    // if the cutting line is not closed, the starting - ending vertices of the
    // polygonal line must come from this list

    std::vector< CartVect > b_pos;
    std::vector< EntityHandle > boundary_nodes;
    std::vector< EntityHandle > splittingNodes;
    Range boundary_mesh_edges;
    if( !closed )
    {
        rval = boundary_nodes_on_face( face, boundary_nodes );
        MBERRORR( rval, "Failed to get boundary nodes." );
        b_pos.resize( boundary_nodes.size() );
        rval = _mbImpl->get_coords( &( boundary_nodes[0] ), boundary_nodes.size(), (double*)( &b_pos[0][0] ) );
        MBERRORR( rval, "Failed to get coordinates for boundary nodes." );
        rval = boundary_mesh_edges_on_face( face, boundary_mesh_edges );
        MBERRORR( rval, "Failed to get mesh boundary edges for face." );
    }
    //
    int i = 0;
    CartVect dirct( direction );
    dirct.normalize();  // maybe an overkill?
    for( ; i < numIniPoints; i++ )
    {

        const double point[] = { xyz[3 * i], xyz[3 * i + 1], xyz[3 * i + 2] };  // or even point( &(xyz[3*i]) ); //
        CartVect p1( point );
        if( !closed && ( ( 0 == i ) || ( numIniPoints - 1 == i ) ) )
        {

            // find the intersection point between a plane and boundary mesh edges
            // this will be the closest point on the boundary of face
            /// the segment is the first or last segment in the polyline
            int i1 = i + 1;
            if( i == numIniPoints - 1 ) i1 = i - 1;  // previous point if the last
            // the direction is from point to point1
            const double point1[] = { xyz[3 * i1], xyz[3 * i1 + 1], xyz[3 * i1 + 2] };
            CartVect p2( point1 );
            CartVect normPlane = ( p2 - p1 ) * dirct;
            normPlane.normalize();
            //(roughly, from p1 to p2, perpendicular to dirct, in the "xy" plane
            // if the intx point is "outside" p1 - p2, skip if the intx point is closer to p2
            CartVect perpDir    = dirct * normPlane;
            Range::iterator ite = boundary_mesh_edges.begin();
            // do a linear search for the best intersection point position (on a boundary edge)
            if( debug_splits )
            {
                std::cout << " p1:" << p1 << "\n";
                std::cout << " p2:" << p2 << "\n";
                std::cout << " perpDir:" << perpDir << "\n";
                std::cout << " boundary edges size:" << boundary_mesh_edges.size() << "\n";
            }
            for( ; ite != boundary_mesh_edges.end(); ++ite )
            {
                EntityHandle candidateEdge = *ite;
                const EntityHandle* conn2;
                int nno;
                rval = _mbImpl->get_connectivity( candidateEdge, conn2, nno );
                MBERRORR( rval, "Failed to get conn for boundary edge" );
                CartVect pts[2];
                rval = _mbImpl->get_coords( conn2, 2, &( pts[0][0] ) );
                MBERRORR( rval, "Failed to get coords of nodes for boundary edge" );
                CartVect intx_point;
                double parPos;
                bool intersect =
                    intersect_segment_and_plane_slice( pts[0], pts[1], p1, p2, dirct, normPlane, intx_point, parPos );
                if( debug_splits )
                {
                    std::cout << "   Edge:" << _mbImpl->id_from_handle( candidateEdge ) << "\n";
                    std::cout << "   Node 1:" << _mbImpl->id_from_handle( conn2[0] ) << pts[0] << "\n";
                    std::cout << "   Node 2:" << _mbImpl->id_from_handle( conn2[1] ) << pts[1] << "\n";
                    std::cout << "    Intersect bool:" << intersect << "\n";
                }
                if( intersect )
                {
                    double proj1 = ( intx_point - p1 ) % perpDir;
                    double proj2 = ( intx_point - p2 ) % perpDir;
                    if( ( fabs( proj1 ) > fabs( proj2 ) )  // this means it is closer to p2 than p1
                    )
                        continue;  // basically, this means the intersection point is with a
                                   //  boundary edge on the other side, closer to p2 than p1, so we
                                   //  skip it
                    if( parPos == 0 )
                    {
                        // close to vertex 1, nothing to do
                        nodes.push_back( conn2[0] );
                        splittingNodes.push_back( conn2[0] );
                    }
                    else if( parPos == 1. )
                    {
                        // close to vertex 2, nothing to do
                        nodes.push_back( conn2[1] );
                        splittingNodes.push_back( conn2[1] );
                    }
                    else
                    {
                        // break the edge, create a new node at intersection point (will be smoothed
                        // out)
                        EntityHandle newVertex;
                        rval = _mbImpl->create_vertex( &( intx_point[0] ), newVertex );
                        MBERRORR( rval, "can't create vertex" );
                        nodes.push_back( newVertex );
                        split_internal_edge( candidateEdge, newVertex );
                        splittingNodes.push_back( newVertex );
                        _brokenEdges[newVertex] = candidateEdge;
                        _piercedEdges.insert( candidateEdge );
                    }
                    break;  // break from the loop over boundary edges, we are interested in the
                            // first
                            //      split (hopefully, the only split)
                }
            }
            if( ite == boundary_mesh_edges.end() )
                MBERRORR( MB_FAILURE, "Failed to find boundary intersection edge. Bail out" );
        }
        else
        {
            std::vector< double > distances_out;
            std::vector< EntityHandle > sets_out;
            std::vector< EntityHandle > facets_out;
            rval = _my_geomTopoTool->obb_tree()->ray_intersect_sets( distances_out, sets_out, facets_out, rootForFace,
                                                                     tolerance, point, dir );
            MBERRORR( rval, "Failed to get ray intersections." );
            if( distances_out.size() < 1 )
                MBERRORR( MB_FAILURE, "Failed to get one intersection point, bad direction." );

            if( distances_out.size() > 1 )
            {
                std::cout << " too many intersection points. Only the first one considered\n";
            }
            std::vector< EntityHandle >::iterator pFace = std::find( sets_out.begin(), sets_out.end(), face );

            if( pFace == sets_out.end() ) MBERRORR( MB_FAILURE, "Failed to intersect given face, bad direction." );
            unsigned int index = pFace - sets_out.begin();
            // get the closest node of the triangle, and modify locally the triangle(s), so the
            // intersection point is at a new vertex, if needed
            CartVect P( point );
            CartVect Dir( dir );
            CartVect newPoint = P + distances_out[index] * Dir;
            // get the triangle coordinates

            double area_coord[3];
            EntityHandle boundary_handle = 0;  // if 0, not on a boundary
            rval = area_coordinates( _mbImpl, facets_out[index], newPoint, area_coord, boundary_handle );
            MBERRORR( rval, "Failed to get area coordinates" );

            if( debug_splits )
            {
                std::cout << " int point:" << newPoint << " area coord " << area_coord[0] << " " << area_coord[1] << " "
                          << area_coord[2] << "\n";
                std::cout << " triangle: " << _mbImpl->id_from_handle( facets_out[index] )
                          << " boundary:" << boundary_handle << "\n";
            }
            EntityType type;
            if( boundary_handle ) type = _mbImpl->type_from_handle( boundary_handle );
            if( boundary_handle && ( type == MBVERTEX ) )
            {
                // nothing to do, we are happy
                nodes.push_back( boundary_handle );
            }
            else
            {
                // for an edge, we will split 2 triangles
                // for interior point, we will create 3 triangles out of one
                // create a new vertex
                EntityHandle newVertex;
                rval = _mbImpl->create_vertex( &( newPoint[0] ), newVertex );
                if( boundary_handle )  // this is edge
                {
                    split_internal_edge( boundary_handle, newVertex );
                    _piercedEdges.insert( boundary_handle );  // to be removed at the end
                }
                else
                    divide_triangle( facets_out[index], newVertex );

                nodes.push_back( newVertex );
            }
        }
    }
    // now, we have to find more intersection points, either interior to triangles, or on edges, or
    // on vertices use the same tolerance as before starting from 2 points on 2 triangles, and
    // having the direction, get more intersection points between the plane formed by direction and
    // those 2 points, and edges from triangulation (the triangles involved will be part of the same
    // gentity , original face ( moab set)
    // int closed = 1;// closed = 0 if the polyline is not closed

    CartVect Dir( direction );
    std::vector< EntityHandle > chainedEdges;

    for( i = 0; i < numIniPoints - 1 + closed; i++ )
    {
        int nextIndex = ( i + 1 ) % numIniPoints;
        std::vector< EntityHandle > trianglesAlong;
        std::vector< CartVect > points;
        // otherwise to edges or even nodes
        std::vector< EntityHandle > entities;
        // start with initial points, intersect along the direction, find the facets
        rval = compute_intersection_points( face, nodes[i], nodes[nextIndex], Dir, points, entities, trianglesAlong );
        MBERRORR( rval, "can't get intersection points along a line" );
        std::vector< EntityHandle > nodesAlongPolyline;
        // refactor code; move some edge creation for each 2 intersection points
        nodesAlongPolyline.push_back( entities[0] );  // it is for sure a node
        int num_points = (int)points.size();          // it should be num_triangles + 1
        for( int j = 0; j < num_points - 1; j++ )
        {
            EntityHandle tri = trianglesAlong[j];  // this is happening in trianglesAlong i
            EntityHandle e1  = entities[j];
            EntityHandle e2  = entities[j + 1];
            EntityType et1   = _mbImpl->type_from_handle( e1 );
            // EntityHandle vertex1 = nodesAlongPolyline[i];// irrespective of the entity type i,
            // we already have the vertex there
            EntityType et2 = _mbImpl->type_from_handle( e2 );
            if( et2 == MBVERTEX )
            {
                nodesAlongPolyline.push_back( e2 );
            }
            else  // if (et2==MBEDGE)
            {
                CartVect coord_vert = points[j + 1];
                EntityHandle newVertex;
                rval = _mbImpl->create_vertex( (double*)&coord_vert, newVertex );
                MBERRORR( rval, "can't create vertex" );
                nodesAlongPolyline.push_back( newVertex );
            }
            // if vertices, do not split anything, just get the edge for polyline
            if( et2 == MBVERTEX && et1 == MBVERTEX )
            {
                // nothing to do, just continue;
                continue;  // continue the for loop
            }

            if( debug_splits )
            {
                std::cout << "tri: type: " << _mbImpl->type_from_handle( tri )
                          << " id:" << _mbImpl->id_from_handle( tri ) << "\n    e1:" << e1
                          << " id:" << _mbImpl->id_from_handle( e1 ) << "   e2:" << e2
                          << " id:" << _mbImpl->id_from_handle( e2 ) << "\n";
            }
            // here, at least one is an edge
            rval = BreakTriangle2( tri, e1, e2, nodesAlongPolyline[j], nodesAlongPolyline[j + 1] );
            MBERRORR( rval, "can't break triangle 2" );
            if( et2 == MBEDGE ) _piercedEdges.insert( e2 );
            _piercedTriangles.insert( tri );
        }
        // nodesAlongPolyline will define a new geometric edge
        if( debug_splits )
        {
            std::cout << "nodesAlongPolyline: " << nodesAlongPolyline.size() << "\n";
            std::cout << "points: " << num_points << "\n";
        }
        // if needed, create edges along polyline, or revert the existing ones, to
        // put them in a new edge set
        EntityHandle new_geo_edge;
        rval = create_new_gedge( nodesAlongPolyline, new_geo_edge );
        MBERRORR( rval, "can't create a new edge" );
        chainedEdges.push_back( new_geo_edge );
        // end copy
    }
    // the segment between point_i and point_i+1 is in trianglesAlong_i
    // points_i is on entities_i
    // all these edges are oriented correctly
    rval = split_surface( face, chainedEdges, splittingNodes, oNewFace );
    MBERRORR( rval, "can't split surface" );
    //
    rval = chain_edges( min_dot );  // acos(0.8)~= 36 degrees
    MBERRORR( rval, "can't chain edges" );
    return MB_SUCCESS;
}
ErrorCode moab::FBEngine::weave_lateral_face_from_edges ( EntityHandle  bEdge,
EntityHandle  tEdge,
double *  direction,
EntityHandle newLatFace 
)

Definition at line 3283 of file FBEngine.cpp.

References _mbImpl, _my_geomTopoTool, moab::Interface::add_entities(), moab::GeomTopoTool::add_geo_set(), moab::Interface::add_parent_child(), moab::Interface::contains_entities(), moab::Interface::create_element(), moab::Interface::create_meshset(), moab::Range::end(), ErrorCode, moab::Range::find(), moab::Interface::get_child_meshsets(), moab::Interface::get_coords(), get_nodes_from_edge(), moab::Interface::get_parent_meshsets(), getEgVtxSense(), moab::CartVect::length_squared(), MB_SUCCESS, MBEDGE, MBERRORR, MBTRI, MESHSET_SET, moab::CartVect::normalize(), and moab::GeomTopoTool::set_sense().

Referenced by create_volume_with_direction().

{
    // in weird cases might need to create new vertices in the interior;
    // in most cases, it is OK

    ErrorCode rval = _mbImpl->create_meshset( MESHSET_SET, newLatFace );
    MBERRORR( rval, "can't create new lateral face" );

    EntityHandle v[4];  // vertex sets
    // bot edge will be v1->v2
    // top edge will be v3->v4
    // we need to create edges from v1 to v3 and from v2 to v4
    std::vector< EntityHandle > adj;
    rval = _mbImpl->get_child_meshsets( bEdge, adj );
    MBERRORR( rval, "can't get children nodes" );
    bool periodic = false;
    if( adj.size() == 1 )
    {
        v[0] = v[1] = adj[0];
        periodic    = true;
    }
    else
    {
        v[0] = adj[0];
        v[1] = adj[1];
    }
    int senseB;
    rval = getEgVtxSense( bEdge, v[0], v[1], senseB );
    MBERRORR( rval, "can't get bottom edge sense" );
    if( -1 == senseB )
    {
        v[1] = adj[0];  // so , bEdge will be oriented from v[0] to v[1], and will start at
                        // nodes1, coords1..
        v[0] = adj[1];
    }
    adj.clear();
    rval = _mbImpl->get_child_meshsets( tEdge, adj );
    MBERRORR( rval, "can't get children nodes" );
    if( adj.size() == 1 )
    {
        v[2] = v[3] = adj[0];
        if( !periodic ) MBERRORR( MB_FAILURE, "top edge is periodic, but bottom edge is not" );
    }
    else
    {
        v[2] = adj[0];
        v[3] = adj[1];
        if( periodic ) MBERRORR( MB_FAILURE, "top edge is not periodic, but bottom edge is" );
    }

    // now, using nodes on bottom edge and top edge, create triangles, oriented outwards the
    //  volume (sense positive on bottom edge)
    std::vector< EntityHandle > nodes1;
    rval = get_nodes_from_edge( bEdge, nodes1 );
    MBERRORR( rval, "can't get nodes from bott edge" );

    std::vector< EntityHandle > nodes2;
    rval = get_nodes_from_edge( tEdge, nodes2 );
    MBERRORR( rval, "can't get nodes from top edge" );

    std::vector< CartVect > coords1, coords2;
    coords1.resize( nodes1.size() );
    coords2.resize( nodes2.size() );

    int N1 = (int)nodes1.size();
    int N2 = (int)nodes2.size();

    rval = _mbImpl->get_coords( &( nodes1[0] ), nodes1.size(), (double*)&( coords1[0] ) );
    MBERRORR( rval, "can't get coords of nodes from bott edge" );

    rval = _mbImpl->get_coords( &( nodes2[0] ), nodes2.size(), (double*)&( coords2[0] ) );
    MBERRORR( rval, "can't get coords of nodes from top edge" );
    CartVect up( direction );

    // see if the start and end coordinates are matching, if not, reverse edge 2 nodes and
    // coordinates
    CartVect v1 = ( coords1[0] - coords2[0] ) * up;
    CartVect v2 = ( coords1[0] - coords2[N2 - 1] ) * up;
    if( v1.length_squared() > v2.length_squared() )
    {
        // we need to reverse coords2 and node2, as nodes are not above each other
        // the edges need to be found between v0 and v3, v1 and v2!
        for( unsigned int k = 0; k < nodes2.size() / 2; k++ )
        {
            EntityHandle tmp    = nodes2[k];
            nodes2[k]           = nodes2[N2 - 1 - k];
            nodes2[N2 - 1 - k]  = tmp;
            CartVect tv         = coords2[k];
            coords2[k]          = coords2[N2 - 1 - k];
            coords2[N2 - 1 - k] = tv;
        }
    }
    // make sure v[2] has nodes2[0], if not, reverse v[2] and v[3]
    if( !_mbImpl->contains_entities( v[2], &( nodes2[0] ), 1 ) )
    {
        // reverse v[2] and v[3], so v[2] will be above v[0]
        EntityHandle tmp = v[2];
        v[2]             = v[3];
        v[3]             = tmp;
    }
    // now we know that v[0]--v[3] will be vertex sets in the order we want
    EntityHandle nd[4] = { nodes1[0], nodes1[N1 - 1], nodes2[0], nodes2[N2 - 1] };

    adj.clear();
    EntityHandle e1, e2;
    // find edge 1 between v[0] and v[2], and e2 between v[1] and v[3]
    rval = _mbImpl->get_parent_meshsets( v[0], adj );
    MBERRORR( rval, "can't get edges connected to vertex set 1" );
    bool found = false;
    for( unsigned int j = 0; j < adj.size(); j++ )
    {
        EntityHandle ed = adj[j];
        Range vertices;
        rval = _mbImpl->get_child_meshsets( ed, vertices );
        if( vertices.find( v[2] ) != vertices.end() )
        {
            found = true;
            e1    = ed;
            break;
        }
    }
    if( !found )
    {
        // create an edge from v[0] to v[2]
        rval = _mbImpl->create_meshset( MESHSET_SET, e1 );
        MBERRORR( rval, "can't create edge 1" );

        rval = _mbImpl->add_parent_child( e1, v[0] );
        MBERRORR( rval, "can't add parent - child relation" );

        rval = _mbImpl->add_parent_child( e1, v[2] );
        MBERRORR( rval, "can't add parent - child relation" );

        EntityHandle nn2[2] = { nd[0], nd[2] };
        EntityHandle medge;
        rval = _mbImpl->create_element( MBEDGE, nn2, 2, medge );
        MBERRORR( rval, "can't create mesh edge" );

        rval = _mbImpl->add_entities( e1, &medge, 1 );
        MBERRORR( rval, "can't add mesh edge to geo edge" );

        rval = this->_my_geomTopoTool->add_geo_set( e1, 1 );
        MBERRORR( rval, "can't add edge to gtt" );
    }

    // find the edge from v2 to v4 (if not, create one)
    rval = _mbImpl->get_parent_meshsets( v[1], adj );
    MBERRORR( rval, "can't get edges connected to vertex set 2" );
    found = false;
    for( unsigned int i = 0; i < adj.size(); i++ )
    {
        EntityHandle ed = adj[i];
        Range vertices;
        rval = _mbImpl->get_child_meshsets( ed, vertices );
        if( vertices.find( v[3] ) != vertices.end() )
        {
            found = true;
            e2    = ed;
            break;
        }
    }
    if( !found )
    {
        // create an edge from v2 to v4
        rval = _mbImpl->create_meshset( MESHSET_SET, e2 );
        MBERRORR( rval, "can't create edge 1" );

        rval = _mbImpl->add_parent_child( e2, v[1] );
        MBERRORR( rval, "can't add parent - child relation" );

        rval = _mbImpl->add_parent_child( e2, v[3] );
        MBERRORR( rval, "can't add parent - child relation" );

        EntityHandle nn2[2] = { nd[1], nd[3] };
        EntityHandle medge;
        rval = _mbImpl->create_element( MBEDGE, nn2, 2, medge );
        MBERRORR( rval, "can't create mesh edge" );

        rval = _mbImpl->add_entities( e2, &medge, 1 );
        MBERRORR( rval, "can't add mesh edge to geo edge" );

        rval = _my_geomTopoTool->add_geo_set( e2, 1 );
        MBERRORR( rval, "can't add edge to gtt" );
    }

    // now we have the four edges, add them to the face, as children

    // add children to face
    rval = _mbImpl->add_parent_child( newLatFace, bEdge );
    MBERRORR( rval, "can't add parent - child relation" );

    rval = _mbImpl->add_parent_child( newLatFace, tEdge );
    MBERRORR( rval, "can't add parent - child relation" );

    rval = _mbImpl->add_parent_child( newLatFace, e1 );
    MBERRORR( rval, "can't add parent - child relation" );

    rval = _mbImpl->add_parent_child( newLatFace, e2 );
    MBERRORR( rval, "can't add parent - child relation" );

    rval = _my_geomTopoTool->add_geo_set( newLatFace, 2 );
    MBERRORR( rval, "can't add face to gtt" );
    // add senses
    //
    rval = _my_geomTopoTool->set_sense( bEdge, newLatFace, 1 );
    MBERRORR( rval, "can't set bottom edge sense to the lateral face" );

    int Tsense;
    rval = getEgVtxSense( tEdge, v[3], v[2], Tsense );
    MBERRORR( rval, "can't get top edge sense" );
    // we need to see what sense has topEdge in face
    rval = _my_geomTopoTool->set_sense( tEdge, newLatFace, Tsense );
    MBERRORR( rval, "can't set top edge sense to the lateral face" );

    rval = _my_geomTopoTool->set_sense( e1, newLatFace, -1 );
    MBERRORR( rval, "can't set first vert edge sense" );

    rval = _my_geomTopoTool->set_sense( e2, newLatFace, 1 );
    MBERRORR( rval, "can't set second edge sense to the lateral face" );
    // first, create edges along direction, for the

    int indexB = 0, indexT = 0;  // indices of the current nodes in the weaving process
    // weaving is either up or down; the triangles are oriented positively either way
    // up is nodes1[indexB], nodes2[indexT+1], nodes2[indexT]
    // down is nodes1[indexB], nodes1[indexB+1], nodes2[indexT]
    // the decision to weave up or down is based on which one is closer to the direction normal
    /*
     *
     *     --------*------*-----------*                           ^
     *            /   .    \ .      .           ------> dir1      |  up
     *           /.         \   .  .                              |
     *     -----*------------*----*
     *
     */
    // we have to change the logic to account for cases when the curve in xy plane is not straight
    // (for example, when merging curves with a min_dot < -0.5, which means that the edges
    // could be even closed (periodic), with one vertex
    // the top and bottom curves should follow the same path in the "xy" plane (better to say
    // the plane perpendicular to the direction of weaving)
    // in this logic, the vector dir1 varies along the curve !!!
    CartVect dir1 = coords1[1] - coords1[0];  // we should have at least 2 nodes, N1>=2!!

    CartVect planeNormal = dir1 * up;
    dir1                 = up * planeNormal;
    dir1.normalize();
    // this direction will be updated constantly with the direction of last edge added
    bool weaveDown = true;

    CartVect startP = coords1[0];  // used for origin of comparisons
    while( 1 )
    {
        if( ( indexB == N1 - 1 ) && ( indexT == N2 - 1 ) ) break;  // we cannot advance anymore
        if( indexB == N1 - 1 )
        {
            weaveDown = false;
        }
        else if( indexT == N2 - 1 )
        {
            weaveDown = true;
        }
        else
        {
            // none are at the end, we have to decide which way to go, based on which  index + 1 is
            // closer
            double proj1 = ( coords1[indexB + 1] - startP ) % dir1;
            double proj2 = ( coords2[indexT + 1] - startP ) % dir1;
            if( proj1 < proj2 )
                weaveDown = true;
            else
                weaveDown = false;
        }
        EntityHandle nTri[3] = { nodes1[indexB], 0, nodes2[indexT] };
        if( weaveDown )
        {
            nTri[1] = nodes1[indexB + 1];
            nTri[2] = nodes2[indexT];
            indexB++;
        }
        else
        {
            nTri[1] = nodes2[indexT + 1];
            indexT++;
        }
        EntityHandle triangle;
        rval = _mbImpl->create_element( MBTRI, nTri, 3, triangle );
        MBERRORR( rval, "can't create triangle" );

        rval = _mbImpl->add_entities( newLatFace, &triangle, 1 );
        MBERRORR( rval, "can't add triangle to face set" );
        if( weaveDown )
        {
            // increase was from nodes1[indexB-1] to nodes1[indexb]
            dir1 = coords1[indexB] - coords1[indexB - 1];  // we should have at least 2 nodes, N1>=2!!
        }
        else
        {
            dir1 = coords2[indexT] - coords2[indexT - 1];
        }
        dir1 = up * ( dir1 * up );
        dir1.normalize();
    }
    // we do not check yet if the triangles are inverted
    // if yes, we should correct them. HOW?
    // we probably need a better meshing strategy, one that can overcome really bad meshes.
    // again, these faces are not what we should use for geometry, maybe we should use the
    //  extruded quads, identified AFTER hexa are created.
    // right now, I see only a cosmetic use of these lateral faces
    // the whole idea of volume maybe is overrated
    // volume should be just quads extruded from bottom ?
    //
    return MB_SUCCESS;
}

Member Data Documentation

Definition at line 315 of file FBEngine.hpp.

Referenced by split_boundary(), and split_surface_with_direction().

Definition at line 298 of file FBEngine.hpp.

Referenced by Init().

Definition at line 310 of file FBEngine.hpp.

Referenced by clean(), and initializeSmoothing().

Definition at line 309 of file FBEngine.hpp.

Referenced by clean(), delete_smooth_tags(), and initializeSmoothing().

Definition at line 296 of file FBEngine.hpp.

Referenced by clean(), FBEngine(), and split_quads().

List of all members.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines