Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
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

◆ FBEngine()

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

Definition at line 192 of file FBEngine.cpp.

193  : _mbImpl( impl ), _my_geomTopoTool( topoTool ), _t_created( false ), _smooth( smooth ), _initialized( false ),
194  _smthFace( NULL ), _smthCurve( NULL )
195 {
196  if( !_my_geomTopoTool )
197  {
198  _my_geomTopoTool = new GeomTopoTool( _mbImpl );
199  _t_created = true;
200  }
201  // should this be part of the constructor or not?
202  // Init();
203 }

References _mbImpl, _my_geomTopoTool, and _t_created.

◆ ~FBEngine()

moab::FBEngine::~FBEngine ( )

Definition at line 204 of file FBEngine.cpp.

205 {
206  clean();
207  _smooth = false;
208 }

References _smooth, and clean().

Member Function Documentation

◆ addEntArrToSet()

ErrorCode moab::FBEngine::addEntArrToSet ( const Range entities,
EntityHandle  set 
)

Definition at line 488 of file FBEngine.cpp.

489 {
490  return MBI->add_entities( set, entities );
491 }

References MBI.

◆ addEntSet()

ErrorCode moab::FBEngine::addEntSet ( EntityHandle  entity_set_to_add,
EntityHandle  entity_set_handle 
)

Definition at line 493 of file FBEngine.cpp.

494 {
495  return MBI->add_entities( entity_set_handle, &entity_set_to_add, 1 );
496 }

References MBI.

◆ boundary_mesh_edges_on_face()

ErrorCode moab::FBEngine::boundary_mesh_edges_on_face ( EntityHandle  face,
Range boundary_mesh_edges 
)
private

Definition at line 3043 of file FBEngine.cpp.

3044 {
3045  // this list is used only for finding the intersecting mesh edge for starting the
3046  // polygonal cutting line at boundary (if !closed)
3047  Range bound_edges;
3048  ErrorCode rval = getAdjacentEntities( face, 1, bound_edges );
3049  MBERRORR( rval, " can't get boundary edges" );
3050  for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
3051  {
3052  EntityHandle b_edge = *it;
3053  // get all edges in range
3054  // Range mesh_edges;
3055  rval = _mbImpl->get_entities_by_dimension( b_edge, 1, boundary_mesh_edges );
3056  MBERRORR( rval, " can't get mesh edges" );
3057  }
3058  return MB_SUCCESS;
3059 }

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().

◆ boundary_nodes_on_face()

ErrorCode moab::FBEngine::boundary_nodes_on_face ( EntityHandle  face,
std::vector< EntityHandle > &  boundary_nodes 
)
private

Definition at line 3060 of file FBEngine.cpp.

3061 {
3062  // even if we repeat some nodes, it is OK
3063  // this list is used only for finding the closest boundary node for starting the
3064  // polygonal cutting line at boundary (if !closed)
3065  Range bound_edges;
3066  ErrorCode rval = getAdjacentEntities( face, 1, bound_edges );
3067  MBERRORR( rval, " can't get boundary edges" );
3068  Range b_nodes;
3069  for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
3070  {
3071  EntityHandle b_edge = *it;
3072  // get all edges in range
3073  Range mesh_edges;
3074  rval = _mbImpl->get_entities_by_dimension( b_edge, 1, mesh_edges );
3075  MBERRORR( rval, " can't get mesh edges" );
3076  rval = _mbImpl->get_connectivity( mesh_edges, b_nodes );
3077  MBERRORR( rval, " can't get nodes from mesh edges" );
3078  }
3079  // create now a vector based on Range of nodes
3080  boundary_nodes.resize( b_nodes.size() );
3081  std::copy( b_nodes.begin(), b_nodes.end(), boundary_nodes.begin() );
3082  return MB_SUCCESS;
3083 }

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().

◆ BreakTriangle()

ErrorCode moab::FBEngine::BreakTriangle ( EntityHandle  tri,
EntityHandle  e1,
EntityHandle  e3,
EntityHandle  n1,
EntityHandle  n2,
EntityHandle  n3 
)
private

Definition at line 1985 of file FBEngine.cpp.

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

◆ BreakTriangle2()

ErrorCode moab::FBEngine::BreakTriangle2 ( EntityHandle  tri,
EntityHandle  e1,
EntityHandle  e2,
EntityHandle  n1,
EntityHandle  n2 
)
private

Definition at line 1991 of file FBEngine.cpp.

1996 {
1997  // we have the nodes, we just need to reconnect to form new triangles
1998  ErrorCode rval;
1999  const EntityHandle* conn3 = NULL;
2000  int nnodes = 0;
2001  rval = _mbImpl->get_connectivity( tri, conn3, nnodes );
2002  MBERRORR( rval, "Failed to get connectivity" );
2003  assert( 3 == nnodes );
2004 
2005  EntityType et1 = _mbImpl->type_from_handle( e1 );
2006  EntityType et2 = _mbImpl->type_from_handle( e2 );
2007 
2008  if( MBVERTEX == et1 )
2009  {
2010  // find the vertex in conn3, and form 2 other triangles
2011  int index = -1;
2012  for( index = 0; index < 3; index++ )
2013  {
2014  if( conn3[index] == e1 ) // also n1
2015  break;
2016  }
2017  if( index == 3 ) return MB_FAILURE;
2018  // 2 triangles: n1, index+1, n2, and n1, n2, index+2
2019  EntityHandle conn[6] = { n1, conn3[( index + 1 ) % 3], n2, n1, n2, conn3[( index + 2 ) % 3] };
2020  EntityHandle newTriangle;
2021  rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle );
2022  MBERRORR( rval, "Failed to create a new triangle" );
2023  _newTriangles.insert( newTriangle );
2024  if( debug_splits ) print_debug_triangle( newTriangle );
2025  rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle ); // the second triangle
2026  MBERRORR( rval, "Failed to create a new triangle" );
2027  _newTriangles.insert( newTriangle );
2028  if( debug_splits ) print_debug_triangle( newTriangle );
2029  }
2030  else if( MBVERTEX == et2 )
2031  {
2032  int index = -1;
2033  for( index = 0; index < 3; index++ )
2034  {
2035  if( conn3[index] == e2 ) // also n2
2036  break;
2037  }
2038  if( index == 3 ) return MB_FAILURE;
2039  // 2 triangles: n1, index+1, n2, and n1, n2, index+2
2040  EntityHandle conn[6] = { n2, conn3[( index + 1 ) % 3], n1, n2, n1, conn3[( index + 2 ) % 3] };
2041  EntityHandle newTriangle;
2042  rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle );
2043  MBERRORR( rval, "Failed to create a new triangle" );
2044  _newTriangles.insert( newTriangle );
2045  if( debug_splits ) print_debug_triangle( newTriangle );
2046  rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle ); // the second triangle
2047  MBERRORR( rval, "Failed to create a new triangle" );
2048  _newTriangles.insert( newTriangle );
2049  if( debug_splits ) print_debug_triangle( newTriangle );
2050  }
2051  else
2052  {
2053  // both are edges adjacent to triangle tri
2054  // there are several configurations possible for n1, n2, between conn3 nodes.
2055  int num1, num2, sense, offset;
2056  rval = _mbImpl->side_number( tri, e1, num1, sense, offset );
2057  MBERRORR( rval, "edge not adjacent" );
2058 
2059  rval = _mbImpl->side_number( tri, e2, num2, sense, offset );
2060  MBERRORR( rval, "edge not adjacent" );
2061 
2062  const EntityHandle* conn12; // connectivity for edge 1
2063  const EntityHandle* conn22; // connectivity for edge 2
2064  // int nnodes;
2065  rval = _mbImpl->get_connectivity( e1, conn12, nnodes );
2066  MBERRORR( rval, "Failed to get connectivity of edge 1" );
2067  assert( 2 == nnodes );
2068  rval = _mbImpl->get_connectivity( e2, conn22, nnodes );
2069  MBERRORR( rval, "Failed to get connectivity of edge 2" );
2070  assert( 2 == nnodes );
2071  // now, having the side number, work through
2072  if( debug_splits )
2073  {
2074  std::cout << "tri conn3:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n";
2075  std::cout << " edge1: conn12:" << conn12[0] << " " << conn12[1] << " side: " << num1 << "\n";
2076  std::cout << " edge2: conn22:" << conn22[0] << " " << conn22[1] << " side: " << num2 << "\n";
2077  }
2078  int unaffectedSide = 3 - num1 - num2;
2079  int i3 = ( unaffectedSide + 2 ) % 3; // to 0 is 2, to 1 is 0, to 2 is 1
2080  // triangles will be formed with triVertexIndex , n1, n2 (in what order?)
2081  EntityHandle v1, v2; // to hold the 2 nodes on edges
2082  if( num1 == i3 )
2083  {
2084  v1 = n1;
2085  v2 = n2;
2086  }
2087  else // if (num2==i3)
2088  {
2089  v1 = n2;
2090  v2 = n1;
2091  }
2092  // three triangles are formed
2093  int i1 = ( i3 + 1 ) % 3;
2094  int i2 = ( i3 + 2 ) % 3;
2095  // we could break the surface differently
2096  EntityHandle conn[9] = { conn3[i3], v1, v2, v1, conn3[i1], conn3[i2], v2, v1, conn3[i2] };
2097  EntityHandle newTriangle;
2098  if( debug_splits ) std::cout << "Split 2 edges :\n";
2099  rval = _mbImpl->create_element( MBTRI, conn, 3, newTriangle );
2100  MBERRORR( rval, "Failed to create a new triangle" );
2101  _newTriangles.insert( newTriangle );
2102  if( debug_splits ) print_debug_triangle( newTriangle );
2103  rval = _mbImpl->create_element( MBTRI, conn + 3, 3, newTriangle ); // the second triangle
2104  MBERRORR( rval, "Failed to create a new triangle" );
2105  _newTriangles.insert( newTriangle );
2106  if( debug_splits ) print_debug_triangle( newTriangle );
2107  rval = _mbImpl->create_element( MBTRI, conn + 6, 3, newTriangle ); // the second triangle
2108  MBERRORR( rval, "Failed to create a new triangle" );
2109  _newTriangles.insert( newTriangle );
2110  if( debug_splits ) print_debug_triangle( newTriangle );
2111  }
2112 
2113  return MB_SUCCESS;
2114 }

References _mbImpl, _newTriangles, moab::Interface::create_element(), moab::debug_splits, ErrorCode, moab::Interface::get_connectivity(), moab::index, 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().

◆ chain_able_edge()

ErrorCode moab::FBEngine::chain_able_edge ( EntityHandle  edge,
double  min_dot,
EntityHandle next_edge,
bool &  chainable 
)
private

Definition at line 3661 of file FBEngine.cpp.

3662 {
3663  // get the end, then get all parents of end
3664  // see if some are the starts of
3665  chainable = false;
3666  EntityHandle v1, v2;
3667  ErrorCode rval = get_vert_edges( edge, v1, v2 );
3668  MBERRORR( rval, "can't get vertices" );
3669  if( v1 == v2 ) return MB_SUCCESS; // it is periodic, can't chain it with another edge!
3670 
3671  // v2 is a set, get its parents, which should be edges
3672  Range edges;
3673  rval = _mbImpl->get_parent_meshsets( v2, edges );
3674  MBERRORR( rval, "can't get parents of vertex set" );
3675  // get parents of current edge (faces)
3676  Range faces;
3677  rval = _mbImpl->get_parent_meshsets( edge, faces );
3678  MBERRORR( rval, "can't get parents of edge set" );
3679  // get also the last edge "tangent" at the vertex
3680  std::vector< EntityHandle > mesh_edges;
3681  rval = _mbImpl->get_entities_by_type( edge, MBEDGE, mesh_edges );
3682  MBERRORR( rval, "can't get mesh edges from edge set" );
3683  EntityHandle lastMeshEdge = mesh_edges[mesh_edges.size() - 1];
3684  const EntityHandle* conn2 = NULL;
3685  int len = 0;
3686  rval = _mbImpl->get_connectivity( lastMeshEdge, conn2, len );
3687  MBERRORR( rval, "can't connectivity of last mesh edge" );
3688  // get the coordinates of last edge
3689  if( len != 2 ) MBERRORR( MB_FAILURE, "bad number of vertices" );
3690  CartVect P[2];
3691  rval = _mbImpl->get_coords( conn2, len, (double*)&P[0] );
3692  MBERRORR( rval, "Failed to get coordinates" );
3693 
3694  CartVect vec1( P[1] - P[0] );
3695  vec1.normalize();
3696  for( Range::iterator edgeIter = edges.begin(); edgeIter != edges.end(); ++edgeIter )
3697  {
3698  EntityHandle otherEdge = *edgeIter;
3699  if( edge == otherEdge ) continue;
3700  // get faces parents of this edge
3701  Range faces2;
3702  rval = _mbImpl->get_parent_meshsets( otherEdge, faces2 );
3703  MBERRORR( rval, "can't get parents of other edge set" );
3704  if( faces != faces2 ) continue;
3705  // now, if the first mesh edge is within given angle, we can go on
3706  std::vector< EntityHandle > mesh_edges2;
3707  rval = _mbImpl->get_entities_by_type( otherEdge, MBEDGE, mesh_edges2 );
3708  MBERRORR( rval, "can't get mesh edges from other edge set" );
3709  EntityHandle firstMeshEdge = mesh_edges2[0];
3710  const EntityHandle* conn22;
3711  int len2;
3712  rval = _mbImpl->get_connectivity( firstMeshEdge, conn22, len2 );
3713  MBERRORR( rval, "can't connectivity of first mesh edge" );
3714  if( len2 != 2 ) MBERRORR( MB_FAILURE, "bad number of vertices" );
3715  if( conn2[1] != conn22[0] ) continue; // the mesh edges are not one after the other
3716  // get the coordinates of first edge
3717 
3718  // CartVect P2[2];
3719  rval = _mbImpl->get_coords( conn22, len, (double*)&P[0] );
3720  CartVect vec2( P[1] - P[0] );
3721  vec2.normalize();
3722  if( vec1 % vec2 < min_dot ) continue;
3723  // we found our edge, victory! we can get out
3724  next_edge = otherEdge;
3725  chainable = true;
3726  return MB_SUCCESS;
3727  }
3728 
3729  return MB_SUCCESS; // in general, hard to come by chain-able edges
3730 }

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().

◆ chain_edges()

ErrorCode moab::FBEngine::chain_edges ( double  min_dot)

Definition at line 3623 of file FBEngine.cpp.

3624 {
3625  Range sets[5];
3626  ErrorCode rval;
3627  while( 1 ) // break out only if no edges are chained
3628  {
3629  rval = _my_geomTopoTool->find_geomsets( sets );
3630  MBERRORR( rval, "can't get geo sets" );
3631  // these gentities are "always" current, while the ones in this-> _my_gsets[0:4] are
3632  // the "originals" before FBEngine modifications
3633  int nedges = (int)sets[1].size();
3634  // as long as we have chainable edges, continue;
3635  bool chain = false;
3636  for( int i = 0; i < nedges; i++ )
3637  {
3638  EntityHandle edge = sets[1][i];
3639  EntityHandle next_edge;
3640  bool chainable = false;
3641  rval = chain_able_edge( edge, min_dot, next_edge, chainable );
3642  MBERRORR( rval, "can't determine chain-ability" );
3643  if( chainable )
3644  {
3645  rval = chain_two_edges( edge, next_edge );
3646  MBERRORR( rval, "can't chain 2 edges" );
3647  chain = true;
3648  break; // interrupt for loop
3649  }
3650  }
3651  if( !chain )
3652  {
3653  break; // break out of while loop
3654  }
3655  }
3656  return MB_SUCCESS;
3657 }

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

Referenced by split_surface_with_direction().

◆ chain_two_edges()

ErrorCode moab::FBEngine::chain_two_edges ( EntityHandle  edge,
EntityHandle  next_edge 
)

Definition at line 3731 of file FBEngine.cpp.

3732 {
3733  // the biggest thing is to see the sense tags; or maybe not...
3734  // they should be correct :)
3735  // get the vertex sets
3736  EntityHandle v11, v12, v21, v22;
3737  ErrorCode rval = get_vert_edges( edge, v11, v12 );
3738  MBERRORR( rval, "can't get vert sets" );
3739  rval = get_vert_edges( next_edge, v21, v22 );
3740  MBERRORR( rval, "can't get vert sets" );
3741  assert( v12 == v21 );
3742  std::vector< EntityHandle > mesh_edges;
3743  rval = MBI->get_entities_by_type( next_edge, MBEDGE, mesh_edges );
3744  MBERRORR( rval, "can't get mesh edges" );
3745 
3746  rval = _mbImpl->add_entities( edge, &mesh_edges[0], (int)mesh_edges.size() );
3747  MBERRORR( rval, "can't add new mesh edges" );
3748  // remove the child - parent relation for second vertex of first edge
3749  rval = _mbImpl->remove_parent_child( edge, v12 );
3750  MBERRORR( rval, "can't remove parent - child relation between first edge and middle vertex" );
3751 
3752  if( v22 != v11 ) // the edge would become periodic, do not add again the relationship
3753  {
3754  rval = _mbImpl->add_parent_child( edge, v22 );
3755  MBERRORR( rval, "can't add second vertex to edge " );
3756  }
3757  // we can now safely eliminate next_edge
3758  rval = _mbImpl->remove_parent_child( next_edge, v21 );
3759  MBERRORR( rval, "can't remove child - parent relation " );
3760 
3761  rval = _mbImpl->remove_parent_child( next_edge, v22 );
3762  MBERRORR( rval, "can't remove child - parent relation " );
3763 
3764  // remove the next_edge relation to the faces
3765  Range faces;
3766  rval = _mbImpl->get_parent_meshsets( next_edge, faces );
3767  MBERRORR( rval, "can't get parent faces " );
3768 
3769  for( Range::iterator it = faces.begin(); it != faces.end(); ++it )
3770  {
3771  EntityHandle ff = *it;
3772  rval = _mbImpl->remove_parent_child( ff, next_edge );
3773  MBERRORR( rval, "can't remove parent-edge rel " );
3774  }
3775 
3776  rval = _mbImpl->delete_entities( &next_edge, 1 );
3777  MBERRORR( rval, "can't remove edge set " );
3778 
3779  // delete the vertex set that is idle now (v12 = v21)
3780  rval = _mbImpl->delete_entities( &v12, 1 );
3781  MBERRORR( rval, "can't remove edge set " );
3782  return MB_SUCCESS;
3783 }

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().

◆ clean()

void moab::FBEngine::clean ( )

Definition at line 210 of file FBEngine.cpp.

211 {
212  if( _smooth )
213  {
214  _faces.clear();
215  _edges.clear();
216  int size1 = _my_gsets[1].size();
217  int i = 0;
218  for( i = 0; i < size1; i++ )
219  delete _smthCurve[i];
220  delete[] _smthCurve;
221  _smthCurve = NULL;
222  size1 = _my_gsets[2].size();
223  for( i = 0; i < size1; i++ )
224  delete _smthFace[i];
225  delete[] _smthFace;
226  _smthFace = NULL;
227  //_smooth = false;
228  }
229 
230  for( int j = 0; j < 5; j++ )
231  _my_gsets[j].clear();
232  if( _t_created ) delete _my_geomTopoTool;
233  _my_geomTopoTool = NULL;
234  _t_created = false;
235 }

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

Referenced by ~FBEngine().

◆ compute_intersection_points()

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 2120 of file FBEngine.cpp.

2127 {
2128  // keep a stack of triangles to process, and do not add those already processed
2129  // either mark them, or maybe keep them in a local set?
2130  // from and to are now nodes, start from them
2131  CartVect p1, p2; // the position of from and to
2132  ErrorCode rval = _mbImpl->get_coords( &from, 1, (double*)&p1 );
2133  MBERRORR( rval, "failed to get 'from' coordinates" );
2134  rval = _mbImpl->get_coords( &to, 1, (double*)&p2 );
2135  MBERRORR( rval, "failed to get 'from' coordinates" );
2136 
2137  CartVect vect( p2 - p1 );
2138  double dist2 = vect.length();
2139  if( dist2 < tolerance_segment )
2140  {
2141  // we are done, return
2142  return MB_SUCCESS;
2143  }
2144  CartVect normPlane = Dir * vect;
2145  normPlane.normalize();
2146  std::set< EntityHandle > visitedTriangles;
2147  CartVect currentPoint = p1;
2148  // push the current point if it is empty
2149  if( points.size() == 0 )
2150  {
2151  points.push_back( p1 );
2152  entities.push_back( from ); // this is a node now
2153  }
2154 
2155  // these will be used in many loops
2156  CartVect intx = p1; // somewhere to start
2157  double param = -1.;
2158 
2159  // first intersection
2160  EntityHandle currentBoundary = from; // it is a node, in the beginning
2161 
2162  vect = p2 - currentPoint;
2163  while( vect.length() > 0. )
2164  {
2165  // advance towards "to" node, from boundary handle
2166  EntityType etype = _mbImpl->type_from_handle( currentBoundary );
2167  // if vertex, look for other triangles connected which intersect our plane (defined by p1,
2168  // p2, dir)
2169  std::vector< EntityHandle > adj_tri;
2170  rval = _mbImpl->get_adjacencies( &currentBoundary, 1, 2, false, adj_tri );
2171  unsigned int j = 0;
2172  EntityHandle tri;
2173  for( ; j < adj_tri.size(); j++ )
2174  {
2175  tri = adj_tri[j];
2176  if( visitedTriangles.find( tri ) != visitedTriangles.end() )
2177  continue; // get another triangle, this one was already visited
2178  // check if it is one of the triangles that was pierced already
2179  if( _piercedTriangles.find( tri ) != _piercedTriangles.end() ) continue;
2180  // if vertex, look for opposite edge
2181  // if edge, look for 2 opposite edges
2182  // get vertices
2183  int nnodes;
2184  const EntityHandle* conn3;
2185  rval = _mbImpl->get_connectivity( tri, conn3, nnodes );
2186  MBERRORR( rval, "Failed to get connectivity" );
2187  // if one of the nodes is to, stop right there
2188  {
2189  if( conn3[0] == to || conn3[1] == to || conn3[2] == to )
2190  {
2191  visitedTriangles.insert( tri );
2192  triangles.push_back( tri );
2193  currentPoint = p2;
2194  points.push_back( p2 );
2195  entities.push_back( to ); // we are done
2196  break; // this is break from for loop, we still need to get out of while
2197  // we will get out, because vect will become 0, (p2-p2)
2198  }
2199  }
2200  EntityHandle nn2[2];
2201  if( MBVERTEX == etype )
2202  {
2203  nn2[0] = conn3[0];
2204  nn2[1] = conn3[1];
2205  if( nn2[0] == currentBoundary ) nn2[0] = conn3[2];
2206  if( nn2[1] == currentBoundary ) nn2[1] = conn3[2];
2207  // get coordinates
2208  CartVect Pt[2];
2209 
2210  rval = _mbImpl->get_coords( nn2, 2, (double*)&Pt[0] );
2211  MBERRORR( rval, "Failed to get coordinates" );
2212  // see the intersection
2213  if( intersect_segment_and_plane_slice( Pt[0], Pt[1], currentPoint, p2, Dir, normPlane, intx, param ) )
2214  {
2215  // we should stop for loop, and decide if it is edge or vertex
2216  if( param == 0. )
2217  currentBoundary = nn2[0];
2218  else
2219  {
2220  if( param == 1. )
2221  currentBoundary = nn2[1];
2222  else // param between 0 and 1, so edge
2223  {
2224  // find the edge between vertices
2225  std::vector< EntityHandle > edges1;
2226  // change the create flag to true, because that edge must exist in
2227  // current triangle if we want to advance; nn2 are 2 nodes in current
2228  // triangle!!
2229  rval = _mbImpl->get_adjacencies( nn2, 2, 1, true, edges1, Interface::INTERSECT );
2230  MBERRORR( rval, "Failed to get edges" );
2231  if( edges1.size() != 1 ) MBERRORR( MB_FAILURE, "Failed to get adjacent edges to 2 nodes" );
2232  currentBoundary = edges1[0];
2233  }
2234  }
2235  visitedTriangles.insert( tri );
2236  currentPoint = intx;
2237  points.push_back( intx );
2238  entities.push_back( currentBoundary );
2239  triangles.push_back( tri );
2240  if( debug_splits )
2241  std::cout << "vtx new tri : " << _mbImpl->id_from_handle( tri )
2242  << " type bdy:" << _mbImpl->type_from_handle( currentBoundary ) << "\n";
2243  break; // out of for loop over triangles
2244  }
2245  }
2246  else // this is MBEDGE, we have the other triangle to try out
2247  {
2248  // first find the nodes from existing boundary
2249  int nnodes2;
2250  const EntityHandle* conn2;
2251  rval = _mbImpl->get_connectivity( currentBoundary, conn2, nnodes2 );
2252  MBERRORR( rval, "Failed to get connectivity" );
2253  int thirdIndex = -1;
2254  for( int tj = 0; tj < 3; tj++ )
2255  {
2256  if( ( conn3[tj] != conn2[0] ) && ( conn3[tj] != conn2[1] ) )
2257  {
2258  thirdIndex = tj;
2259  break;
2260  }
2261  }
2262  if( -1 == thirdIndex ) MBERRORR( MB_FAILURE, " can't get third node" );
2263  CartVect Pt[3];
2264  rval = _mbImpl->get_coords( conn3, 3, (double*)&Pt[0] );
2265  MBERRORR( rval, "Failed to get coords" );
2266  int indexFirst = ( thirdIndex + 1 ) % 3;
2267  int indexSecond = ( thirdIndex + 2 ) % 3;
2268  int index[2] = { indexFirst, indexSecond };
2269  for( int k = 0; k < 2; k++ )
2270  {
2271  nn2[0] = conn3[index[k]], nn2[1] = conn3[thirdIndex];
2272  if( intersect_segment_and_plane_slice( Pt[index[k]], Pt[thirdIndex], currentPoint, p2, Dir,
2273  normPlane, intx, param ) )
2274  {
2275  // we should stop for loop, and decide if it is edge or vertex
2276  if( param == 0. )
2277  currentBoundary = conn3[index[k]]; // it is not really possible
2278  else
2279  {
2280  if( param == 1. )
2281  currentBoundary = conn3[thirdIndex];
2282  else // param between 0 and 1, so edge is fine
2283  {
2284  // find the edge between vertices
2285  std::vector< EntityHandle > edges1;
2286  // change the create flag to true, because that edge must exist in
2287  // current triangle if we want to advance; nn2 are 2 nodes in
2288  // current triangle!!
2289  rval = _mbImpl->get_adjacencies( nn2, 2, 1, true, edges1, Interface::INTERSECT );
2290  MBERRORR( rval, "Failed to get edges" );
2291  if( edges1.size() != 1 )
2292  MBERRORR( MB_FAILURE, "Failed to get adjacent edges to 2 nodes" );
2293  currentBoundary = edges1[0];
2294  }
2295  }
2296  visitedTriangles.insert( tri );
2297  currentPoint = intx;
2298  points.push_back( intx );
2299  entities.push_back( currentBoundary );
2300  triangles.push_back( tri );
2301  if( debug_splits )
2302  {
2303  std::cout << "edge new tri : " << _mbImpl->id_from_handle( tri )
2304  << " type bdy: " << _mbImpl->type_from_handle( currentBoundary )
2305  << " id: " << _mbImpl->id_from_handle( currentBoundary ) << "\n";
2306  _mbImpl->list_entity( currentBoundary );
2307  }
2308  break; // out of for loop over triangles
2309  }
2310  // we should not reach here
2311  }
2312  }
2313  }
2314  /*if (j==adj_tri.size())
2315  {
2316  MBERRORR(MB_FAILURE, "did not advance");
2317  }*/
2318  vect = p2 - currentPoint;
2319  }
2320 
2321  if( debug_splits )
2322  std::cout << "nb entities: " << entities.size() << " triangles:" << triangles.size()
2323  << " points.size(): " << points.size() << "\n";
2324 
2325  return MB_SUCCESS;
2326 }

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::index, 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().

◆ create_new_gedge()

ErrorCode moab::FBEngine::create_new_gedge ( std::vector< EntityHandle > &  nodesAlongPolyline,
EntityHandle new_geo_edge 
)
private

Definition at line 1868 of file FBEngine.cpp.

1869 {
1870 
1871  ErrorCode rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_geo_edge );
1872  MBERRORR( rval, "can't create geo edge" );
1873 
1874  // now, get the edges, or create if not existing
1875  std::vector< EntityHandle > mesh_edges;
1876  for( unsigned int i = 0; i < nodesAlongPolyline.size() - 1; i++ )
1877  {
1878  EntityHandle n1 = nodesAlongPolyline[i], n2 = nodesAlongPolyline[i + 1];
1879 
1880  EntityHandle nn2[2];
1881  nn2[0] = n1;
1882  nn2[1] = n2;
1883 
1884  std::vector< EntityHandle > adjacent;
1885  rval = _mbImpl->get_adjacencies( nn2, 2, 1, false, adjacent, Interface::INTERSECT );
1886  // see if there is an edge between those 2 already, and if it is oriented as we like
1887  bool new_edge = true;
1888  if( adjacent.size() >= 1 )
1889  {
1890  // check the orientation
1891  const EntityHandle* conn2 = NULL;
1892  int nnodes = 0;
1893  rval = _mbImpl->get_connectivity( adjacent[0], conn2, nnodes );
1894  MBERRORR( rval, "can't get connectivity" );
1895  if( conn2[0] == nn2[0] && conn2[1] == nn2[1] )
1896  {
1897  // everything is fine
1898  mesh_edges.push_back( adjacent[0] );
1899  new_edge = false; // we found one that's good, no need to create a new one
1900  }
1901  else
1902  {
1903  _piercedEdges.insert( adjacent[0] ); // we want to remove this one, it will be not needed
1904  }
1905  }
1906  if( new_edge )
1907  {
1908  // there is no edge between n1 and n2, create one
1909  EntityHandle mesh_edge;
1910  rval = _mbImpl->create_element( MBEDGE, nn2, 2, mesh_edge );
1911  MBERRORR( rval, "Failed to create a new edge" );
1912  mesh_edges.push_back( mesh_edge );
1913  }
1914  }
1915 
1916  // add loops edges to the edge set
1917  rval = _mbImpl->add_entities( new_geo_edge, &mesh_edges[0],
1918  mesh_edges.size() ); // only one edge
1919  MBERRORR( rval, "can't add edges to new_geo_set" );
1920  // check vertex sets for vertex 1 and vertex 2?
1921  // get all sets of dimension 0 from database, and see if our ends are here or not
1922 
1923  Range ends_geo_edge;
1924  ends_geo_edge.insert( nodesAlongPolyline[0] );
1925  ends_geo_edge.insert( nodesAlongPolyline[nodesAlongPolyline.size() - 1] );
1926 
1927  for( unsigned int k = 0; k < ends_geo_edge.size(); k++ )
1928  {
1929  EntityHandle node = ends_geo_edge[k];
1930  EntityHandle nodeSet;
1931  bool found = find_vertex_set_for_node( node, nodeSet );
1932 
1933  if( !found )
1934  {
1935  // create a node set and add the node
1936 
1937  rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet );
1938  MBERRORR( rval, "Failed to create a new vertex set" );
1939 
1940  rval = _mbImpl->add_entities( nodeSet, &node, 1 );
1941  MBERRORR( rval, "Failed to add the node to the set" );
1942 
1943  rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 ); //
1944  MBERRORR( rval, "Failed to commit the node set" );
1945 
1946  if( debug_splits )
1947  {
1948  std::cout << " create a vertex set " << _mbImpl->id_from_handle( nodeSet )
1949  << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << " for node " << node
1950  << "\n";
1951  }
1952  }
1953 
1954  rval = _mbImpl->add_parent_child( new_geo_edge, nodeSet );
1955  MBERRORR( rval, "Failed to add parent child relation" );
1956  }
1957  // finally, put the edge in the range of edges
1958  MB_CHK_ERR( _my_geomTopoTool->add_geo_set( new_geo_edge, 1 ) );
1959 
1960  return rval;
1961 }

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().

◆ create_volume_with_direction()

ErrorCode moab::FBEngine::create_volume_with_direction ( EntityHandle  newFace1,
EntityHandle  newFace2,
double *  direction,
EntityHandle volume 
)

Definition at line 3187 of file FBEngine.cpp.

3191 {
3192 
3193  // MESHSET
3194  // ErrorCode rval = _mbImpl->create_meshset(MESHSET_ORDERED, new_geo_edge);
3195  ErrorCode rval = _mbImpl->create_meshset( MESHSET_SET, volume );
3196  MBERRORR( rval, "can't create volume" );
3197 
3198  int volumeMatId = 1; // just give a mat id, for debugging, mostly
3199  Tag matTag;
3201  MBERRORR( rval, "can't get material tag" );
3202 
3203  rval = _mbImpl->tag_set_data( matTag, &volume, 1, &volumeMatId );
3204  MBERRORR( rval, "can't set material tag value on volume" );
3205 
3206  // get the edges of those 2 faces, and get the vertices of those edges
3207  // in order, they should be created in the same order (?); check for that, anyway
3208  rval = _mbImpl->add_parent_child( volume, newFace1 );
3209  MBERRORR( rval, "can't add first face to volume" );
3210 
3211  rval = _mbImpl->add_parent_child( volume, newFace2 );
3212  MBERRORR( rval, "can't add second face to volume" );
3213 
3214  // first is bottom, so it is negatively oriented
3215  rval = _my_geomTopoTool->add_geo_set( volume, 3 );
3216  MBERRORR( rval, "can't add volume to the gtt" );
3217 
3218  // set senses
3219  // bottom face is negatively oriented, its normal is toward interior of the volume
3220  rval = _my_geomTopoTool->set_sense( newFace1, volume, -1 );
3221  MBERRORR( rval, "can't set bottom face sense to the volume" );
3222 
3223  // the top face is positively oriented
3224  rval = _my_geomTopoTool->set_sense( newFace2, volume, 1 );
3225  MBERRORR( rval, "can't set top face sense to the volume" );
3226 
3227  // the children should be in the same direction
3228  // get the side edges of each face, and form lateral faces, along direction
3229  std::vector< EntityHandle > edges1;
3230  std::vector< EntityHandle > edges2;
3231 
3232  rval = _mbImpl->get_child_meshsets( newFace1, edges1 ); // no hops
3233  MBERRORR( rval, "can't get children edges or first face, bottom" );
3234 
3235  rval = _mbImpl->get_child_meshsets( newFace2, edges2 ); // no hops
3236  MBERRORR( rval, "can't get children edges for second face, top" );
3237 
3238  if( edges1.size() != edges2.size() ) MBERRORR( MB_FAILURE, "wrong correspondence " );
3239 
3240  for( unsigned int i = 0; i < edges1.size(); ++i )
3241  {
3242  EntityHandle newLatFace;
3243  rval = weave_lateral_face_from_edges( edges1[i], edges2[i], direction, newLatFace );
3244  MBERRORR( rval, "can't weave lateral face" );
3245  rval = _mbImpl->add_parent_child( volume, newLatFace );
3246  MBERRORR( rval, "can't add lateral face to volume" );
3247 
3248  // set sense as positive
3249  rval = _my_geomTopoTool->set_sense( newLatFace, volume, 1 );
3250  MBERRORR( rval, "can't set lateral face sense to the volume" );
3251  }
3252 
3253  rval = set_default_neumann_tags();
3254  MBERRORR( rval, "can't set new neumann tags" );
3255 
3256  return MB_SUCCESS;
3257 }

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().

◆ createEntSet()

ErrorCode moab::FBEngine::createEntSet ( int  isList,
EntityHandle pSet 
)

Definition at line 448 of file FBEngine.cpp.

449 {
450  ErrorCode rval;
451 
452  if( isList )
453  rval = MBI->create_meshset( MESHSET_ORDERED, *pSet );
454  else
455  rval = MBI->create_meshset( MESHSET_SET, *pSet );
456 
457  return rval;
458 }

References ErrorCode, MBI, and MESHSET_SET.

◆ createTag()

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

Definition at line 916 of file FBEngine.cpp.

917 {
918  // Create a tag with the specified name, size, and type
919  // This function provides a simplified interface for tag creation
920 
922  MB_TYPE_HANDLE };
923  moab::TagType storage = MB_TAG_SPARSE;
924  ErrorCode result;
925 
926  result =
927  MBI->tag_get_handle( tag_name, tag_size, mb_data_type_table2[tag_type], tag_handle_out, storage | MB_TAG_EXCL );
928 
929  if( MB_SUCCESS != result )
930  {
931  std::string msg( "createTag: " );
932  if( MB_ALREADY_ALLOCATED == result )
933  {
934  msg += "Tag already exists with name: \"";
935  msg += tag_name;
936  std::cout << msg << "\n";
937  }
938  else
939  {
940  std::cout << "Failed to create tag with name: " << tag_name << "\n";
941  return MB_FAILURE;
942  }
943  }
944 
945  return MB_SUCCESS;
946 }

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.

◆ delete_smooth_tags()

void moab::FBEngine::delete_smooth_tags ( )

Definition at line 389 of file FBEngine.cpp.

390 {
391  // get all tags from database that are created for smooth data, and
392  // delete them; it will delete all data associated with them
393  // first tags from faces, edges:
394  std::vector< Tag > smoothTags;
395  int size1 = (int)_my_gsets[2].size();
396 
397  for( int i = 0; i < size1; i++ )
398  {
399  // these 2 will append gradient tag and plane tag
400  _smthFace[i]->append_smooth_tags( smoothTags );
401  }
402  // then , get other tags:
403  // "TANGENTS", "MARKER", "CONTROLEDGE", "CONTROLFACE", "CONTROLEDGEFACE"
404  Tag tag_handle;
405  ErrorCode rval = _mbImpl->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tag_handle );
406  if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
407 
408  rval = _mbImpl->tag_get_handle( "MARKER", 1, MB_TYPE_BIT, tag_handle );
409  if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
410 
411  rval = _mbImpl->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, tag_handle );
412  if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
413 
414  rval = _mbImpl->tag_get_handle( "CONTROLFACE", 18, MB_TYPE_DOUBLE, tag_handle );
415  if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
416 
417  rval = _mbImpl->tag_get_handle( "CONTROLEDGEFACE", 27, MB_TYPE_DOUBLE, tag_handle );
418  if( rval != MB_TAG_NOT_FOUND ) smoothTags.push_back( tag_handle );
419 
420  // a lot of tags, delete them
421  for( unsigned int k = 0; k < smoothTags.size(); k++ )
422  {
423  // could be a lot of data
424  _mbImpl->tag_delete( smoothTags[k] );
425  }
426 }

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

◆ divide_triangle()

ErrorCode moab::FBEngine::divide_triangle ( EntityHandle  triangle,
EntityHandle newVertex 
)
private

Definition at line 3141 of file FBEngine.cpp.

3142 {
3143  //
3144  _piercedTriangles.insert( triangle );
3145  int nnodes = 0;
3146  const EntityHandle* conn3 = NULL;
3147  ErrorCode rval = _mbImpl->get_connectivity( triangle, conn3, nnodes );
3148  MBERRORR( rval, "can't get nodes" );
3149  EntityHandle t1[] = { conn3[0], conn3[1], newVertex };
3150  EntityHandle t2[] = { conn3[1], conn3[2], newVertex };
3151  EntityHandle t3[] = { conn3[2], conn3[0], newVertex };
3152  EntityHandle newTriangle, newTriangle2, newTriangle3;
3153  rval = _mbImpl->create_element( MBTRI, t1, 3, newTriangle );
3154  MBERRORR( rval, "can't create triangle" );
3155  _newTriangles.insert( newTriangle );
3156  rval = _mbImpl->create_element( MBTRI, t2, 3, newTriangle3 );
3157  MBERRORR( rval, "can't create triangle" );
3158  _newTriangles.insert( newTriangle3 );
3159  rval = _mbImpl->create_element( MBTRI, t3, 3, newTriangle2 );
3160  MBERRORR( rval, "can't create triangle" );
3161  _newTriangles.insert( newTriangle2 );
3162 
3163  // create all edges
3164  std::vector< EntityHandle > edges0;
3165  rval = _mbImpl->get_adjacencies( &newTriangle, 1, 1, true, edges0 );
3166  MBERRORR( rval, "can't get new edges" );
3167  edges0.clear();
3168  rval = _mbImpl->get_adjacencies( &newTriangle2, 1, 1, true, edges0 );
3169  MBERRORR( rval, "can't get new edges" );
3170  if( debug_splits )
3171  {
3172  std::cout << "3 triangles formed:\n";
3173  _mbImpl->list_entity( newTriangle );
3174  print_debug_triangle( newTriangle );
3175  _mbImpl->list_entity( newTriangle3 );
3176  print_debug_triangle( newTriangle3 );
3177  _mbImpl->list_entity( newTriangle2 );
3178  print_debug_triangle( newTriangle2 );
3179  std::cout << "original nodes in tri:\n";
3180  _mbImpl->list_entity( conn3[0] );
3181  _mbImpl->list_entity( conn3[1] );
3182  _mbImpl->list_entity( conn3[2] );
3183  }
3184  return MB_SUCCESS;
3185 }

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().

