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