ElemEvaluator.cpp

Go to the documentation of this file.

#include <limits>
 
#include "moab/ElemEvaluator.hpp"
#include "moab/CartVect.hpp"
#include "moab/Matrix3.hpp"
 
// need to include eval set types here to support get_eval_set; alternative would be to have some
// type of registration, but we'd still need static registration for the built-in types
#include "moab/LocalDiscretization/LinearTri.hpp"
#include "moab/LocalDiscretization/LinearQuad.hpp"
#include "moab/LocalDiscretization/LinearTet.hpp"
#include "moab/LocalDiscretization/LinearHex.hpp"
#include "moab/LocalDiscretization/QuadraticHex.hpp"
//#include "moab/SpectralQuad.hpp"
//#include "moab/SpectralHex.hpp"
 
namespace moab
{
ErrorCode EvalSet::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 );
 
 // initialize to center of element
 *cvparams = CartVect( -.4 );
 
 CartVect new_pos;
 // evaluate that first guess to get a new position
 ErrorCode rval = ( *eval )( cvparams->array(), verts, ndim,
 3, // hardwire to num_tuples to 3 since the field is coords
 work, new_pos.array() );
 if( MB_SUCCESS != rval ) return rval;
 
 // residual is diff between old and new pos; need to minimize that
 CartVect res = new_pos - *cvposn;
 Matrix3 J;
 int dum, *tmp_inside = ( inside ? inside : &dum );
 
 int iters = 0;
 // while |res| larger than tol
 while( res % res > error_tol_sqr )
 {
 if( ++iters > 10 )
 {
 // if we haven't converged but we're outside, that's defined as success
 *tmp_inside = ( *inside_f )( params, ndim, inside_tol );
 if( !( *tmp_inside ) )
 return MB_SUCCESS;
 else
 return MB_FAILURE;
 }
 
 // get jacobian at current params
 rval = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J.array() );
 double det = J.determinant();
 if( det < std::numeric_limits< double >::epsilon() )
 {
 *tmp_inside = ( *inside_f )( params, ndim, inside_tol );
 if( !( *tmp_inside ) )
 return MB_SUCCESS;
 else
 return MB_INDEX_OUT_OF_RANGE;
 }
 
 // new params tries to eliminate residual
 *cvparams -= J.inverse() * res;
 
 // get the new forward-evaluated position, and its difference from the target pt
 rval = ( *eval )( params, verts, ndim,
 3, // hardwire to num_tuples to 3 since the field is coords
 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()
 
int EvalSet::inside_function( const double* params, const int ndims, const double tol )
{
 if( params[0] >= -1 - tol && params[0] <= 1 + tol &&
 ( ndims < 2 || ( params[1] >= -1 - tol && params[1] <= 1 + tol ) ) &&
 ( ndims < 3 || ( params[2] >= -1 - tol && params[2] <= 1 + tol ) ) )
 return true;
 else
 return false;
}
 
/** \brief Given type & #vertices, get an appropriate eval set */
ErrorCode EvalSet::get_eval_set( EntityType tp, unsigned int num_vertices, EvalSet& eval_set )
{
 switch( tp )
 {
 case MBEDGE:
 break;
 case MBTRI:
 if( LinearTri::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
 break;
 case MBQUAD:
 if( LinearQuad::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
 // if (SpectralQuad::compatible(tp, num_vertices, eval_set)) return
 // MB_SUCCESS;
 break;
 case MBTET:
 if( LinearTet::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
 break;
 case MBHEX:
 if( LinearHex::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
 if( QuadraticHex::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
 // if (SpectralHex::compatible(tp, num_vertices, eval_set)) return
 // MB_SUCCESS;
 break;
 default:
 break;
 }
 
 return MB_NOT_IMPLEMENTED;
}
 
ErrorCode ElemEvaluator::find_containing_entity( Range& entities,
 const double* point,
 const double iter_tol,
 const double inside_tol,
 EntityHandle& containing_ent,
 double* params,
 unsigned int* num_evals )
{
 int is_inside;
 ErrorCode rval = MB_SUCCESS;
 unsigned int nevals = 0;
 Range::iterator i;
 for( i = entities.begin(); i != entities.end(); ++i )
 {
 nevals++;
 set_ent_handle( *i );
 rval = reverse_eval( point, iter_tol, inside_tol, params, &is_inside );
 if( MB_SUCCESS != rval ) return rval;
 if( is_inside ) break;
 }
 containing_ent = ( i == entities.end() ? 0 : *i );
 if( num_evals ) *num_evals += nevals;
 return MB_SUCCESS;
}
} // namespace moab