SmoothCurve.cpp

Go to the documentation of this file.

/*
 * SmoothCurve.cpp
 *
 */
 
#include "SmoothCurve.hpp"
//#include "SmoothVertex.hpp"
#include "SmoothFace.hpp"
#include <cassert>
#include <iostream>
#include "moab/GeomTopoTool.hpp"
 
namespace moab
{
 
SmoothCurve::SmoothCurve( Interface* mb, EntityHandle curve, GeomTopoTool* gTool )
 : _mb( mb ), _set( curve ), _gtt( gTool )
{
 //_mbOut->create_meshset(MESHSET_ORDERED, _oSet);
 /*_cmlEdgeMesher = new CMLEdgeMesher (this, CML::STANDARD);
 _cmlEdgeMesher->set_sizing_function(CML::LINEAR_SIZING);*/
 _leng = 0; // not initialized
 _edgeTag = 0; // not initialized
}
SmoothCurve::~SmoothCurve()
{
 // TODO Auto-generated destructor stub
}
 
double SmoothCurve::arc_length()
{
 
 return _leng;
}
 
//! \brief Get the parametric status of the curve.
//!
//! \return \a true if curve is parametric, \a false otherwise.
bool SmoothCurve::is_parametric()
{
 return true;
}
 
//! \brief Get the periodic status of the curve.
//!
//! \param period The period of the curve if periodic.
//!
//! \return \a true if curve is periodic, \a false otherwise.
bool SmoothCurve::is_periodic( double& period )
{
 // assert(_ref_edge);
 // return _ref_edge->is_periodic( period);
 Range vsets;
 _mb->get_child_meshsets( _set, vsets ); // num_hops =1
 if( vsets.size() == 1 )
 {
 period = _leng;
 return true; // true , especially for ice sheet data
 }
 return false;
}
 
//! \brief Get the parameter range of the curve.
//!
//! \param u_start The beginning curve parameter
//! \param u_end The ending curve parameter
//!
//! \note The numerical value of \a u_start may be greater
//! than the numerical value of \a u_end.
void SmoothCurve::get_param_range( double& u_start, double& u_end )
{
 // assert(_ref_edge);
 u_start = 0;
 u_end = 1.;
 
 return;
}
 
//! Compute the parameter value at a specified distance along the curve.
//!
//! \param u_root The start parameter from which to compute the distance
//! along the curve.
//! \param arc_length The distance to move along the curve.
//!
//! \note For positive values of \a arc_length the distance will be
//! computed in the direction of increasing parameter value along the
//! curve. For negative values of \a arc_length the distance will be
//! computed in the direction of decreasing parameter value along the
//! curve.
//!
//! \return The parametric coordinate u along the curve
double SmoothCurve::u_from_arc_length( double u_root, double arc_leng )
{
 
 if( _leng <= 0 ) return 0;
 return u_root + arc_leng / _leng;
}
 
//! \brief Evaluate the curve at a specified parameter value.
//!
//! \param u The parameter at which to evaluate the curve
//! \param x The x coordinate of the evaluated point
//! \param y The y coordinate of the evaluated point
//! \param z The z coordinate of the evaluated point
bool SmoothCurve::position_from_u( double u, double& x, double& y, double& z, double* tg )
{
 
 // _fractions are increasing, so find the
 double* ptr = std::lower_bound( &_fractions[0], ( &_fractions[0] ) + _fractions.size(), u );
 int index = ptr - &_fractions[0];
 double nextFraction = _fractions[index];
 double prevFraction = 0;
 if( index > 0 )
 {
 prevFraction = _fractions[index - 1];
 }
 double t = ( u - prevFraction ) / ( nextFraction - prevFraction );
 
 EntityHandle edge = _entities[index];
 
 CartVect position, tangent;
 ErrorCode rval = evaluate_smooth_edge( edge, t, position, tangent );
 if( MB_SUCCESS != rval ) return false;
 assert( rval == MB_SUCCESS );
 x = position[0];
 y = position[1];
 z = position[2];
 if( tg )
 {
 // we need to do some scaling,
 double dtdu = 1 / ( nextFraction - prevFraction );
 tg[0] = tangent[0] * dtdu;
 tg[1] = tangent[1] * dtdu;
 tg[2] = tangent[2] * dtdu;
 }
 
 return true;
}
//! \brief Move a point near the curve to the closest point on the curve.
//!
//! \param x The x coordinate of the point
//! \param y The y coordinate of the point
//! \param z The z coordinate of the point
void SmoothCurve::move_to_curve( double& x, double& y, double& z )
{
 
 // find closest point to the curve, and the parametric position
 // must be close by, but how close ???
 EntityHandle v;
 int edgeIndex;
 double u = u_from_position( x, y, z, v, edgeIndex );
 position_from_u( u, x, y, z );
 
 return;
}
 
//! Get the u parameter value on the curve closest to x,y,z
//! and the point on the curve.
//!
//! \param x The x coordinate of the point
//! \param y The y coordinate of the point
//! \param z The z coordinate of the point
//!
//! \return The parametric coordinate u on the curve
double SmoothCurve::u_from_position( double x, double y, double z, EntityHandle& v, int& edgeIndex )
{
 // this is an iterative process, expensive usually
 // get first all nodes , and their positions
 // find the closest node (and edge), and from there do some
 // iterations up to a point
 // do not exaggerate with convergence criteria
 
 v = 0; // we do not have a close by vertex yet
 CartVect initialPos( x, y, z );
 double u = 0;
 int nbNodes = (int)_entities.size() * 2; // the mesh edges are stored
 std::vector< EntityHandle > nodesConnec;
 nodesConnec.resize( nbNodes );
 ErrorCode rval = this->_mb->get_connectivity( &( _entities[0] ), nbNodes / 2, nodesConnec );
 if( MB_SUCCESS != rval )
 {
 std::cout << "error in getting connectivity\n";
 return 0;
 }
 // collapse nodesConnec, nodes should be in order
 for( int k = 0; k < nbNodes / 2; k++ )
 {
 nodesConnec[k + 1] = nodesConnec[2 * k + 1];
 }
 int numNodes = nbNodes / 2 + 1;
 std::vector< CartVect > coordNodes;
 coordNodes.resize( numNodes );
 
 rval = _mb->get_coords( &( nodesConnec[0] ), numNodes, (double*)&( coordNodes[0] ) );
 if( MB_SUCCESS != rval )
 {
 std::cout << "error in getting node positions\n";
 return 0;
 }
 // find the closest node, then find the closest edge, based on closest node
 
 int indexNode = 0;
 double minDist = 1.e30;
 // expensive linear search
 for( int i = 0; i < numNodes; i++ )
 {
 double d1 = ( initialPos - coordNodes[i] ).length();
 if( d1 < minDist )
 {
 indexNode = i;
 minDist = d1;
 }
 }
 double tolerance = 0.00001; // what is the unit?
 // something reasonable
 if( minDist < tolerance )
 {
 v = nodesConnec[indexNode];
 // we are done, just return the proper u (from fractions)
 if( indexNode == 0 )
 {
 return 0; // first node has u = 0
 }
 else
 return _fractions[indexNode - 1]; // fractions[0] > 0!!)
 }
 // find the mesh edge; could be previous or next edge
 edgeIndex = indexNode; // could be the previous one, though!!
 if( edgeIndex == numNodes - 1 )
 edgeIndex--; // we have one less edge, and do not worry about next edge
 else
 {
 if( edgeIndex > 0 )
 {
 // could be the previous; decide based on distance to the other
 // nodes of the 2 connected edges
 CartVect prevNodePos = coordNodes[edgeIndex - 1];
 CartVect nextNodePos = coordNodes[edgeIndex + 1];
 if( ( prevNodePos - initialPos ).length_squared() < ( nextNodePos - initialPos ).length_squared() )
 {
 edgeIndex--;
 }
 }
 }
 // now, we know for sure that the closest point is somewhere on edgeIndex edge
 //
 
 // do newton iteration for local t between 0 and 1
 
 // copy from evaluation method
 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;
 
 P[0] = coordNodes[edgeIndex];
 P[1] = coordNodes[edgeIndex + 1];
 
 if( 0 == _edgeTag )
 {
 rval = _mb->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, _edgeTag );
 if( rval != MB_SUCCESS ) return 0;
 }
 rval = _mb->tag_get_data( _edgeTag, &( _entities[edgeIndex] ), 1, (double*)&controlPoints[0] );
 if( rval != MB_SUCCESS ) return rval;
 
