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