◆ find_vertex_set_for_node()

bool moab::FBEngine::find_vertex_set_for_node ( EntityHandle  iNode,
EntityHandle oVertexSet 
)
private

Definition at line 2843 of file FBEngine.cpp.

2844 {
2845  bool found = false;
2846  Range vertex_sets;
2847 
2848  const int zero = 0;
2849  const void* const zero_val[] = { &zero };
2850  Tag geom_tag;
2851  ErrorCode rval = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag );
2852  if( MB_SUCCESS != rval ) return false;
2853  rval = _mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom_tag, zero_val, 1, vertex_sets );
2854  if( MB_SUCCESS != rval ) return false;
2855  // local _gsets, as we might have not updated the local lists
2856  // see if ends of geo edge generated is in a geo set 0 or not
2857 
2858  for( Range::iterator vsit = vertex_sets.begin(); vsit != vertex_sets.end(); ++vsit )
2859  {
2860  EntityHandle vset = *vsit;
2861  // is the node part of this set?
2862  if( _mbImpl->contains_entities( vset, &iNode, 1 ) )
2863  {
2864 
2865  found = true;
2866  oVertexSet = vset;
2867  break;
2868  }
2869  }
2870  return found;
2871 }

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().

◆ get_gtt()

GeomTopoTool* moab::FBEngine::get_gtt ( )
inline

Definition at line 183 of file FBEngine.hpp.

184  {
185  return this->_my_geomTopoTool;
186  }

References _my_geomTopoTool.

◆ get_nodes_from_edge()

ErrorCode moab::FBEngine::get_nodes_from_edge ( EntityHandle  gedge,
std::vector< EntityHandle > &  nodes 
)

Definition at line 3259 of file FBEngine.cpp.

3260 {
3261  std::vector< EntityHandle > ents;
3262  ErrorCode rval = _mbImpl->get_entities_by_type( gedge, MBEDGE, ents );
3263  if( MB_SUCCESS != rval ) return rval;
3264  if( ents.size() < 1 ) return MB_FAILURE;
3265 
3266  nodes.resize( ents.size() + 1 );
3267  const EntityHandle* conn = NULL;
3268  int len;
3269  for( unsigned int i = 0; i < ents.size(); ++i )
3270  {
3271  rval = _mbImpl->get_connectivity( ents[i], conn, len );
3272  MBERRORR( rval, "can't get edge connectivity" );
3273  nodes[i] = conn[0];
3274  }
3275  // the last one is conn[1]
3276  nodes[ents.size()] = conn[1];
3277  return MB_SUCCESS;
3278 }

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().

◆ get_vert_edges()

ErrorCode moab::FBEngine::get_vert_edges ( EntityHandle  edge,
EntityHandle v1,
EntityHandle v2 
)

Definition at line 3784 of file FBEngine.cpp.

3785 {
3786  // need to decide first or second vertex
3787  // important for moab
3788 
3789  Range children;
3790  // EntityHandle v1, v2;
3791  ErrorCode rval = _mbImpl->get_child_meshsets( edge, children );
3792  MBERRORR( rval, "can't get child meshsets" );
3793  if( children.size() == 1 )
3794  {
3795  // this is periodic edge, get out early
3796  v1 = children[0];
3797  v2 = v1;
3798  return MB_SUCCESS;
3799  }
3800  else if( children.size() > 2 )
3801  MBERRORR( MB_FAILURE, "too many vertices in one edge" );
3802  // edge: get one vertex as part of the vertex set
3803  Range entities;
3804  rval = MBI->get_entities_by_type( children[0], MBVERTEX, entities );
3805  MBERRORR( rval, "can't get entities from vertex set" );
3806  if( entities.size() < 1 ) MBERRORR( MB_FAILURE, "no mesh nodes in vertex set" );
3807  EntityHandle node0 = entities[0]; // the first vertex
3808  entities.clear();
3809 
3810  // now get the edges, and get the first node and the last node in sequence of edges
3811  // the order is important...
3812  // these are ordered sets !!
3813  std::vector< EntityHandle > ents;
3814  rval = MBI->get_entities_by_type( edge, MBEDGE, ents );
3815  MBERRORR( rval, "can't get mesh edges" );
3816  if( ents.size() < 1 ) MBERRORR( MB_FAILURE, "no mesh edges in edge set" );
3817 
3818  const EntityHandle* conn = NULL;
3819  int len;
3820  rval = MBI->get_connectivity( ents[0], conn, len );
3821  MBERRORR( rval, "can't connectivity of first mesh edge" );
3822 
3823  if( conn[0] == node0 )
3824  {
3825  v1 = children[0];
3826  v2 = children[1];
3827  }
3828  else // the other way around, although we should check (we are paranoid)
3829  {
3830  v2 = children[0];
3831  v1 = children[1];
3832  }
3833 
3834  return MB_SUCCESS;
3835 }

References _mbImpl, moab::Range::clear(), 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().

◆ getAdjacentEntities()

ErrorCode moab::FBEngine::getAdjacentEntities ( const EntityHandle  from,
const int  to_dim,
Range adj_ents 
)
private

Definition at line 850 of file FBEngine.cpp.

851 {
852  int this_dim = -1;
853  for( int i = 0; i < 4; i++ )
854  {
855  if( _my_geomTopoTool->geoRanges()[i].find( from ) != _my_geomTopoTool->geoRanges()[i].end() )
856  {
857  this_dim = i;
858  break;
859  }
860  }
861 
862  // check target dimension
863  if( -1 == this_dim )
864  {
865  // Entity not a geometry entity
866  return MB_FAILURE;
867  }
868  else if( 0 > to_dim || 3 < to_dim )
869  {
870  // To dimension must be between 0 and 3
871  return MB_FAILURE;
872  }
873  else if( to_dim == this_dim )
874  {
875  // To dimension must be different from entity dimension
876  return MB_FAILURE;
877  }
878 
879  ErrorCode rval = MB_SUCCESS;
880  adjs.clear();
881  if( to_dim > this_dim )
882  {
883  int diffDim = to_dim - this_dim;
884  rval = MBI->get_parent_meshsets( from, adjs, diffDim );
885  if( MB_SUCCESS != rval ) return rval;
886  if( diffDim > 1 )
887  {
888  // subtract the parents that come with diffDim-1 hops
889  Range extra;
890  rval = MBI->get_parent_meshsets( from, extra, diffDim - 1 );
891  if( MB_SUCCESS != rval ) return rval;
892  adjs = subtract( adjs, extra );
893  }
894  }
895  else
896  {
897  int diffDim = this_dim - to_dim;
898  rval = MBI->get_child_meshsets( from, adjs, diffDim );
899  if( MB_SUCCESS != rval ) return rval;
900  if( diffDim > 1 )
901  {
902  // subtract the children that come with diffDim-1 hops
903  Range extra;
904  rval = MBI->get_child_meshsets( from, extra, diffDim - 1 );
905  if( MB_SUCCESS != rval ) return rval;
906  adjs = subtract( adjs, extra );
907  }
908  }
909 
910  return rval;
911 }

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().

◆ getArrData()

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

Definition at line 948 of file FBEngine.cpp.

952 {
953  // responsibility of the user to have tag_values_out properly allocated
954  // only some types of Tags are possible (double, int, etc)
955  return MBI->tag_get_data( tag_handle, entity_handles, entity_handles_size, tag_values_out );
956 }

References MBI.

◆ getEgEvalXYZ()

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 1059 of file FBEngine.cpp.

1072 {
1073  return MB_NOT_IMPLEMENTED; // not implemented
1074 }

References MB_NOT_IMPLEMENTED.

◆ getEgFcSense()

ErrorCode moab::FBEngine::getEgFcSense ( EntityHandle  mbedge,
EntityHandle  mbface,
int &  sense 
)

Definition at line 973 of file FBEngine.cpp.

974 {
975 
976  // this one is important, for establishing the orientation of the edges in faces
977  // use senses
978  std::vector< EntityHandle > faces;
979  std::vector< int > senses; // 0 is forward and 1 is backward
980  ErrorCode rval = _my_geomTopoTool->get_senses( mbedge, faces, senses );
981  if( MB_SUCCESS != rval ) return rval;
982 
983  for( unsigned int i = 0; i < faces.size(); i++ )
984  {
985  if( faces[i] == mbface )
986  {
987  sense_out = senses[i];
988  return MB_SUCCESS;
989  }
990  }
991  return MB_FAILURE;
992 }

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

◆ getEgVtxSense()

ErrorCode moab::FBEngine::getEgVtxSense ( EntityHandle  edge,
EntityHandle  vtx1,
EntityHandle  vtx2,
int &  sense 
)

Definition at line 1095 of file FBEngine.cpp.

1096 {
1097  // need to decide first or second vertex
1098  // important for moab
1099  int type;
1100 
1101  EntityHandle v1, v2;
1102  ErrorCode rval = getEntType( vtx1, &type );
1103  if( MB_SUCCESS != rval || type != 0 ) return MB_FAILURE;
1104  // edge: get one vertex as part of the vertex set
1105  Range entities;
1106  rval = MBI->get_entities_by_type( vtx1, MBVERTEX, entities );
1107  if( MB_SUCCESS != rval ) return rval;
1108  if( entities.size() < 1 ) return MB_FAILURE;
1109  v1 = entities[0]; // the first vertex
1110  entities.clear();
1111  rval = getEntType( vtx2, &type );
1112  if( MB_SUCCESS != rval || type != 0 ) return MB_FAILURE;
1113  rval = MBI->get_entities_by_type( vtx2, MBVERTEX, entities );
1114  if( MB_SUCCESS != rval ) return rval;
1115  if( entities.size() < 1 ) return MB_FAILURE;
1116  v2 = entities[0]; // the first vertex
1117  entities.clear();
1118  // now get the edges, and get the first node and the last node in sequence of edges
1119  // the order is important...
1120  // these are ordered sets !!
1121  std::vector< EntityHandle > ents;
1122  rval = MBI->get_entities_by_type( edge, MBEDGE, ents );
1123  if( MB_SUCCESS != rval ) return rval;
1124  if( ents.size() < 1 ) return MB_FAILURE;
1125 
1126  const EntityHandle* conn = NULL;
1127  int len;
1128  EntityHandle startNode, endNode;
1129  rval = MBI->get_connectivity( ents[0], conn, len );
1130  if( MB_SUCCESS != rval ) return rval;
1131  startNode = conn[0];
1132  rval = MBI->get_connectivity( ents[ents.size() - 1], conn, len );
1133  if( MB_SUCCESS != rval ) return rval;
1134 
1135  endNode = conn[1];
1136  sense = 1; //
1137  if( ( startNode == endNode ) && ( v1 == startNode ) )
1138  {
1139  sense = 0; // periodic
1140  }
1141  if( ( startNode == v1 ) && ( endNode == v2 ) )
1142  {
1143  sense = 1; // forward
1144  }
1145  if( ( startNode == v2 ) && ( endNode == v1 ) )
1146  {
1147  sense = -1; // reverse
1148  }
1149  return MB_SUCCESS;
1150 }

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

Referenced by weave_lateral_face_from_edges().

◆ getEntAdj()

ErrorCode moab::FBEngine::getEntAdj ( EntityHandle  handle,
int  type_requested,
Range adjEnts 
)

Definition at line 968 of file FBEngine.cpp.

969 {
970  return getAdjacentEntities( handle, type_requested, adjEnts );
971 }

References getAdjacentEntities().

◆ getEntBoundBox()

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.

556 {
557  ErrorCode rval;
558  int type;
559  rval = getEntType( gent, &type );
560  MBERRORR( rval, "Failed to get entity type." );
561 
562  if( type == 0 )
563  {
564  rval = getVtxCoord( gent, min_x, min_y, min_z );
565  MBERRORR( rval, "Failed to get vertex coordinates." );
566  max_x = min_x;
567  max_y = min_y;
568  max_z = min_z;
569  }
570  else if( type == 1 )
571  {
572  MBERRORR( MB_FAILURE, "iGeom_getEntBoundBox is not supported for Edge entity type." );
573  }
574  else if( type == 2 || type == 3 )
575  {
576 
577  EntityHandle root;
578  CartVect center, axis[3];
579  rval = _my_geomTopoTool->get_root( gent, root );
580  MBERRORR( rval, "Failed to get tree root in iGeom_getEntBoundBox." );
581  rval = _my_geomTopoTool->obb_tree()->box( root, center.array(), axis[0].array(), axis[1].array(),
582  axis[2].array() );
583  MBERRORR( rval, "Failed to get closest point in iGeom_getEntBoundBox." );
584 
585  CartVect absv[3];
586  for( int i = 0; i < 3; i++ )
587  {
588  absv[i] = CartVect( fabs( axis[i][0] ), fabs( axis[i][1] ), fabs( axis[i][2] ) );
589  }
590  CartVect min, max;
591  min = center - absv[0] - absv[1] - absv[2];
592  max = center + absv[0] + absv[1] + absv[2];
593  *min_x = min[0];
594  *min_y = min[1];
595  *min_z = min[2];
596  *max_x = max[0];
597  *max_y = max[1];
598  *max_z = max[2];
599  }
600  else
601  return MB_FAILURE;
602 
603  return MB_SUCCESS;
604 }

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().

◆ getEntClosestPt()

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.

612 {
613  ErrorCode rval;
614  int type;
615  rval = getEntType( this_gent, &type );
616  MBERRORR( rval, "Failed to get entity type." );
617 
618  if( type == 0 )
619  {
620  rval = getVtxCoord( this_gent, on_x, on_y, on_z );
621  MBERRORR( rval, "Failed to get vertex coordinates." );
622  }
623  else if( _smooth && type == 1 )
624  {
625  *on_x = near_x;
626  *on_y = near_y;
627  *on_z = near_z;
628  SmoothCurve* smthcurve = _edges[this_gent];
629  // call the new method from smooth edge
630  smthcurve->move_to_curve( *on_x, *on_y, *on_z );
631  }
632  else if( type == 2 || type == 3 )
633  {
634  double point[3] = { near_x, near_y, near_z };
635  double point_out[3];
636  EntityHandle root, facet_out;
637  if( _smooth && 2 == type )
638  {
639  SmoothFace* smthFace = _faces[this_gent];
640  *on_x = near_x;
641  *on_y = near_y;
642  *on_z = near_z;
643  smthFace->move_to_surface( *on_x, *on_y, *on_z );
644  }
645  else
646  {
647  rval = _my_geomTopoTool->get_root( this_gent, root );
648  MBERRORR( rval, "Failed to get tree root in iGeom_getEntClosestPt." );
649  rval = _my_geomTopoTool->obb_tree()->closest_to_location( point, root, point_out, facet_out );
650  MBERRORR( rval, "Failed to get closest point in iGeom_getEntClosestPt." );
651 
652  *on_x = point_out[0];
653  *on_y = point_out[1];
654  *on_z = point_out[2];
655  }
656  }
657  else
658  return MB_TYPE_OUT_OF_RANGE;
659 
660  return MB_SUCCESS;
661 }

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().

◆ getEntities()

ErrorCode moab::FBEngine::getEntities ( EntityHandle  root_set,
int  ent_type,
Range gentities 
)

Definition at line 460 of file FBEngine.cpp.

461 {
462  int i;
463  if( 0 > entity_type || 4 < entity_type )
464  {
465  return MB_FAILURE;
466  }
467  else if( entity_type < 4 )
468  { // 4 means all entities
469  gentities = _my_geomTopoTool->geoRanges()[entity_type]; // all from root set!
470  }
471  else
472  {
473  gentities.clear();
474  for( i = 0; i < 4; i++ )
475  {
476  gentities.merge( _my_geomTopoTool->geoRanges()[i] );
477  }
478  }
479  Range sets;
480  // see now if they are in the set passed as input or not
481  ErrorCode rval = MBI->get_entities_by_type( set_handle, MBENTITYSET, sets );
482  MBERRORR( rval, "can't get sets in the initial set" );
483  gentities = intersect( gentities, sets );
484 
485  return MB_SUCCESS;
486 }

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

◆ getEntNrmlSense()

ErrorCode moab::FBEngine::getEntNrmlSense ( EntityHandle  face,
EntityHandle  region,
int &  sense 
)

Definition at line 1054 of file FBEngine.cpp.

1055 {
1056  return MB_NOT_IMPLEMENTED; // not implemented
1057 }

References MB_NOT_IMPLEMENTED.

◆ getEntNrmlXYZ()

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.

720 {
721  // just do for surface and volume
722  int type;
723  ErrorCode rval = getEntType( entity_handle, &type );
724  MBERRORR( rval, "Failed to get entity type in iGeom_getEntNrmlXYZ." );
725 
726  if( type != 2 && type != 3 )
727  {
728  MBERRORR( MB_FAILURE, "Entities passed into gentityNormal must be face or volume." );
729  }
730 
731  if( _smooth && 2 == type )
732  {
733  SmoothFace* smthFace = _faces[entity_handle];
734  //*on_x = near_x; *on_y = near_y; *on_z = near_z;
735  smthFace->normal_at( x, y, z, *nrml_i, *nrml_j, *nrml_k );
736  }
737  else
738  {
739  // get closest location and facet
740  double point[3] = { x, y, z };
741  double point_out[3];
742  EntityHandle root, facet_out;
743  _my_geomTopoTool->get_root( entity_handle, root );
744  rval = _my_geomTopoTool->obb_tree()->closest_to_location( point, root, point_out, facet_out );
745  MBERRORR( rval, "Failed to get closest location in iGeom_getEntNrmlXYZ." );
746 
747  // get facet normal
748  const EntityHandle* conn;
749  int len;
750  CartVect coords[3], normal;
751  rval = MBI->get_connectivity( facet_out, conn, len );
752  MBERRORR( rval, "Failed to get triangle connectivity in iGeom_getEntNrmlXYZ." );
753  if( len != 3 ) MBERRORR( MB_FAILURE, " not a triangle, error " );
754 
755  rval = MBI->get_coords( conn, len, coords[0].array() );
756  MBERRORR( rval, "Failed to get triangle coordinates in iGeom_getEntNrmlXYZ." );
757 
758  coords[1] -= coords[0];
759  coords[2] -= coords[0];
760  normal = coords[1] * coords[2];
761  normal.normalize();
762  *nrml_i = normal[0];
763  *nrml_j = normal[1];
764  *nrml_k = normal[2];
765  }
766  return MB_SUCCESS;
767 }

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().

◆ getEntTgntU()

ErrorCode moab::FBEngine::getEntTgntU ( EntityHandle  edge,
double  u,
double &  i,
double &  j,
double &  k 
)

Definition at line 1168 of file FBEngine.cpp.

1169 {
1170  SmoothCurve* smoothCurve = _edges[edge]; // this is a map
1171  // now, call smoothCurve methods
1172  double tg[3];
1173  double x, y, z;
1174  smoothCurve->position_from_u( u, x, y, z, tg );
1175  i = tg[0];
1176  j = tg[1];
1177  k = tg[2];
1178  return MB_SUCCESS;
1179 }

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

◆ getEntType()

ErrorCode moab::FBEngine::getEntType ( EntityHandle  gent,
int *  type 
)

Definition at line 536 of file FBEngine.cpp.

537 {
538  for( int i = 0; i < 4; i++ )
539  {
540  if( _my_geomTopoTool->geoRanges()[i].find( gent ) != _my_geomTopoTool->geoRanges()[i].end() )
541  {
542  *type = i;
543  return MB_SUCCESS;
544  }
545  }
546  *type = -1; // failure
547  return MB_FAILURE;
548 }

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

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

◆ getEntURange()

ErrorCode moab::FBEngine::getEntURange ( EntityHandle  edge,
double &  u_min,
double &  u_max 
)

Definition at line 1152 of file FBEngine.cpp.

1153 {
1154  SmoothCurve* smoothCurve = _edges[edge]; // this is a map
1155  // now, call smoothCurve methods
1156  smoothCurve->get_param_range( u_min, u_max );
1157  return MB_SUCCESS;
1158 }

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

◆ getEntUtoXYZ()

ErrorCode moab::FBEngine::getEntUtoXYZ ( EntityHandle  edge,
double  u,
double &  x,
double &  y,
double &  z 
)

Definition at line 1160 of file FBEngine.cpp.

1161 {
1162  SmoothCurve* smoothCurve = _edges[edge]; // this is a map
1163  // now, call smoothCurve methods
1164  smoothCurve->position_from_u( u, x, y, z );
1165  return MB_SUCCESS;
1166 }

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

◆ getFcEvalXYZ()

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 1075 of file FBEngine.cpp.

1091 {
1092  return MB_NOT_IMPLEMENTED; // not implemented
1093 }

References MB_NOT_IMPLEMENTED.

◆ getNumEntSets()

ErrorCode moab::FBEngine::getNumEntSets ( EntityHandle  set,
int  num_hops,
int *  all_sets 
)

Definition at line 442 of file FBEngine.cpp.

443 {
444  ErrorCode rval = MBI->num_contained_meshsets( set, all_sets, num_hops + 1 );
445  return rval;
446 }

References ErrorCode, and MBI.

◆ getNumOfType()

ErrorCode moab::FBEngine::getNumOfType ( EntityHandle  set,
int  ent_type,
int *  pNum 
)

Definition at line 498 of file FBEngine.cpp.

499 {
500  if( 0 > ent_type || 4 < ent_type )
501  {
502  std::cout << "Invalid type\n";
503  return MB_FAILURE;
504  }
505  // get sets of geom dimension tag from here, and intersect with the gentities from geo
506  // ranges
507 
508  // get the geom dimensions sets in the set (AKA gentities)
509  Range geom_sets;
510  Tag geom_tag;
511  ErrorCode rval =
513  MBERRORR( rval, "Failed to get geom tag." );
514  rval = _mbImpl->get_entities_by_type_and_tag( set, MBENTITYSET, &geom_tag, NULL, 1, geom_sets, Interface::UNION );
515  MBERRORR( rval, "Failed to get gentities from set" );
516 
517  if( ent_type == 4 )
518  {
519  *pNum = 0;
520  for( int k = 0; k <= 3; k++ )
521  {
522  Range gEntsOfTypeK = intersect( geom_sets, _my_geomTopoTool->geoRanges()[k] );
523  *pNum += (int)gEntsOfTypeK.size();
524  }
525  }
526  else
527  {
528  Range gEntsOfType = intersect( geom_sets, _my_geomTopoTool->geoRanges()[ent_type] );
529  *pNum = (int)gEntsOfType.size();
530  }
531  // we do not really check if it is in the set or not;
532  // _my_gsets[i].find(gent) != _my_gsets[i].end()
533  return MB_SUCCESS;
534 }

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.

◆ getPntRayIntsct()

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.

779 {
780  // this is pretty cool
781  // we will return only surfaces (gentities )
782  //
783  ErrorCode rval;
784 
785  unsigned int numfaces = _my_gsets[2].size();
786  // do ray fire
787  const double point[] = { x, y, z };
788  const double dir[] = { dir_x, dir_y, dir_z };
789  CartVect P( point );
790  CartVect V( dir );
791 
792  // std::vector<double> distances;
793  std::vector< EntityHandle > facets;
794  // std::vector<EntityHandle> sets;
795  unsigned int i;
796  for( i = 0; i < numfaces; i++ )
797  {
798  EntityHandle face = _my_gsets[2][i];
799  EntityHandle rootForFace;
800  rval = _my_geomTopoTool->get_root( face, rootForFace );
801  MBERRORR( rval, "Failed to get root of face." );
802  std::vector< double > distances_out;
803  std::vector< EntityHandle > sets_out;
804  std::vector< EntityHandle > facets_out;
805  rval = _my_geomTopoTool->obb_tree()->ray_intersect_sets( distances_out, sets_out, facets_out, rootForFace,
806  tolerance, point, dir );
807  unsigned int j;
808  for( j = 0; j < distances_out.size(); j++ )
809  param_coords.push_back( distances_out[j] );
810  for( j = 0; j < sets_out.size(); j++ )
811  intersect_entity_handles.push_back( sets_out[j] );
812  for( j = 0; j < facets_out.size(); j++ )
813  facets.push_back( facets_out[j] );
814 
815  MBERRORR( rval, "Failed to get ray intersections." );
816  }
817  // facets.size == distances.size()!!
818  for( i = 0; i < param_coords.size(); i++ )
819  {
820  CartVect intx = P + param_coords[i] * V;
821  for( int j = 0; j < 3; j++ )
822  intersect_coords.push_back( intx[j] );
823  }
824  if( _smooth )
825  {
826  // correct the intersection point and the distance for smooth surfaces
827  for( i = 0; i < intersect_entity_handles.size(); i++ )
828  {
829  // EntityHandle geoSet = MBH_cast(sets[i]);
830  SmoothFace* sFace = _faces[intersect_entity_handles[i]];
831  // correct coordinates and distance from point
832  /*moab::ErrorCode ray_intersection_correct(moab::EntityHandle facet, // (IN) the facet
833  where the patch is defined moab::CartVect &pt, // (IN) shoot from moab::CartVect &ray,
834  // (IN) ray direction moab::CartVect &eval_pt, // (INOUT) The intersection point double
835  & distance, // (IN OUT) the new distance bool &outside);*/
836  CartVect pos( &( intersect_coords[3 * i] ) );
837  double dist = param_coords[i];
838  bool outside = false;
839  rval = sFace->ray_intersection_correct( facets[i], P, V, pos, dist, outside );
840  MBERRORR( rval, "Failed to get better point on ray." );
841  param_coords[i] = dist;
842 
843  for( int j = 0; j < 3; j++ )
844  intersect_coords[3 * i + j] = pos[j];
845  }
846  }
847  return MB_SUCCESS;
848 }

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.

◆ getRootSet()

ErrorCode moab::FBEngine::getRootSet ( EntityHandle root_set)

Definition at line 436 of file FBEngine.cpp.

437 {
438  *root_set = _my_geomTopoTool->get_root_model_set();
439  return MB_SUCCESS;
440 }

References _my_geomTopoTool, moab::GeomTopoTool::get_root_model_set(), and MB_SUCCESS.

◆ getVtxCoord()

ErrorCode moab::FBEngine::getVtxCoord ( EntityHandle  this_gent,
double *  x0,
double *  y0,
double *  z0 
)

Definition at line 663 of file FBEngine.cpp.

664 {
665  int type;
666  ErrorCode rval = getEntType( vertex_handle, &type );
667  MBERRORR( rval, "Failed to get entity type in getVtxCoord." );
668 
669  if( type != 0 )
670  {
671  MBERRORR( MB_FAILURE, "Entity is not a vertex type." );
672  }
673 
674  Range entities;
675  rval = MBI->get_entities_by_type( vertex_handle, MBVERTEX, entities );
676  MBERRORR( rval, "can't get nodes in vertex set." );
677 
678  if( entities.size() != 1 )
679  {
680  MBERRORR( MB_FAILURE, "Vertex has multiple points." );
681  }
682  double coords[3];
683  EntityHandle node = entities[0];
684  rval = MBI->get_coords( &node, 1, coords );
685  MBERRORR( rval, "can't get coordinates." );
686  *x0 = coords[0];
687  *y0 = coords[1];
688  *z0 = coords[2];
689 
690  return MB_SUCCESS;
691 }

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

Referenced by getEntBoundBox(), and getEntClosestPt().

◆ gsubtract()

ErrorCode moab::FBEngine::gsubtract ( EntityHandle  entity_set_1,
EntityHandle  entity_set_2,
EntityHandle  result_entity_set 
)

Definition at line 693 of file FBEngine.cpp.

694 {
695  /*result_entity_set = subtract(entity_set_1, entity_set_2);*/
696  Range ents1, ents2;
697  ErrorCode rval = MBI->get_entities_by_type( entity_set_1, MBENTITYSET, ents1 );
698  MBERRORR( rval, "can't get entities from set 1." );
699 
700  rval = MBI->get_entities_by_type( entity_set_2, MBENTITYSET, ents2 );
701  MBERRORR( rval, "can't get entities from set 2." );
702 
703  ents1 = subtract( ents1, ents2 );
704  rval = MBI->clear_meshset( &result_entity_set, 1 );
705  MBERRORR( rval, "can't empty set." );
706 
707  rval = MBI->add_entities( result_entity_set, ents1 );
708  MBERRORR( rval, "can't add result to set." );
709 
710  return rval;
711 }

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

◆ Init()

ErrorCode moab::FBEngine::Init ( )

Definition at line 237 of file FBEngine.cpp.

238 {
239  if( !_initialized )
240  {
241  if( !_my_geomTopoTool ) return MB_FAILURE;
242 
244  assert( rval == MB_SUCCESS );
245  if( MB_SUCCESS != rval )
246  {
247  return rval;
248  }
249 
250  rval = split_quads();
251  assert( rval == MB_SUCCESS );
252 
254  assert( rval == MB_SUCCESS );
255 
256  if( _smooth ) rval = initializeSmoothing();
257  assert( rval == MB_SUCCESS );
258 
259  _initialized = true;
260  }
261  return MB_SUCCESS;
262 }

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

◆ initializeSmoothing()

ErrorCode moab::FBEngine::initializeSmoothing ( )
private

Definition at line 263 of file FBEngine.cpp.

264 {
265  //
266  /*ErrorCode rval = Init();
267  MBERRORR(rval, "failed initialize");*/
268  // first of all, we need to retrieve all the surfaces from the (root) set
269  // in icesheet_test we use iGeom, but maybe that is a stretch
270  // get directly the sets with geom dim 2, and from there create the SmoothFace
271  Tag geom_tag;
272  ErrorCode rval = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag );
273  MBERRORR( rval, "can't get geom tag" );
274 
275  int numSurfaces = _my_gsets[2].size();
276  // SmoothFace ** smthFace = new SmoothFace *[numSurfaces];
277  _smthFace = new SmoothFace*[numSurfaces];
278 
279  // there should also be a map from surfaces to evaluators
280  // std::map<MBEntityHandle, SmoothFace*> mapSurfaces;
281 
282  int i = 0;
283  Range::iterator it;
284  for( it = _my_gsets[2].begin(); it != _my_gsets[2].end(); ++it, i++ )
285  {
286  EntityHandle face = *it;
287  _smthFace[i] = new SmoothFace( MBI, face, _my_geomTopoTool ); // geom topo tool will be used for searching,
288  // among other things; also for senses in edge sets...
289  _faces[face] = _smthFace[i];
290  }
291 
292  int numCurves = _my_gsets[1].size(); // csets.size();
293  // SmoothCurve ** smthCurve = new SmoothCurve *[numCurves];
294  _smthCurve = new SmoothCurve*[numCurves];
295  // there should also be a map from surfaces to evaluators
296  // std::map<MBEntityHandle, SmoothCurve*> mapCurves;
297 
298  i = 0;
299  for( it = _my_gsets[1].begin(); it != _my_gsets[1].end(); ++it, i++ )
300  {
301  EntityHandle curve = *it;
302  _smthCurve[i] = new SmoothCurve( MBI, curve, _my_geomTopoTool );
303  _edges[curve] = _smthCurve[i];
304  }
305 
306  for( i = 0; i < numSurfaces; i++ )
307  {
308  _smthFace[i]->init_gradient(); // this will also retrieve the triangles in each surface
309  _smthFace[i]->compute_tangents_for_each_edge(); // this one will consider all edges
310  // internal, so the
311  // tangents are all in the direction of the edge; a little bit of waste, as we store
312  // one tangent for each edge node , even though they are equal here...
313  // no loops are considered
314  }
315 
316  // this will be used to mark boundary edges, so for them the control points are computed earlier
317  unsigned char value = 0; // default value is "not used"=0 for the tag
318  // unsigned char def_data_bit = 1;// valid by default
319  // rval = mb->tag_create("valid", 1, MB_TAG_BIT, validTag, &def_data_bit);
320  Tag markTag;
321  rval = MBI->tag_get_handle( "MARKER", 1, MB_TYPE_BIT, markTag, MB_TAG_EXCL | MB_TAG_BIT,
322  &value ); // default value : 0 = not computed yet
323  // each feature edge will need to have a way to retrieve at every moment the surfaces it belongs
324  // to from edge sets, using the sense tag, we can get faces, and from each face, using the map,
325  // we can get the SmoothFace (surface evaluator), that has everything, including the normals!!!
326  assert( rval == MB_SUCCESS );
327 
328  // create the tag also for control points on the edges
329  double defCtrlPoints[9] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. };
330  Tag edgeCtrlTag;
331  rval = MBI->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, edgeCtrlTag, MB_TAG_DENSE | MB_TAG_CREAT,
332  &defCtrlPoints );
333  assert( rval == MB_SUCCESS );
334 
335  Tag facetCtrlTag;
336  double defControls[18] = { 0. };
337  rval = MBI->tag_get_handle( "CONTROLFACE", 18, MB_TYPE_DOUBLE, facetCtrlTag, MB_TAG_CREAT | MB_TAG_DENSE,
338  &defControls );
339  assert( rval == MB_SUCCESS );
340 
341  Tag facetEdgeCtrlTag;
342  double defControls2[27] = { 0. }; // corresponding to 9 control points on edges, in order
343  // from edge 0, 1, 2 ( 1-2, 2-0, 0-1 )
344  rval = MBI->tag_get_handle( "CONTROLEDGEFACE", 27, MB_TYPE_DOUBLE, facetEdgeCtrlTag, MB_TAG_CREAT | MB_TAG_DENSE,
345  &defControls2 );
346  assert( rval == MB_SUCCESS );
347  // if the
348  double min_dot = -1.0; // depends on _angle, but now we ignore it, for the time being
349  for( i = 0; i < numCurves; i++ )
350  {
351  _smthCurve[i]->compute_tangents_for_each_edge(); // do we need surfaces now? or just the chains?
352  // the computed edges will be marked accordingly; later one, only internal edges to surfaces
353  // are left
354  _smthCurve[i]->compute_control_points_on_boundary_edges( min_dot, _faces, edgeCtrlTag, markTag );
355  }
356 
357  // when done with boundary edges, compute the control points on all edges in the surfaces
358 
359  for( i = 0; i < numSurfaces; i++ )
360  {
361  // also pass the tags for
362  _smthFace[i]->compute_control_points_on_edges( min_dot, edgeCtrlTag, markTag );
363  }
364 
365  // now we should be able to compute the control points for the facets
366 
367  for( i = 0; i < numSurfaces; i++ )
368  {
369  // also pass the tags for edge and facet control points
370  _smthFace[i]->compute_internal_control_points_on_facets( min_dot, facetCtrlTag, facetEdgeCtrlTag );
371  }
372  // we will need to compute the tangents for all edges in the model
373  // they will be needed for control points for each edge
374  // the boundary edges and the feature edges are more complicated
375  // the boundary edges need to consider their loops, but feature edges need to consider loops and
376  // the normals on each connected surface
377 
378  // some control points
379  if( Debug_surf_eval )
380  for( i = 0; i < numSurfaces; i++ )
381  _smthFace[i]->DumpModelControlPoints();
382 
383  return MB_SUCCESS;
384 }

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::SmoothCurve::compute_tangents_for_each_edge(), moab::SmoothFace::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().

◆ isEntAdj()

ErrorCode moab::FBEngine::isEntAdj ( EntityHandle  entity1,
EntityHandle  entity2,
bool &  adjacent_out 
)

Definition at line 1180 of file FBEngine.cpp.

1181 {
1182  int type1, type2;
1183  ErrorCode rval = getEntType( entity1, &type1 );
1184  if( MB_SUCCESS != rval ) return rval;
1185  rval = getEntType( entity2, &type2 );
1186  if( MB_SUCCESS != rval ) return rval;
1187 
1188  Range adjs;
1189  if( type1 < type2 )
1190  {
1191  rval = MBI->get_parent_meshsets( entity1, adjs, type2 - type1 );
1192  if( MB_SUCCESS != rval ) return rval; // MBERRORR("Failed to get parent meshsets in iGeom_isEntAdj.");
1193  }
1194  else
1195  {
1196  // note: if they ave the same type, they will not be adjacent, in our definition
1197  rval = MBI->get_child_meshsets( entity1, adjs, type1 - type2 );
1198  if( MB_SUCCESS != rval ) return rval; // MBERRORR("Failed to get child meshsets in iGeom_isEntAdj.");
1199  }
1200 
1201  // adjacent_out = adjs.find(entity2) != _my_gsets[type2].end();
1202  // hmmm, possible bug here; is this called?
1203  adjacent_out = adjs.find( entity2 ) != adjs.end();
1204 
1205  return MB_SUCCESS;
1206 }

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

◆ measure()

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

Definition at line 994 of file FBEngine.cpp.

995 {
996  ErrorCode rval;
997  for( int i = 0; i < entities_size; i++ )
998  {
999  measures[i] = 0.;
1000 
1001  int type;
1002  EntityHandle gset = moab_entities[i];
1003  rval = getEntType( gset, &type );
1004  if( MB_SUCCESS != rval ) return rval;
1005  if( type == 1 )
1006  { // edge: get all edges part of the edge set
1007  Range entities;
1008  rval = MBI->get_entities_by_type( gset, MBEDGE, entities );
1009  if( MB_SUCCESS != rval ) return rval;
1010 
1011  for( Range::iterator it = entities.begin(); it != entities.end(); ++it )
1012  {
1013  EntityHandle edge = *it;
1014  CartVect vv[2];
1015  const EntityHandle* conn2 = NULL;
1016  int num_nodes;
1017  rval = MBI->get_connectivity( edge, conn2, num_nodes );
1018  if( MB_SUCCESS != rval || num_nodes != 2 ) return MB_FAILURE;
1019  rval = MBI->get_coords( conn2, 2, (double*)&( vv[0][0] ) );
1020  if( MB_SUCCESS != rval ) return rval;
1021 
1022  vv[0] = vv[1] - vv[0];
1023  measures[i] += vv[0].length();
1024  }
1025  }
1026  if( type == 2 )
1027  { // surface
1028  // get triangles in surface; TODO: quads!
1029  Range entities;
1030  rval = MBI->get_entities_by_type( gset, MBTRI, entities );
1031  if( MB_SUCCESS != rval ) return rval;
1032 
1033  for( Range::iterator it = entities.begin(); it != entities.end(); ++it )
1034  {
1035  EntityHandle tri = *it;
1036  CartVect vv[3];
1037  const EntityHandle* conn3 = NULL;
1038  int num_nodes;
1039  rval = MBI->get_connectivity( tri, conn3, num_nodes );
1040  if( MB_SUCCESS != rval || num_nodes != 3 ) return MB_FAILURE;
1041  rval = MBI->get_coords( conn3, 3, (double*)&( vv[0][0] ) );
1042  if( MB_SUCCESS != rval ) return rval;
1043 
1044  vv[1] = vv[1] - vv[0];
1045  vv[2] = vv[2] - vv[0];
1046  vv[0] = vv[1] * vv[2];
1047  measures[i] += vv[0].length() / 2; // area of triangle
1048  }
1049  }
1050  }
1051  return MB_SUCCESS;
1052 }

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

◆ moab_instance()

Interface* moab::FBEngine::moab_instance ( )
inline

Definition at line 91 of file FBEngine.hpp.

92  {
93  return _mbImpl;
94  }

References _mbImpl.

◆ print_debug_triangle()

void moab::FBEngine::print_debug_triangle ( EntityHandle  triangle)
private

Definition at line 1963 of file FBEngine.cpp.

1964 {
1965  std::cout << " triangle id:" << _mbImpl->id_from_handle( t ) << "\n";
1966  const EntityHandle* conn3 = NULL;
1967  int nnodes = 0;
1968  _mbImpl->get_connectivity( t, conn3, nnodes );
1969  // get coords
1970  CartVect P[3];
1971  _mbImpl->get_coords( conn3, 3, (double*)&P[0] );
1972  std::cout << " nodes:" << conn3[0] << " " << conn3[1] << " " << conn3[2] << "\n";
1973  CartVect PP[3];
1974  PP[0] = P[1] - P[0];
1975  PP[1] = P[2] - P[1];
1976  PP[2] = P[0] - P[2];
1977 
1978  std::cout << " pos:" << P[0] << " " << P[1] << " " << P[2] << "\n";
1979  std::cout << " x,y diffs " << PP[0][0] << " " << PP[0][1] << ", " << PP[1][0] << " " << PP[1][1] << ", "
1980  << PP[2][0] << " " << PP[2][1] << "\n";
1981  return;
1982 }

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().

◆ redistribute_boundary_edges_to_faces()

ErrorCode moab::FBEngine::redistribute_boundary_edges_to_faces ( EntityHandle  face,
EntityHandle  newFace,
std::vector< EntityHandle > &  chainedEdges 
)
private

