MOAB: Mesh Oriented datABase  (version 5.3.0)
ElemEvaluator.cpp
Go to the documentation of this file.
00001 #include <limits>
00002 
00003 #include "moab/ElemEvaluator.hpp"
00004 #include "moab/CartVect.hpp"
00005 #include "moab/Matrix3.hpp"
00006 
00007 // need to include eval set types here to support get_eval_set; alternative would be to have some
00008 // type of registration, but we'd still need static registration for the built-in types
00009 #include "moab/LocalDiscretization/LinearTri.hpp"
00010 #include "moab/LocalDiscretization/LinearQuad.hpp"
00011 #include "moab/LocalDiscretization/LinearTet.hpp"
00012 #include "moab/LocalDiscretization/LinearHex.hpp"
00013 #include "moab/LocalDiscretization/QuadraticHex.hpp"
00014 //#include "moab/SpectralQuad.hpp"
00015 //#include "moab/SpectralHex.hpp"
00016 
00017 namespace moab
00018 {
00019 ErrorCode EvalSet::evaluate_reverse( EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f, const double* posn,
00020                                      const double* verts, const int nverts, const int ndim, const double iter_tol,
00021                                      const double inside_tol, double* work, double* params, int* inside )
00022 {
00023     // TODO: should differentiate between epsilons used for
00024     // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
00025     // right now, fix the tolerance used for NR
00026     const double error_tol_sqr = iter_tol * iter_tol;
00027     CartVect* cvparams         = reinterpret_cast< CartVect* >( params );
00028     const CartVect* cvposn     = reinterpret_cast< const CartVect* >( posn );
00029 
00030     // initialize to center of element
00031     *cvparams = CartVect( -.4 );
00032 
00033     CartVect new_pos;
00034     // evaluate that first guess to get a new position
00035     ErrorCode rval = ( *eval )( cvparams->array(), verts, ndim,
00036                                 3,  // hardwire to num_tuples to 3 since the field is coords
00037                                 work, new_pos.array() );
00038     if( MB_SUCCESS != rval ) return rval;
00039 
00040     // residual is diff between old and new pos; need to minimize that
00041     CartVect res = new_pos - *cvposn;
00042     Matrix3 J;
00043     int dum, *tmp_inside = ( inside ? inside : &dum );
00044 
00045     int iters = 0;
00046     // while |res| larger than tol
00047     while( res % res > error_tol_sqr )
00048     {
00049         if( ++iters > 10 )
00050         {
00051             // if we haven't converged but we're outside, that's defined as success
00052             *tmp_inside = ( *inside_f )( params, ndim, inside_tol );
00053             if( !( *tmp_inside ) )
00054                 return MB_SUCCESS;
00055             else
00056                 return MB_FAILURE;
00057         }
00058 
00059         // get jacobian at current params
00060         rval       = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J.array() );
00061         double det = J.determinant();
00062         if( det < std::numeric_limits< double >::epsilon() )
00063         {
00064             *tmp_inside = ( *inside_f )( params, ndim, inside_tol );
00065             if( !( *tmp_inside ) )
00066                 return MB_SUCCESS;
00067             else
00068                 return MB_INDEX_OUT_OF_RANGE;
00069         }
00070 
00071         // new params tries to eliminate residual
00072         *cvparams -= J.inverse() * res;
00073 
00074         // get the new forward-evaluated position, and its difference from the target pt
00075         rval = ( *eval )( params, verts, ndim,
00076                           3,  // hardwire to num_tuples to 3 since the field is coords
00077                           work, new_pos.array() );
00078         if( MB_SUCCESS != rval ) return rval;
00079         res = new_pos - *cvposn;
00080     }
00081 
00082     if( inside ) *inside = ( *inside_f )( params, ndim, inside_tol );
00083 
00084     return MB_SUCCESS;
00085 }  // Map::evaluate_reverse()
00086 
00087 int EvalSet::inside_function( const double* params, const int ndims, const double tol )
00088 {
00089     if( params[0] >= -1 - tol && params[0] <= 1 + tol &&
00090         ( ndims < 2 || ( params[1] >= -1 - tol && params[1] <= 1 + tol ) ) &&
00091         ( ndims < 3 || ( params[2] >= -1 - tol && params[2] <= 1 + tol ) ) )
00092         return true;
00093     else
00094         return false;
00095 }
00096 
00097 /** \brief Given type & #vertices, get an appropriate eval set */
00098 ErrorCode EvalSet::get_eval_set( EntityType tp, unsigned int num_vertices, EvalSet& eval_set )
00099 {
00100     switch( tp )
00101     {
00102         case MBEDGE:
00103             break;
00104         case MBTRI:
00105             if( LinearTri::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
00106             break;
00107         case MBQUAD:
00108             if( LinearQuad::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
00109             //            if (SpectralQuad::compatible(tp, num_vertices, eval_set)) return
00110             //            MB_SUCCESS;
00111             break;
00112         case MBTET:
00113             if( LinearTet::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
00114             break;
00115         case MBHEX:
00116             if( LinearHex::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
00117             if( QuadraticHex::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
00118             //            if (SpectralHex::compatible(tp, num_vertices, eval_set)) return
00119             //            MB_SUCCESS;
00120             break;
00121         default:
00122             break;
00123     }
00124 
00125     return MB_NOT_IMPLEMENTED;
00126 }
00127 
00128 ErrorCode ElemEvaluator::find_containing_entity( Range& entities, const double* point, const double iter_tol,
00129                                                  const double inside_tol, EntityHandle& containing_ent, double* params,
00130                                                  unsigned int* num_evals )
00131 {
00132     int is_inside;
00133     ErrorCode rval      = MB_SUCCESS;
00134     unsigned int nevals = 0;
00135     Range::iterator i;
00136     for( i = entities.begin(); i != entities.end(); ++i )
00137     {
00138         nevals++;
00139         set_ent_handle( *i );
00140         rval = reverse_eval( point, iter_tol, inside_tol, params, &is_inside );
00141         if( MB_SUCCESS != rval ) return rval;
00142         if( is_inside ) break;
00143     }
00144     containing_ent = ( i == entities.end() ? 0 : *i );
00145     if( num_evals ) *num_evals += nevals;
00146     return MB_SUCCESS;
00147 }
00148 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines