moab::SmoothFace Class Reference

Implement CAMAL geometry callbacks using smooth iMesh. More...

#include <SmoothFace.hpp>

Collaboration diagram for moab::SmoothFace:

Public Member Functions
	SmoothFace (Interface *mb, EntityHandle surface_set, GeomTopoTool *gTool)
virtual	~SmoothFace ()
virtual double	area ()
virtual void	bounding_box (double box_min[3], double box_max[3])
virtual void	move_to_surface (double &x, double &y, double &z)
virtual bool	normal_at (double x, double y, double z, double &nx, double &ny, double &nz)
int	init_gradient ()
ErrorCode	evaluate_smooth_edge (EntityHandle eh, double &t, CartVect &outv)
ErrorCode	eval_bezier_patch (EntityHandle tri, CartVect &areacoord, CartVect &pt)
void	project_to_facet_plane (EntityHandle tri, CartVect &pt, CartVect &point_on_plane, double &dist_to_plane)
void	facet_area_coordinate (EntityHandle facet, CartVect &pt_on_plane, CartVect &areacoord)
ErrorCode	project_to_facets_main (CartVect &this_point, bool trim, bool &outside, CartVect *closest_point_ptr=NULL, CartVect *normal_ptr=NULL)
ErrorCode	project_to_facets (std::vector< EntityHandle > &facet_list, EntityHandle &lastFacet, int interpOrder, double compareTol, CartVect &this_point, bool trim, bool &outside, CartVect *closest_point_ptr, CartVect *normal_ptr)
ErrorCode	project_to_facet (EntityHandle facet, CartVect &pt, CartVect &areacoord, CartVect &close_point, bool &outside_facet, double compare_tol)
bool	is_at_vertex (EntityHandle facet, CartVect &pt, CartVect &ac, double compare_tol, CartVect &eval_pt, CartVect *eval_norm_ptr)
ErrorCode	project_to_patch (EntityHandle facet, CartVect &ac, CartVect &pt, CartVect &eval_pt, CartVect *eval_norm, bool &outside, double compare_tol, int edge_id)
ErrorCode	eval_bezier_patch_normal (EntityHandle facet, CartVect &areacoord, CartVect &normal)
ErrorCode	compute_tangents_for_each_edge ()
ErrorCode	get_normals_for_vertices (const EntityHandle *conn2, CartVect N[2])
ErrorCode	init_edge_control_points (CartVect &P0, CartVect &P3, CartVect &N0, CartVect &N3, CartVect &T0, CartVect &T3, CartVect *Pi)
ErrorCode	compute_control_points_on_edges (double min_dot, Tag edgeCtrlTag, Tag markTag)
ErrorCode	compute_internal_control_points_on_facets (double min_dot, Tag facetCtrlTag, Tag facetEdgeCtrlTag)
void	DumpModelControlPoints ()
int	eval_counter ()
ErrorCode	ray_intersection_correct (EntityHandle facet, CartVect &pt, CartVect &ray, CartVect &eval_pt, double &distance, bool &outside)
void	append_smooth_tags (std::vector< Tag > &smoothTags)

Private Member Functions
bool	move_ac_inside (CartVect &ac, double tol)
void	ac_at_edge (CartVect &fac, CartVect &eac, int edge_id)
ErrorCode	init_bezier_edge (EntityHandle edge, double min_dot)
ErrorCode	find_edges_orientations (EntityHandle edges[3], const EntityHandle *conn3, int orient[3])
ErrorCode	init_facet_control_points (CartVect N[6], CartVect P[3][5], CartVect G[6])
void	adjust_bounding_box (CartVect &vect)

Private Attributes
CartVect	_minim
CartVect	_maxim
Range	_triangles
Range	_edges
Range	_nodes
Tag	_markTag
Tag	_gradientTag
Tag	_tangentsTag
Tag	_edgeCtrlTag
Tag	_facetCtrlTag
Tag	_facetEdgeCtrlTag
Tag	_planeTag
Interface *	_mb
EntityHandle	_set
GeomTopoTool *	_my_geomTopoTool
EntityHandle	_obb_root
long	_evaluationsCounter

Detailed Description

Implement CAMAL geometry callbacks using smooth iMesh.

Definition at line 29 of file SmoothFace.hpp.

Constructor & Destructor Documentation

SmoothFace()

moab::SmoothFace::SmoothFace	(	Interface *	mb,
		EntityHandle	surface_set,
		GeomTopoTool *	gTool
	)

Definition at line 41 of file SmoothFace.cpp.

 : _markTag( 0 ), _gradientTag( 0 ), _tangentsTag( 0 ), _edgeCtrlTag( 0 ), _facetCtrlTag( 0 ),
 _facetEdgeCtrlTag( 0 ), _planeTag( 0 ), _mb( mb ), _set( surface_set ), _my_geomTopoTool( gTool ), _obb_root( 0 ),
 _evaluationsCounter( 0 )
{
 //_smooth_face = NULL;
 //_mbOut->create_meshset(MESHSET_SET, _oSet); //will contain the
 // get also the obb_root
 if( _my_geomTopoTool )
 {
 _my_geomTopoTool->get_root( this->_set, _obb_root );
 if( debug_surf_eval1 ) _my_geomTopoTool->obb_tree()->stats( _obb_root, std::cout );
 }
}

References _my_geomTopoTool, _obb_root, _set, moab::GeomTopoTool::get_root(), moab::GeomTopoTool::obb_tree(), and moab::OrientedBoxTreeTool::stats().

~SmoothFace()

moab::SmoothFace::~SmoothFace ( )

virtual

Definition at line 56 of file SmoothFace.cpp.

{}

Member Function Documentation

ac_at_edge()

void moab::SmoothFace::ac_at_edge ( CartVect &  fac, CartVect &  eac, int  edge_id  )

private

Definition at line 1504 of file SmoothFace.cpp.

{
 double u, v, w;
 switch( edge_id )
 {
 case 0:
 u = 0.0;
 v = fac[1] / ( fac[1] + fac[2] ); // v = fac.y() / (fac.y() + fac.z());
 w = 1.0 - v;
 break;
 case 1:
 u = fac[0] / ( fac[0] + fac[2] ); // u = fac.x() / (fac.x() + fac.z());
 v = 0.0;
 w = 1.0 - u;
 break;
 case 2:
 u = fac[0] / ( fac[0] + fac[1] ); // u = fac.x() / (fac.x() + fac.y());
 v = 1.0 - u;
 w = 0.0;
 break;
 default:
 assert( 0 );
 u = -1; // needed to eliminate warnings about used before set
 v = -1; // needed to eliminate warnings about used before set
 w = -1; // needed to eliminate warnings about used before set
 break;
 }
 eac[0] = u;
 eac[1] = v;
 eac[2] = w; //= CartVect(u, v, w);
}

Referenced by project_to_patch().

adjust_bounding_box()

void moab::SmoothFace::adjust_bounding_box ( CartVect &  vect )

private

Definition at line 565 of file SmoothFace.cpp.

{
 // _minim, _maxim
 for( int j = 0; j < 3; j++ )
 {
 if( _minim[j] > vect[j] ) _minim[j] = vect[j];
 if( _maxim[j] < vect[j] ) _maxim[j] = vect[j];
 }
}

References _maxim, and _minim.

Referenced by compute_internal_control_points_on_facets().

append_smooth_tags()

void moab::SmoothFace::append_smooth_tags

(

std::vector< Tag > &

smoothTags

)

Definition at line 87 of file SmoothFace.cpp.

{
 // these are created locally, for each smooth face
 smoothTags.push_back( _gradientTag );
 smoothTags.push_back( _planeTag );
}

References _gradientTag, and _planeTag.

Referenced by moab::FBEngine::delete_smooth_tags().

area()

double moab::SmoothFace::area ( )

virtual

Definition at line 58 of file SmoothFace.cpp.

{
 // find the area of this entity
 // assert(_smooth_face);
 // double area1 = _smooth_face->area();
 double totArea = 0.;
 for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
 {
 EntityHandle tria = *it;
 const EntityHandle* conn3;
 int nnodes;
 _mb->get_connectivity( tria, conn3, nnodes );
 //
 // double coords[9]; // store the coordinates for the nodes
 //_mb->get_coords(conn3, 3, coords);
 CartVect p[3];
 _mb->get_coords( conn3, 3, (double*)&p[0] );
 // need to compute the angles
 // compute angles and the normal
 // CartVect p1(&coords[0]), p2(&coords[3]), p3(&coords[6]);
 
 CartVect AB( p[1] - p[0] ); //(p2 - p1);
 CartVect BC( p[2] - p[1] ); //(p3 - p2);
 CartVect normal = AB * BC;
 totArea += normal.length() / 2;
 }
 return totArea;
}

References _mb, _triangles, moab::Range::begin(), moab::Range::end(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), and moab::CartVect::length().

bounding_box()

void moab::SmoothFace::bounding_box ( double  box_min[3], double  box_max[3]  )

virtual

Definition at line 93 of file SmoothFace.cpp.

{
 
 for( int i = 0; i < 3; i++ )
 {
 box_min[i] = _minim[i];
 box_max[i] = _maxim[i];
 }
 // _minim, _maxim
}

References _maxim, _minim, box_max(), and box_min().

compute_control_points_on_edges()

ErrorCode moab::SmoothFace::compute_control_points_on_edges	(	double	min_dot,
		Tag	edgeCtrlTag,
		Tag	markTag
	)

Definition at line 147 of file SmoothFace.cpp.

{
 
 _edgeCtrlTag = edgeCtrlTag;
 _markTag = markTag;
 
 // now, compute control points for all edges that are not marked already (they are no on the
 // boundary!)
 for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
 {
 EntityHandle edg = *it;
 // is the edge marked? already computed
 unsigned char tagVal = 0;
 _mb->tag_get_data( _markTag, &edg, 1, &tagVal );
 if( tagVal ) continue;
 // double min_dot;
 init_bezier_edge( edg, min_dot );
 tagVal = 1;
 _mb->tag_set_data( _markTag, &edg, 1, &tagVal );
 }
 return MB_SUCCESS;
}

References _edgeCtrlTag, _edges, _markTag, _mb, moab::Range::begin(), moab::Range::end(), init_bezier_edge(), MB_SUCCESS, moab::Interface::tag_get_data(), and moab::Interface::tag_set_data().

Referenced by moab::FBEngine::initializeSmoothing().

compute_internal_control_points_on_facets()

ErrorCode moab::SmoothFace::compute_internal_control_points_on_facets	(	double	min_dot,
		Tag	facetCtrlTag,
		Tag	facetEdgeCtrlTag
	)

Definition at line 470 of file SmoothFace.cpp.