Definition at line 2872 of file FBEngine.cpp.

2875 {
2876 
2877  // so far, original boundary edges are all parent/child relations for face
2878  // we should get them all, and see if they are truly adjacent to face or newFace
2879  // decide also on the orientation/sense with respect to the triangles
2880  Range r1; // range in old face
2881  Range r2; // range of tris in new face
2882  ErrorCode rval = _mbImpl->get_entities_by_dimension( face, 2, r1 );
2883  MBERRORR( rval, " can't get triangles from old face" );
2884  rval = _mbImpl->get_entities_by_dimension( newFace, 2, r2 );
2885  MBERRORR( rval, " can't get triangles from new face" );
2886  // right now, all boundary edges are children of face
2887  // we need to get them all, and verify if they indeed are adjacent to face
2888  Range children;
2889  rval = _mbImpl->get_child_meshsets( face, children ); // only direct children are of interest
2890  MBERRORR( rval, " can't get children sets from face" );
2891 
2892  for( Range::iterator it = children.begin(); it != children.end(); ++it )
2893  {
2894  EntityHandle edge = *it;
2895  if( std::find( chainedEdges.begin(), chainedEdges.end(), edge ) != chainedEdges.end() )
2896  continue; // we already set this one fine
2897  // get one mesh edge from the edge; we have to get all of them!!
2898  if( _my_geomTopoTool->dimension( edge ) != 1 ) continue; // not an edge
2899  Range mesh_edges;
2900  rval = _mbImpl->get_entities_by_handle( edge, mesh_edges );
2901  MBERRORR( rval, " can't get mesh edges from edge" );
2902  if( mesh_edges.empty() ) MBERRORR( MB_FAILURE, " no mesh edges" );
2903  EntityHandle mesh_edge = mesh_edges[0]; // first one is enough
2904  // get its triangles; see which one is in r1 or r2; it should not be in both
2905  Range triangles;
2906  rval = _mbImpl->get_adjacencies( &mesh_edge, 1, 2, false, triangles );
2907  MBERRORR( rval, " can't get adjacent triangles" );
2908  Range intx1 = intersect( triangles, r1 );
2909  Range intx2 = intersect( triangles, r2 );
2910  if( !intx1.empty() && !intx2.empty() ) MBERRORR( MB_FAILURE, " at least one should be empty" );
2911 
2912  if( intx2.empty() )
2913  {
2914  // it means it is in the range r1; the sense should be fine, no need to reset
2915  // the sense should have been fine, also
2916  continue;
2917  }
2918  // so it must be a triangle in r2;
2919  EntityHandle triangle = intx2[0]; // one triangle only
2920  // remove the edge from boundary of face, and add it to the boundary of newFace
2921  // remove_parent_child( EntityHandle parent, EntityHandle child )
2922  rval = _mbImpl->remove_parent_child( face, edge );
2923  MBERRORR( rval, " can't remove parent child relation for edge" );
2924  // add to the new face
2925  rval = _mbImpl->add_parent_child( newFace, edge );
2926  MBERRORR( rval, " can't add parent child relation for edge" );
2927 
2928  // set some sense, based on the sense of the mesh_edge in triangle
2929  int num1, sense, offset;
2930  rval = _mbImpl->side_number( triangle, mesh_edge, num1, sense, offset );
2931  MBERRORR( rval, "mesh edge not adjacent to triangle" );
2932 
2933  rval = this->_my_geomTopoTool->set_sense( edge, newFace, sense );
2934  MBERRORR( rval, "can't set proper sense of edge in face" );
2935  }
2936 
2937  return MB_SUCCESS;
2938 }

References _mbImpl, _my_geomTopoTool, moab::Interface::add_parent_child(), moab::Range::begin(), 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().

◆ separate()

ErrorCode moab::FBEngine::separate ( EntityHandle  face,
std::vector< EntityHandle > &  chainedEdges,
Range first,
Range second 
)
private

Definition at line 1676 of file FBEngine.cpp.

1680 {
1681  // Range unaffectedTriangles = subtract(iniTriangles, _piercedTriangles);
1682  // insert in each
1683  // start with a new triangle, and flood to get the first range; what is left is the
1684  // second range
1685  // flood fill is considering edges adjacent to seed triangles; if there is
1686  // an edge in the new_geo_edge, it is skipped; triangles in the
1687  // triangles to delete are not added
1688  // first, create all edges of the new triangles
1689 
1690  //
1691  // new face will have the new edge oriented positively
1692  // get mesh edges from geo edge (splitting gedge);
1693 
1694  Range mesh_edges;
1695  ErrorCode rval;
1696  // mesh_edges
1697  for( unsigned int j = 0; j < chainedEdges.size(); j++ )
1698  {
1699  // this will keep adding edges to the mesh_edges range
1700  // when we get out, the mesh_edges will be in this range, but not ordered
1701  rval = _mbImpl->get_entities_by_type( chainedEdges[j], MBEDGE, mesh_edges );
1702  MBERRORR( rval, "can't get new polyline edges" );
1703  if( debug_splits )
1704  {
1705  std::cout << " At chained edge " << j << " " << _mbImpl->id_from_handle( chainedEdges[j] )
1706  << " mesh_edges Range size:" << mesh_edges.size() << "\n";
1707  }
1708  }
1709 
1710  // get a positive triangle adjacent to mesh_edge[0]
1711  // add to first triangles to the left, second triangles to the right of the mesh_edges ;
1712 
1713  // create a temp tag, and when done, delete it
1714  // default value: 0
1715  // 3 to be deleted, pierced
1716  // 1 first set
1717  // 2 second set
1718  // create the tag also for control points on the edges
1719  int defVal = 0;
1720  Tag separateTag;
1721  rval = MBI->tag_get_handle( "SEPARATE_TAG", 1, MB_TYPE_INTEGER, separateTag, MB_TAG_DENSE | MB_TAG_CREAT, &defVal );
1722  MBERRORR( rval, "can't create temp tag for separation" );
1723  // the deleted triangles will get a value 3, from start
1724  int delVal = 3;
1725  for( Range::iterator it1 = this->_piercedTriangles.begin(); it1 != _piercedTriangles.end(); ++it1 )
1726  {
1727  EntityHandle trToDelete = *it1;
1728  rval = _mbImpl->tag_set_data( separateTag, &trToDelete, 1, &delVal );
1729  MBERRORR( rval, "can't set delete tag value" );
1730  }
1731 
1732  // find a triangle that will be in the first range, positively oriented about the splitting edge
1733  EntityHandle seed1 = 0;
1734  for( Range::iterator it = mesh_edges.begin(); it != mesh_edges.end() && !seed1; ++it )
1735  {
1736  EntityHandle meshEdge = *it;
1737  Range adj_tri;
1738  rval = _mbImpl->get_adjacencies( &meshEdge, 1, 2, false, adj_tri );
1739  MBERRORR( rval, "can't get adj_tris to mesh edge" );
1740 
1741  for( Range::iterator it2 = adj_tri.begin(); it2 != adj_tri.end(); ++it2 )
1742  {
1743  EntityHandle tr = *it2;
1744  if( _piercedTriangles.find( tr ) != _piercedTriangles.end() )
1745  continue; // do not attach pierced triangles, they are not good
1746  int num1, sense, offset;
1747  rval = _mbImpl->side_number( tr, meshEdge, num1, sense, offset );
1748  MBERRORR( rval, "edge not adjacent" );
1749  if( sense == 1 )
1750  {
1751  // firstSet.insert(tr);
1752  if( !seed1 )
1753  {
1754  seed1 = tr;
1755  break;
1756  }
1757  }
1758  }
1759  }
1760 
1761  // flood fill first set, the rest will be in second set
1762  // the edges from new_geo_edge will not be crossed
1763 
1764  // get edges of face (adjacencies)
1765  // also get the old boundary edges, from face; they will be edges to not cross, too
1766  Range bound_edges;
1767  rval = getAdjacentEntities( face, 1, bound_edges );
1768  MBERRORR( rval, "can't get boundary edges" );
1769 
1770  // add to the do not cross edges range, all edges from initial boundary
1771  Range initialBoundaryEdges;
1772  for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
1773  {
1774  EntityHandle bound_edge = *it;
1775  rval = _mbImpl->get_entities_by_dimension( bound_edge, 1, initialBoundaryEdges );
1776  }
1777 
1778  Range doNotCrossEdges = unite( initialBoundaryEdges, mesh_edges ); // add the splitting edges !
1779 
1780  // use a second method, with tags
1781  //
1782  std::queue< EntityHandle > queue1;
1783  queue1.push( seed1 );
1784  std::vector< EntityHandle > arr1;
1785  while( !queue1.empty() )
1786  {
1787  // start copy
1788  EntityHandle currentTriangle = queue1.front();
1789  queue1.pop();
1790  arr1.push_back( currentTriangle );
1791  // add new triangles that share an edge
1792  Range currentEdges;
1793  rval = _mbImpl->get_adjacencies( &currentTriangle, 1, 1, true, currentEdges, Interface::UNION );
1794  MBERRORR( rval, "can't get adjacencies" );
1795  for( Range::iterator it = currentEdges.begin(); it != currentEdges.end(); ++it )
1796  {
1797  EntityHandle frontEdge = *it;
1798  if( doNotCrossEdges.find( frontEdge ) == doNotCrossEdges.end() )
1799  {
1800  // this is an edge that can be crossed
1801  Range adj_tri;
1802  rval = _mbImpl->get_adjacencies( &frontEdge, 1, 2, false, adj_tri, Interface::UNION );
1803  MBERRORR( rval, "can't get adj_tris" );
1804  // if the triangle is not in first range, add it to the queue
1805  for( Range::iterator it2 = adj_tri.begin(); it2 != adj_tri.end(); ++it2 )
1806  {
1807  EntityHandle tri2 = *it2;
1808  int val = 0;
1809  rval = _mbImpl->tag_get_data( separateTag, &tri2, 1, &val );
1810  MBERRORR( rval, "can't get tag value" );
1811  if( val ) continue;
1812  // else, set it to 1
1813  val = 1;
1814  rval = _mbImpl->tag_set_data( separateTag, &tri2, 1, &val );
1815  MBERRORR( rval, "can't get tag value" );
1816 
1817  queue1.push( tri2 );
1818  }
1819  } // end edge do not cross
1820  } // end while
1821  }
1822 
1823  std::sort( arr1.begin(), arr1.end() );
1824  // Range first1;
1825  std::copy( arr1.rbegin(), arr1.rend(), range_inserter( first ) );
1826 
1827  // std::cout<< "\n first1.size() " << first1.size() << " first.size(): " << first.size() <<
1828  // "\n";
1829  if( debug_splits )
1830  {
1831  EntityHandle tmpSet;
1832  _mbImpl->create_meshset( MESHSET_SET, tmpSet );
1833  _mbImpl->add_entities( tmpSet, first );
1834  _mbImpl->write_file( "dbg1.vtk", "vtk", 0, &tmpSet, 1 );
1835  }
1836  // now, decide the set 2:
1837  // first, get all ini tris
1838  Range initr;
1839  rval = _mbImpl->get_entities_by_type( face, MBTRI, initr );
1840  MBERRORR( rval, "can't get tris " );
1841  second = unite( initr, _newTriangles );
1842  Range second2 = subtract( second, _piercedTriangles );
1843  second = subtract( second2, first );
1844  _newTriangles.clear();
1845  if( debug_splits )
1846  {
1847  std::cout << "\n second.size() " << second.size() << " first.size(): " << first.size() << "\n";
1848  // debugging code
1849  EntityHandle tmpSet2;
1850  _mbImpl->create_meshset( MESHSET_SET, tmpSet2 );
1851  _mbImpl->add_entities( tmpSet2, second );
1852  _mbImpl->write_file( "dbg2.vtk", "vtk", 0, &tmpSet2, 1 );
1853  }
1854  /*Range intex = intersect(first, second);
1855  if (!intex.empty() && debug_splits)
1856  {
1857  std::cout << "error, the sets should be disjoint\n";
1858  for (Range::iterator it1=intex.begin(); it1!=intex.end(); ++it1)
1859  {
1860  std::cout<<_mbImpl->id_from_handle(*it1) << "\n";
1861  }
1862  }*/
1863  rval = _mbImpl->tag_delete( separateTag );
1864  MBERRORR( rval, "can't delete tag " );
1865  return MB_SUCCESS;
1866 }

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::GeomUtil::first(), 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().

◆ set_default_neumann_tags()

ErrorCode moab::FBEngine::set_default_neumann_tags ( )
private

Definition at line 3598 of file FBEngine.cpp.

3599 {
3600  // these are for debugging purposes only
3601  // check the initial tag, then
3602  Tag ntag;
3604  MBERRORR( rval, "can't get tag handle" );
3605  // get all surfaces in the model now
3606  Range sets[5];
3607  rval = _my_geomTopoTool->find_geomsets( sets );
3608  MBERRORR( rval, "can't get geo sets" );
3609  int nfaces = (int)sets[2].size();
3610  int* vals = new int[nfaces];
3611  for( int i = 0; i < nfaces; i++ )
3612  {
3613  vals[i] = i + 1;
3614  }
3615  rval = _mbImpl->tag_set_data( ntag, sets[2], (void*)vals );
3616  MBERRORR( rval, "can't set tag values for neumann sets" );
3617 
3618  delete[] vals;
3619 
3620  return MB_SUCCESS;
3621 }

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

Referenced by create_volume_with_direction().

◆ set_neumann_tags()

ErrorCode moab::FBEngine::set_neumann_tags ( EntityHandle  face,
EntityHandle  newFace 
)
private

Definition at line 2940 of file FBEngine.cpp.

2941 {
2942  // these are for debugging purposes only
2943  // check the initial tag, then
2944  Tag ntag;
2946  MBERRORR( rval, "can't get tag handle" );
2947  // check the value for face
2948  int nval;
2949  rval = _mbImpl->tag_get_data( ntag, &face, 1, &nval );
2950  if( MB_SUCCESS == rval )
2951  {
2952  nval++;
2953  }
2954  else
2955  {
2956  nval = 1;
2957  rval = _mbImpl->tag_set_data( ntag, &face, 1, &nval );
2958  MBERRORR( rval, "can't set tag" );
2959  nval = 2;
2960  }
2961  rval = _mbImpl->tag_set_data( ntag, &newFace, 1, &nval );
2962  MBERRORR( rval, "can't set tag" );
2963 
2964  return MB_SUCCESS;
2965 }

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().

◆ set_smooth()

void moab::FBEngine::set_smooth ( )
inline

Definition at line 218 of file FBEngine.hpp.

219  {
220  _smooth = true;
221  }

References _smooth.

◆ setArrData()

ErrorCode moab::FBEngine::setArrData ( const EntityHandle entity_handles,
int  entity_handles_size,
Tag  tag_handle,
const void *  tag_values 
)

Definition at line 958 of file FBEngine.cpp.

962 {
963  // responsibility of the user to have tag_values_out properly allocated
964  // only some types of Tags are possible (double, int, etc)
965  return MBI->tag_set_data( tag_handle, entity_handles, entity_handles_size, tag_values );
966 }

References MBI.

◆ smooth_new_intx_points()

ErrorCode moab::FBEngine::smooth_new_intx_points ( EntityHandle  face,
std::vector< EntityHandle > &  chainedEdges 
)
private

Definition at line 1622 of file FBEngine.cpp.

1623 {
1624 
1625  // do not move nodes from the original face
1626  // first get all triangles, and then all nodes from those triangles
1627 
1628  Range tris;
1629  ErrorCode rval = _mbImpl->get_entities_by_type( face, MBTRI, tris );
1630  MBERRORR( rval, "can't get triangles" );
1631 
1632  Range ini_nodes;
1633  rval = _mbImpl->get_connectivity( tris, ini_nodes );
1634  MBERRORR( rval, "can't get connectivities" );
1635 
1636  SmoothFace* smthFace = _faces[face];
1637 
1638  // get all nodes from chained edges
1639  Range mesh_edges;
1640  for( unsigned int j = 0; j < chainedEdges.size(); j++ )
1641  {
1642  // keep adding to the range of mesh edges
1643  rval = _mbImpl->get_entities_by_dimension( chainedEdges[j], 1, mesh_edges );
1644  MBERRORR( rval, "can't get mesh edges" );
1645  }
1646  // nodes along polyline
1647  Range nodes_on_polyline;
1648  rval = _mbImpl->get_connectivity( mesh_edges, nodes_on_polyline, true ); // corners only
1649  MBERRORR( rval, "can't get nodes on the polyline" );
1650 
1651  Range new_intx_nodes = subtract( nodes_on_polyline, ini_nodes );
1652 
1653  std::vector< double > ini_coords;
1654  int num_points = (int)new_intx_nodes.size();
1655  ini_coords.resize( 3 * num_points );
1656  rval = _mbImpl->get_coords( new_intx_nodes, &( ini_coords[0] ) );
1657  MBERRORR( rval, "can't get coordinates" );
1658 
1659  int i = 0;
1660  for( Range::iterator it = new_intx_nodes.begin(); it != new_intx_nodes.end(); ++it )
1661  {
1662  /*EntityHandle node = *it;*/
1663  int i3 = 3 * i;
1664  smthFace->move_to_surface( ini_coords[i3], ini_coords[i3 + 1], ini_coords[i3 + 2] );
1665  // reset the coordinates of this node
1666  ++i;
1667  }
1668  rval = _mbImpl->set_coords( new_intx_nodes, &( ini_coords[0] ) );
1669  MBERRORR( rval, "can't set smoothed coordinates for the new nodes" );
1670 
1671  return MB_SUCCESS;
1672 }

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().

◆ split_bedge_at_new_mesh_node()

ErrorCode moab::FBEngine::split_bedge_at_new_mesh_node ( EntityHandle  b_edge,
EntityHandle  atNode,
EntityHandle  brokenEdge,
EntityHandle new_edge 
)

Definition at line 2561 of file FBEngine.cpp.

