LinearTri.cpp

Go to the documentation of this file.

#include "moab/LocalDiscretization/LinearTri.hpp"
#include "moab/Forward.hpp"
#include <algorithm>
#include <cmath>
#include <limits>
 
namespace moab
{
 
const double LinearTri::corner[3][2] = { { 0, 0 }, { 1, 0 }, { 0, 1 } };
 
ErrorCode LinearTri::initFcn( const double* verts, const int nverts, double*& work )
{
 // 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;
}
 
ErrorCode LinearTri::evalFcn( const double* params,
 const double* field,
 const int /*ndim*/,
 const int num_tuples,
 double* /*work*/,
 double* result )
{
 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 LinearTri::integrateFcn( const double* field,
 const double* /*verts*/,
 const int nverts,
 const int /*ndim*/,
 const int num_tuples,
 double* work,
 double* result )
{
 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 LinearTri::jacobianFcn( const double*, const double*, const int, const int, double* work, double* result )
{
 // jacobian is cached in work array
 assert( work );
 std::copy( work, work + 9, result );
 return MB_SUCCESS;
}
 
ErrorCode 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 )
{
 assert( posn && verts );
 return evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
 is_inside );
}
 
int LinearTri::insideFcn( const double* params, const int, const double tol )
{
 return ( params[0] >= -1.0 - tol && params[1] >= -1.0 - tol && params[0] + params[1] <= 1.0 + tol );
}
 
ErrorCode 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 )
{
 // 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 LinearTri::get_normal( int facet, double *work, double *normal)
 {
 ErrorCode error;
 //Get the local vertex ids of local edge
 int id1 = ledges[facet][0];
 int id2 = ledges[facet][1];
 
 //Find the normal to the face
 double face_normal[3];
 
 
 }*/
 
ErrorCode LinearTri::normalFcn( const int ientDim,
 const int facet,
 const int nverts,
 const double* verts,
 double normal[3] )
{
 // 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;
}
 
} // namespace moab