{
 // collect from each triangle the control points in order
 //
 
 _facetCtrlTag = facetCtrlTag;
 _facetEdgeCtrlTag = facetEdgeCtrlTag;
 
 for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
 {
 EntityHandle tri = *it;
 // first get connectivity, and the edges
 // we need a fast method to retrieve the adjacent edges to each triangle
 const EntityHandle* conn3;
 int nnodes;
 ErrorCode rval = _mb->get_connectivity( tri, conn3, nnodes );
 assert( MB_SUCCESS == rval );
 if( MB_SUCCESS != rval ) return rval;
 assert( 3 == nnodes );
 
 // would it be easier to do
 CartVect vNode[3]; // position at nodes
 rval = _mb->get_coords( conn3, 3, (double*)&vNode[0] );
 assert( MB_SUCCESS == rval );
 if( MB_SUCCESS != rval ) return rval;
 
 // get gradients (normal) at each node of triangle
 CartVect NN[3];
 rval = _mb->tag_get_data( _gradientTag, conn3, 3, &NN[0] );
 assert( MB_SUCCESS == rval );
 if( MB_SUCCESS != rval ) return rval;
 
 EntityHandle edges[3];
 int orient[3]; // + 1 or -1, if the edge is positive or negative within the face
 rval = find_edges_orientations( edges, conn3, orient ); // maybe we will set it?
 assert( MB_SUCCESS == rval );
 if( MB_SUCCESS != rval ) return rval;
 // maybe we will store some tags with edges and their orientation with respect to
 // a triangle;
 CartVect P[3][5];
 CartVect N[6], G[6];
 // create the linear array for control points on edges, for storage (expensive !!!)
 CartVect CP[9];
 int index = 0;
 // maybe store a tag / entity handle for edges?
 for( int i = 0; i < 3; i++ )
 {
 // populate P and N with the right vectors
 int i1 = ( i + 1 ) % 3; // the first node of the edge
 int i2 = ( i + 2 ) % 3; // the second node of the edge
 N[2 * i] = NN[i1];
 N[2 * i + 1] = NN[i2];
 P[i][0] = vNode[i1];
 rval = _mb->tag_get_data( _edgeCtrlTag, &edges[i], 1, &( P[i][1] ) );
 // if sense is -1, swap 1 and 3 control points
 if( orient[i] == -1 )
 {
 CartVect tmp;
 tmp = P[i][1];
 P[i][1] = P[i][3];
 P[i][3] = tmp;
 }
 P[i][4] = vNode[i2];
 for( int j = 1; j < 4; j++ )
 CP[index++] = P[i][j];
 
 // the first edge control points
 }
 
 // stat = facet->get_edge_control_points( P );
 init_facet_control_points( N, P, G );
 // what do we need to store in the tag control points?
 rval = _mb->tag_set_data( _facetCtrlTag, &tri, 1, &G[0] );
 assert( MB_SUCCESS == rval );
 if( MB_SUCCESS != rval ) return rval;
 
 // store here again the 9 control points on the edges
 rval = _mb->tag_set_data( _facetEdgeCtrlTag, &tri, 1, &CP[0] );
 assert( MB_SUCCESS == rval );
 if( MB_SUCCESS != rval ) return rval;
 // look at what we retrieve later
 
 // adjust the bounding box
 int j = 0;
 for( j = 0; j < 3; j++ )
 adjust_bounding_box( vNode[j] );
 // edge control points
 for( j = 0; j < 9; j++ )
 adjust_bounding_box( CP[j] );
 // internal facet control points
 for( j = 0; j < 6; j++ )
 adjust_bounding_box( G[j] );
 }
 return MB_SUCCESS;
}

References _edgeCtrlTag, _facetCtrlTag, _facetEdgeCtrlTag, _gradientTag, _mb, _triangles, adjust_bounding_box(), moab::Range::begin(), moab::Range::end(), ErrorCode, find_edges_orientations(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), init_facet_control_points(), MB_SUCCESS, moab::Interface::tag_get_data(), and moab::Interface::tag_set_data().

Referenced by moab::FBEngine::initializeSmoothing().

compute_tangents_for_each_edge()

ErrorCode moab::SmoothFace::compute_tangents_for_each_edge

(

)

Definition at line 362 of file SmoothFace.cpp.

{
 double defTangents[6] = { 0., 0., 0., 0., 0., 0. };
 ErrorCode rval =
 _mb->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, _tangentsTag, MB_TAG_DENSE | MB_TAG_CREAT, &defTangents );
 if( MB_SUCCESS != rval ) return MB_FAILURE;
 
 // now, compute Tangents for all edges that are not on boundary, so they are not marked
 for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
 {
 EntityHandle edg = *it;
 
 int nnodes;
 const EntityHandle* conn2; //
 _mb->get_connectivity( edg, conn2, nnodes );
 assert( nnodes == 2 );
 CartVect P[2]; // store the coordinates for the nodes
 rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
 if( MB_SUCCESS != rval ) return rval;
 assert( rval == MB_SUCCESS );
 CartVect T[2];
 T[0] = P[1] - P[0];
 T[0].normalize();
 T[1] = T[0]; //
 _mb->tag_set_data( _tangentsTag, &edg, 1,
 (double*)&T[0] ); // set the tangents computed at every edge
 }
 return MB_SUCCESS;
}

References _edges, _mb, _tangentsTag, moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, T, moab::Interface::tag_get_handle(), and moab::Interface::tag_set_data().

Referenced by moab::FBEngine::initializeSmoothing().

DumpModelControlPoints()

void moab::SmoothFace::DumpModelControlPoints

(

)

Definition at line 622 of file SmoothFace.cpp.

{
 // here, we will dump all control points from edges and facets (6 control points for each facet)
 // we may also create some edges; maybe later...
 // create a point3D file
 // output a Point3D file (special visit format)
 unsigned long setId = _mb->id_from_handle( _set );
 char name[50] = { 0 };
 sprintf( name, "%lucontrol.Point3D", setId ); // name should be something 2control.Point3D
 std::ofstream point3DFile;
 point3DFile.open( name ); //("control.Point3D");
 point3DFile << "# x y z \n";
 std::ofstream point3DEdgeFile;
 sprintf( name, "%lucontrolEdge.Point3D", setId ); //
 point3DEdgeFile.open( name ); //("controlEdge.Point3D");
 point3DEdgeFile << "# x y z \n";
 std::ofstream smoothPoints;
 sprintf( name, "%lusmooth.Point3D", setId ); //
 smoothPoints.open( name ); //("smooth.Point3D");
 smoothPoints << "# x y z \n";
 CartVect controlPoints[3]; // edge control points
 for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
 {
 EntityHandle edge = *it;
 
 _mb->tag_get_data( _edgeCtrlTag, &edge, 1, (double*)&controlPoints[0] );
 for( int i = 0; i < 3; i++ )
 {
 CartVect& c = controlPoints[i];
 point3DEdgeFile << std::setprecision( 11 ) << c[0] << " " << c[1] << " " << c[2] << " \n";
 }
 }
 CartVect controlTriPoints[6]; // triangle control points
 CartVect P_facet[3]; // result in 3 "mid" control points
 for( Range::iterator it2 = _triangles.begin(); it2 != _triangles.end(); ++it2 )
 {
 EntityHandle tri = *it2;
 
 _mb->tag_get_data( _facetCtrlTag, &tri, 1, (double*)&controlTriPoints[0] );
 
 // draw a line of points between pairs of control points
 int numPoints = 7;
 for( int n = 0; n < numPoints; n++ )
 {
 double a = 1. * n / ( numPoints - 1 );
 double b = 1.0 - a;
 
 P_facet[0] = a * controlTriPoints[3] + b * controlTriPoints[4];
 // 1,2,1
 P_facet[1] = a * controlTriPoints[0] + b * controlTriPoints[5];
 // 1,1,2
 P_facet[2] = a * controlTriPoints[1] + b * controlTriPoints[2];
 for( int i = 0; i < 3; i++ )
 {
 CartVect& c = P_facet[i];
 point3DFile << std::setprecision( 11 ) << c[0] << " " << c[1] << " " << c[2] << " \n";
 }
 }
 
 // evaluate for each triangle a lattice of points
 int N = 40;
 for( int k = 0; k <= N; k++ )
 {
 for( int m = 0; m <= N - k; m++ )
 {
 int n = N - m - k;
 CartVect areacoord( 1. * k / N, 1. * m / N, 1. * n / N );
 CartVect pt;
 eval_bezier_patch( tri, areacoord, pt );
 smoothPoints << std::setprecision( 11 ) << pt[0] << " " << pt[1] << " " << pt[2] << " \n";
 }
 }
 }
 point3DFile.close();
 smoothPoints.close();
 point3DEdgeFile.close();
 return;
}

References _edgeCtrlTag, _edges, _facetCtrlTag, _mb, _set, _triangles, moab::Range::begin(), moab::Range::end(), eval_bezier_patch(), moab::Interface::id_from_handle(), and moab::Interface::tag_get_data().

eval_bezier_patch()

ErrorCode moab::SmoothFace::eval_bezier_patch	(	EntityHandle	tri,
		CartVect &	areacoord,
		CartVect &	pt
	)

Definition at line 748 of file SmoothFace.cpp.

{
 //
 // interpolate internal control points
 
 CartVect gctrl_pts[6];
 // get the control points facet->get_control_points( gctrl_pts );
 // init_facet_control_points( N, P, G) ;
 // what do we need to store in the tag control points?
 ErrorCode rval = _mb->tag_get_data( _facetCtrlTag, &tri, 1, &gctrl_pts[0] ); // get all 6 control points
 assert( MB_SUCCESS == rval );
 if( MB_SUCCESS != rval ) return rval;
 const EntityHandle* conn3 = NULL;
 int nnodes = 0;
 rval = _mb->get_connectivity( tri, conn3, nnodes );
 assert( MB_SUCCESS == rval );
 
 CartVect vN[3];
 _mb->get_coords( conn3, 3, (double*)&vN[0] ); // fill the coordinates of the vertices
 
 if( fabs( areacoord[1] + areacoord[2] ) < 1.0e-6 )
 {
 pt = vN[0];
 return MB_SUCCESS;
 }
 if( fabs( areacoord[0] + areacoord[2] ) < 1.0e-6 )
 {
 pt = vN[0];
 return MB_SUCCESS;
 }
 if( fabs( areacoord[0] + areacoord[1] ) < 1.0e-6 )
 {
 pt = vN[0];
 return MB_SUCCESS;
 }
 
 CartVect P_facet[3];
 // 2,1,1
 P_facet[0] =
 ( 1.0e0 / ( areacoord[1] + areacoord[2] ) ) * ( areacoord[1] * gctrl_pts[3] + areacoord[2] * gctrl_pts[4] );
 // 1,2,1
 P_facet[1] =
 ( 1.0e0 / ( areacoord[0] + areacoord[2] ) ) * ( areacoord[0] * gctrl_pts[0] + areacoord[2] * gctrl_pts[5] );
 // 1,1,2
 P_facet[2] =
 ( 1.0e0 / ( areacoord[0] + areacoord[1] ) ) * ( areacoord[0] * gctrl_pts[1] + areacoord[1] * gctrl_pts[2] );
 
 // sum the contribution from each of the control points
 
 pt = CartVect( 0. ); // set all to 0, we start adding / accumulating different parts
 // first edge is from node 0 to 1, index 2 in
 
 // retrieve the points, in order, and the control points on edges
 
 // store here again the 9 control points on the edges
 CartVect CP[9];
 rval = _mb->tag_get_data( _facetEdgeCtrlTag, &tri, 1, &CP[0] );
 assert( MB_SUCCESS == rval );
 
 // CubitFacetEdge *edge;
 // edge = facet->edge(2);! start with edge 2, from 0-1
 int k = 0;
 CartVect ctrl_pts[5];
 // edge->control_points(facet, ctrl_pts);
 ctrl_pts[0] = vN[0]; //
 for( k = 1; k < 4; k++ )
 ctrl_pts[k] = CP[k + 5]; // for edge index 2
 ctrl_pts[4] = vN[1]; //
 
 // i=4; j=0; k=0;
 double B = mbquart( areacoord[0] );
 pt += B * ctrl_pts[0];
 
 // i=3; j=1; k=0;
 B = 4.0 * mbcube( areacoord[0] ) * areacoord[1];
 pt += B * ctrl_pts[1];
 
 // i=2; j=2; k=0;
 B = 6.0 * mbsqr( areacoord[0] ) * mbsqr( areacoord[1] );
 pt += B * ctrl_pts[2];
 
 // i=1; j=3; k=0;
 B = 4.0 * areacoord[0] * mbcube( areacoord[1] );
 pt += B * ctrl_pts[3];
 
 // edge = facet->edge(0);
 // edge->control_points(facet, ctrl_pts);
 // edge index 0, from 1 to 2
 ctrl_pts[0] = vN[1]; //
 for( k = 1; k < 4; k++ )
 ctrl_pts[k] = CP[k - 1]; // for edge index 0
 ctrl_pts[4] = vN[2]; //
 
 // i=0; j=4; k=0;
 B = mbquart( areacoord[1] );
 pt += B * ctrl_pts[0];
 
 // i=0; j=3; k=1;
 B = 4.0 * mbcube( areacoord[1] ) * areacoord[2];
 pt += B * ctrl_pts[1];
 
 // i=0; j=2; k=2;
 B = 6.0 * mbsqr( areacoord[1] ) * mbsqr( areacoord[2] );
 pt += B * ctrl_pts[2];
 
 // i=0; j=1; k=3;
 B = 4.0 * areacoord[1] * mbcube( areacoord[2] );
 pt += B * ctrl_pts[3];
 
 // edge = facet->edge(1);
 // edge->control_points(facet, ctrl_pts);
 // edge index 1, from 2 to 0
 ctrl_pts[0] = vN[2]; //
 for( k = 1; k < 4; k++ )
 ctrl_pts[k] = CP[k + 2]; // for edge index 0
 ctrl_pts[4] = vN[0]; //
 
 // i=0; j=0; k=4;
 B = mbquart( areacoord[2] );
 pt += B * ctrl_pts[0];
 
 // i=1; j=0; k=3;
 B = 4.0 * areacoord[0] * mbcube( areacoord[2] );
 pt += B * ctrl_pts[1];
 
 // i=2; j=0; k=2;
 B = 6.0 * mbsqr( areacoord[0] ) * mbsqr( areacoord[2] );
 pt += B * ctrl_pts[2];
 
 // i=3; j=0; k=1;
 B = 4.0 * mbcube( areacoord[0] ) * areacoord[2];
 pt += B * ctrl_pts[3];
 
 // i=2; j=1; k=1;
 B = 12.0 * mbsqr( areacoord[0] ) * areacoord[1] * areacoord[2];
 pt += B * P_facet[0];
 
 // i=1; j=2; k=1;
 B = 12.0 * areacoord[0] * mbsqr( areacoord[1] ) * areacoord[2];
 pt += B * P_facet[1];
 
 // i=1; j=1; k=2;
 B = 12.0 * areacoord[0] * areacoord[1] * mbsqr( areacoord[2] );
 pt += B * P_facet[2];
 
 return MB_SUCCESS;
}

References _facetCtrlTag, _facetEdgeCtrlTag, _mb, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), MB_SUCCESS, mbcube, mbquart, mbsqr, and moab::Interface::tag_get_data().

Referenced by DumpModelControlPoints(), and project_to_patch().

eval_bezier_patch_normal()

ErrorCode moab::SmoothFace::eval_bezier_patch_normal	(	EntityHandle	facet,
		CartVect &	areacoord,
		CartVect &	normal
	)

Definition at line 1710 of file SmoothFace.cpp.

{
 // interpolate internal control points
 
 CartVect gctrl_pts[6];
 // facet->get_control_points( gctrl_pts );
 ErrorCode rval = _mb->tag_get_data( _facetCtrlTag, &facet, 1, &gctrl_pts[0] );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 // _gradientTag
 // get normals at points
 const EntityHandle* conn3 = NULL;
 int nnodes = 0;
 rval = _mb->get_connectivity( facet, conn3, nnodes );
 if( MB_SUCCESS != rval ) return rval;
 
 CartVect NN[3];
 rval = _mb->tag_get_data( _gradientTag, conn3, 3, &NN[0] );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 if( fabs( areacoord[1] + areacoord[2] ) < 1.0e-6 )
 {
 normal = NN[0];
 return MB_SUCCESS;
 }
 if( fabs( areacoord[0] + areacoord[2] ) < 1.0e-6 )
 {
 normal = NN[1]; // facet->point(1)->normal(facet);
 return MB_SUCCESS;
 }
 if( fabs( areacoord[0] + areacoord[1] ) < 1.0e-6 )
 {
 normal = NN[2]; // facet->point(2)->normal(facet);
 return MB_SUCCESS;
 }
 
 // compute the hodograph of the quartic Gregory patch
 
 CartVect Nijk[10];
 // hodograph(facet,areacoord,Nijk);
 // start copy from hodograph
 // CubitVector gctrl_pts[6];
 // facet->get_control_points( gctrl_pts );
 CartVect P_facet[3];
 
 // 2,1,1
 /*P_facet[0] = (1.0e0 / (areacoord.y() + areacoord.z())) *
 (areacoord.y() * gctrl_pts[3] +
 areacoord.z() * gctrl_pts[4]);*/
 P_facet[0] =
 ( 1.0e0 / ( areacoord[1] + areacoord[2] ) ) * ( areacoord[1] * gctrl_pts[3] + areacoord[2] * gctrl_pts[4] );
 // 1,2,1
 /*P_facet[1] = (1.0e0 / (areacoord.x() + areacoord.z())) *
 (areacoord.x() * gctrl_pts[0] +
 areacoord.z() * gctrl_pts[5]);*/
 P_facet[1] =
 ( 1.0e0 / ( areacoord[0] + areacoord[2] ) ) * ( areacoord[0] * gctrl_pts[0] + areacoord[2] * gctrl_pts[5] );
 // 1,1,2
 /*P_facet[2] = (1.0e0 / (areacoord.x() + areacoord.y())) *
 (areacoord.x() * gctrl_pts[1] +
 areacoord.y() * gctrl_pts[2]);*/
 P_facet[2] =
 ( 1.0e0 / ( areacoord[0] + areacoord[1] ) ) * ( areacoord[0] * gctrl_pts[1] + areacoord[1] * gctrl_pts[2] );
 
 // corner control points are just the normals at the points
 
 // 3, 0, 0
 Nijk[0] = NN[0];
 // 0, 3, 0
 Nijk[3] = NN[1];
 // 0, 0, 3
 Nijk[9] = NN[2]; // facet->point(2)->normal(facet);
 
 // fill in the boundary control points. Define as the normal to the local
 // triangle formed by the quartic control point lattice
 
 // store here again the 9 control points on the edges
 CartVect CP[9]; // 9 control points on the edges,
 rval = _mb->tag_get_data( _facetEdgeCtrlTag, &facet, 1, &CP[0] );
 if( MB_SUCCESS != rval ) return rval;
 // there are 3 CP for each edge, 0, 1, 2; first edge is 1-2
 // CubitFacetEdge *edge;
 // edge = facet->edge( 2 );
 // CubitVector ctrl_pts[5];
 // edge->control_points(facet, ctrl_pts);
 
 // 2, 1, 0
 // Nijk[1] = (ctrl_pts[2] - ctrl_pts[1]) * (P_facet[0] - ctrl_pts[1]);
 Nijk[1] = ( CP[7] - CP[6] ) * ( P_facet[0] - CP[6] );
 Nijk[1].normalize();
 
 // 1, 2, 0
 // Nijk[2] = (ctrl_pts[3] - ctrl_pts[2]) * (P_facet[1] - ctrl_pts[2]);
 Nijk[2] = ( CP[8] - CP[7] ) * ( P_facet[1] - CP[7] );
 Nijk[2].normalize();
 
 // edge = facet->edge( 0 );
 // edge->control_points(facet, ctrl_pts);
 
 // 0, 2, 1
 // Nijk[6] = (ctrl_pts[1] - P_facet[1]) * (ctrl_pts[2] - P_facet[1]);
 Nijk[6] = ( CP[0] - P_facet[1] ) * ( CP[1] - P_facet[1] );
 Nijk[6].normalize();
 
 // 0, 1, 2
 // Nijk[8] = (ctrl_pts[2] - P_facet[2]) * (ctrl_pts[3] - P_facet[2]);
 Nijk[8] = ( CP[1] - P_facet[2] ) * ( CP[2] - P_facet[2] );
 Nijk[8].normalize();
 
 // edge = facet->edge( 1 );
 // edge->control_points(facet, ctrl_pts);
 
 // 1, 0, 2
 // Nijk[7] = (P_facet[2] - ctrl_pts[2]) * (ctrl_pts[1] - ctrl_pts[2]);
 Nijk[7] = ( P_facet[2] - CP[4] ) * ( CP[3] - CP[4] );
 Nijk[7].normalize();
 
 // 2, 0, 1
 // Nijk[4] = (P_facet[0] - ctrl_pts[3]) * (ctrl_pts[2] - ctrl_pts[3]);
 Nijk[4] = ( P_facet[0] - CP[5] ) * ( CP[4] - CP[5] );
 Nijk[4].normalize();
 
 // 1, 1, 1
 Nijk[5] = ( P_facet[1] - P_facet[0] ) * ( P_facet[2] - P_facet[0] );
 Nijk[5].normalize();
 // end copy from hodograph
 
 // sum the contribution from each of the control points
 
 normal = CartVect( 0.0e0, 0.0e0, 0.0e0 );
 
 // i=3; j=0; k=0;
 // double Bsum = 0.0;
 double B = mbcube( areacoord[0] );
 // Bsum += B;
 normal += B * Nijk[0];
 
 // i=2; j=1; k=0;
 B = 3.0 * mbsqr( areacoord[0] ) * areacoord[1];
 // Bsum += B;
 normal += B * Nijk[1];
 
 // i=1; j=2; k=0;
 B = 3.0 * areacoord[0] * mbsqr( areacoord[1] );
 // Bsum += B;
 normal += B * Nijk[2];
 
 // i=0; j=3; k=0;
 B = mbcube( areacoord[1] );
 // Bsum += B;
 normal += B * Nijk[3];
 
 // i=2; j=0; k=1;
 B = 3.0 * mbsqr( areacoord[0] ) * areacoord[2];
 // Bsum += B;
 normal += B * Nijk[4];
 
 // i=1; j=1; k=1;
 B = 6.0 * areacoord[0] * areacoord[1] * areacoord[2];
 // Bsum += B;
 normal += B * Nijk[5];
 
 // i=0; j=2; k=1;
 B = 3.0 * mbsqr( areacoord[1] ) * areacoord[2];
 // Bsum += B;
 normal += B * Nijk[6];
 
 // i=1; j=0; k=2;
 B = 3.0 * areacoord[0] * mbsqr( areacoord[2] );
 // Bsum += B;
 normal += B * Nijk[7];
 
 // i=0; j=1; k=2;
 B = 3.0 * areacoord[1] * mbsqr( areacoord[2] );
 // Bsum += B;
 normal += B * Nijk[8];
 
 // i=0; j=0; k=3;
 B = mbcube( areacoord[2] );
 // Bsum += B;
 normal += B * Nijk[9];
 
 // assert(fabs(Bsum - 1.0) < 1e-9);
 
 normal.normalize();
 
 return MB_SUCCESS;
}

References _facetCtrlTag, _facetEdgeCtrlTag, _gradientTag, _mb, ErrorCode, moab::Interface::get_connectivity(), MB_SUCCESS, mbcube, mbsqr, moab::CartVect::normalize(), and moab::Interface::tag_get_data().

Referenced by project_to_facets().

eval_counter()

int moab::SmoothFace::eval_counter ( )

inline

Definition at line 125 of file SmoothFace.hpp.

 {
 return _evaluationsCounter;
 }

References _evaluationsCounter.

evaluate_smooth_edge()

ErrorCode moab::SmoothFace::evaluate_smooth_edge	(	EntityHandle	eh,
		double &	t,
		CartVect &	outv
	)

Definition at line 708 of file SmoothFace.cpp.

{
 CartVect P[2]; // P0 and P1
 CartVect controlPoints[3]; // edge control points
 double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
 
 // project the position to the linear edge
 // t is from 0 to 1 only!!
 // double tt = (t + 1) * 0.5;
 if( tt <= 0.0 ) tt = 0.0;
 if( tt >= 1.0 ) tt = 1.0;
 
 int nnodes = 0;
 const EntityHandle* conn2 = NULL;
 ErrorCode rval = _mb->get_connectivity( eh, conn2, nnodes );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 rval = _mb->tag_get_data( _edgeCtrlTag, &eh, 1, (double*)&controlPoints[0] );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 t2 = tt * tt;
 t3 = t2 * tt;
 t4 = t3 * tt;
 one_minus_t = 1. - tt;
 one_minus_t2 = one_minus_t * one_minus_t;
 one_minus_t3 = one_minus_t2 * one_minus_t;
 one_minus_t4 = one_minus_t3 * one_minus_t;
 
 outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] + 6. * one_minus_t2 * t2 * controlPoints[1] +
 4. * one_minus_t * t3 * controlPoints[2] + t4 * P[1];
 
 return MB_SUCCESS;
}

References _edgeCtrlTag, _mb, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), MB_SUCCESS, and moab::Interface::tag_get_data().

facet_area_coordinate()

void moab::SmoothFace::facet_area_coordinate	(	EntityHandle	facet,
		CartVect &	pt_on_plane,
		CartVect &	areacoord
	)

Definition at line 929 of file SmoothFace.cpp.

{
 
 const EntityHandle* conn3 = NULL;
 int nnodes = 0;
 ErrorCode rval = _mb->get_connectivity( facet, conn3, nnodes );
 assert( MB_SUCCESS == rval );
 if( rval )
 {
 } // empty statement to prevent compiler warning
 
 // double coords[9]; // store the coordinates for the nodes
 //_mb->get_coords(conn3, 3, coords);
 CartVect p[3];
 rval = _mb->get_coords( conn3, 3, (double*)&p[0] );
 assert( MB_SUCCESS == rval );
 if( rval )
 {
 } // empty statement to prevent compiler warning
 double plane[4];
 
 rval = _mb->tag_get_data( _planeTag, &facet, 1, plane );
 assert( rval == MB_SUCCESS );
 if( rval )
 {
 } // empty statement to prevent compiler warning
 CartVect normal( &plane[0] ); // just first 3 components are used
 
 double area2;
 
 double tol = GEOMETRY_RESABS * 1.e-5; // 1.e-11;
 
 CartVect v1( p[1] - p[0] );
 CartVect v2( p[2] - p[0] );
 
 area2 = ( v1 * v2 ).length_squared(); // the same for CartVect
 if( area2 < 100 * tol )
 {
 tol = .01 * area2;
 }
 CartVect absnorm( fabs( normal[0] ), fabs( normal[1] ), fabs( normal[2] ) );
 
 // project to the closest coordinate plane so we only have to do this in 2D
 
 if( absnorm[0] >= absnorm[1] && absnorm[0] >= absnorm[2] )
 {
 area2 = determ3( p[0][1], p[0][2], p[1][1], p[1][2], p[2][1], p[2][2] );
 if( fabs( area2 ) < tol )
 {
 areacoord = CartVect( -std::numeric_limits< double >::min() ); // .set(
 // -std::numeric_limits<double>::min(),
 // -std::numeric_limits<double>::min(),
 // -std::numeric_limits<double>::min() );
 }
 else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
 {
 areacoord = CartVect( 1., 0., 0. ); //.set( 1.0, 0.0, 0.0 );
 }
 else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
 {
 areacoord = CartVect( 0., 1., 0. ); //.set( 0.0, 1.0, 0.0 );
 }
 else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
 {
 areacoord = CartVect( 0., 0., 1. ); //.set( 0.0, 0.0, 1.0 );
 }
 else
 {
 
 areacoord[0] = determ3( pt_on_plane[1], pt_on_plane[2], p[1][1], p[1][2], p[2][1], p[2][2] ) / area2;
 
 areacoord[1] = determ3( p[0][1], p[0][2], pt_on_plane[1], pt_on_plane[2], p[2][1], p[2][2] ) / area2;
 
 areacoord[2] = determ3( p[0][1], p[0][2], p[1][1], p[1][2], pt_on_plane[1], pt_on_plane[2] ) / area2;
 }
 }
 else if( absnorm[1] >= absnorm[0] && absnorm[1] >= absnorm[2] )
 {
 
 area2 = determ3( p[0][0], p[0][2], p[1][0], p[1][2], p[2][0], p[2][2] );
 if( fabs( area2 ) < tol )
 {
 areacoord = CartVect( -std::numeric_limits< double >::min() ); //.set(
 //-std::numeric_limits<double>::min(),
 //-std::numeric_limits<double>::min(),
 //-std::numeric_limits<double>::min() );
 }
 else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
 {
 areacoord = CartVect( 1., 0., 0. ); //.set( 1.0, 0.0, 0.0 );
 }
 else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
 {
 areacoord = CartVect( 0., 1., 0. ); //.set( 0.0, 1.0, 0.0 );
 }
 else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
 {
 areacoord = CartVect( 0., 0., 1. ); //.set( 0.0, 0.0, 1.0 );
 }
 else
 {
 
 areacoord[0] = determ3( pt_on_plane[0], pt_on_plane[2], p[1][0], p[1][2], p[2][0], p[2][2] ) / area2;
 
 areacoord[1] = determ3( p[0][0], p[0][2], pt_on_plane[0], pt_on_plane[2], p[2][0], p[2][2] ) / area2;
 
 areacoord[2] = determ3( p[0][0], p[0][2], p[1][0], p[1][2], pt_on_plane[0], pt_on_plane[2] ) / area2;
 }
 }
 else
 {
 /*area2 = determ3(pt0->x(), pt0->y(),
 pt1->x(), pt1->y(),
 pt2->x(), pt2->y());*/
 area2 = determ3( p[0][0], p[0][1], p[1][0], p[1][1], p[2][0], p[2][1] );
 if( fabs( area2 ) < tol )
 {
 areacoord = CartVect( -std::numeric_limits< double >::min() ); //.set(
 //-std::numeric_limits<double>::min(),
 //-std::numeric_limits<double>::min(),
 //-std::numeric_limits<double>::min() );
 }
 else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
 {
 areacoord = CartVect( 1., 0., 0. ); //.set( 1.0, 0.0, 0.0 );
 }
 else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
 {
 areacoord = CartVect( 0., 1., 0. ); //.set( 0.0, 1.0, 0.0 );
 }
 else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
 {
 areacoord = CartVect( 0., 0., 1. ); //.set( 0.0, 0.0, 1.0 );
 }
 else
 {
 
 areacoord[0] = determ3( pt_on_plane[0], pt_on_plane[1], p[1][0], p[1][1], p[2][0], p[2][1] ) / area2;
 
 areacoord[1] = determ3( p[0][0], p[0][1], pt_on_plane[0], pt_on_plane[1], p[2][0], p[2][1] ) / area2;
 
 areacoord[2] = determ3( p[0][0], p[0][1], p[1][0], p[1][1], pt_on_plane[0], pt_on_plane[1] ) / area2;
 }
 }
}

References _mb, _planeTag, determ3, ErrorCode, GEOMETRY_RESABS, moab::Interface::get_connectivity(), moab::Interface::get_coords(), length_squared(), MB_SUCCESS, moab::Interface::tag_get_data(), and moab::within_tolerance().

Referenced by project_to_facets().

find_edges_orientations()

ErrorCode moab::SmoothFace::find_edges_orientations ( EntityHandle  edges[3], const EntityHandle *  conn3, int  orient[3]  )

private

Definition at line 436 of file SmoothFace.cpp.

{
 // find the edge that is adjacent to 2 vertices at a time
 for( int i = 0; i < 3; i++ )
 {
 // edge 0 is 1-2, 1 is 3-1, 2 is 0-1
 EntityHandle v[2];
 v[0] = conn3[( i + 1 ) % 3];
 v[1] = conn3[( i + 2 ) % 3];
 std::vector< EntityHandle > adjacencies;
 // generate all edges for these two hexes
 ErrorCode rval = _mb->get_adjacencies( v, 2, 1, false, adjacencies, Interface::INTERSECT );
 if( MB_SUCCESS != rval ) return rval;
 
 // find the edge connected to both vertices, and then see its orientation
 assert( adjacencies.size() == 1 );
 const EntityHandle* conn2 = NULL;
 int nnodes = 0;
 rval = _mb->get_connectivity( adjacencies[0], conn2, nnodes );
 assert( rval == MB_SUCCESS );
 assert( 2 == nnodes );
 edges[i] = adjacencies[0];
 // what is the story morning glory?
 if( conn2[0] == v[0] && conn2[1] == v[1] )
 orient[i] = 1;
 else if( conn2[0] == v[1] && conn2[1] == v[0] )
 orient[i] = -1;
 else
 return MB_FAILURE;
 }
 return MB_SUCCESS;
}

References _mb, ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::INTERSECT, and MB_SUCCESS.

Referenced by compute_internal_control_points_on_facets().

get_normals_for_vertices()

ErrorCode moab::SmoothFace::get_normals_for_vertices	(	const EntityHandle *	conn2,
		CartVect	N[2]
	)

Definition at line 1900 of file SmoothFace.cpp.

{
 // CartVect N[2];
 //_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
 ErrorCode rval = _mb->tag_get_data( _gradientTag, conn2, 2, (double*)&N[0] );
 return rval;
}

References _gradientTag, _mb, ErrorCode, and moab::Interface::tag_get_data().

init_bezier_edge()

ErrorCode moab::SmoothFace::init_bezier_edge ( EntityHandle  edge, double  min_dot  )

private

Definition at line 305 of file SmoothFace.cpp.

{
 // min dot was used for angle here
 // int stat = 0; // CUBIT_SUCCESS;
 // all boundaries will be simple, initially
 // we may complicate them afterwards
 
 CartVect ctrl_pts[3];
 int nnodes = 0;
 const EntityHandle* conn2 = NULL;
 ErrorCode rval = _mb->get_connectivity( edge, conn2, nnodes );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 assert( 2 == nnodes );
 // double coords[6]; // store the coordinates for the nodes
 CartVect P[2];
 // ErrorCode rval = _mb->get_coords(conn2, 2, coords);
 rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 // CartVect P0(&coords[0]);
 // CartVect P3(&coords[3]);
 
 // double normalVec[6];
 CartVect N[2];
 //_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
 rval = _mb->tag_get_data( _gradientTag, conn2, 2, (double*)&N[0] );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 CartVect T[2]; // T0, T3
 
 rval = _mb->tag_get_data( _tangentsTag, &edge, 1, &T[0] );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 rval = init_edge_control_points( P[0], P[1], N[0], N[1], T[0], T[1], ctrl_pts );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 rval = _mb->tag_set_data( _edgeCtrlTag, &edge, 1, &ctrl_pts[0] );
 assert( rval == MB_SUCCESS );
 if( MB_SUCCESS != rval ) return rval;
 
 if( debug_surf_eval1 )
 {
 std::cout << "edge: " << _mb->id_from_handle( edge ) << " tangents: " << T[0] << T[1] << std::endl;
 std::cout << " points: " << P[0] << " " << P[1] << std::endl;
 std::cout << " normals: " << N[0] << " " << N[1] << std::endl;
 std::cout << " Control points " << ctrl_pts[0] << " " << ctrl_pts[1] << " " << ctrl_pts[2] << std::endl;
 }
 
 return MB_SUCCESS;
}

References _edgeCtrlTag, _gradientTag, _mb, _tangentsTag, moab::debug_surf_eval1, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::id_from_handle(), init_edge_control_points(), MB_SUCCESS, T, moab::Interface::tag_get_data(), and moab::Interface::tag_set_data().

Referenced by compute_control_points_on_edges().

init_edge_control_points()

ErrorCode moab::SmoothFace::init_edge_control_points	(	CartVect &	P0,
		CartVect &	P3,
		CartVect &	N0,
		CartVect &	N3,
		CartVect &	T0,
		CartVect &	T3,
		CartVect *	Pi
	)

Definition at line 399 of file SmoothFace.cpp.

{
 CartVect Vi[4];
 Vi[0] = P0;
 Vi[3] = P3;
 CartVect P03( P3 - P0 );
 double di = P03.length();
 double ai = N0 % N3; // this is the dot operator, the same as in cgm for CubitVector
 double ai0 = N0 % T0;
 double ai3 = N3 % T3;
 double denom = 4 - ai * ai;
 if( fabs( denom ) < 1e-20 )
 {
 return MB_FAILURE; // CUBIT_FAILURE;
 }
 double row = 6.0e0 * ( 2.0e0 * ai0 + ai * ai3 ) / denom;
 double omega = 6.0e0 * ( 2.0e0 * ai3 + ai * ai0 ) / denom;
 Vi[1] = Vi[0] + ( di * ( ( ( 6.0e0 * T0 ) - ( ( 2.0e0 * row ) * N0 ) + ( omega * N3 ) ) / 18.0e0 ) );
 Vi[2] = Vi[3] - ( di * ( ( ( 6.0e0 * T3 ) + ( row * N0 ) - ( ( 2.0e0 * omega ) * N3 ) ) / 18.0e0 ) );
 // CartVect Wi[3];
 // Wi[0] = Vi[1] - Vi[0];
 // Wi[1] = Vi[2] - Vi[1];
 // Wi[2] = Vi[3] - Vi[2];
 
 Pi[0] = 0.25 * Vi[0] + 0.75 * Vi[1];
 Pi[1] = 0.50 * Vi[1] + 0.50 * Vi[2];
 Pi[2] = 0.75 * Vi[2] + 0.25 * Vi[3];
 
 return MB_SUCCESS;
}

References moab::CartVect::length(), and MB_SUCCESS.

Referenced by init_bezier_edge().

init_facet_control_points()

ErrorCode moab::SmoothFace::init_facet_control_points ( CartVect  N[6], CartVect  P[3][5], CartVect  G[6]  )

private

Definition at line 580 of file SmoothFace.cpp.

{
 CartVect Di[4], Ai[3], N0, N3, Vi[4], Wi[3];
 double denom;
 double lambda[2], mu[2];
 
 ErrorCode rval = MB_SUCCESS;
 
 for( int i = 0; i < 3; i++ )
 {
 N0 = N[i * 2];
 N3 = N[i * 2 + 1];
 Vi[0] = P[i][0];
 Vi[1] = ( P[i][1] - 0.25 * P[i][0] ) / 0.75;
 Vi[2] = ( P[i][3] - 0.25 * P[i][4] ) / 0.75;
 Vi[3] = P[i][4];
 Wi[0] = Vi[1] - Vi[0];
 Wi[1] = Vi[2] - Vi[1];
 Wi[2] = Vi[3] - Vi[2];
 Di[0] = P[( i + 2 ) % 3][3] - 0.5 * ( P[i][1] + P[i][0] );
 Di[3] = P[( i + 1 ) % 3][1] - 0.5 * ( P[i][4] + P[i][3] );
 Ai[0] = ( N0 * Wi[0] ) / Wi[0].length();
 Ai[2] = ( N3 * Wi[2] ) / Wi[2].length();
 Ai[1] = Ai[0] + Ai[2];
 denom = Ai[1].length();
 Ai[1] /= denom;
 lambda[0] = ( Di[0] % Wi[0] ) / ( Wi[0] % Wi[0] );
 lambda[1] = ( Di[3] % Wi[2] ) / ( Wi[2] % Wi[2] );
 mu[0] = ( Di[0] % Ai[0] );
 mu[1] = ( Di[3] % Ai[2] );
 G[i * 2] = 0.5 * ( P[i][1] + P[i][2] ) + 0.66666666666666 * lambda[0] * Wi[1] +
 0.33333333333333 * lambda[1] * Wi[0] + 0.66666666666666 * mu[0] * Ai[1] +
 0.33333333333333 * mu[1] * Ai[0];
 G[i * 2 + 1] = 0.5 * ( P[i][2] + P[i][3] ) + 0.33333333333333 * lambda[0] * Wi[2] +
 0.66666666666666 * lambda[1] * Wi[1] + 0.33333333333333 * mu[0] * Ai[2] +
 0.66666666666666 * mu[1] * Ai[1];
 }
 return rval;
}

References ErrorCode, moab::CartVect::length(), length(), and MB_SUCCESS.

Referenced by compute_internal_control_points_on_facets().

init_gradient()

int moab::SmoothFace::init_gradient

(

)

Definition at line 170 of file SmoothFace.cpp.

{
 // first, create a Tag for gradient (or normal)
 // loop over all triangles in set, and modify the normals according to the angle as weight
 // get all the edges from this subset
 if( NULL == _mb ) return 1; // fail
 _triangles.clear();
 ErrorCode rval = _mb->get_entities_by_type( _set, MBTRI, _triangles );
 if( MB_SUCCESS != rval ) return 1;
 // get a new range of edges, and decide the loops from here
 _edges.clear();
 rval = _mb->get_adjacencies( _triangles, 1, true, _edges, Interface::UNION );
 assert( rval == MB_SUCCESS );
 
 rval = _mb->get_adjacencies( _triangles, 0, false, _nodes, Interface::UNION );
 assert( rval == MB_SUCCESS );
 
 // initialize bounding box limits
 CartVect vert1;
 EntityHandle v1 = _nodes[0]; // first vertex
 _mb->get_coords( &v1, 1, (double*)&vert1 );
 _minim = vert1;
 _maxim = vert1;
 
 double defNormal[] = { 0., 0., 0. };
 // look for a tag name here that is definitely unique. We do not want the tags to interfere with
 // each other this normal will be for each node in a face some nodes have multiple normals, if
 // they are at the feature edges
 unsigned long setId = _mb->id_from_handle( _set );
 char name[50] = { 0 };
 sprintf( name, "GRADIENT%lu",
 setId ); // name should be something like GRADIENT29, where 29 is the set ID of the face
 rval = _mb->tag_get_handle( name, 3, MB_TYPE_DOUBLE, _gradientTag, MB_TAG_DENSE | MB_TAG_CREAT, &defNormal );
 assert( rval == MB_SUCCESS );
 
 double defPlane[4] = { 0., 0., 1., 0. };
 // also define a plane tag ; this will be for each triangle
 char namePlaneTag[50] = { 0 };
 sprintf( namePlaneTag, "PLANE%lu", setId );
 rval = _mb->tag_get_handle( "PLANE", 4, MB_TYPE_DOUBLE, _planeTag, MB_TAG_DENSE | MB_TAG_CREAT, &defPlane );
 assert( rval == MB_SUCCESS );
 // the fourth double is for weight, accumulated at each vertex so far
 // maybe not needed in the end
 for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
 {
 EntityHandle tria = *it;
 const EntityHandle* conn3;
 int nnodes;
 _mb->get_connectivity( tria, conn3, nnodes );
 if( nnodes != 3 ) return 1; // error
 // double coords[9]; // store the coordinates for the nodes
 //_mb->get_coords(conn3, 3, coords);
 CartVect p[3];
 _mb->get_coords( conn3, 3, (double*)&p[0] );
 // need to compute the angles
 // compute angles and the normal
 // CartVect p1(&coords[0]), p2(&coords[3]), p3(&coords[6]);
 
 CartVect AB( p[1] - p[0] ); //(p2 - p1);
 CartVect BC( p[2] - p[1] ); //(p3 - p2);
 CartVect CA( p[0] - p[2] ); //(p1 - p3);
 double a[3];
 a[1] = angle( AB, -BC ); // angle at B (p2), etc.
 a[2] = angle( BC, -CA );
 a[0] = angle( CA, -AB );
 CartVect normal = -AB * CA;
 normal.normalize();
 double plane[4];
 
 const double* coordNormal = normal.array();
 
 plane[0] = coordNormal[0];
 plane[1] = coordNormal[1];
 plane[2] = coordNormal[2];
 plane[3] = -normal % p[0]; // dot product
 // set the plane
 rval = _mb->tag_set_data( _planeTag, &tria, 1, plane );
 assert( rval == MB_SUCCESS );
 
 // add to each vertex the tag value the normal multiplied by the angle
 double values[9];
 
 _mb->tag_get_data( _gradientTag, conn3, 3, values );
 for( int i = 0; i < 3; i++ )
 {
 // values[4*i]+=a[i]; // first is the weight, which we do not really need
 values[3 * i + 0] += a[i] * coordNormal[0];
 values[3 * i + 1] += a[i] * coordNormal[1];
 values[3 * i + 2] += a[i] * coordNormal[2];
 }
 
 // reset those values
 _mb->tag_set_data( _gradientTag, conn3, 3, values );
 }
 // normalize the gradients at each node; maybe not needed here?
 // no, we will do it, it is important
 int numNodes = _nodes.size();
 double* normalVal = new double[3 * numNodes];
 _mb->tag_get_data( _gradientTag, _nodes,
 normalVal ); // get all the normal values at the _nodes
 for( int i = 0; i < numNodes; i++ )
 {
 CartVect p1( &normalVal[3 * i] );
 p1.normalize();
 p1.get( &normalVal[3 * i] );
 }
 
 // reset the normal values after normalization
 _mb->tag_set_data( _gradientTag, _nodes, normalVal );
 // print the loops size and some other stuff
 if( debug_surf_eval1 )
 {
 std::cout << " normals at " << numNodes << " nodes" << std::endl;
 int i = 0;
 for( Range::iterator it = _nodes.begin(); it != _nodes.end(); ++it, i++ )
 {
 EntityHandle node = *it;
 std::cout << " Node id " << _mb->id_from_handle( node ) << " " << normalVal[3 * i] << " "
 << normalVal[3 * i + 1] << " " << normalVal[3 * i + 2] << std::endl;
 }
 }
 
 delete[] normalVal;
 
 return 0;
}

References _edges, _gradientTag, _maxim, _mb, _minim, _nodes, _planeTag, _set, _triangles, moab::angle(), moab::CartVect::array(), moab::Range::begin(), moab::Range::clear(), moab::debug_surf_eval1, moab::Range::end(), ErrorCode, moab::CartVect::get(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), moab::Interface::id_from_handle(), MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MBTRI, moab::CartVect::normalize(), moab::Range::size(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), and moab::Interface::UNION.

Referenced by moab::FBEngine::initializeSmoothing().

is_at_vertex()

bool moab::SmoothFace::is_at_vertex	(	EntityHandle	facet,
		CartVect &	pt,
		CartVect &	ac,
		double	compare_tol,
		CartVect &	eval_pt,
		CartVect *	eval_norm_ptr
	)

Definition at line 1576 of file SmoothFace.cpp.

{
 double dist;
 CartVect vert_loc;
 const double actol = 0.1;
 // get coordinates get_coords
 const EntityHandle* conn3 = NULL;
 int nnodes = 0;
 _mb->get_connectivity( facet, conn3, nnodes );
 //
 // double coords[9]; // store the coordinates for the nodes
 //_mb->get_coords(conn3, 3, coords);
 CartVect p[3];
 _mb->get_coords( conn3, 3, (double*)&p[0] );
 // also get the normals at nodes
 CartVect NN[3];
 _mb->tag_get_data( _gradientTag, conn3, 3, (double*)&NN[0] );
 if( fabs( ac[0] ) < actol && fabs( ac[1] ) < actol )
 {
 vert_loc = p[2];
 dist = ( pt - vert_loc ).length(); // pt.distance_between( vert_loc );
 if( dist <= compare_tol )
 {
 eval_pt = vert_loc;
 if( eval_norm_ptr )
 {
 *eval_norm_ptr = NN[2];
 }
 return true;
 }
 }
 
 if( fabs( ac[0] ) < actol && fabs( ac[2] ) < actol )
 {
 vert_loc = p[1];
 dist = ( pt - vert_loc ).length(); // pt.distance_between( vert_loc );
 if( dist <= compare_tol )
 {
 eval_pt = vert_loc;
 if( eval_norm_ptr )
 {
 *eval_norm_ptr = NN[1]; // facet->point(1)->normal( facet );
 }
 return true;
 }
 }
 
 if( fabs( ac[1] ) < actol && fabs( ac[2] ) < actol )
 {
 vert_loc = p[0];
 dist = ( pt - vert_loc ).length(); // pt.distance_between( vert_loc );
 if( dist <= compare_tol )
 {
 eval_pt = vert_loc;
 if( eval_norm_ptr )
 {
 *eval_norm_ptr = NN[0];
 }
 return true;
 }
 }
 
 return false;
}

References _gradientTag, _mb, moab::Interface::get_connectivity(), moab::Interface::get_coords(), length(), and moab::Interface::tag_get_data().

Referenced by project_to_patch().

move_ac_inside()

bool moab::SmoothFace::move_ac_inside ( CartVect &  ac, double  tol  )

private

Definition at line 1654 of file SmoothFace.cpp.

{
 int nout = 0;
 if( ac[0] < -tol )
 {
 ac[0] = 0.0;
 ac[1] = ac[1] / ( ac[1] + ac[2] ); //( ac.y() / (ac.y() + ac.z()) ;
 ac[2] = 1. - ac[1]; // ac.z( 1.0 - ac.y() );
 nout++;
 }
 if( ac[1] < -tol )
 {
 ac[1] = 0.; // ac.y( 0.0 );
 ac[0] = ac[0] / ( ac[0] + ac[2] ); // ac.x( ac.x() / (ac.x() + ac.z()) );
 ac[2] = 1. - ac[0]; // ac.z( 1.0 - ac.x() );
 nout++;
 }
 if( ac[2] < -tol )
 {
 ac[2] = 0.; // ac.z( 0.0 );
 ac[0] = ac[0] / ( ac[0] + ac[1] ); // ac.x( ac.x() / (ac.x() + ac.y()) );
 ac[1] = 1. - ac[0]; // ac.y( 1.0 - ac.x() );
 nout++;
 }
 return ( nout > 0 ) ? true : false;
}

Referenced by project_to_patch().

move_to_surface()

void moab::SmoothFace::move_to_surface ( double &  x, double &  y, double &  z  )

virtual

Definition at line 104 of file SmoothFace.cpp.

{
 
 CartVect loc2( x, y, z );
 bool trim = false; // is it needed?
 bool outside = true;
 CartVect closestPoint;
 
 ErrorCode rval = project_to_facets_main( loc2, trim, outside, &closestPoint, NULL );
 if( MB_SUCCESS != rval ) return;
 assert( rval == MB_SUCCESS );
 x = closestPoint[0];
 y = closestPoint[1];
 z = closestPoint[2];
}

References ErrorCode, MB_SUCCESS, and project_to_facets_main().

Referenced by moab::FBEngine::getEntClosestPt(), and moab::FBEngine::smooth_new_intx_points().

normal_at()

bool moab::SmoothFace::normal_at ( double  x, double  y, double  z, double &  nx, double &  ny, double &  nz  )

virtual

Definition at line 127 of file SmoothFace.cpp.

{
 
 CartVect loc2( x, y, z );
 
 bool trim = false; // is it needed?
 bool outside = true;
 // CartVect closestPoint;// not needed
 // not interested in normal
 CartVect normal;
 ErrorCode rval = project_to_facets_main( loc2, trim, outside, NULL, &normal );
 if( MB_SUCCESS != rval ) return false;
 assert( rval == MB_SUCCESS );
 nx = normal[0];
 ny = normal[1];
 nz = normal[2];
 
 return true;
}

References ErrorCode, MB_SUCCESS, and project_to_facets_main().

Referenced by moab::FBEngine::getEntNrmlXYZ().

project_to_facet()

ErrorCode moab::SmoothFace::project_to_facet	(	EntityHandle	facet,
		CartVect &	pt,
		CartVect &	areacoord,
		CartVect &	close_point,
		bool &	outside_facet,
		double	compare_tol
	)

Definition at line 1545 of file SmoothFace.cpp.

{
 const EntityHandle* conn3 = NULL;
 int nnodes = 0;
 _mb->get_connectivity( facet, conn3, nnodes );
 //
 // double coords[9]; // store the coordinates for the nodes
 //_mb->get_coords(conn3, 3, coords);
 CartVect p[3];
 _mb->get_coords( conn3, 3, (double*)&p[0] );
 
 int edge_id = -1;
 ErrorCode stat = project_to_patch( facet, areacoord, pt, close_point, NULL, outside_facet, compare_tol, edge_id );
 /* }
 break;
 }*/
 
 return stat;
}

References _mb, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), and project_to_patch().

Referenced by project_to_facets().

project_to_facet_plane()

void moab::SmoothFace::project_to_facet_plane	(	EntityHandle	tri,
		CartVect &	pt,
		CartVect &	point_on_plane,
		double &	dist_to_plane
	)

Definition at line 902 of file SmoothFace.cpp.

{
 double plane[4];
 
 ErrorCode rval = _mb->tag_get_data( _planeTag, &tri, 1, plane );
 if( MB_SUCCESS != rval ) return;
 assert( rval == MB_SUCCESS );
 // _planeTag
 CartVect normal( &plane[0] ); // just first 3 components are used
 
 double dist = normal % pt + plane[3]; // coeff d is saved!!!
 dist_to_plane = fabs( dist );
 point_on_plane = pt - dist * normal;
 
 return;
}

References _mb, _planeTag, ErrorCode, MB_SUCCESS, and moab::Interface::tag_get_data().

Referenced by project_to_facets().

project_to_facets()

ErrorCode moab::SmoothFace::project_to_facets	(	std::vector< EntityHandle > &	facet_list,
		EntityHandle &	lastFacet,
		int	interpOrder,
		double	compareTol,
		CartVect &	this_point,
		bool	trim,
		bool &	outside,
		CartVect *	closest_point_ptr,
		CartVect *	normal_ptr
	)

best_facet

Definition at line 1101 of file SmoothFace.cpp.

{
 
 bool outside_facet = false;
 bool best_outside_facet = true;
 double mindist = 1.e20;
 CartVect close_point, best_point( mindist, mindist, mindist ), best_areacoord;
 EntityHandle best_facet = 0L; // no best facet found yet
 EntityHandle facet;
 assert( facet_list.size() > 0 );
 
 double big_dist = compareTol * 1.0e3;
 
 // from the list of close facets, determine the closest point
 for( size_t i = 0; i < facet_list.size(); i++ )
 {
 facet = facet_list[i];
 CartVect pt_on_plane;
 double dist_to_plane;
 project_to_facet_plane( facet, this_point, pt_on_plane, dist_to_plane );
 
 CartVect areacoord;
 // CartVect close_point;
 facet_area_coordinate( facet, pt_on_plane, areacoord );
 if( interpOrder != 0 )
 {
 
 // modify the areacoord - project to the bezier patch- snaps to the
 // edge of the patch if necessary
 
 if( project_to_facet( facet, this_point, areacoord, close_point, outside_facet, compareTol ) != MB_SUCCESS )
 {
 return MB_FAILURE;
 }
 // if (closest_point_ptr)
 //*closest_point_ptr = close_point;
 }
 // keep track of the minimum distance
 
 double dist = ( close_point - this_point ).length(); // close_point.distance_between(this_point);
 if( ( best_outside_facet == outside_facet && dist < mindist ) ||
 ( best_outside_facet && !outside_facet && ( dist < big_dist || best_facet == 0L /*!best_facet*/ ) ) )
 {
 mindist = dist;
 best_point = close_point;
 best_facet = facet;
 best_areacoord = areacoord;
 best_outside_facet = outside_facet;
 
 if( dist < compareTol )
 {
 break;
 }
 big_dist = 10.0 * mindist;
 }
 // facet->marked(1);
 // used_facet_list.append(facet);
 }
 
 if( normal_ptr )
 {
 CartVect normal;
 if( eval_bezier_patch_normal( best_facet, best_areacoord, normal ) != MB_SUCCESS )
 {
 return MB_FAILURE;
 }
 *normal_ptr = normal;
 }
 
 if( closest_point_ptr )
 {
 *closest_point_ptr = best_point;
 }
 
 outside = best_outside_facet;
 lastFacet = best_facet;
 
 return MB_SUCCESS;
 // end copy
}

References eval_bezier_patch_normal(), facet_area_coordinate(), length(), MB_SUCCESS, project_to_facet(), and project_to_facet_plane().

Referenced by project_to_facets_main().

project_to_facets_main()

ErrorCode moab::SmoothFace::project_to_facets_main	(	CartVect &	this_point,
		bool	trim,
		bool &	outside,
		CartVect *	closest_point_ptr = NULL,
		CartVect *	normal_ptr = NULL
	)

Definition at line 1075 of file SmoothFace.cpp.

{
 
 // if there are a lot of facets on this surface - use the OBB search first
 // to narrow the selection
 
 _evaluationsCounter++;
 double tolerance = 1.e-5;
 std::vector< EntityHandle > facets_out;
 // we will start with a list of facets anyway, the best among them wins
 ErrorCode rval =
 _my_geomTopoTool->obb_tree()->closest_to_location( (double*)&this_point, _obb_root, tolerance, facets_out );
 if( MB_SUCCESS != rval ) return rval;
 
 int interpOrder = 4;
 double compareTol = 1.e-5;
 EntityHandle lastFacet = facets_out.front();
 rval = project_to_facets( facets_out, lastFacet, interpOrder, compareTol, this_point, trim, outside,
 closest_point_ptr, normal_ptr );
 
 return rval;
}

References _evaluationsCounter, _my_geomTopoTool, _obb_root, moab::OrientedBoxTreeTool::closest_to_location(), ErrorCode, MB_SUCCESS, moab::GeomTopoTool::obb_tree(), project_to_facets(), and moab::tolerance.

Referenced by move_to_surface(), normal_at(), and ray_intersection_correct().

project_to_patch()

ErrorCode moab::SmoothFace::project_to_patch	(	EntityHandle	facet,
		CartVect &	ac,
		CartVect &	pt,
		CartVect &	eval_pt,
		CartVect *	eval_norm,
		bool &	outside,
		double	compare_tol,
		int	edge_id
	)

Definition at line 1199 of file SmoothFace.cpp.

{
 ErrorCode status = MB_SUCCESS;
 
 // see if we are at a vertex
 
#define INCR 0.01
 const double tol = compare_tol;
 
 if( is_at_vertex( facet, pt, ac, compare_tol, eval_pt, eval_norm ) )
 {
 outside = false;
 return MB_SUCCESS;
 }
 
 // check if the start ac is inside the patch -if not, then move it there
 
 int nout = 0;
 const double atol = 0.001;
 if( move_ac_inside( ac, atol ) ) nout++;
 
 int diverge = 0;
 int iter = 0;
 CartVect newpt;
 eval_bezier_patch( facet, ac, newpt );
 CartVect move = pt - newpt;
 double lastdist = move.length();
 double bestdist = lastdist;
 CartVect bestac = ac;
 CartVect bestpt = newpt;
 CartVect bestnorm( 0, 0, 0 );
 
 // If we are already close enough, then return now
 
 if( lastdist <= tol && !eval_norm && nout == 0 )
 {
 eval_pt = pt;
 outside = false;
 return status;
 }
 
 double ratio, mag, umove, vmove, det, distnew, movedist;
 CartVect lastpt = newpt;
 CartVect lastac = ac;
 CartVect norm;
 CartVect xpt, ypt, zpt, xac, yac, zac, xvec, yvec, zvec;
 CartVect du, dv, newac;
 bool done = false;
 while( !done )
 {
 
 // We will be locating the projected point within the u,v,w coordinate
 // system of the triangular bezier patch. Since u+v+w=1, the components
 // are not linearly independent. We will choose only two of the
 // coordinates to use and compute the third.
 
 int system;
 if( lastac[0] >= lastac[1] && lastac[0] >= lastac[2] )
 {
 system = 0;
 }
 else if( lastac[1] >= lastac[2] )
 {
 system = 1;
 }
 else
 {
 system = 2;
 }
 
 // compute the surface derivatives with respect to each
 // of the barycentric coordinates
 
 if( system == 1 || system == 2 )
 {
 xac[0] = lastac[0] + INCR; // xac.x( lastac.x() + INCR );
 if( lastac[1] + lastac[2] == 0.0 ) return MB_FAILURE;
 ratio = lastac[2] / ( lastac[1] + lastac[2] ); // ratio = lastac.z() / (lastac.y() + lastac.z());
 xac[1] = ( 1.0 - xac[0] ) * ( 1.0 - ratio ); // xac.y( (1.0 - xac.x()) * (1.0 - ratio) );
 xac[2] = 1.0 - xac[0] - xac[1]; // xac.z( 1.0 - xac.x() - xac.y() );
 eval_bezier_patch( facet, xac, xpt );
 xvec = xpt - lastpt;
 xvec /= INCR;
 }
 if( system == 0 || system == 2 )
 {
 yac[1] = ( lastac[1] + INCR ); // yac.y( lastac.y() + INCR );
 if( lastac[0] + lastac[2] == 0.0 ) // if (lastac.x() + lastac.z() == 0.0)
 return MB_FAILURE;
 ratio = lastac[2] / ( lastac[0] + lastac[2] ); // ratio = lastac.z() / (lastac.x() + lastac.z());
 yac[0] = ( ( 1.0 - yac[1] ) * ( 1.0 - ratio ) ); // yac.x( (1.0 - yac.y()) * (1.0 - ratio) );
 yac[2] = ( 1.0 - yac[0] - yac[1] ); // yac.z( 1.0 - yac.x() - yac.y() );
 eval_bezier_patch( facet, yac, ypt );
 yvec = ypt - lastpt;
 yvec /= INCR;
 }
 if( system == 0 || system == 1 )
 {
 zac[2] = ( lastac[2] + INCR ); // zac.z( lastac.z() + INCR );
 if( lastac[0] + lastac[1] == 0.0 ) // if (lastac.x() + lastac.y() == 0.0)
 return MB_FAILURE;
 ratio = lastac[1] / ( lastac[0] + lastac[1] ); // ratio = lastac.y() / (lastac.x() + lastac.y());
 zac[0] = ( ( 1.0 - zac[2] ) * ( 1.0 - ratio ) ); // zac.x( (1.0 - zac.z()) * (1.0 - ratio) );
 zac[1] = ( 1.0 - zac[0] - zac[2] ); // zac.y( 1.0 - zac.x() - zac.z() );
 eval_bezier_patch( facet, zac, zpt );
 zvec = zpt - lastpt;
 zvec /= INCR;
 }
 
 // compute the surface normal
 
 switch( system )
 {
 case 0:
 du = yvec;
 dv = zvec;
 break;
 case 1:
 du = zvec;
 dv = xvec;
 break;
 case 2:
 du = xvec;
 dv = yvec;
 break;
 }
 norm = du * dv;
 mag = norm.length();
 if( mag < DBL_EPSILON )
 {
 return MB_FAILURE;
 // do something else here (it is likely a flat triangle -
 // so try evaluating just an edge of the bezier patch)
 }
 norm /= mag;
 if( iter == 0 ) bestnorm = norm;
 
 // project the move vector to the tangent plane
 
 move = ( norm * move ) * norm;
 
 // compute an equivalent u-v-w vector
 
 CartVect absnorm( fabs( norm[0] ), fabs( norm[1] ), fabs( norm[2] ) );
 if( absnorm[2] >= absnorm[1] && absnorm[2] >= absnorm[0] )
 {
 det = du[0] * dv[1] - dv[0] * du[1];
 if( fabs( det ) <= DBL_EPSILON )
 {
 return MB_FAILURE; // do something else here
 }
 umove = ( move[0] * dv[1] - dv[0] * move[1] ) / det;
 vmove = ( du[0] * move[1] - move[0] * du[1] ) / det;
 }
 else if( absnorm[1] >= absnorm[2] && absnorm[1] >= absnorm[0] )
 {
 det = du[0] * dv[2] - dv[0] * du[2];
 if( fabs( det ) <= DBL_EPSILON )
 {
 return MB_FAILURE;
 }
 umove = ( move[0] * dv[2] - dv[0] * move[2] ) / det;
 vmove = ( du[0] * move[2] - move[0] * du[2] ) / det;
 }
 else
 {
 det = du[1] * dv[2] - dv[1] * du[2];
 if( fabs( det ) <= DBL_EPSILON )
 {
 return MB_FAILURE;
 }
 umove = ( move[1] * dv[2] - dv[1] * move[2] ) / det;
 vmove = ( du[1] * move[2] - move[1] * du[2] ) / det;
 }
 
 /* === compute the new u-v coords and evaluate surface at new location */
 
 switch( system )
 {
 case 0:
 newac[1] = ( lastac[1] + umove ); // newac.y( lastac.y() + umove );
 newac[2] = ( lastac[2] + vmove ); // newac.z( lastac.z() + vmove );
 newac[0] = ( 1.0 - newac[1] - newac[2] ); // newac.x( 1.0 - newac.y() - newac.z() );
 break;
 case 1:
 newac[2] = ( lastac[2] + umove ); // newac.z( lastac.z() + umove );
 newac[0] = ( lastac[0] + vmove ); // newac.x( lastac.x() + vmove );
 newac[1] = ( 1.0 - newac[2] - newac[0] ); // newac.y( 1.0 - newac.z() - newac.x() );
 break;
 case 2:
 newac[0] = ( lastac[0] + umove ); // newac.x( lastac.x() + umove );
 newac[1] = ( lastac[1] + vmove ); // newac.y( lastac.y() + vmove );
 newac[2] = ( 1.0 - newac[0] - newac[1] ); // newac.z( 1.0 - newac.x() - newac.y() );
 break;
 }
 
 // Keep it inside the patch
 
 if( newac[0] >= -atol && newac[1] >= -atol && newac[2] >= -atol )
 {
 nout = 0;
 }
 else
 {
 if( move_ac_inside( newac, atol ) ) nout++;
 }
 
 // Evaluate at the new location
 
 if( edge_id != -1 ) ac_at_edge( newac, newac, edge_id ); // move to edge first
 eval_bezier_patch( facet, newac, newpt );
 
 // Check for convergence
 
 distnew = ( pt - newpt ).length(); // pt.distance_between(newpt);
 move = newpt - lastpt;
 movedist = move.length();
 if( movedist < tol || distnew < tol )
 {
 done = true;
 if( distnew < bestdist )
 {
 bestdist = distnew;
 bestac = newac;
 bestpt = newpt;
 bestnorm = norm;
 }
 }
 else
 {
 
 // don't allow more than 30 iterations
 
 iter++;
 if( iter > 30 )
 {
 // if (movedist > tol * 100.0) nout=1;
 done = true;
 }
 
 // Check for divergence - don't allow more than 5 divergent
 // iterations
 
 if( distnew > lastdist )
 {
 diverge++;
 if( diverge > 10 )
 {
 done = true;
 // if (movedist > tol * 100.0) nout=1;
 }
 }
 
 // Check if we are continuing to project outside the facet.
 // If so, then stop now
 
 if( nout > 3 )
 {
 done = true;
 }
 
 // set up for next iteration
 
 if( !done )
 {
 if( distnew < bestdist )
 {
 bestdist = distnew;
 bestac = newac;
 bestpt = newpt;
 bestnorm = norm;
 }
 lastdist = distnew;
 lastpt = newpt;
 lastac = newac;
 move = pt - lastpt;
 }
 }
 }
 
 eval_pt = bestpt;
 if( eval_norm )
 {
 *eval_norm = bestnorm;
 }
 outside = ( nout > 0 ) ? true : false;
 ac = bestac;
 
 return status;
}

References ac_at_edge(), ErrorCode, eval_bezier_patch(), INCR, is_at_vertex(), moab::CartVect::length(), length(), MB_SUCCESS, and move_ac_inside().

Referenced by project_to_facet().

ray_intersection_correct()

ErrorCode moab::SmoothFace::ray_intersection_correct	(	EntityHandle	facet,
		CartVect &	pt,
		CartVect &	ray,
		CartVect &	eval_pt,
		double &	distance,
		bool &	outside
	)

Definition at line 1910 of file SmoothFace.cpp.

{
 // find a point on the smooth surface
 CartVect currentPoint = eval_pt;
 int numIter = 0;
 double improvement = 1.e20;
 CartVect diff;
 while( numIter++ < 5 && improvement > 0.01 )
 {
 CartVect newPos;
 
 bool trim = false; // is it needed?
 outside = true;
 CartVect closestPoint;
 CartVect normal;
 
 ErrorCode rval = project_to_facets_main( currentPoint, trim, outside, &newPos, &normal );
 if( MB_SUCCESS != rval ) return rval;
 assert( rval == MB_SUCCESS );
 diff = newPos - currentPoint;
 improvement = diff.length();
 // ( pt + t * ray - closest ) % normal = 0;
 // intersect tangent plane that goes through closest point with the direction
 // t = normal%(closest-pt) / normal%ray;
 double dot = normal % ray; // if it is 0, get out while we can
 if( dot < 0.00001 )
 {
 // bad convergence, get out, do not modify anything
 return MB_SUCCESS;
 }
 double t = ( ( newPos - pt ) % normal ) / ( dot );
 currentPoint = pt + t * ray;
 }
 eval_pt = currentPoint;
 diff = currentPoint - pt;
 distance = diff.length();
 return MB_SUCCESS;
}

References moab::dot(), ErrorCode, moab::CartVect::length(), MB_SUCCESS, project_to_facets_main(), and t.

Referenced by moab::FBEngine::getPntRayIntsct().

Member Data Documentation

_edgeCtrlTag

Tag moab::SmoothFace::_edgeCtrlTag

private

Definition at line 211 of file SmoothFace.hpp.

Referenced by compute_control_points_on_edges(), compute_internal_control_points_on_facets(), DumpModelControlPoints(), evaluate_smooth_edge(), and init_bezier_edge().

_edges

Range moab::SmoothFace::_edges

private

Definition at line 180 of file SmoothFace.hpp.

Referenced by compute_control_points_on_edges(), compute_tangents_for_each_edge(), DumpModelControlPoints(), and init_gradient().

_evaluationsCounter

long moab::SmoothFace::_evaluationsCounter

private

Definition at line 237 of file SmoothFace.hpp.

Referenced by eval_counter(), and project_to_facets_main().

_facetCtrlTag

Tag moab::SmoothFace::_facetCtrlTag

private

Definition at line 217 of file SmoothFace.hpp.

Referenced by compute_internal_control_points_on_facets(), DumpModelControlPoints(), eval_bezier_patch(), and eval_bezier_patch_normal().

_facetEdgeCtrlTag

Tag moab::SmoothFace::_facetEdgeCtrlTag

private

Definition at line 225 of file SmoothFace.hpp.

Referenced by compute_internal_control_points_on_facets(), eval_bezier_patch(), and eval_bezier_patch_normal().

_gradientTag

Tag moab::SmoothFace::_gradientTag

private

Definition at line 199 of file SmoothFace.hpp.

Referenced by append_smooth_tags(), compute_internal_control_points_on_facets(), eval_bezier_patch_normal(), get_normals_for_vertices(), init_bezier_edge(), init_gradient(), and is_at_vertex().

_markTag

Tag moab::SmoothFace::_markTag

private

Definition at line 194 of file SmoothFace.hpp.

Referenced by compute_control_points_on_edges().

_maxim

CartVect moab::SmoothFace::_maxim

private

Definition at line 177 of file SmoothFace.hpp.

Referenced by adjust_bounding_box(), bounding_box(), and init_gradient().

_mb

Interface* moab::SmoothFace::_mb

private

Definition at line 232 of file SmoothFace.hpp.

Referenced by area(), compute_control_points_on_edges(), compute_internal_control_points_on_facets(), compute_tangents_for_each_edge(), DumpModelControlPoints(), eval_bezier_patch(), eval_bezier_patch_normal(), evaluate_smooth_edge(), facet_area_coordinate(), find_edges_orientations(), get_normals_for_vertices(), init_bezier_edge(), init_gradient(), is_at_vertex(), project_to_facet(), and project_to_facet_plane().

_minim

CartVect moab::SmoothFace::_minim

private

Definition at line 176 of file SmoothFace.hpp.

Referenced by adjust_bounding_box(), bounding_box(), and init_gradient().

_my_geomTopoTool

GeomTopoTool* moab::SmoothFace::_my_geomTopoTool

private

Definition at line 234 of file SmoothFace.hpp.

Referenced by project_to_facets_main(), and SmoothFace().

_nodes

Range moab::SmoothFace::_nodes

private

Definition at line 181 of file SmoothFace.hpp.

Referenced by init_gradient().

_obb_root

EntityHandle moab::SmoothFace::_obb_root

private

Definition at line 235 of file SmoothFace.hpp.

Referenced by project_to_facets_main(), and SmoothFace().

_planeTag

Tag moab::SmoothFace::_planeTag

private

Definition at line 230 of file SmoothFace.hpp.

Referenced by append_smooth_tags(), facet_area_coordinate(), init_gradient(), and project_to_facet_plane().

_set

EntityHandle moab::SmoothFace::_set

private

Definition at line 233 of file SmoothFace.hpp.

Referenced by DumpModelControlPoints(), init_gradient(), and SmoothFace().

_tangentsTag

Tag moab::SmoothFace::_tangentsTag

private

Definition at line 204 of file SmoothFace.hpp.

Referenced by compute_tangents_for_each_edge(), and init_bezier_edge().

_triangles

Range moab::SmoothFace::_triangles

private

Definition at line 179 of file SmoothFace.hpp.

Referenced by area(), compute_internal_control_points_on_facets(), DumpModelControlPoints(), and init_gradient().

---

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

• SmoothFace.hpp
• SmoothFace.cpp