2565 {
2566  // the node should be in the list of nodes
2567 
2568  int dim = _my_geomTopoTool->dimension( edge );
2569  if( dim != 1 ) return MB_FAILURE; // not an edge
2570 
2571  if( debug_splits )
2572  {
2573  std::cout << "Split edge " << _mbImpl->id_from_handle( edge )
2574  << " with global id: " << _my_geomTopoTool->global_id( edge )
2575  << " at new node:" << _mbImpl->id_from_handle( node ) << "breaking mesh edge"
2576  << _mbImpl->id_from_handle( brokenEdge ) << "\n";
2577  }
2578 
2579  // now get the edges
2580  // the order is important...
2581  // these are ordered sets !!
2582  std::vector< EntityHandle > ents;
2583  ErrorCode rval = _mbImpl->get_entities_by_type( edge, MBEDGE, ents );
2584  if( MB_SUCCESS != rval ) return rval;
2585  if( ents.size() < 1 ) return MB_FAILURE; // no edges
2586 
2587  const EntityHandle* conn = NULL;
2588  int len;
2589  // the edge connected to the splitting node is brokenEdge
2590  // find the small edges it is broken into, which are connected to
2591  // new vertex (node) and its ends; also, if necessary, reorient these small edges
2592  // for proper orientation
2593  rval = _mbImpl->get_connectivity( brokenEdge, conn, len );
2594  MBERRORR( rval, "Failed to get connectivity of broken edge" );
2595 
2596  // we already created the new edges, make sure they are oriented fine; if not, revert them
2597  EntityHandle conn02[] = { conn[0], node }; // first node and new node
2598  // default, intersect
2599  std::vector< EntityHandle > adj_edges;
2600  rval = _mbImpl->get_adjacencies( conn02, 2, 1, false, adj_edges );
2601  if( adj_edges.size() < 1 || rval != MB_SUCCESS ) MBERRORR( MB_FAILURE, " Can't find small split edge" );
2602 
2603  // get this edge connectivity;
2604  EntityHandle firstEdge = adj_edges[0];
2605  const EntityHandle* connActual = NULL;
2606  rval = _mbImpl->get_connectivity( firstEdge, connActual, len );
2607  MBERRORR( rval, "Failed to get connectivity of first split edge" );
2608  // if it is the same as conn02, we are happy
2609  if( conn02[0] != connActual[0] )
2610  {
2611  // reset connectivity of edge
2612  rval = _mbImpl->set_connectivity( firstEdge, conn02, 2 );
2613  MBERRORR( rval, "Failed to reset connectivity of first split edge" );
2614  }
2615  // now treat the second edge
2616  adj_edges.clear(); //
2617  EntityHandle conn21[] = { node, conn[1] }; // new node and second node
2618  rval = _mbImpl->get_adjacencies( conn21, 2, 1, false, adj_edges );
2619  if( adj_edges.size() < 1 || rval != MB_SUCCESS ) MBERRORR( MB_FAILURE, " Can't find second small split edge" );
2620 
2621  // get this edge connectivity;
2622  EntityHandle secondEdge = adj_edges[0];
2623  rval = _mbImpl->get_connectivity( firstEdge, connActual, len );
2624  MBERRORR( rval, "Failed to get connectivity of first split edge" );
2625  // if it is the same as conn21, we are happy
2626  if( conn21[0] != connActual[0] )
2627  {
2628  // reset connectivity of edge
2629  rval = _mbImpl->set_connectivity( secondEdge, conn21, 2 );
2630  MBERRORR( rval, "Failed to reset connectivity of first split edge" );
2631  }
2632 
2633  int num_mesh_edges = (int)ents.size();
2634  int index_edge; // this is the index of the edge that will be removed (because it is split)
2635  // the rest of edges will be put in order in the (remaining) edge and new edge
2636 
2637  for( index_edge = 0; index_edge < num_mesh_edges; index_edge++ )
2638  if( brokenEdge == ents[index_edge] ) break;
2639  if( index_edge >= num_mesh_edges ) MBERRORR( MB_FAILURE, "can't find the broken edge" );
2640 
2641  // so the edges up to index_edge and firstEdge, will form the "old" edge
2642  // the edges secondEdge and from index_edge+1 to end will form the new_edge
2643 
2644  // here, we have 0 ... index_edge edges in the first set,
2645  // create a vertex set and add the node to it
2646 
2647  if( debug_splits )
2648  {
2649  std::cout << "Split edge with " << num_mesh_edges << " mesh edges, at index (0 based) " << index_edge << "\n";
2650  }
2651 
2652  // at this moment, the node set should have been already created in new_geo_edge
2653  EntityHandle nodeSet; // the node set that has the node (find it!)
2654  bool found = find_vertex_set_for_node( node, nodeSet );
2655 
2656  if( !found )
2657  {
2658  // create a node set and add the node
2659 
2660  // must be an error, but create one nevertheless
2661  rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet );
2662  MBERRORR( rval, "Failed to create a new vertex set" );
2663 
2664  rval = _mbImpl->add_entities( nodeSet, &node, 1 );
2665  MBERRORR( rval, "Failed to add the node to the set" );
2666 
2667  rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 ); //
2668  MBERRORR( rval, "Failed to commit the node set" );
2669 
2670  if( debug_splits )
2671  {
2672  std::cout << " create a vertex set (this should have been created before!)"
2673  << _mbImpl->id_from_handle( nodeSet )
2674  << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << "\n";
2675  }
2676  }
2677 
2678  // we need to remove the remaining mesh edges from first set, and add it
2679  // to the second set, in order
2680 
2681  rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_edge );
2682  MBERRORR( rval, "can't create geo edge" );
2683 
2684  int remaining = num_mesh_edges - 1 - index_edge;
2685 
2686  // add edges to the new edge set
2687  rval = _mbImpl->add_entities( new_edge, &secondEdge, 1 ); // add the second split edge to new edge
2688  MBERRORR( rval, "can't add second split edge to the new edge" );
2689  // then add the rest
2690  rval = _mbImpl->add_entities( new_edge, &ents[index_edge + 1], remaining );
2691  MBERRORR( rval, "can't add edges to the new edge" );
2692 
2693  // also, remove the second node set from old edge
2694  // remove the edges index_edge and up
2695 
2696  rval = _mbImpl->remove_entities( edge, &ents[index_edge], remaining + 1 ); // include the
2697  MBERRORR( rval, "can't remove edges from the old edge" );
2698  // add the firstEdge too
2699  rval = _mbImpl->add_entities( edge, &firstEdge, 1 ); // add the second split edge to new edge
2700  MBERRORR( rval, "can't add first split edge to the old edge" );
2701 
2702  // need to find the adjacent vertex sets
2703  Range vertexRange;
2704  rval = getAdjacentEntities( edge, 0, vertexRange );
2705 
2706  EntityHandle secondSet;
2707  if( vertexRange.size() == 1 )
2708  {
2709  // initially a periodic edge, OK to add the new set to both edges, and the
2710  // second set
2711  secondSet = vertexRange[0];
2712  }
2713  else
2714  {
2715  if( vertexRange.size() > 2 ) return MB_FAILURE; // something must be wrong with too many vertices
2716  // find first node
2717  EntityHandle firstNode;
2718 
2719  rval = MBI->get_connectivity( ents[0], conn, len );
2720  if( MB_SUCCESS != rval ) return rval;
2721  firstNode = conn[0]; // this is the first node of the initial edge
2722  // we will use it to identify the vertex set associated with it
2723  int k;
2724  for( k = 0; k < 2; k++ )
2725  {
2726  Range verts;
2727  rval = _mbImpl->get_entities_by_type( vertexRange[k], MBVERTEX, verts );
2728 
2729  MBERRORR( rval, "can't get vertices from vertex set" );
2730  if( verts.size() != 1 ) MBERRORR( MB_FAILURE, " node set not defined well" );
2731  if( firstNode == verts[0] )
2732  {
2733  secondSet = vertexRange[1 - k]; // the other set; it is 1 or 0
2734  break;
2735  }
2736  }
2737  if( k >= 2 )
2738  {
2739  // it means we didn't find the right node set
2740  MBERRORR( MB_FAILURE, " can't find the right vertex" );
2741  }
2742  // remove the second set from the connectivity with the
2743  // edge (parent-child relation)
2744  // remove_parent_child( EntityHandle parent,
2745  // EntityHandle child )
2746  rval = _mbImpl->remove_parent_child( edge, secondSet );
2747  MBERRORR( rval, " can't remove the second vertex from edge" );
2748  }
2749  // at this point, we just need to add to both edges the new node set (vertex)
2750  rval = _mbImpl->add_parent_child( edge, nodeSet );
2751  MBERRORR( rval, " can't add new vertex to old edge" );
2752 
2753  rval = _mbImpl->add_parent_child( new_edge, nodeSet );
2754  MBERRORR( rval, " can't add new vertex to new edge" );
2755 
2756  // one thing that I forgot: add the secondSet as a child to new edge!!!
2757  // (so far, the new edge has only one end vertex!)
2758  rval = _mbImpl->add_parent_child( new_edge, secondSet );
2759  MBERRORR( rval, " can't add second vertex to new edge" );
2760 
2761  // now, add the edge and vertex to geo tool
2762 
2763  rval = _my_geomTopoTool->add_geo_set( new_edge, 1 );
2764  MBERRORR( rval, " can't add new edge" );
2765 
2766  // next, get the adjacent faces to initial edge, and add them as parents
2767  // to the new edge
2768 
2769  // need to find the adjacent face sets
2770  Range faceRange;
2771  rval = getAdjacentEntities( edge, 2, faceRange );
2772 
2773  // these faces will be adjacent to the new edge too!
2774  // need to set the senses of new edge within faces
2775 
2776  for( Range::iterator it = faceRange.begin(); it != faceRange.end(); ++it )
2777  {
2778  EntityHandle face = *it;
2779  rval = _mbImpl->add_parent_child( face, new_edge );
2780  MBERRORR( rval, " can't add new edge - face parent relation" );
2781  int sense;
2782  rval = _my_geomTopoTool->get_sense( edge, face, sense );
2783  MBERRORR( rval, " can't get initial sense of edge in the adjacent face" );
2784  // keep the same sense for the new edge within the faces
2785  rval = _my_geomTopoTool->set_sense( new_edge, face, sense );
2786  MBERRORR( rval, " can't set sense of new edge in the adjacent face" );
2787  }
2788 
2789  return MB_SUCCESS;
2790 }

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, 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().

◆ split_boundary()

ErrorCode moab::FBEngine::split_boundary ( EntityHandle  face,
EntityHandle  atNode 
)
private

Definition at line 2792 of file FBEngine.cpp.

2793 {
2794  // find the boundary edges, and split the one that we find it is a part of
2795  if( debug_splits )
2796  {
2797  std::cout << "Split face " << _mbImpl->id_from_handle( face )
2798  << " at node:" << _mbImpl->id_from_handle( atNode ) << "\n";
2799  }
2800  Range bound_edges;
2801  ErrorCode rval = getAdjacentEntities( face, 1, bound_edges );
2802  MBERRORR( rval, " can't get boundary edges" );
2803  bool brokEdge = _brokenEdges.find( atNode ) != _brokenEdges.end();
2804 
2805  for( Range::iterator it = bound_edges.begin(); it != bound_edges.end(); ++it )
2806  {
2807  EntityHandle b_edge = *it;
2808  // get all edges in range
2809  Range mesh_edges;
2810  rval = _mbImpl->get_entities_by_dimension( b_edge, 1, mesh_edges );
2811  MBERRORR( rval, " can't get mesh edges" );
2812  if( brokEdge )
2813  {
2814  EntityHandle brokenEdge = _brokenEdges[atNode];
2815  if( mesh_edges.find( brokenEdge ) != mesh_edges.end() )
2816  {
2817  EntityHandle new_edge;
2818  rval = split_bedge_at_new_mesh_node( b_edge, atNode, brokenEdge, new_edge );
2819  return rval;
2820  }
2821  }
2822  else
2823  {
2824  Range nodes;
2825  rval = _mbImpl->get_connectivity( mesh_edges, nodes );
2826  MBERRORR( rval, " can't get nodes from mesh edges" );
2827 
2828  if( nodes.find( atNode ) != nodes.end() )
2829  {
2830  // we found our boundary edge candidate
2831  EntityHandle new_edge;
2832  rval = split_edge_at_mesh_node( b_edge, atNode, new_edge );
2833  return rval;
2834  }
2835  }
2836  }
2837  // if the node was not found in any "current" boundary, it broke an existing
2838  // boundary edge
2839  MBERRORR( MB_FAILURE, " we did not find an appropriate boundary edge" );
2840  ; //
2841 }

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().

◆ split_edge_at_mesh_node()

ErrorCode moab::FBEngine::split_edge_at_mesh_node ( EntityHandle  edge,
EntityHandle  node,
EntityHandle new_edge 
)

Definition at line 2369 of file FBEngine.cpp.

2370 {
2371  // the node should be in the list of nodes
2372 
2373  int dim = _my_geomTopoTool->dimension( edge );
2374  if( dim != 1 ) return MB_FAILURE; // not an edge
2375 
2376  if( debug_splits )
2377  {
2378  std::cout << "Split edge " << _mbImpl->id_from_handle( edge )
2379  << " with global id: " << _my_geomTopoTool->global_id( edge )
2380  << " at node:" << _mbImpl->id_from_handle( node ) << "\n";
2381  }
2382 
2383  // now get the edges
2384  // the order is important...
2385  // these are ordered sets !!
2386  std::vector< EntityHandle > ents;
2387  ErrorCode rval = _mbImpl->get_entities_by_type( edge, MBEDGE, ents );
2388  if( MB_SUCCESS != rval ) return rval;
2389  if( ents.size() < 1 ) return MB_FAILURE; // no edges
2390 
2391  const EntityHandle* conn = NULL;
2392  int len;
2393  // find the edge connected to the splitting node
2394  int num_mesh_edges = (int)ents.size();
2395  int index_edge;
2396  EntityHandle firstNode = 0;
2397  for( index_edge = 0; index_edge < num_mesh_edges; index_edge++ )
2398  {
2399  rval = MBI->get_connectivity( ents[index_edge], conn, len );
2400  if( MB_SUCCESS != rval ) return rval;
2401  if( index_edge == 0 ) firstNode = conn[0]; // will be used to decide vertex sets adjacencies
2402  if( conn[0] == node )
2403  {
2404  if( index_edge == 0 )
2405  {
2406  new_edge = 0; // no need to split, it is the first node
2407  return MB_SUCCESS; // no split
2408  }
2409  else
2410  return MB_FAILURE; // we should have found the node already , wrong
2411  // connectivity
2412  }
2413  else if( conn[1] == node )
2414  {
2415  // we found the index of interest
2416  break;
2417  }
2418  }
2419  if( index_edge == num_mesh_edges - 1 )
2420  {
2421  new_edge = 0; // no need to split, it is the last node
2422  return MB_SUCCESS; // no split
2423  }
2424 
2425  // here, we have 0 ... index_edge edges in the first set,
2426  // create a vertex set and add the node to it
2427 
2428  if( debug_splits )
2429  {
2430  std::cout << "Split edge with " << num_mesh_edges << " mesh edges, at index (0 based) " << index_edge << "\n";
2431  }
2432 
2433  // at this moment, the node set should have been already created in new_geo_edge
2434  EntityHandle nodeSet; // the node set that has the node (find it!)
2435  bool found = find_vertex_set_for_node( node, nodeSet );
2436 
2437  if( !found )
2438  {
2439  // create a node set and add the node
2440 
2441  // must be an error, but create one nevertheless
2442  rval = _mbImpl->create_meshset( MESHSET_SET, nodeSet );
2443  MBERRORR( rval, "Failed to create a new vertex set" );
2444 
2445  rval = _mbImpl->add_entities( nodeSet, &node, 1 );
2446  MBERRORR( rval, "Failed to add the node to the set" );
2447 
2448  rval = _my_geomTopoTool->add_geo_set( nodeSet, 0 ); //
2449  MBERRORR( rval, "Failed to commit the node set" );
2450 
2451  if( debug_splits )
2452  {
2453  std::cout << " create a vertex set (this should have been created before!)"
2454  << _mbImpl->id_from_handle( nodeSet )
2455  << " global id:" << this->_my_geomTopoTool->global_id( nodeSet ) << "\n";
2456  }
2457  }
2458 
2459  // we need to remove the remaining mesh edges from first set, and add it
2460  // to the second set, in order
2461 
2462  rval = _mbImpl->create_meshset( MESHSET_ORDERED, new_edge );
2463  MBERRORR( rval, "can't create geo edge" );
2464 
2465  int remaining = num_mesh_edges - 1 - index_edge;
2466 
2467  // add edges to the edge set
2468  rval = _mbImpl->add_entities( new_edge, &ents[index_edge + 1], remaining );
2469  MBERRORR( rval, "can't add edges to the new edge" );
2470 
2471  // also, remove the second node set from old edge
2472  // remove the edges index_edge+1 and up
2473 
2474  rval = _mbImpl->remove_entities( edge, &ents[index_edge + 1], remaining );
2475  MBERRORR( rval, "can't remove edges from the old edge" );
2476 
2477  // need to find the adjacent vertex sets
2478  Range vertexRange;
2479  rval = getAdjacentEntities( edge, 0, vertexRange );
2480 
2481  EntityHandle secondSet;
2482  if( vertexRange.size() == 1 )
2483  {
2484  // initially a periodic edge, OK to add the new set to both edges, and the
2485  // second set
2486  secondSet = vertexRange[0];
2487  }
2488  else
2489  {
2490  if( vertexRange.size() > 2 ) return MB_FAILURE; // something must be wrong with too many vertices
2491  // find first node
2492  int k;
2493  for( k = 0; k < 2; k++ )
2494  {
2495  Range verts;
2496  rval = _mbImpl->get_entities_by_type( vertexRange[k], MBVERTEX, verts );
2497 
2498  MBERRORR( rval, "can't get vertices from vertex set" );
2499  if( verts.size() != 1 ) MBERRORR( MB_FAILURE, " node set not defined well" );
2500  if( firstNode == verts[0] )
2501  {
2502  secondSet = vertexRange[1 - k]; // the other set; it is 1 or 0
2503  break;
2504  }
2505  }
2506  if( k >= 2 )
2507  {
2508  // it means we didn't find the right node set
2509  MBERRORR( MB_FAILURE, " can't find the right vertex" );
2510  }
2511  // remove the second set from the connectivity with the
2512  // edge (parent-child relation)
2513  // remove_parent_child( EntityHandle parent,
2514  // EntityHandle child )
2515  rval = _mbImpl->remove_parent_child( edge, secondSet );
2516  MBERRORR( rval, " can't remove the second vertex from edge" );
2517  }
2518  // at this point, we just need to add to both edges the new node set (vertex)
2519  rval = _mbImpl->add_parent_child( edge, nodeSet );
2520  MBERRORR( rval, " can't add new vertex to old edge" );
2521 
2522  rval = _mbImpl->add_parent_child( new_edge, nodeSet );
2523  MBERRORR( rval, " can't add new vertex to new edge" );
2524 
2525  // one thing that I forgot: add the secondSet as a child to new edge!!!
2526  // (so far, the new edge has only one end vertex!)
2527  rval = _mbImpl->add_parent_child( new_edge, secondSet );
2528  MBERRORR( rval, " can't add second vertex to new edge" );
2529 
2530  // now, add the edge and vertex to geo tool
2531 
2532  rval = _my_geomTopoTool->add_geo_set( new_edge, 1 );
2533  MBERRORR( rval, " can't add new edge" );
2534 
2535  // next, get the adjacent faces to initial edge, and add them as parents
2536  // to the new edge
2537 
2538  // need to find the adjacent face sets
2539  Range faceRange;
2540  rval = getAdjacentEntities( edge, 2, faceRange );
2541 
2542  // these faces will be adjacent to the new edge too!
2543  // need to set the senses of new edge within faces
2544 
2545  for( Range::iterator it = faceRange.begin(); it != faceRange.end(); ++it )
2546  {
2547  EntityHandle face = *it;
2548  rval = _mbImpl->add_parent_child( face, new_edge );
2549  MBERRORR( rval, " can't add new edge - face parent relation" );
2550  int sense;
2551  rval = _my_geomTopoTool->get_sense( edge, face, sense );
2552  MBERRORR( rval, " can't get initial sense of edge in the adjacent face" );
2553  // keep the same sense for the new edge within the faces
2554  rval = _my_geomTopoTool->set_sense( new_edge, face, sense );
2555  MBERRORR( rval, " can't set sense of new edge in the adjacent face" );
2556  }
2557 
2558  return MB_SUCCESS;
2559 }

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, 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().

◆ split_edge_at_point()

ErrorCode moab::FBEngine::split_edge_at_point ( EntityHandle  edge,
CartVect point,
EntityHandle new_edge 
)

Definition at line 2328 of file FBEngine.cpp.

2329 {
2330  // return MB_NOT_IMPLEMENTED;
2331  // first, we need to find the closest point on the smooth edge, then
2332  // split the smooth edge, then call the split_edge_at_mesh_node
2333  // or maybe just find the closest node??
2334  // first of all, we need to find a point on the smooth edge, close by
2335  // then split the mesh edge (if needed)
2336  // this would be quite a drama, as the new edge has to be inserted in
2337  // order for proper geo edge definition
2338 
2339  // first of all, find a closest point
2340  // usually called for
2341  if( debug_splits )
2342  {
2343  std::cout << "Split edge " << _mbImpl->id_from_handle( edge ) << " at point:" << point << "\n";
2344  }
2345  int dim = _my_geomTopoTool->dimension( edge );
2346  if( dim != 1 ) return MB_FAILURE;
2347  if( !_smooth ) return MB_FAILURE; // call it only for smooth option...
2348  // maybe we should do something for linear too
2349 
2350  SmoothCurve* curve = this->_edges[edge];
2351  EntityHandle closeNode;
2352  int edgeIndex;
2353  double u = curve->u_from_position( point[0], point[1], point[2], closeNode, edgeIndex );
2354  if( 0 == closeNode )
2355  {
2356  // we really need to split an existing edge
2357  // do not do that yet
2358  // evaluate from u:
2359  /*double pos[3];
2360  curve->position_from_u(u, pos[0], pos[1], pos[2] );*/
2361  // create a new node here, and split one edge
2362  // change also connectivity, create new triangles on both sides ...
2363  std::cout << "not found a close node, u is: " << u << " edge index: " << edgeIndex << "\n";
2364  return MB_FAILURE; // not implemented yet
2365  }
2366  return split_edge_at_mesh_node( edge, closeNode, new_edge );
2367 }

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

◆ split_internal_edge()

ErrorCode moab::FBEngine::split_internal_edge ( EntityHandle edge,
EntityHandle newVertex 
)
private

Definition at line 3085 of file FBEngine.cpp.

3086 {
3087  // split the edge, and form 4 new triangles and 2 edges
3088  // get 2 triangles adjacent to edge
3089  Range adj_tri;
3090  ErrorCode rval = _mbImpl->get_adjacencies( &edge, 1, 2, false, adj_tri );
3091  MBERRORR( rval, "can't get adj_tris" );
3092  adj_tri = subtract( adj_tri, _piercedTriangles );
3093  if( adj_tri.size() >= 3 )
3094  {
3095  std::cout << "WARNING: non manifold geometry. Are you sure?";
3096  }
3097  for( Range::iterator it = adj_tri.begin(); it != adj_tri.end(); ++it )
3098  {
3099  EntityHandle tri = *it;
3100  _piercedTriangles.insert( tri );
3101  const EntityHandle* conn3;
3102  int nnodes;
3103  rval = _mbImpl->get_connectivity( tri, conn3, nnodes );
3104  MBERRORR( rval, "can't get nodes" );
3105  int num1, sense, offset;
3106  rval = _mbImpl->side_number( tri, edge, num1, sense, offset );
3107  MBERRORR( rval, "can't get side number" );
3108  // after we know the side number, we can split in 2 triangles
3109  // node i is opposite to edge i
3110  int num2 = ( num1 + 1 ) % 3;
3111  int num3 = ( num2 + 1 ) % 3;
3112  // the edge from num1 to num2 is split into 2 edges
3113  EntityHandle t1[] = { conn3[num2], conn3[num3], newVertex };
3114  EntityHandle t2[] = { conn3[num1], newVertex, conn3[num3] };
3115  EntityHandle newTriangle, newTriangle2;
3116  rval = _mbImpl->create_element( MBTRI, t1, 3, newTriangle );
3117  MBERRORR( rval, "can't create triangle" );
3118  _newTriangles.insert( newTriangle );
3119  rval = _mbImpl->create_element( MBTRI, t2, 3, newTriangle2 );
3120  MBERRORR( rval, "can't create triangle" );
3121  _newTriangles.insert( newTriangle2 );
3122  // create edges with this, indirectly
3123  std::vector< EntityHandle > edges0;
3124  rval = _mbImpl->get_adjacencies( &newTriangle, 1, 1, true, edges0 );
3125  MBERRORR( rval, "can't get new edges" );
3126  edges0.clear();
3127  rval = _mbImpl->get_adjacencies( &newTriangle2, 1, 1, true, edges0 );
3128  MBERRORR( rval, "can't get new edges" );
3129  if( debug_splits )
3130  {
3131  std::cout << "2 (out of 4) triangles formed:\n";
3132  _mbImpl->list_entity( newTriangle );
3133  print_debug_triangle( newTriangle );
3134  _mbImpl->list_entity( newTriangle2 );
3135  print_debug_triangle( newTriangle2 );
3136  }
3137  }
3138  return MB_SUCCESS;
3139 }

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_quads()

ErrorCode moab::FBEngine::split_quads ( )
private

Definition at line 2969 of file FBEngine.cpp.

2970 {
2971  // first see if we have any quads in the 2d faces
2972  // _my_gsets[2] is a range of surfaces (moab sets of dimension 2)
2973  int num_quads = 0;
2974  for( Range::iterator it = _my_gsets[2].begin(); it != _my_gsets[2].end(); ++it )
2975  {
2976  EntityHandle surface = *it;
2977  int num_q = 0;
2978  _mbImpl->get_number_entities_by_type( surface, MBQUAD, num_q );
2979  num_quads += num_q;
2980  }
2981 
2982  if( num_quads == 0 ) return MB_SUCCESS; // nothing to do
2983 
2984  GeomTopoTool* new_gtt = NULL;
2985  ErrorCode rval = _my_geomTopoTool->duplicate_model( new_gtt );
2986  MBERRORR( rval, "can't duplicate model" );
2987  if( this->_t_created ) delete _my_geomTopoTool;
2988 
2989  _t_created = true; // this one is for sure created here, even if the original gtt was not
2990 
2991  // if we were using smart pointers, we would decrease the reference to the _my_geomTopoTool, at
2992  // least
2993  _my_geomTopoTool = new_gtt;
2994 
2995  // replace the _my_gsets with the new ones, from the new set
2997 
2998  // we have duplicated now the model, and the quads are part of the new _my_gsets[2]
2999  // we will split them now, by creating triangles along the smallest diagonal
3000  // maybe we will come up with a better criteria, but for the time being, it should be fine.
3001  // what we will do: we will remove all quads from the surface sets, and split them
3002 
3003  for( Range::iterator it2 = _my_gsets[2].begin(); it2 != _my_gsets[2].end(); ++it2 )
3004  {
3005  EntityHandle surface = *it2;
3006  Range quads;
3007  rval = _mbImpl->get_entities_by_type( surface, MBQUAD, quads );
3008  MBERRORR( rval, "can't get quads from the surface set" );
3009  rval = _mbImpl->remove_entities( surface, quads );
3010  MBERRORR( rval, "can't remove quads from the surface set" ); // they are not deleted, just
3011  // removed from the set
3012  for( Range::iterator it = quads.begin(); it != quads.end(); ++it )
3013  {
3014  EntityHandle quad = *it;
3015  int nnodes;
3016  const EntityHandle* conn;
3017  rval = _mbImpl->get_connectivity( quad, conn, nnodes );
3018  MBERRORR( rval, "can't get quad connectivity" );
3019  // get all vertices position, to see which one is the shortest diagonal
3020  CartVect pos[4];
3021  rval = _mbImpl->get_coords( conn, 4, (double*)&pos[0] );
3022  MBERRORR( rval, "can't get node coordinates" );
3023  bool diag1 = ( ( pos[2] - pos[0] ).length_squared() < ( pos[3] - pos[1] ).length_squared() );
3024  EntityHandle newTris[2];
3025  EntityHandle tri1[3] = { conn[0], conn[1], conn[2] };
3026  EntityHandle tri2[3] = { conn[0], conn[2], conn[3] };
3027  if( !diag1 )
3028  {
3029  tri1[2] = conn[3];
3030  tri2[0] = conn[1];
3031  }
3032  rval = _mbImpl->create_element( MBTRI, tri1, 3, newTris[0] );
3033  MBERRORR( rval, "can't create triangle 1" );
3034  rval = _mbImpl->create_element( MBTRI, tri2, 3, newTris[1] );
3035  MBERRORR( rval, "can't create triangle 2" );
3036  rval = _mbImpl->add_entities( surface, newTris, 2 );
3037  MBERRORR( rval, "can't add triangles to the set" );
3038  }
3039  //
3040  }
3041  return MB_SUCCESS;
3042 }

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().

◆ split_surface()

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 1509 of file FBEngine.cpp.

1513 {
1514  // use now the chained edges to create a new face (loop or clean split)
1515  // use a fill to determine the new sets, up to the polyline
1516  // points on the polyline will be moved to the closest point location, with some constraints
1517  // then the sets will be reset, geometry recomputed. new vertices, new edges, etc.
1518 
1519  Range iniTris;
1520  ErrorCode rval;
1521  rval = _mbImpl->get_entities_by_type( face, MBTRI, iniTris );
1522  MBERRORR( rval, "can't get initial triangles" );
1523 
1524  // start from a triangle that is not in the triangles to delete
1525  // flood fill
1526 
1527  bool closed = splittingNodes.size() == 0;
1528  if( !closed )
1529  {
1530  //
1531  if( splittingNodes.size() != 2 ) MBERRORR( MB_FAILURE, "need to have exactly 2 nodes for splitting" );
1532  // we will have to split the boundary edges
1533  // first, find the actual boundary, and try to split with the 2 end points (nodes)
1534  // get the adjacent edges, and see which one has the end nodes
1535  rval = split_boundary( face, splittingNodes[0] );
1536  MBERRORR( rval, "can't split with first node" );
1537  rval = split_boundary( face, splittingNodes[1] );
1538  MBERRORR( rval, "can't split with second node)" );
1539  }
1540  // we will separate triangles to delete, unaffected, new_triangles,
1541  // nodesAlongPolyline,
1542  Range first, second;
1543  rval = separate( face, chainedEdges, first, second );
1544 
1545  // now, we are done with the computations;
1546  // we need to put the new nodes on the smooth surface
1547  if( this->_smooth )
1548  {
1549  rval = smooth_new_intx_points( face, chainedEdges );
1550  MBERRORR( rval, "can't smooth new points" );
1551  }
1552 
1553  // create the new set
1554  rval = _mbImpl->create_meshset( MESHSET_SET, newFace );
1555  MBERRORR( rval, "can't create a new face" );
1556 
1557  _my_geomTopoTool->add_geo_set( newFace, 2 );
1558 
1559  // the new face will have the first set (positive sense triangles, to the left)
1560  rval = _mbImpl->add_entities( newFace, first );
1561  MBERRORR( rval, "can't add first range triangles to new face" );
1562 
1563  for( unsigned int j = 0; j < chainedEdges.size(); j++ )
1564  {
1565  EntityHandle new_geo_edge = chainedEdges[j];
1566  // both faces will have the edge now
1567  rval = _mbImpl->add_parent_child( face, new_geo_edge );
1568  MBERRORR( rval, "can't add parent child relations for new edge" );
1569 
1570  rval = _mbImpl->add_parent_child( newFace, new_geo_edge );
1571  MBERRORR( rval, "can't add parent child relations for new edge" );
1572  // add senses
1573  // sense newFace is 1, old face is -1
1574  rval = _my_geomTopoTool->set_sense( new_geo_edge, newFace, 1 );
1575  MBERRORR( rval, "can't set sense for new edge" );
1576 
1577  rval = _my_geomTopoTool->set_sense( new_geo_edge, face, -1 );
1578  MBERRORR( rval, "can't set sense for new edge in original face" );
1579  }
1580 
1581  rval = set_neumann_tags( face, newFace );
1582  MBERRORR( rval, "can't set NEUMANN set tags" );
1583 
1584  // now, we should remove from the original set all tris, and put the "second" range
1585  rval = _mbImpl->remove_entities( face, iniTris );
1586  MBERRORR( rval, "can't remove original tris from initial face set" );
1587 
1588  rval = _mbImpl->add_entities( face, second );
1589  MBERRORR( rval, "can't add second range to the original set" );
1590 
1591  if( !closed )
1592  {
1593  rval = redistribute_boundary_edges_to_faces( face, newFace, chainedEdges );
1594  MBERRORR( rval, "fail to reset the proper boundary faces" );
1595  }
1596 
1597  /*if (_smooth)
1598  delete_smooth_tags();// they need to be recomputed, anyway
1599  // this will remove the extra smooth faces and edges
1600  clean();*/
1601  // also, these nodes need to be moved to the smooth surface, sometimes before deleting the old
1602  // triangles
1603  // remove the triangles from the set, then delete triangles (also some edges need to be
1604  // deleted!)
1606 
1607  MBERRORR( rval, "can't delete triangles" );
1609  // delete edges that are broke up in 2
1611  MBERRORR( rval, "can't delete edges" );
1612  _piercedEdges.clear();
1613 
1614  if( debug_splits )
1615  {
1616  _mbImpl->write_file( "newFace.vtk", "vtk", 0, &newFace, 1 );
1617  _mbImpl->write_file( "leftoverFace.vtk", "vtk", 0, &face, 1 );
1618  }
1619  return MB_SUCCESS;
1620 }

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().

◆ split_surface_with_direction()

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 1208 of file FBEngine.cpp.

1214 {
1215 
1216  // first of all, find all intersection points (piercing in the face along the direction)
1217  // assume it is robust; what if it is not sufficiently robust?
1218  // if the polyline is open, find the intersection with the boundary edges, of the
1219  // polyline extruded at ends
1220 
1221  ErrorCode rval;
1222 
1223  // then find the position
1224  int numIniPoints = (int)xyz.size() / 3;
1225  if( ( closed && numIniPoints < 3 ) || ( !closed && numIniPoints < 2 ) )
1226  MBERRORR( MB_FAILURE, "not enough polyline points " );
1227  EntityHandle rootForFace;
1228 
1229  rval = _my_geomTopoTool->get_root( face, rootForFace );
1230  MBERRORR( rval, "Failed to get root of face." );
1231 
1232  const double dir[] = { direction[0], direction[1], direction[2] };
1233  std::vector< EntityHandle > nodes; // get the nodes closest to the ray traces of interest
1234 
1235  // these are nodes on the boundary of original face;
1236  // if the cutting line is not closed, the starting - ending vertices of the
1237  // polygonal line must come from this list
1238 
1239  std::vector< CartVect > b_pos;
1240  std::vector< EntityHandle > boundary_nodes;
1241  std::vector< EntityHandle > splittingNodes;
1242  Range boundary_mesh_edges;
1243  if( !closed )
1244  {
1245  rval = boundary_nodes_on_face( face, boundary_nodes );
1246  MBERRORR( rval, "Failed to get boundary nodes." );
1247  b_pos.resize( boundary_nodes.size() );
1248  rval = _mbImpl->get_coords( &( boundary_nodes[0] ), boundary_nodes.size(), (double*)( &b_pos[0][0] ) );
1249  MBERRORR( rval, "Failed to get coordinates for boundary nodes." );
1250  rval = boundary_mesh_edges_on_face( face, boundary_mesh_edges );
1251  MBERRORR( rval, "Failed to get mesh boundary edges for face." );
1252  }
1253  //
1254  int i = 0;
1255  CartVect dirct( direction );
1256  dirct.normalize(); // maybe an overkill?
1257  for( ; i < numIniPoints; i++ )
1258  {
1259 
1260  const double point[] = { xyz[3 * i], xyz[3 * i + 1], xyz[3 * i + 2] }; // or even point( &(xyz[3*i]) ); //
1261  CartVect p1( point );
1262  if( !closed && ( ( 0 == i ) || ( numIniPoints - 1 == i ) ) )
1263  {
1264 
1265  // find the intersection point between a plane and boundary mesh edges
1266  // this will be the closest point on the boundary of face
1267  /// the segment is the first or last segment in the polyline
1268  int i1 = i + 1;
1269  if( i == numIniPoints - 1 ) i1 = i - 1; // previous point if the last
1270  // the direction is from point to point1
1271  const double point1[] = { xyz[3 * i1], xyz[3 * i1 + 1], xyz[3 * i1 + 2] };
1272  CartVect p2( point1 );
1273  CartVect normPlane = ( p2 - p1 ) * dirct;
1274  normPlane.normalize();
1275  //(roughly, from p1 to p2, perpendicular to dirct, in the "xy" plane
1276  // if the intx point is "outside" p1 - p2, skip if the intx point is closer to p2
1277  CartVect perpDir = dirct * normPlane;
1278  Range::iterator ite = boundary_mesh_edges.begin();
1279  // do a linear search for the best intersection point position (on a boundary edge)
1280  if( debug_splits )
1281  {
1282  std::cout << " p1:" << p1 << "\n";
1283  std::cout << " p2:" << p2 << "\n";
1284  std::cout << " perpDir:" << perpDir << "\n";
1285  std::cout << " boundary edges size:" << boundary_mesh_edges.size() << "\n";
1286  }
1287  for( ; ite != boundary_mesh_edges.end(); ++ite )
1288  {
1289  EntityHandle candidateEdge = *ite;
1290  const EntityHandle* conn2;
1291  int nno;
1292  rval = _mbImpl->get_connectivity( candidateEdge, conn2, nno );
1293  MBERRORR( rval, "Failed to get conn for boundary edge" );
1294  CartVect pts[2];
1295  rval = _mbImpl->get_coords( conn2, 2, &( pts[0][0] ) );
1296  MBERRORR( rval, "Failed to get coords of nodes for boundary edge" );
1297  CartVect intx_point;
1298  double parPos;
1299  bool intersect =
1300  intersect_segment_and_plane_slice( pts[0], pts[1], p1, p2, dirct, normPlane, intx_point, parPos );
1301  if( debug_splits )
1302  {
1303  std::cout << " Edge:" << _mbImpl->id_from_handle( candidateEdge ) << "\n";
1304  std::cout << " Node 1:" << _mbImpl->id_from_handle( conn2[0] ) << pts[0] << "\n";
1305  std::cout << " Node 2:" << _mbImpl->id_from_handle( conn2[1] ) << pts[1] << "\n";
1306  std::cout << " Intersect bool:" << intersect << "\n";
1307  }
1308  if( intersect )
1309  {
1310  double proj1 = ( intx_point - p1 ) % perpDir;
1311  double proj2 = ( intx_point - p2 ) % perpDir;
1312  if( ( fabs( proj1 ) > fabs( proj2 ) ) // this means it is closer to p2 than p1
1313  )
1314  continue; // basically, this means the intersection point is with a
1315  // boundary edge on the other side, closer to p2 than p1, so we
1316  // skip it
1317  if( parPos == 0 )
1318  {
1319  // close to vertex 1, nothing to do
1320  nodes.push_back( conn2[0] );
1321  splittingNodes.push_back( conn2[0] );
1322  }
1323  else if( parPos == 1. )
1324  {
1325  // close to vertex 2, nothing to do
1326  nodes.push_back( conn2[1] );
1327  splittingNodes.push_back( conn2[1] );
1328  }
1329  else
1330  {
1331  // break the edge, create a new node at intersection point (will be smoothed
1332  // out)
1333  EntityHandle newVertex;
1334  rval = _mbImpl->create_vertex( &( intx_point[0] ), newVertex );
1335  MBERRORR( rval, "can't create vertex" );
1336  nodes.push_back( newVertex );
1337  split_internal_edge( candidateEdge, newVertex );
1338  splittingNodes.push_back( newVertex );
1339  _brokenEdges[newVertex] = candidateEdge;
1340  _piercedEdges.insert( candidateEdge );
1341  }
1342  break; // break from the loop over boundary edges, we are interested in the
1343  // first
1344  // split (hopefully, the only split)
1345  }
1346  }
1347  if( ite == boundary_mesh_edges.end() )
1348  MBERRORR( MB_FAILURE, "Failed to find boundary intersection edge. Bail out" );
1349  }
1350  else
1351  {
1352  std::vector< double > distances_out;
1353  std::vector< EntityHandle > sets_out;
1354  std::vector< EntityHandle > facets_out;
1355  rval = _my_geomTopoTool->obb_tree()->ray_intersect_sets( distances_out, sets_out, facets_out, rootForFace,
1356  tolerance, point, dir );
1357  MBERRORR( rval, "Failed to get ray intersections." );
1358  if( distances_out.size() < 1 )
1359  MBERRORR( MB_FAILURE, "Failed to get one intersection point, bad direction." );
1360 
1361  if( distances_out.size() > 1 )
1362  {
1363  std::cout << " too many intersection points. Only the first one considered\n";
1364  }
1365  std::vector< EntityHandle >::iterator pFace = std::find( sets_out.begin(), sets_out.end(), face );
1366 
1367  if( pFace == sets_out.end() ) MBERRORR( MB_FAILURE, "Failed to intersect given face, bad direction." );
1368  unsigned int index = pFace - sets_out.begin();
1369  // get the closest node of the triangle, and modify locally the triangle(s), so the
1370  // intersection point is at a new vertex, if needed
1371  CartVect P( point );
1372  CartVect Dir( dir );
1373  CartVect newPoint = P + distances_out[index] * Dir;
1374  // get the triangle coordinates
1375 
1376  double area_coord[3];
1377  EntityHandle boundary_handle = 0; // if 0, not on a boundary
1378  rval = area_coordinates( _mbImpl, facets_out[index], newPoint, area_coord, boundary_handle );
1379  MBERRORR( rval, "Failed to get area coordinates" );
1380 
1381  if( debug_splits )
1382  {
1383  std::cout << " int point:" << newPoint << " area coord " << area_coord[0] << " " << area_coord[1] << " "
1384  << area_coord[2] << "\n";
1385  std::cout << " triangle: " << _mbImpl->id_from_handle( facets_out[index] )
1386  << " boundary:" << boundary_handle << "\n";
1387  }
1388  EntityType type;
1389  if( boundary_handle ) type = _mbImpl->type_from_handle( boundary_handle );
1390  if( boundary_handle && ( type == MBVERTEX ) )
1391  {
1392  // nothing to do, we are happy
1393  nodes.push_back( boundary_handle );
1394  }
1395  else
1396  {
1397  // for an edge, we will split 2 triangles
1398  // for interior point, we will create 3 triangles out of one
1399  // create a new vertex
1400  EntityHandle newVertex;
1401  rval = _mbImpl->create_vertex( &( newPoint[0] ), newVertex );
1402  if( boundary_handle ) // this is edge
1403  {
1404  split_internal_edge( boundary_handle, newVertex );
1405  _piercedEdges.insert( boundary_handle ); // to be removed at the end
1406  }
1407  else
1408  divide_triangle( facets_out[index], newVertex );
1409 
1410  nodes.push_back( newVertex );
1411  }
1412  }
1413  }
1414  // now, we have to find more intersection points, either interior to triangles, or on edges, or
1415  // on vertices use the same tolerance as before starting from 2 points on 2 triangles, and
1416  // having the direction, get more intersection points between the plane formed by direction and
1417  // those 2 points, and edges from triangulation (the triangles involved will be part of the same
1418  // gentity , original face ( moab set)
1419  // int closed = 1;// closed = 0 if the polyline is not closed
1420 
1421  CartVect Dir( direction );
1422  std::vector< EntityHandle > chainedEdges;
1423 
1424  for( i = 0; i < numIniPoints - 1 + closed; i++ )
1425  {
1426  int nextIndex = ( i + 1 ) % numIniPoints;
1427  std::vector< EntityHandle > trianglesAlong;
1428  std::vector< CartVect > points;
1429  // otherwise to edges or even nodes
1430  std::vector< EntityHandle > entities;
1431  // start with initial points, intersect along the direction, find the facets
1432  rval = compute_intersection_points( face, nodes[i], nodes[nextIndex], Dir, points, entities, trianglesAlong );
1433  MBERRORR( rval, "can't get intersection points along a line" );
1434  std::vector< EntityHandle > nodesAlongPolyline;
1435  // refactor code; move some edge creation for each 2 intersection points
1436  nodesAlongPolyline.push_back( entities[0] ); // it is for sure a node
1437  int num_points = (int)points.size(); // it should be num_triangles + 1
1438  for( int j = 0; j < num_points - 1; j++ )
1439  {
1440  EntityHandle tri = trianglesAlong[j]; // this is happening in trianglesAlong i
1441  EntityHandle e1 = entities[j];
1442  EntityHandle e2 = entities[j + 1];
1443  EntityType et1 = _mbImpl->type_from_handle( e1 );
1444  // EntityHandle vertex1 = nodesAlongPolyline[i];// irrespective of the entity type i,
1445  // we already have the vertex there
1446  EntityType et2 = _mbImpl->type_from_handle( e2 );
1447  if( et2 == MBVERTEX )
1448  {
1449  nodesAlongPolyline.push_back( e2 );
1450  }
1451  else // if (et2==MBEDGE)
1452  {
1453  CartVect coord_vert = points[j + 1];
1454  EntityHandle newVertex;
1455  rval = _mbImpl->create_vertex( (double*)&coord_vert, newVertex );
1456  MBERRORR( rval, "can't create vertex" );
1457  nodesAlongPolyline.push_back( newVertex );
1458  }
1459  // if vertices, do not split anything, just get the edge for polyline
1460  if( et2 == MBVERTEX && et1 == MBVERTEX )
1461  {
1462  // nothing to do, just continue;
1463  continue; // continue the for loop
1464  }
1465 
1466  if( debug_splits )
1467  {
1468  std::cout << "tri: type: " << _mbImpl->type_from_handle( tri )
1469  << " id:" << _mbImpl->id_from_handle( tri ) << "\n e1:" << e1
1470  << " id:" << _mbImpl->id_from_handle( e1 ) << " e2:" << e2
1471  << " id:" << _mbImpl->id_from_handle( e2 ) << "\n";
1472  }
1473  // here, at least one is an edge
1474  rval = BreakTriangle2( tri, e1, e2, nodesAlongPolyline[j], nodesAlongPolyline[j + 1] );
1475  MBERRORR( rval, "can't break triangle 2" );
1476  if( et2 == MBEDGE ) _piercedEdges.insert( e2 );
1477  _piercedTriangles.insert( tri );
1478  }
1479  // nodesAlongPolyline will define a new geometric edge
1480  if( debug_splits )
1481  {
1482  std::cout << "nodesAlongPolyline: " << nodesAlongPolyline.size() << "\n";
1483  std::cout << "points: " << num_points << "\n";
1484  }
1485  // if needed, create edges along polyline, or revert the existing ones, to
1486  // put them in a new edge set
1487  EntityHandle new_geo_edge;
1488  rval = create_new_gedge( nodesAlongPolyline, new_geo_edge );
1489  MBERRORR( rval, "can't create a new edge" );
1490  chainedEdges.push_back( new_geo_edge );
1491  // end copy
1492  }
1493  // the segment between point_i and point_i+1 is in trianglesAlong_i
1494  // points_i is on entities_i
1495  // all these edges are oriented correctly
1496  rval = split_surface( face, chainedEdges, splittingNodes, oNewFace );
1497  MBERRORR( rval, "can't split surface" );
1498  //
1499  rval = chain_edges( min_dot ); // acos(0.8)~= 36 degrees
1500  MBERRORR( rval, "can't chain edges" );
1501  return MB_SUCCESS;
1502 }

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(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::GeomTopoTool::get_root(), moab::Interface::id_from_handle(), moab::index, 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().

◆ weave_lateral_face_from_edges()

ErrorCode moab::FBEngine::weave_lateral_face_from_edges ( EntityHandle  bEdge,
EntityHandle  tEdge,
double *  direction,
EntityHandle newLatFace 
)

Definition at line 3279 of file FBEngine.cpp.

3283 {
3284  // in weird cases might need to create new vertices in the interior;
3285  // in most cases, it is OK
3286 
3287  ErrorCode rval = _mbImpl->create_meshset( MESHSET_SET, newLatFace );
3288  MBERRORR( rval, "can't create new lateral face" );
3289 
3290  EntityHandle v[4]; // vertex sets
3291  // bot edge will be v1->v2
3292  // top edge will be v3->v4
3293  // we need to create edges from v1 to v3 and from v2 to v4
3294  std::vector< EntityHandle > adj;
3295  rval = _mbImpl->get_child_meshsets( bEdge, adj );
3296  MBERRORR( rval, "can't get children nodes" );
3297  bool periodic = false;
3298  if( adj.size() == 1 )
3299  {
3300  v[0] = v[1] = adj[0];
3301  periodic = true;
3302  }
3303  else
3304  {
3305  v[0] = adj[0];
3306  v[1] = adj[1];
3307  }
3308  int senseB;
3309  rval = getEgVtxSense( bEdge, v[0], v[1], senseB );
3310  MBERRORR( rval, "can't get bottom edge sense" );
3311  if( -1 == senseB )
3312  {
3313  v[1] = adj[0]; // so , bEdge will be oriented from v[0] to v[1], and will start at
3314  // nodes1, coords1..
3315  v[0] = adj[1];
3316  }
3317  adj.clear();
3318  rval = _mbImpl->get_child_meshsets( tEdge, adj );
3319  MBERRORR( rval, "can't get children nodes" );
3320  if( adj.size() == 1 )
3321  {
3322  v[2] = v[3] = adj[0];
3323  if( !periodic ) MBERRORR( MB_FAILURE, "top edge is periodic, but bottom edge is not" );
3324  }
3325  else
3326  {
3327  v[2] = adj[0];
3328  v[3] = adj[1];
3329  if( periodic ) MBERRORR( MB_FAILURE, "top edge is not periodic, but bottom edge is" );
3330  }
3331 
3332  // now, using nodes on bottom edge and top edge, create triangles, oriented outwards the
3333  // volume (sense positive on bottom edge)
3334  std::vector< EntityHandle > nodes1;
3335  rval = get_nodes_from_edge( bEdge, nodes1 );
3336  MBERRORR( rval, "can't get nodes from bott edge" );
3337 
3338  std::vector< EntityHandle > nodes2;
3339  rval = get_nodes_from_edge( tEdge, nodes2 );
3340  MBERRORR( rval, "can't get nodes from top edge" );
3341 
3342  std::vector< CartVect > coords1, coords2;
3343  coords1.resize( nodes1.size() );
3344  coords2.resize( nodes2.size() );
3345 
3346  int N1 = (int)nodes1.size();
3347  int N2 = (int)nodes2.size();
3348 
3349  rval = _mbImpl->get_coords( &( nodes1[0] ), nodes1.size(), (double*)&( coords1[0] ) );
3350  MBERRORR( rval, "can't get coords of nodes from bott edge" );
3351 
3352  rval = _mbImpl->get_coords( &( nodes2[0] ), nodes2.size(), (double*)&( coords2[0] ) );
3353  MBERRORR( rval, "can't get coords of nodes from top edge" );
3354  CartVect up( direction );
3355 
3356  // see if the start and end coordinates are matching, if not, reverse edge 2 nodes and
3357  // coordinates
3358  CartVect v1 = ( coords1[0] - coords2[0] ) * up;
3359  CartVect v2 = ( coords1[0] - coords2[N2 - 1] ) * up;
3360  if( v1.length_squared() > v2.length_squared() )
3361  {
3362  // we need to reverse coords2 and node2, as nodes are not above each other
3363  // the edges need to be found between v0 and v3, v1 and v2!
3364  for( unsigned int k = 0; k < nodes2.size() / 2; k++ )
3365  {
3366  EntityHandle tmp = nodes2[k];
3367  nodes2[k] = nodes2[N2 - 1 - k];
3368  nodes2[N2 - 1 - k] = tmp;
3369  CartVect tv = coords2[k];
3370  coords2[k] = coords2[N2 - 1 - k];
3371  coords2[N2 - 1 - k] = tv;
3372  }
3373  }
3374  // make sure v[2] has nodes2[0], if not, reverse v[2] and v[3]
3375  if( !_mbImpl->contains_entities( v[2], &( nodes2[0] ), 1 ) )
3376  {
3377  // reverse v[2] and v[3], so v[2] will be above v[0]
3378  EntityHandle tmp = v[2];
3379  v[2] = v[3];
3380  v[3] = tmp;
3381  }
3382  // now we know that v[0]--v[3] will be vertex sets in the order we want
3383  EntityHandle nd[4] = { nodes1[0], nodes1[N1 - 1], nodes2[0], nodes2[N2 - 1] };
3384 
3385  adj.clear();
3386  EntityHandle e1, e2;
3387  // find edge 1 between v[0] and v[2], and e2 between v[1] and v[3]
3388  rval = _mbImpl->get_parent_meshsets( v[0], adj );
3389  MBERRORR( rval, "can't get edges connected to vertex set 1" );
3390  bool found = false;
3391  for( unsigned int j = 0; j < adj.size(); j++ )
3392  {
3393  EntityHandle ed = adj[j];
3394  Range vertices;
3395  rval = _mbImpl->get_child_meshsets( ed, vertices );
3396  if( vertices.find( v[2] ) != vertices.end() )
3397  {
3398  found = true;
3399  e1 = ed;
3400  break;
3401  }
3402  }
3403  if( !found )
3404  {
3405  // create an edge from v[0] to v[2]
3406  rval = _mbImpl->create_meshset( MESHSET_SET, e1 );
3407  MBERRORR( rval, "can't create edge 1" );
3408 
3409  rval = _mbImpl->add_parent_child( e1, v[0] );
3410  MBERRORR( rval, "can't add parent - child relation" );
3411 
3412  rval = _mbImpl->add_parent_child( e1, v[2] );
3413  MBERRORR( rval, "can't add parent - child relation" );
3414 
3415  EntityHandle nn2[2] = { nd[0], nd[2] };
3416  EntityHandle medge;
3417  rval = _mbImpl->create_element( MBEDGE, nn2, 2, medge );
3418  MBERRORR( rval, "can't create mesh edge" );
3419 
3420  rval = _mbImpl->add_entities( e1, &medge, 1 );
3421  MBERRORR( rval, "can't add mesh edge to geo edge" );
3422 
3423  rval = this->_my_geomTopoTool->add_geo_set( e1, 1 );
3424  MBERRORR( rval, "can't add edge to gtt" );
3425  }
3426 
3427  // find the edge from v2 to v4 (if not, create one)
3428  rval = _mbImpl->get_parent_meshsets( v[1], adj );
3429  MBERRORR( rval, "can't get edges connected to vertex set 2" );
3430  found = false;
3431  for( unsigned int i = 0; i < adj.size(); i++ )
3432  {
3433  EntityHandle ed = adj[i];
3434  Range vertices;
3435  rval = _mbImpl->get_child_meshsets( ed, vertices );
3436  if( vertices.find( v[3] ) != vertices.end() )
3437  {
3438  found = true;
3439  e2 = ed;
3440  break;
3441  }
3442  }
3443  if( !found )
3444  {
3445  // create an edge from v2 to v4
3446  rval = _mbImpl->create_meshset( MESHSET_SET, e2 );
3447  MBERRORR( rval, "can't create edge 1" );
3448 
3449  rval = _mbImpl->add_parent_child( e2, v[1] );
3450  MBERRORR( rval, "can't add parent - child relation" );
3451 
3452  rval = _mbImpl->add_parent_child( e2, v[3] );
3453  MBERRORR( rval, "can't add parent - child relation" );
3454 
3455  EntityHandle nn2[2] = { nd[1], nd[3] };
3456  EntityHandle medge;
3457  rval = _mbImpl->create_element( MBEDGE, nn2, 2, medge );
3458  MBERRORR( rval, "can't create mesh edge" );
3459 
3460  rval = _mbImpl->add_entities( e2, &medge, 1 );
3461  MBERRORR( rval, "can't add mesh edge to geo edge" );
3462 
3463  rval = _my_geomTopoTool->add_geo_set( e2, 1 );
3464  MBERRORR( rval, "can't add edge to gtt" );
3465  }
3466 
3467  // now we have the four edges, add them to the face, as children
3468 
3469  // add children to face
3470  rval = _mbImpl->add_parent_child( newLatFace, bEdge );
3471  MBERRORR( rval, "can't add parent - child relation" );
3472 
3473  rval = _mbImpl->add_parent_child( newLatFace, tEdge );
3474  MBERRORR( rval, "can't add parent - child relation" );
3475 
3476  rval = _mbImpl->add_parent_child( newLatFace, e1 );
3477  MBERRORR( rval, "can't add parent - child relation" );
3478 
3479  rval = _mbImpl->add_parent_child( newLatFace, e2 );
3480  MBERRORR( rval, "can't add parent - child relation" );
3481 
3482  rval = _my_geomTopoTool->add_geo_set( newLatFace, 2 );
3483  MBERRORR( rval, "can't add face to gtt" );
3484  // add senses
3485  //
3486  rval = _my_geomTopoTool->set_sense( bEdge, newLatFace, 1 );
3487  MBERRORR( rval, "can't set bottom edge sense to the lateral face" );
3488 
3489  int Tsense;
3490  rval = getEgVtxSense( tEdge, v[3], v[2], Tsense );
3491  MBERRORR( rval, "can't get top edge sense" );
3492  // we need to see what sense has topEdge in face
3493  rval = _my_geomTopoTool->set_sense( tEdge, newLatFace, Tsense );
3494  MBERRORR( rval, "can't set top edge sense to the lateral face" );
3495 
3496  rval = _my_geomTopoTool->set_sense( e1, newLatFace, -1 );
3497  MBERRORR( rval, "can't set first vert edge sense" );
3498 
3499  rval = _my_geomTopoTool->set_sense( e2, newLatFace, 1 );
3500  MBERRORR( rval, "can't set second edge sense to the lateral face" );
3501  // first, create edges along direction, for the
3502 
3503  int indexB = 0, indexT = 0; // indices of the current nodes in the weaving process
3504  // weaving is either up or down; the triangles are oriented positively either way
3505  // up is nodes1[indexB], nodes2[indexT+1], nodes2[indexT]
3506  // down is nodes1[indexB], nodes1[indexB+1], nodes2[indexT]
3507  // the decision to weave up or down is based on which one is closer to the direction normal
3508  /*
3509  *
3510  * --------*------*-----------* ^
3511  * / . \ . . ------> dir1 | up
3512  * /. \ . . |
3513  * -----*------------*----*
3514  *
3515  */
3516  // we have to change the logic to account for cases when the curve in xy plane is not straight
3517  // (for example, when merging curves with a min_dot < -0.5, which means that the edges
3518  // could be even closed (periodic), with one vertex
3519  // the top and bottom curves should follow the same path in the "xy" plane (better to say
3520  // the plane perpendicular to the direction of weaving)
3521  // in this logic, the vector dir1 varies along the curve !!!
3522  CartVect dir1 = coords1[1] - coords1[0]; // we should have at least 2 nodes, N1>=2!!
3523 
3524  CartVect planeNormal = dir1 * up;
3525  dir1 = up * planeNormal;
3526  dir1.normalize();
3527  // this direction will be updated constantly with the direction of last edge added
3528  bool weaveDown = true;
3529 
3530  CartVect startP = coords1[0]; // used for origin of comparisons
3531  while( 1 )
3532  {
3533  if( ( indexB == N1 - 1 ) && ( indexT == N2 - 1 ) ) break; // we cannot advance anymore
3534  if( indexB == N1 - 1 )
3535  {
3536  weaveDown = false;
3537  }
3538  else if( indexT == N2 - 1 )
3539  {
3540  weaveDown = true;
3541  }
3542  else
3543  {
3544  // none are at the end, we have to decide which way to go, based on which index + 1 is
3545  // closer
3546  double proj1 = ( coords1[indexB + 1] - startP ) % dir1;
3547  double proj2 = ( coords2[indexT + 1] - startP ) % dir1;
3548  if( proj1 < proj2 )
3549  weaveDown = true;
3550  else
3551  weaveDown = false;
3552  }
3553  EntityHandle nTri[3] = { nodes1[indexB], 0, nodes2[indexT] };
3554  if( weaveDown )
3555  {
3556  nTri[1] = nodes1[indexB + 1];
3557  nTri[2] = nodes2[indexT];
3558  indexB++;
3559  }
3560  else
3561  {
3562  nTri[1] = nodes2[indexT + 1];
3563  indexT++;
3564  }
3565  EntityHandle triangle;
3566  rval = _mbImpl->create_element( MBTRI, nTri, 3, triangle );
3567  MBERRORR( rval, "can't create triangle" );
3568 
3569  rval = _mbImpl->add_entities( newLatFace, &triangle, 1 );
3570  MBERRORR( rval, "can't add triangle to face set" );
3571  if( weaveDown )
3572  {
3573  // increase was from nodes1[indexB-1] to nodes1[indexb]
3574  dir1 = coords1[indexB] - coords1[indexB - 1]; // we should have at least 2 nodes, N1>=2!!
3575  }
3576  else
3577  {
3578  dir1 = coords2[indexT] - coords2[indexT - 1];
3579  }
3580  dir1 = up * ( dir1 * up );
3581  dir1.normalize();
3582  }
3583  // we do not check yet if the triangles are inverted
3584  // if yes, we should correct them. HOW?
3585  // we probably need a better meshing strategy, one that can overcome really bad meshes.
3586  // again, these faces are not what we should use for geometry, maybe we should use the
3587  // extruded quads, identified AFTER hexa are created.
3588  // right now, I see only a cosmetic use of these lateral faces
3589  // the whole idea of volume maybe is overrated
3590  // volume should be just quads extruded from bottom ?
3591  //
3592  return MB_SUCCESS;
3593 }

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().

Member Data Documentation

◆ _brokenEdges

std::map< EntityHandle, EntityHandle > moab::FBEngine::_brokenEdges
private

Definition at line 315 of file FBEngine.hpp.

Referenced by split_boundary(), and split_surface_with_direction().

◆ _edges

std::map< EntityHandle, SmoothCurve* > moab::FBEngine::_edges
private

◆ _faces

std::map< EntityHandle, SmoothFace* > moab::FBEngine::_faces
private

◆ _initialized

bool moab::FBEngine::_initialized
private

Definition at line 298 of file FBEngine.hpp.

Referenced by Init().

◆ _mbImpl

◆ _my_geomTopoTool

◆ _my_gsets

Range moab::FBEngine::_my_gsets[5]
private

◆ _newTriangles

Range moab::FBEngine::_newTriangles
private

Definition at line 313 of file FBEngine.hpp.

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

◆ _piercedEdges

Range moab::FBEngine::_piercedEdges
private

Definition at line 314 of file FBEngine.hpp.

Referenced by create_new_gedge(), split_surface(), and split_surface_with_direction().

◆ _piercedTriangles

Range moab::FBEngine::_piercedTriangles
private

◆ _smooth

bool moab::FBEngine::_smooth
private

◆ _smthCurve

SmoothCurve** moab::FBEngine::_smthCurve
private

Definition at line 310 of file FBEngine.hpp.

Referenced by clean(), and initializeSmoothing().

◆ _smthFace

SmoothFace** moab::FBEngine::_smthFace
private

Definition at line 309 of file FBEngine.hpp.

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

◆ _t_created

bool moab::FBEngine::_t_created
private

Definition at line 296 of file FBEngine.hpp.

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


The documentation for this class was generated from the following files: