#include <LinearTri.hpp>

Static Public Member Functions
static ErrorCode	evalFcn (const double params, const double field, const int ndim, const int num_tuples, double work, double result)
	Forward-evaluation of field at parametric coordinates.
static ErrorCode	reverseEvalFcn (EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double posn, const double verts, const int nverts, const int ndim, const double iter_tol, const double inside_tol, double work, double params, int *is_inside)
	Reverse-evaluation of parametric coordinates at physical space position.
static ErrorCode	normalFcn (const int ientDim, const int facet, const int nverts, const double *verts, double normal[])
	Evaluate the normal at a specified facet.
static ErrorCode	jacobianFcn (const double params, const double verts, const int nverts, const int ndim, double work, double result)
	Evaluate the jacobian at a specified parametric position.
static ErrorCode	integrateFcn (const double field, const double verts, const int nverts, const int ndim, const int num_tuples, double work, double result)
	Forward-evaluation of field at parametric coordinates.
static ErrorCode	initFcn (const double verts, const int nverts, double &work)
	Initialize this EvalSet.
static int	insideFcn (const double *params, const int ndim, const double tol)
	Function that returns whether or not the parameters are inside the natural space of the element.
static ErrorCode	evaluate_reverse (EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f, const double posn, const double verts, const int nverts, const int ndim, const double iter_tol, const double inside_tol, double work, double params, int *inside)
static EvalSet	eval_set ()
static bool	compatible (EntityType tp, int numv, EvalSet &eset)
Static Protected Attributes
static const double	corner [3][2] = { { 0, 0 }, { 1, 0 }, { 0, 1 } }

Detailed Description

Definition at line 12 of file LinearTri.hpp.

Member Function Documentation

static bool moab::LinearTri::compatible	(	EntityType	tp,
		int	numv,
		EvalSet &	eset
	)		`[inline, static]`

Definition at line 86 of file LinearTri.hpp.

References eval_set(), and MBTRI.

Referenced by moab::EvalSet::get_eval_set().

    {
        if( tp == MBTRI && numv >= 3 )
        {
            eset = eval_set();
            return true;
        }
        else
            return false;
    }

static EvalSet moab::LinearTri::eval_set ( ) [inline, static]

Definition at line 81 of file LinearTri.hpp.

Referenced by compatible().

    {
        return EvalSet( evalFcn, reverseEvalFcn, normalFcn, jacobianFcn, integrateFcn, initFcn, insideFcn );
    }

ErrorCode moab::LinearTri::evalFcn	(	const double *	params,
		const double *	field,
		const int	ndim,
		const int	num_tuples,
		double *	work,
		double *	result
	)		`[static]`

Forward-evaluation of field at parametric coordinates.

Definition at line 43 of file LinearTri.cpp.

References MB_SUCCESS.

{
    assert( params && field && num_tuples > 0 );
    // convert to [0,1]
    double p1 = 0.5 * ( 1.0 + params[0] ), p2 = 0.5 * ( 1.0 + params[1] ), p0 = 1.0 - p1 - p2;

    for( int j = 0; j < num_tuples; j++ )
        result[j] = p0 * field[0 * num_tuples + j] + p1 * field[1 * num_tuples + j] + p2 * field[2 * num_tuples + j];

    return MB_SUCCESS;
}

ErrorCode moab::LinearTri::evaluate_reverse	(	EvalFcn	eval,
		JacobianFcn	jacob,
		InsideFcn	inside_f,
		const double *	posn,
		const double *	verts,
		const int	nverts,
		const int	ndim,
		const double	iter_tol,
		const double	inside_tol,
		double *	work,
		double *	params,
		int *	inside
	)		`[static]`

Definition at line 113 of file LinearTri.cpp.

References moab::CartVect::array(), moab::Matrix3::determinant(), ErrorCode, moab::Matrix3::inverse(), length_squared(), and MB_SUCCESS.

Referenced by reverseEvalFcn().

{
    // TODO: should differentiate between epsilons used for
    // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
    // right now, fix the tolerance used for NR
    const double error_tol_sqr = iter_tol * iter_tol;
    CartVect* cvparams         = reinterpret_cast< CartVect* >( params );
    const CartVect* cvposn     = reinterpret_cast< const CartVect* >( posn );

    // find best initial guess to improve convergence
    CartVect tmp_params[] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( -1, 1, -1 ) };
    double resl           = std::numeric_limits< double >::max();
    CartVect new_pos, tmp_pos;
    ErrorCode rval;
    for( unsigned int i = 0; i < 3; i++ )
    {
        rval = ( *eval )( tmp_params[i].array(), verts, ndim, 3, work, tmp_pos.array() );
        if( MB_SUCCESS != rval ) return rval;
        double tmp_resl = ( tmp_pos - *cvposn ).length_squared();
        if( tmp_resl < resl )
        {
            *cvparams = tmp_params[i];
            new_pos   = tmp_pos;
            resl      = tmp_resl;
        }
    }

    // residual is diff between old and new pos; need to minimize that
    CartVect res = new_pos - *cvposn;
    Matrix3 J;
    rval = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J[0] );
#ifndef NDEBUG
    double det = J.determinant();
    assert( det > std::numeric_limits< double >::epsilon() );
#endif
    Matrix3 Ji = J.inverse();

    int iters = 0;
    // while |res| larger than tol
    while( res % res > error_tol_sqr )
    {
        if( ++iters > 25 ) return MB_FAILURE;

        // new params tries to eliminate residual
        *cvparams -= Ji * res;

        // get the new forward-evaluated position, and its difference from the target pt
        rval = ( *eval )( params, verts, ndim, 3, work, new_pos.array() );
        if( MB_SUCCESS != rval ) return rval;
        res = new_pos - *cvposn;
    }

    if( inside ) *inside = ( *inside_f )( params, ndim, inside_tol );

    return MB_SUCCESS;
}  // Map::evaluate_reverse()

ErrorCode moab::LinearTri::initFcn	(	const double *	verts,
		const int	nverts,
		double *&	work
	)		`[static]`

Initialize this EvalSet.

Definition at line 12 of file LinearTri.cpp.

References moab::Matrix3::copyto(), moab::Matrix3::determinant(), moab::Matrix3::inverse(), MB_SUCCESS, and moab::Matrix3::size.

{
    // allocate work array as:
    // work[0..8] = T
    // work[9..17] = Tinv
    // work[18] = detT
    // work[19] = detTinv
    if( nverts != 3 )
    {
        std::cout << "Invalid Triangle. Expected 3 vertices.\n";
        return MB_FAILURE;
    }

    assert( verts );

    Matrix3 J( verts[1 * 3 + 0] - verts[0 * 3 + 0], verts[2 * 3 + 0] - verts[0 * 3 + 0], 0.0,
               verts[1 * 3 + 1] - verts[0 * 3 + 1], verts[2 * 3 + 1] - verts[0 * 3 + 1], 0.0,
               verts[1 * 3 + 2] - verts[0 * 3 + 2], verts[2 * 3 + 2] - verts[0 * 3 + 2], 1.0 );
    J *= 0.5;

    // Update the work array
    if( !work ) work = new double[20];

    J.copyto( work );
    J.inverse().copyto( work + Matrix3::size );
    work[18] = J.determinant();
    work[19] = ( work[18] < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / work[18] );

    return MB_SUCCESS;
}

int moab::LinearTri::insideFcn	(	const double *	params,
		const int	ndim,
		const double	tol
	)		`[static]`

Function that returns whether or not the parameters are inside the natural space of the element.

Definition at line 108 of file LinearTri.cpp.

{
    return ( params[0] >= -1.0 - tol && params[1] >= -1.0 - tol && params[0] + params[1] <= 1.0 + tol );
}

ErrorCode moab::LinearTri::integrateFcn	(	const double *	field,
		const double *	verts,
		const int	nverts,
		const int	ndim,
		const int	num_tuples,
		double *	work,
		double *	result
	)		`[static]`

Forward-evaluation of field at parametric coordinates.

Definition at line 60 of file LinearTri.cpp.

References MB_SUCCESS.

{
    assert( field && num_tuples > 0 );
    std::fill( result, result + num_tuples, 0.0 );
    for( int i = 0; i < nverts; ++i )
    {
        for( int j = 0; j < num_tuples; j++ )
            result[j] += field[i * num_tuples + j];
    }
    double tmp = work[18] / 6.0;
    for( int i = 0; i < num_tuples; i++ )
        result[i] *= tmp;

    return MB_SUCCESS;
}

ErrorCode moab::LinearTri::jacobianFcn	(	const double *	params,
		const double *	verts,
		const int	nverts,
		const int	ndim,
		double *	work,
		double *	result
	)		`[static]`

Evaluate the jacobian at a specified parametric position.

Definition at line 82 of file LinearTri.cpp.

References MB_SUCCESS.

{
    // jacobian is cached in work array
    assert( work );
    std::copy( work, work + 9, result );
    return MB_SUCCESS;
}

ErrorCode moab::LinearTri::normalFcn	(	const int	ientDim,
		const int	facet,
		const int	nverts,
		const double *	verts,
		double	normal[]
	)		`[static]`

Evaluate the normal at a specified facet.

Definition at line 194 of file LinearTri.cpp.

References moab::CN::ConnMap::conn, MB_SET_ERR, MB_SUCCESS, MBTRI, and moab::CN::mConnectivityMap.

{
    // assert(facet < 3 && ientDim == 1 && nverts==3);
    if( nverts != 3 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed triangle :: expected value = 3 " );
    if( ientDim != 1 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 1 " );
    if( facet > 3 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local edge id :: expected value = one of 0-2" );

    // Get the local vertex ids of  local edge
    int id0 = CN::mConnectivityMap[MBTRI][ientDim - 1].conn[facet][0];
    int id1 = CN::mConnectivityMap[MBTRI][ientDim - 1].conn[facet][1];

    // Find a vector along the edge
    double edge[3];
    for( int i = 0; i < 3; i++ )
    {
        edge[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
    }
    // Find the normal of the face
    double x0[3], x1[3], fnrm[3];
    for( int i = 0; i < 3; i++ )
    {
        x0[i] = verts[3 * 1 + i] - verts[3 * 0 + i];
        x1[i] = verts[3 * 2 + i] - verts[3 * 0 + i];
    }
    fnrm[0] = x0[1] * x1[2] - x1[1] * x0[2];
    fnrm[1] = x1[0] * x0[2] - x0[0] * x1[2];
    fnrm[2] = x0[0] * x1[1] - x1[0] * x0[1];

    // Find the normal of the edge as the cross product of edge and face normal

    double a   = edge[1] * fnrm[2] - fnrm[1] * edge[2];
    double b   = edge[2] * fnrm[0] - fnrm[2] * edge[0];
    double c   = edge[0] * fnrm[1] - fnrm[0] * edge[1];
    double nrm = sqrt( a * a + b * b + c * c );

    if( nrm > std::numeric_limits< double >::epsilon() )
    {
        normal[0] = a / nrm;
        normal[1] = b / nrm;
        normal[2] = c / nrm;
    }
    return MB_SUCCESS;
}

ErrorCode moab::LinearTri::reverseEvalFcn	(	EvalFcn	eval,
		JacobianFcn	jacob,
		InsideFcn	ins,
		const double *	posn,
		const double *	verts,
		const int	nverts,
		const int	ndim,
		const double	iter_tol,
		const double	inside_tol,
		double *	work,
		double *	params,
		int *	is_inside
	)		`[static]`

Reverse-evaluation of parametric coordinates at physical space position.

Definition at line 90 of file LinearTri.cpp.

References evaluate_reverse().

{
    assert( posn && verts );
    return evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
                             is_inside );
}

Member Data Documentation

const double moab::LinearTri::corner = { { 0, 0 }, { 1, 0 }, { 0, 1 } } [static, protected]

Definition at line 98 of file LinearTri.hpp.

List of all members.

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

Static Public Member Functions

Static Protected Attributes

Detailed Description

Member Function Documentation

Member Data Documentation