 // starting point
 double tt = 0.5; // between the 2 ends of the edge
 int iterations = 0;
 // find iteratively a better point
 int maxIterations = 10; // not too many
 CartVect outv;
 // we will solve minimize F = 0.5 * ( ini - r(t) )^2
 // so solve F'(t) = 0
 // Iteration: t_ -> t - F'(t)/F"(t)
 // F'(t) = r'(t) (ini-r(t) )
 // F"(t) = r"(t) (ini-r(t) ) - (r'(t))^2
 while( iterations < maxIterations )
 //
 {
 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];
 
 CartVect out_tangent = -4. * one_minus_t3 * P[0] +
 4. * ( one_minus_t3 - 3. * tt * one_minus_t2 ) * controlPoints[0] +
 12. * ( tt * one_minus_t2 - t2 * one_minus_t ) * controlPoints[1] +
 4. * ( 3. * t2 * one_minus_t - t3 ) * controlPoints[2] + 4. * t3 * P[1];
 
 CartVect second_deriv =
 12. * one_minus_t2 * P[0] +
 4. * ( -3. * one_minus_t2 - 3. * one_minus_t2 + 6. * tt * one_minus_t ) * controlPoints[0] +
 12. * ( one_minus_t2 - 4 * tt * one_minus_t + t2 ) * controlPoints[1] +
 4. * ( 6. * tt - 12 * t2 ) * controlPoints[2] + 12. * t2 * P[1];
 CartVect diff = outv - initialPos;
 double F_d = out_tangent % diff;
 double F_dd = second_deriv % diff + out_tangent.length_squared();
 
 if( 0 == F_dd ) break; // get out, we found minimum?
 
 double delta_t = -F_d / F_dd;
 
 if( fabs( delta_t ) < 0.000001 ) break;
 tt = tt + delta_t;
 if( tt < 0 )
 {
 tt = 0.;
 v = nodesConnec[edgeIndex]; // we are at end of mesh edge
 break;
 }
 if( tt > 1 )
 {
 tt = 1;
 v = nodesConnec[edgeIndex + 1]; // we are at one end
 break;
 }
 iterations++;
 }
 // so we have t on the segment, convert to u, which should
 // be between _fractions[edgeIndex] numbers
 double prevFraction = 0;
 if( edgeIndex > 0 ) prevFraction = _fractions[edgeIndex - 1];
 
 u = prevFraction + tt * ( _fractions[edgeIndex] - prevFraction );
 return u;
}
 
//! \brief Get the starting point of the curve.
//!
//! \param x The x coordinate of the start point
//! \param y The y coordinate of the start point
//! \param z The z coordinate of the start point
void SmoothCurve::start_coordinates( double& x, double& y, double& z )
{
 
 int nnodes = 0;
 const EntityHandle* conn2 = NULL;
 _mb->get_connectivity( _entities[0], conn2, nnodes );
 double c[3];
 _mb->get_coords( conn2, 1, c );
 
 x = c[0];
 y = c[1];
 z = c[2];
 
 return;
}
 
//! \brief Get the ending point of the curve.
//!
//! \param x The x coordinate of the start point
//! \param y The y coordinate of the start point
//! \param z The z coordinate of the start point
void SmoothCurve::end_coordinates( double& x, double& y, double& z )
{
 
 int nnodes = 0;
 const EntityHandle* conn2 = NULL;
 _mb->get_connectivity( _entities[_entities.size() - 1], conn2, nnodes );
 double c[3];
 // careful, the second node here
 _mb->get_coords( &conn2[1], 1, c );
 
 x = c[0];
 y = c[1];
 z = c[2];
 return;
}
 
// this will recompute the 2 tangents for each edge, considering the geo edge they are into
void SmoothCurve::compute_tangents_for_each_edge()
{
 // will retrieve the edges in each set; they are retrieved in order they were put into, because
 // these sets are "MESHSET_ORDERED"
 // retrieve the tag handle for the tangents; it should have been created already
 // this tangents are computed for the chain of edges that form a geometric edge
 // some checks should be performed on the vertices, but we trust the correctness of the model
 // completely (like the vertices should match in the chain...)
 Tag tangentsTag;
 ErrorCode rval = _mb->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tangentsTag );
 if( rval != MB_SUCCESS ) return; // some error should be thrown
 std::vector< EntityHandle > entities;
 _mb->get_entities_by_type( _set, MBEDGE, entities ); // no recursion!!
 // basically, each tangent at a node will depend on previous tangent
 int nbEdges = entities.size();
 // now, we can advance in the loop
 // the only special problem is if the first node coincides with the last node, then we should
 // consider the closed loop; or maybe we should look at angles in that case too?
 // also, should we look at the 2 semi-circles case? How to decide if we need to continue the
 // "tangents" maybe we can do that later, and we can alter the tangents at the feature nodes, in
 // the directions of the loops again, do we need to decide the "closed" loop or not? Not yet...
 EntityHandle previousEdge = entities[0]; // this is the first edge in the chain
 CartVect TP[2]; // tangents for the previous edge
 rval = _mb->tag_get_data( tangentsTag, &previousEdge, 1, &TP[0] ); // tangents for previous edge
 if( rval != MB_SUCCESS ) return; // some error should be thrown
 CartVect TC[2]; // tangents for the current edge
 EntityHandle currentEdge;
 for( int i = 1; i < nbEdges; i++ )
 {
 // current edge will start after first one
 currentEdge = entities[i];
 rval = _mb->tag_get_data( tangentsTag, &currentEdge, 1, &TC[0] ); //
 if( rval != MB_SUCCESS ) return; // some error should be thrown
 // now compute the new tangent at common vertex; reset tangents for previous edge and
 // current edge a little bit of CPU and memory waste, but this is life
 CartVect T = 0.5 * TC[0] + 0.5 * TP[1]; //
 T.normalize();
 TP[1] = T;
 rval = _mb->tag_set_data( tangentsTag, &previousEdge, 1, &TP[0] ); //
 if( rval != MB_SUCCESS ) return; // some error should be thrown
 TC[0] = T;
 rval = _mb->tag_set_data( tangentsTag, &currentEdge, 1, &TC[0] ); //
 if( rval != MB_SUCCESS ) return; // some error should be thrown
 // now set the next edge
 previousEdge = currentEdge;
 TP[0] = TC[0];
 TP[1] = TC[1];
 }
 return;
}
 
void SmoothCurve::compute_control_points_on_boundary_edges( double,
 std::map< EntityHandle, SmoothFace* >& mapSurfaces,
 Tag controlPointsTag,
 Tag markTag )
{
 // these points really need the surfaces they belong to, because the control points on edges
 // depend on the normals on surfaces
 // the control points are averaged from different surfaces, by simple mean average
 // the surfaces have
 // do we really need min_dot here?
 // first of all, find out the SmoothFace for each surface set that is adjacent here
 // GeomTopoTool gTopoTool(_mb);
 std::vector< EntityHandle > faces;
 std::vector< int > senses;
 ErrorCode rval = _gtt->get_senses( _set, faces, senses );
 if( MB_SUCCESS != rval ) return;
 
 // need to find the smooth face attached
 unsigned int numSurfacesAdjacent = faces.size();
 // get the edges, and then get the
 // std::vector<EntityHandle> entities;
 _mb->get_entities_by_type( _set, MBEDGE, _entities ); // no recursion!!
 // each edge has the tangent computed already
 Tag tangentsTag;
 rval = _mb->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tangentsTag );
 if( rval != MB_SUCCESS ) return; // some error should be thrown
 
 // we do not want to search every time
 std::vector< SmoothFace* > smoothFaceArray;
 unsigned int i = 0;
 for( i = 0; i < numSurfacesAdjacent; i++ )
 {
 SmoothFace* sms = mapSurfaces[faces[i]];
 smoothFaceArray.push_back( sms );
 }
 
 unsigned int e = 0;
 for( e = 0; e < _entities.size(); e++ )
 {
 CartVect zero( 0. );
 CartVect ctrlP[3] = { zero, zero, zero }; // null positions initially
 // the control points are averaged from connected faces
 EntityHandle edge = _entities[e]; // the edge in the chain
 
 int nnodes;
 const EntityHandle* conn2;
 rval = _mb->get_connectivity( edge, conn2, nnodes );
 if( rval != MB_SUCCESS || 2 != nnodes ) return; // or continue or return error
 
 // 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] );
 if( rval != MB_SUCCESS ) return;
 
 CartVect chord = P[1] - P[0];
 _leng += chord.length();
 _fractions.push_back( _leng );
 CartVect N[2];
 
 // MBCartVect N0(&normalVec[0]);
 // MBCartVect N3(&normalVec[3]);
 CartVect T[2]; // T0, T3
 // if (edge->num_adj_facets() <= 1) {
 // stat = compute_curve_tangent(edge, min_dot, T0, T3);
 // if (stat != CUBIT_SUCCESS)
 // return stat;
 //} else {
 //}
 rval = _mb->tag_get_data( tangentsTag, &edge, 1, &T[0] );
 if( rval != MB_SUCCESS ) return;
 
 for( i = 0; i < numSurfacesAdjacent; i++ )
 {
 CartVect controlForEdge[3];
 rval = smoothFaceArray[i]->get_normals_for_vertices( conn2, N );
 if( rval != MB_SUCCESS ) return;
 
 rval = smoothFaceArray[i]->init_edge_control_points( P[0], P[1], N[0], N[1], T[0], T[1], controlForEdge );
 if( rval != MB_SUCCESS ) return;
 
 // accumulate those over faces!!!
 for( int j = 0; j < 3; j++ )
 {
 ctrlP[j] += controlForEdge[j];
 }
 }
 // now divide them for the average position!
 for( int j = 0; j < 3; j++ )
 {
 ctrlP[j] /= numSurfacesAdjacent;
 }
 // we are done, set the control points now!
 // edge->control_points(ctrl_pts, 4);
 rval = _mb->tag_set_data( controlPointsTag, &edge, 1, &ctrlP[0] );
 if( rval != MB_SUCCESS ) return;
 
 this->_edgeTag = controlPointsTag; // this is a tag that will be stored with the edge
 // is that a waste of memory or not...
 // also mark the edge for later on
 unsigned char used = 1;
 _mb->tag_set_data( markTag, &edge, 1, &used );
 }
 // now divide fractions, to make them vary from 0 to 1
 assert( _leng > 0. );
 for( e = 0; e < _entities.size(); e++ )
 _fractions[e] /= _leng;
}
 
ErrorCode SmoothCurve::evaluate_smooth_edge( EntityHandle eh, double& tt, CartVect& outv, CartVect& out_tangent )
{
 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 );
 if( rval != MB_SUCCESS ) return rval;
 
 rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
 if( rval != MB_SUCCESS ) return rval;
 
 if( 0 == _edgeTag )
 {
 rval = _mb->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, _edgeTag );
 if( rval != MB_SUCCESS ) return rval;
 }
 rval = _mb->tag_get_data( _edgeTag, &eh, 1, (double*)&controlPoints[0] );
 if( rval != MB_SUCCESS ) 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];
 
 out_tangent = -4. * one_minus_t3 * P[0] + 4. * ( one_minus_t3 - 3. * tt * one_minus_t2 ) * controlPoints[0] +
 12. * ( tt * one_minus_t2 - t2 * one_minus_t ) * controlPoints[1] +
 4. * ( 3. * t2 * one_minus_t - t3 ) * controlPoints[2] + 4. * t3 * P[1];
 return MB_SUCCESS;
}
} // namespace moab