Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ElemEvaluator.cpp
Go to the documentation of this file.
1 #include <limits>
2 
3 #include "moab/ElemEvaluator.hpp"
4 #include "moab/CartVect.hpp"
5 #include "moab/Matrix3.hpp"
6 
7 // need to include eval set types here to support get_eval_set; alternative would be to have some
8 // type of registration, but we'd still need static registration for the built-in types
14 //#include "moab/SpectralQuad.hpp"
15 //#include "moab/SpectralHex.hpp"
16 
17 namespace moab
18 {
20  JacobianFcn jacob,
21  InsideFcn inside_f,
22  const double* posn,
23  const double* verts,
24  const int nverts,
25  const int ndim,
26  const double iter_tol,
27  const double inside_tol,
28  double* work,
29  double* params,
30  int* inside )
31 {
32  // TODO: should differentiate between epsilons used for
33  // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
34  // right now, fix the tolerance used for NR
35  const double error_tol_sqr = iter_tol * iter_tol;
36  CartVect* cvparams = reinterpret_cast< CartVect* >( params );
37  const CartVect* cvposn = reinterpret_cast< const CartVect* >( posn );
38 
39  // initialize to center of element
40  *cvparams = CartVect( -.4 );
41 
42  CartVect new_pos;
43  // evaluate that first guess to get a new position
44  ErrorCode rval = ( *eval )( cvparams->array(), verts, ndim,
45  3, // hardwire to num_tuples to 3 since the field is coords
46  work, new_pos.array() );
47  if( MB_SUCCESS != rval ) return rval;
48 
49  // residual is diff between old and new pos; need to minimize that
50  CartVect res = new_pos - *cvposn;
51  Matrix3 J;
52  int dum, *tmp_inside = ( inside ? inside : &dum );
53 
54  int iters = 0;
55  // while |res| larger than tol
56  while( res % res > error_tol_sqr )
57  {
58  if( ++iters > 10 )
59  {
60  // if we haven't converged but we're outside, that's defined as success
61  *tmp_inside = ( *inside_f )( params, ndim, inside_tol );
62  if( !( *tmp_inside ) )
63  return MB_SUCCESS;
64  else
65  return MB_FAILURE;
66  }
67 
68  // get jacobian at current params
69  rval = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J.array() );
70  double det = J.determinant();
71  if( det < std::numeric_limits< double >::epsilon() )
72  {
73  *tmp_inside = ( *inside_f )( params, ndim, inside_tol );
74  if( !( *tmp_inside ) )
75  return MB_SUCCESS;
76  else
77  return MB_INDEX_OUT_OF_RANGE;
78  }
79 
80  // new params tries to eliminate residual
81  *cvparams -= J.inverse() * res;
82 
83  // get the new forward-evaluated position, and its difference from the target pt
84  rval = ( *eval )( params, verts, ndim,
85  3, // hardwire to num_tuples to 3 since the field is coords
86  work, new_pos.array() );
87  if( MB_SUCCESS != rval ) return rval;
88  res = new_pos - *cvposn;
89  }
90 
91  if( inside ) *inside = ( *inside_f )( params, ndim, inside_tol );
92 
93  return MB_SUCCESS;
94 } // Map::evaluate_reverse()
95 
96 int EvalSet::inside_function( const double* params, const int ndims, const double tol )
97 {
98  if( params[0] >= -1 - tol && params[0] <= 1 + tol &&
99  ( ndims < 2 || ( params[1] >= -1 - tol && params[1] <= 1 + tol ) ) &&
100  ( ndims < 3 || ( params[2] >= -1 - tol && params[2] <= 1 + tol ) ) )
101  return true;
102  else
103  return false;
104 }
105 
106 /** \brief Given type & #vertices, get an appropriate eval set */
107 ErrorCode EvalSet::get_eval_set( EntityType tp, unsigned int num_vertices, EvalSet& eval_set )
108 {
109  switch( tp )
110  {
111  case MBEDGE:
112  break;
113  case MBTRI:
114  if( LinearTri::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
115  break;
116  case MBQUAD:
117  if( LinearQuad::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
118  // if (SpectralQuad::compatible(tp, num_vertices, eval_set)) return
119  // MB_SUCCESS;
120  break;
121  case MBTET:
122  if( LinearTet::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
123  break;
124  case MBHEX:
125  if( LinearHex::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
126  if( QuadraticHex::compatible( tp, num_vertices, eval_set ) ) return MB_SUCCESS;
127  // if (SpectralHex::compatible(tp, num_vertices, eval_set)) return
128  // MB_SUCCESS;
129  break;
130  default:
131  break;
132  }
133 
134  return MB_NOT_IMPLEMENTED;
135 }
136 
138  const double* point,
139  const double iter_tol,
140  const double inside_tol,
141  EntityHandle& containing_ent,
142  double* params,
143  unsigned int* num_evals )
144 {
145  int is_inside;
146  ErrorCode rval = MB_SUCCESS;
147  unsigned int nevals = 0;
148  Range::iterator i;
149  for( i = entities.begin(); i != entities.end(); ++i )
150  {
151  nevals++;
152  set_ent_handle( *i );
153  rval = reverse_eval( point, iter_tol, inside_tol, params, &is_inside );
154  if( MB_SUCCESS != rval ) return rval;
155  if( is_inside ) break;
156  }
157  containing_ent = ( i == entities.end() ? 0 : *i );
158  if( num_evals ) *num_evals += nevals;
159  return MB_SUCCESS;
160 }
161 } // namespace moab