#include <SmoothCurve.hpp>

Collaboration diagram for moab::SmoothCurve:

Public Member Functions
	SmoothCurve (Interface mb, EntityHandle curve, GeomTopoTool gTool)
virtual	~SmoothCurve ()
virtual double	arc_length ()
virtual bool	is_parametric ()
	Get the parametric status of the curve.
virtual bool	is_periodic (double &period)
	Get the periodic status of the curve.
virtual void	get_param_range (double &u_start, double &u_end)
	Get the parameter range of the curve.
virtual double	u_from_arc_length (double u_root, double arc_length)
	Compute the parameter value at a specified distance along the curve.
virtual bool	position_from_u (double u, double &x, double &y, double &z, double *tg=NULL)
	Evaluate the curve at a specified parameter value.
virtual void	move_to_curve (double &x, double &y, double &z)
	Move a point near the curve to the closest point on the curve.
virtual double	u_from_position (double x, double y, double z, EntityHandle &v, int &indexEdge)
	Get the u parameter value on the curve closest to x,y,z and the point on the curve.
virtual void	start_coordinates (double &x, double &y, double &z)
	Get the starting point of the curve.
virtual void	end_coordinates (double &x, double &y, double &z)
	Get the ending point of the curve.
void	compute_tangents_for_each_edge ()
void	compute_control_points_on_boundary_edges (double min_dot, std::map< EntityHandle, SmoothFace * > &mapSurfaces, Tag controlPointsTag, Tag markTag)
ErrorCode	evaluate_smooth_edge (EntityHandle eh, double &tt, CartVect &outv, CartVect &out_tangent)
Private Attributes
std::vector< EntityHandle >	_entities
double	_leng
std::vector< double >	_fractions
Tag	_edgeTag
Interface *	_mb
EntityHandle	_set
GeomTopoTool *	_gtt

Detailed Description

Definition at line 29 of file SmoothCurve.hpp.

Constructor & Destructor Documentation

moab::SmoothCurve::SmoothCurve	(	Interface *	mb,
		EntityHandle	curve,
		GeomTopoTool *	gTool
	)

Definition at line 16 of file SmoothCurve.cpp.

References _edgeTag, and _leng.

    : _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
}

moab::SmoothCurve::~SmoothCurve ( ) [virtual]

Definition at line 25 of file SmoothCurve.cpp.

{
    // TODO Auto-generated destructor stub
}

Member Function Documentation

double moab::SmoothCurve::arc_length ( ) [virtual]

Definition at line 30 of file SmoothCurve.cpp.

References _leng.

{

    return _leng;
}

void moab::SmoothCurve::compute_control_points_on_boundary_edges	(	double	min_dot,
		std::map< EntityHandle, SmoothFace * > &	mapSurfaces,
		Tag	controlPointsTag,
		Tag	markTag
	)

Definition at line 428 of file SmoothCurve.cpp.

References _edgeTag, _entities, _fractions, _gtt, _leng, _mb, _set, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), moab::GeomTopoTool::get_senses(), moab::CartVect::length(), MB_SUCCESS, MB_TYPE_DOUBLE, MBEDGE, T, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), and moab::Interface::tag_set_data().

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

{
    // 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;
}

void moab::SmoothCurve::compute_tangents_for_each_edge ( )

Definition at line 377 of file SmoothCurve.cpp.

References _mb, _set, entities, ErrorCode, moab::Interface::get_entities_by_type(), MB_SUCCESS, MB_TYPE_DOUBLE, MBEDGE, moab::CartVect::normalize(), T, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), and TC.

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

{
    // 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 moab::SmoothCurve::end_coordinates	(	double &	x,
		double &	y,
		double &	z
	)		`[virtual]`

Get the ending point of the curve.

Parameters:

x	The x coordinate of the start point
y	The y coordinate of the start point
z	The z coordinate of the start point

Definition at line 360 of file SmoothCurve.cpp.

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

{

    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;
}

ErrorCode moab::SmoothCurve::evaluate_smooth_edge	(	EntityHandle	eh,
		double &	tt,
		CartVect &	outv,
		CartVect &	out_tangent
	)

Definition at line 537 of file SmoothCurve.cpp.

References _edgeTag, _mb, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), MB_SUCCESS, MB_TYPE_DOUBLE, moab::Interface::tag_get_data(), and moab::Interface::tag_get_handle().

Referenced by position_from_u().

{
    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;
}

void moab::SmoothCurve::get_param_range	(	double &	u_start,
		double &	u_end
	)		`[virtual]`

Get the parameter range of the curve.

Parameters:

u_start	The beginning curve parameter
u_end	The ending curve parameter

Note:: The numerical value of u_start may be greater than the numerical value of u_end.

Definition at line 70 of file SmoothCurve.cpp.

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

{
    // assert(_ref_edge);
    u_start = 0;
    u_end   = 1.;

    return;
}

bool moab::SmoothCurve::is_parametric ( ) [virtual]

Get the parametric status of the curve.

Returns:: true if curve is parametric, false otherwise.

Definition at line 39 of file SmoothCurve.cpp.

{
    return true;
}

bool moab::SmoothCurve::is_periodic ( double & period ) [virtual]

Get the periodic status of the curve.

Parameters:

period The period of the curve if periodic.

Returns:: true if curve is periodic, false otherwise.

Definition at line 49 of file SmoothCurve.cpp.

References _leng, _mb, _set, moab::Interface::get_child_meshsets(), and moab::Range::size().

{
    // 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;
}

void moab::SmoothCurve::move_to_curve	(	double &	x,
		double &	y,
		double &	z
	)		`[virtual]`

Move a point near the curve to the closest point on the curve.

Parameters:

x	The x coordinate of the point
y	The y coordinate of the point
z	The z coordinate of the point

Definition at line 144 of file SmoothCurve.cpp.

References position_from_u(), and u_from_position().

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

{

    // 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;
}

bool moab::SmoothCurve::position_from_u	(	double	u,
		double &	x,
		double &	y,
		double &	z,
		double *	tg = `NULL`
	)		`[virtual]`

Evaluate the curve at a specified parameter value.

Parameters:

u	The parameter at which to evaluate the curve
x	The x coordinate of the evaluated point
y	The y coordinate of the evaluated point
z	The z coordinate of the evaluated point

Definition at line 105 of file SmoothCurve.cpp.

References _entities, _fractions, ErrorCode, evaluate_smooth_edge(), MB_SUCCESS, and t.

Referenced by moab::FBEngine::getEntTgntU(), moab::FBEngine::getEntUtoXYZ(), and move_to_curve().

{

    // _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;
}

void moab::SmoothCurve::start_coordinates	(	double &	x,
		double &	y,
		double &	z
	)		`[virtual]`

Get the starting point of the curve.

Parameters:

x	The x coordinate of the start point
y	The y coordinate of the start point
z	The z coordinate of the start point

Definition at line 339 of file SmoothCurve.cpp.

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

{

    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;
}

double moab::SmoothCurve::u_from_arc_length	(	double	u_root,
		double	arc_leng
	)		`[virtual]`

Compute the parameter value at a specified distance along the curve.

Parameters:

u_root	The start parameter from which to compute the distance along the curve.
arc_length	The distance to move along the curve.

Note:: For positive values of arc_length the distance will be computed in the direction of increasing parameter value along the curve. For negative values of arc_length the distance will be computed in the direction of decreasing parameter value along the curve.

Returns:: The parametric coordinate u along the curve

Definition at line 92 of file SmoothCurve.cpp.

References _leng.

{

    if( _leng <= 0 ) return 0;
    return u_root + arc_leng / _leng;
}

double moab::SmoothCurve::u_from_position	(	double	x,
		double	y,
		double	z,
		EntityHandle &	v,
		int &	edgeIndex
	)		`[virtual]`

Get the u parameter value on the curve closest to x,y,z and the point on the curve.

Parameters:

x	The x coordinate of the point
y	The y coordinate of the point
z	The z coordinate of the point

Returns:: The parametric coordinate u on the curve

Definition at line 165 of file SmoothCurve.cpp.

References _edgeTag, _entities, _fractions, _mb, delta_t, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), length(), length_squared(), moab::CartVect::length_squared(), MB_SUCCESS, MB_TYPE_DOUBLE, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), and moab::tolerance.

Referenced by move_to_curve(), and moab::FBEngine::split_edge_at_point().

{
    // 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;
}