 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
#include <LinearTet.hpp>
Static Public Member Functions | |
| static ErrorCode | evalFcn (const double *params, const double *field, const int ndim, const int num_tuples, double *work, double *result) |
| Forward-evaluation of field at parametric coordinates. More... | |
| static ErrorCode | reverseEvalFcn (EvalFcn eval, JacobianFcn jacob, InsideFcn ins, 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 *is_inside) |
| Reverse-evaluation of parametric coordinates at physical space position. More... | |
| static ErrorCode | normalFcn (const int ientDim, const int facet, const int nverts, const double *verts, double normal[3]) |
| Evaluate the normal at a specified facet. More... | |
| static ErrorCode | jacobianFcn (const double *params, const double *verts, const int nverts, const int ndim, double *work, double *result) |
| Evaluate the jacobian at a specified parametric position. More... | |
| static ErrorCode | integrateFcn (const double *field, const double *verts, const int nverts, const int ndim, const int num_tuples, double *work, double *result) |
| Forward-evaluation of field at parametric coordinates. More... | |
| static ErrorCode | initFcn (const double *verts, const int nverts, double *&work) |
| Initialize this EvalSet. More... | |
| static int | insideFcn (const double *params, const int ndim, const double tol) |
| Function that returns whether or not the parameters are inside the natural space of the element. More... | |
| static ErrorCode | 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) |
| static EvalSet | eval_set () |
| static bool | compatible (EntityType tp, int numv, EvalSet &eset) |
Static Protected Attributes | |
| static const double | corner [4][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } } |
Definition at line 12 of file LinearTet.hpp.
|
inlinestatic |
Definition at line 86 of file LinearTet.hpp.
87 {
88 if( tp == MBTET && numv >= 4 )
89 {
90 eset = eval_set();
91 return true;
92 }
93 else
94 return false;
95 }
References eval_set(), and MBTET.
Referenced by moab::EvalSet::get_eval_set().
|
inlinestatic |
Definition at line 81 of file LinearTet.hpp.
82 {
83 return EvalSet( evalFcn, reverseEvalFcn, normalFcn, jacobianFcn, integrateFcn, initFcn, insideFcn );
84 }
References evalFcn(), initFcn(), insideFcn(), integrateFcn(), jacobianFcn(), normalFcn(), and reverseEvalFcn().
Referenced by compatible().
|
static |
Forward-evaluation of field at parametric coordinates.
Definition at line 44 of file LinearTet.cpp.
50 {
51 assert( params && field && num_tuples > 0 );
52 std::vector< double > f0( num_tuples );
53 std::copy( field, field + num_tuples, f0.begin() );
54 std::copy( field, field + num_tuples, result );
55
56 for( unsigned i = 1; i < 4; ++i )
57 {
58 double p = 0.5 * ( params[i - 1] + 1 ); // transform from -1 <= p <= 1 to 0 <= p <= 1
59 for( int j = 0; j < num_tuples; j++ )
60 result[j] += ( field[i * num_tuples + j] - f0[j] ) * p;
61 }
62
63 return MB_SUCCESS;
64 }
References MB_SUCCESS.
Referenced by eval_set().
|
static |
Definition at line 120 of file LinearTet.cpp.
132 {
133 // TODO: should differentiate between epsilons used for
134 // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
135 // right now, fix the tolerance used for NR
136 const double error_tol_sqr = iter_tol * iter_tol;
137 CartVect* cvparams = reinterpret_cast< CartVect* >( params );
138 const CartVect* cvposn = reinterpret_cast< const CartVect* >( posn );
139
140 // find best initial guess to improve convergence
141 CartVect tmp_params[] = { CartVect( -1, -1, -1 ), CartVect( 1, -1, -1 ), CartVect( -1, 1, -1 ),
142 CartVect( -1, -1, 1 ) };
143 double resl = std::numeric_limits< double >::max();
144 CartVect new_pos, tmp_pos;
145 ErrorCode rval;
146 for( unsigned int i = 0; i < 4; i++ )
147 {
148 rval = ( *eval )( tmp_params[i].array(), verts, ndim, ndim, work, tmp_pos.array() );
149 if( MB_SUCCESS != rval ) return rval;
150 double tmp_resl = ( tmp_pos - *cvposn ).length_squared();
151 if( tmp_resl < resl )
152 {
153 *cvparams = tmp_params[i];
154 new_pos = tmp_pos;
155 resl = tmp_resl;
156 }
157 }
158
159 // residual is diff between old and new pos; need to minimize that
160 CartVect res = new_pos - *cvposn;
161 Matrix3 J;
162 rval = ( *jacob )( cvparams->array(), verts, nverts, ndim, work, J.array() );
163 #ifndef NDEBUG
164 double det = J.determinant();
165 assert( det > std::numeric_limits< double >::epsilon() );
166 #endif
167 Matrix3 Ji = J.inverse();
168
169 int iters = 0;
170 // while |res| larger than tol
171 int dum, *tmp_inside = ( inside ? inside : &dum );
172 while( res % res > error_tol_sqr )
173 {
174 if( ++iters > 25 )
175 {
176 // if we haven't converged but we're outside, that's defined as success
177 *tmp_inside = ( *inside_f )( params, ndim, inside_tol );
178 if( !( *tmp_inside ) )
179 return MB_SUCCESS;
180 else
181 return MB_INDEX_OUT_OF_RANGE;
182 }
183
184 // new params tries to eliminate residual
185 *cvparams -= Ji * res;
186
187 // get the new forward-evaluated position, and its difference from the target pt
188 rval = ( *eval )( params, verts, ndim, ndim, work, new_pos.array() );
189 if( MB_SUCCESS != rval ) return rval;
190 res = new_pos - *cvposn;
191 }
192
193 if( inside ) *inside = ( *inside_f )( params, ndim, inside_tol );
194
195 return MB_SUCCESS;
196 } // Map::evaluate_reverse()
References moab::CartVect::array(), moab::Matrix3::array(), moab::Matrix3::determinant(), moab::dum, ErrorCode, moab::Matrix3::inverse(), length_squared(), MB_INDEX_OUT_OF_RANGE, and MB_SUCCESS.
Referenced by reverseEvalFcn().
|
static |
Initialize this EvalSet.
Definition at line 12 of file LinearTet.cpp.
13 {
14 // allocate work array as:
15 // work[0..8] = T
16 // work[9..17] = Tinv
17 // work[18] = detT
18 // work[19] = detTinv
19 if( nverts != 4 )
20 {
21 std::cout << "Invalid Tetrahedron. Expected 4 vertices.\n";
22 return MB_FAILURE;
23 }
24
25 assert( verts );
26
27 Matrix3 J( verts[1 * 3 + 0] - verts[0 * 3 + 0], verts[2 * 3 + 0] - verts[0 * 3 + 0],
28 verts[3 * 3 + 0] - verts[0 * 3 + 0], verts[1 * 3 + 1] - verts[0 * 3 + 1],
29 verts[2 * 3 + 1] - verts[0 * 3 + 1], verts[3 * 3 + 1] - verts[0 * 3 + 1],
30 verts[1 * 3 + 2] - verts[0 * 3 + 2], verts[2 * 3 + 2] - verts[0 * 3 + 2],
31 verts[3 * 3 + 2] - verts[0 * 3 + 2] );
32
33 // Update the work array
34 if( !work ) work = new double[20];
35
36 J.copyto( work );
37 J.inverse().copyto( work + Matrix3::size );
38 work[18] = J.determinant();
39 work[19] = ( work[18] < 1e-12 ? std::numeric_limits< double >::max() : 1.0 / work[18] );
40
41 return MB_SUCCESS;
42 }
References moab::Matrix3::copyto(), moab::Matrix3::determinant(), moab::Matrix3::inverse(), MB_SUCCESS, and moab::Matrix3::size.
Referenced by eval_set().
|
static |
Function that returns whether or not the parameters are inside the natural space of the element.
Definition at line 114 of file LinearTet.cpp.
115 {
116 return ( params[0] >= -1.0 - tol && params[1] >= -1.0 - tol && params[2] >= -1.0 - tol &&
117 params[0] + params[1] + params[2] <= 1.0 + tol );
118 }
Referenced by eval_set().
|
static |
Forward-evaluation of field at parametric coordinates.
Definition at line 66 of file LinearTet.cpp.
73 {
74 assert( field && num_tuples > 0 );
75 std::fill( result, result + num_tuples, 0.0 );
76 for( int i = 0; i < nverts; ++i )
77 {
78 for( int j = 0; j < num_tuples; j++ )
79 result[j] += field[i * num_tuples + j];
80 }
81 double tmp = work[18] / 24.0;
82 for( int i = 0; i < num_tuples; i++ )
83 result[i] *= tmp;
84
85 return MB_SUCCESS;
86 }
References MB_SUCCESS.
Referenced by eval_set().
|
static |
Evaluate the jacobian at a specified parametric position.
Definition at line 88 of file LinearTet.cpp.
89 {
90 // jacobian is cached in work array
91 assert( work );
92 std::copy( work, work + 9, result );
93 return MB_SUCCESS;
94 }
References MB_SUCCESS.
Referenced by eval_set().
|
static |
Evaluate the normal at a specified facet.
Definition at line 198 of file LinearTet.cpp.
203 {
204 // assert(facet < 4 && ientDim == 2 && nverts == 4);
205 if( nverts != 4 ) MB_SET_ERR( MB_FAILURE, "Incorrect vertex count for passed tet :: expected value = 4 " );
206 if( ientDim != 2 ) MB_SET_ERR( MB_FAILURE, "Requesting normal for unsupported dimension :: expected value = 2 " );
207 if( facet > 4 || facet < 0 ) MB_SET_ERR( MB_FAILURE, "Incorrect local face id :: expected value = one of 0-3" );
208
209 int id0 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][0];
210 int id1 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][1];
211 int id2 = CN::mConnectivityMap[MBTET][ientDim - 1].conn[facet][2];
212
213 double x0[3], x1[3];
214
215 for( int i = 0; i < 3; i++ )
216 {
217 x0[i] = verts[3 * id1 + i] - verts[3 * id0 + i];
218 x1[i] = verts[3 * id2 + i] - verts[3 * id0 + i];
219 }
220
221 double a = x0[1] * x1[2] - x1[1] * x0[2];
222 double b = x1[0] * x0[2] - x0[0] * x1[2];
223 double c = x0[0] * x1[1] - x1[0] * x0[1];
224 double nrm = sqrt( a * a + b * b + c * c );
225
226 if( nrm > std::numeric_limits< double >::epsilon() )
227 {
228 normal[0] = a / nrm;
229 normal[1] = b / nrm;
230 normal[2] = c / nrm;
231 }
232 return MB_SUCCESS;
233 }
References moab::CN::ConnMap::conn, MB_SET_ERR, MB_SUCCESS, MBTET, and moab::CN::mConnectivityMap.
Referenced by eval_set().
|
static |
Reverse-evaluation of parametric coordinates at physical space position.
Definition at line 96 of file LinearTet.cpp.
108 {
109 assert( posn && verts );
110 return evaluate_reverse( eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work, params,
111 is_inside );
112 }
References evaluate_reverse().
Referenced by eval_set().
|
staticprotected |
Definition at line 98 of file LinearTet.hpp.
1.9.1.