Branch data Line data Source code
1 : : #ifndef ELEM_EVALUATOR_HPP
2 : : #define ELEM_EVALUATOR_HPP
3 : :
4 : : #include "moab/Interface.hpp"
5 : : #include "moab/CartVect.hpp"
6 : : #include "moab/Matrix3.hpp"
7 : : #include "moab/CN.hpp"
8 : :
9 : : #include <vector>
10 : :
11 : : namespace moab
12 : : {
13 : :
14 : : typedef ErrorCode ( *EvalFcn )( const double* params, const double* field, const int ndim, const int num_tuples,
15 : : double* work, double* result );
16 : :
17 : : typedef ErrorCode ( *JacobianFcn )( const double* params, const double* verts, const int nverts, const int ndim,
18 : : double* work, double* result );
19 : :
20 : : typedef ErrorCode ( *IntegrateFcn )( const double* field, const double* verts, const int nverts, const int ndim,
21 : : const int num_tuples, double* work, double* result );
22 : :
23 : : typedef ErrorCode ( *InitFcn )( const double* verts, const int nverts, double*& work );
24 : :
25 : : typedef int ( *InsideFcn )( const double* verts, const int ndims, const double tol );
26 : :
27 : : typedef ErrorCode ( *ReverseEvalFcn )( EvalFcn eval, JacobianFcn jacob, InsideFcn ins, const double* posn,
28 : : const double* verts, const int nverts, const int ndim, const double iter_tol,
29 : : const double inside_tol, double* work, double* params, int* is_inside );
30 : :
31 : : typedef ErrorCode ( *NormalFcn )( const int ientDim, const int facet, const int nverts, const double* verts,
32 : : double normal[3] );
33 : :
34 : : class EvalSet
35 : : {
36 : : public:
37 : : /** \brief Forward-evaluation of field at parametric coordinates */
38 : : EvalFcn evalFcn;
39 : :
40 : : /** \brief Reverse-evaluation of parametric coordinates at physical space position */
41 : : ReverseEvalFcn reverseEvalFcn;
42 : :
43 : : /** \brief Evaluate the normal at a local facet (edge/face for 2D/3D) */
44 : : NormalFcn normalFcn;
45 : :
46 : : /** \brief Evaluate the jacobian at a specified parametric position */
47 : : JacobianFcn jacobianFcn;
48 : :
49 : : /** \brief Forward-evaluation of field at parametric coordinates */
50 : : IntegrateFcn integrateFcn;
51 : :
52 : : /** \brief Initialization function for an element */
53 : : InitFcn initFcn;
54 : :
55 : : /** \brief Function that returns whether or not the parameters are inside the natural space of
56 : : * the element */
57 : : InsideFcn insideFcn;
58 : :
59 : : /** \brief Bare constructor */
60 : 135 : EvalSet()
61 : : : evalFcn( NULL ), reverseEvalFcn( NULL ), normalFcn( NULL ), jacobianFcn( NULL ), integrateFcn( NULL ),
62 : 135 : initFcn( NULL ), insideFcn( NULL )
63 : : {
64 : 135 : }
65 : :
66 : : /** \brief Constructor */
67 : 18 : EvalSet( EvalFcn eval, ReverseEvalFcn rev, NormalFcn normal, JacobianFcn jacob, IntegrateFcn integ, InitFcn initf,
68 : : InsideFcn insidef )
69 : : : evalFcn( eval ), reverseEvalFcn( rev ), normalFcn( normal ), jacobianFcn( jacob ), integrateFcn( integ ),
70 : 18 : initFcn( initf ), insideFcn( insidef )
71 : : {
72 : 18 : }
73 : :
74 : : /** \brief Given an entity handle, get an appropriate eval set, based on type & #vertices */
75 : : static ErrorCode get_eval_set( Interface* mb, EntityHandle eh, EvalSet& eval_set );
76 : :
77 : : /** \brief Given type & #vertices, get an appropriate eval set */
78 : : static ErrorCode get_eval_set( EntityType tp, unsigned int num_vertices, EvalSet& eval_set );
79 : :
80 : : /** \brief Operator= */
81 : 18 : EvalSet& operator=( const EvalSet& eval )
82 : : {
83 : 18 : evalFcn = eval.evalFcn;
84 : 18 : reverseEvalFcn = eval.reverseEvalFcn;
85 : 18 : normalFcn = eval.normalFcn;
86 : 18 : jacobianFcn = eval.jacobianFcn;
87 : 18 : integrateFcn = eval.integrateFcn;
88 : 18 : initFcn = eval.initFcn;
89 : 18 : insideFcn = eval.insideFcn;
90 : 18 : return *this;
91 : : }
92 : :
93 : : /** \brief Common function to do reverse evaluation based on evaluation and jacobian functions
94 : : */
95 : : static ErrorCode evaluate_reverse( EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f, const double* posn,
96 : : const double* verts, const int nverts, const int ndim, const double iter_tol,
97 : : const double inside_tol, double* work, double* params, int* inside );
98 : : /** \brief Common function that returns true if params is in [-1,1]^ndims */
99 : : static int inside_function( const double* params, const int ndims, const double tol );
100 : : };
101 : :
102 : : /** \brief Given an entity handle, get an appropriate eval set, based on type & #vertices */
103 : 3 : inline ErrorCode EvalSet::get_eval_set( Interface* mb, EntityHandle eh, EvalSet& eval_set )
104 : : {
105 : : int nv;
106 [ + - ]: 3 : EntityType tp = mb->type_from_handle( eh );
107 : : const EntityHandle* connect;
108 [ + - ]: 3 : std::vector< EntityHandle > dum_vec;
109 [ + - ]: 3 : ErrorCode rval = mb->get_connectivity( eh, connect, nv, false, &dum_vec );
110 [ - + ]: 3 : if( MB_SUCCESS != rval ) return rval;
111 : :
112 [ + - ]: 3 : return get_eval_set( tp, nv, eval_set );
113 : : }
114 : :
115 : : /**\brief Class facilitating local discretization-related functions
116 : : * \class ElemEvaluator
117 : : * This class implements discretization-related functionality operating
118 : : * on data in MOAB. A member of this class caches certain data about the element
119 : : * it's currently operating on, but is not meant to be instantiated one-per-element,
120 : : * but rather one-per-search (or other operation on a collection of elements).
121 : : *
122 : : * Actual discretization functionality is accessed through function pointers,
123 : : * allowing applications to specialize the implementation of specific functions
124 : : * while still using this class.
125 : : *
126 : : * This class depends on MOAB functionality for gathering entity-based data; the functions
127 : : * it calls through function pointers depend only on POD (plain old data, or intrinsic data types).
128 : : * This allows the use of other packages for serving these functions, without having to modify
129 : : * them to get data through MOAB. This should also promote efficiency, since in many cases they
130 : : * will be able to read data from its native storage locations.
131 : : */
132 : :
133 : : class ElemEvaluator
134 : : {
135 : : public:
136 : : /** \brief Constructor
137 : : * \param impl MOAB instance
138 : : * \param ent Entity handle to cache on the evaluator; if non-zero, calls set_ent_handle, which
139 : : * does some other stuff. \param tag Tag to cache on the evaluator; if non-zero, calls
140 : : * set_tag_handle, which does some other stuff too. \param tagged_ent_dim Dimension of entities
141 : : * to be tagged to cache on the evaluator
142 : : */
143 : : ElemEvaluator( Interface* impl, EntityHandle ent = 0, Tag tag = 0, int tagged_ent_dim = -1 );
144 : : ~ElemEvaluator();
145 : :
146 : : /** \brief Evaluate cached tag at a given parametric location within the cached entity
147 : : * If evaluating coordinates, call set_tag(0, 0), which indicates coords instead of a tag.
148 : : * \param params Parameters at which to evaluate field
149 : : * \param result Result of evaluation
150 : : * \param num_tuples Size of evaluated field, in number of values
151 : : */
152 : : ErrorCode eval( const double* params, double* result, int num_tuples = -1 ) const;
153 : :
154 : : /** \brief Reverse-evaluate the cached entity at a given physical position
155 : : * \param posn Position at which to evaluate parameters
156 : : * \param iter_tol Tolerance of reverse evaluation non-linear iteration, usually 10^-10 or so
157 : : * \param inside_tol Tolerance of is_inside evaluation, usually 10^-6 or so
158 : : * \param params Result of evaluation
159 : : * \param is_inside If non-NULL, returns true of resulting parameters place the point inside the
160 : : * element (in most cases, within [-1]*(dim)
161 : : */
162 : : ErrorCode reverse_eval( const double* posn, double iter_tol, double inside_tol, double* params,
163 : : int* is_inside = NULL ) const;
164 : :
165 : : /**
166 : : * \brief Evaluate the normal to a facet of an entity
167 : : * \param ientDim Dimension of the facet. Should be (d-1) for d-dimensional entities
168 : : * \param facet Local id of the facet w.r.t the entity
169 : : * \param normal Returns the normal.
170 : : */
171 : : ErrorCode get_normal( const int ientDim, const int facet, double normal[3] ) const;
172 : :
173 : : /** \brief Evaluate the jacobian of the cached entity at a given parametric location
174 : : * \param params Parameters at which to evaluate jacobian
175 : : * \param result Result of evaluation, in the form of a 3x3 matrix, stored in column-major order
176 : : */
177 : : ErrorCode jacobian( const double* params, double* result ) const;
178 : :
179 : : /** \brief Integrate the cached tag over the cached entity
180 : : * \param result Result of the integration
181 : : */
182 : : ErrorCode integrate( double* result ) const;
183 : :
184 : : /** \brief Return whether a physical position is inside the cached entity to within a tolerance
185 : : * \param params Parameters at which to query the element
186 : : * \param tol Tolerance, usually 10^-6 or so
187 : : */
188 : : int inside( const double* params, const double tol ) const;
189 : :
190 : : /** \brief Given a list of entities, return the entity the point is in, or none
191 : : * This function reverse-evaluates the entities, returning the first entity containing the
192 : : * point. If no entity contains the point, containing_ent is returned as 0 and params are
193 : : * unchanged. This function returns something other than MB_SUCCESS only when queries go wrong
194 : : * for some reason. num_evals, if non-NULL, is always incremented for each call to reverse_eval.
195 : : * This function calls set_ent_handle for each entity before calling reverse_eval, so the
196 : : * ElemEvaluator object is changed. \param entities Entities tested \param point Point tested,
197 : : * must have 3 dimensions, even for edge and face entities \param iter_tol Tolerance for
198 : : * non-linear reverse evaluation \param inside_tol Tolerance for is_inside test \param
199 : : * containing_ent Entity containing the point, returned 0 if no entity \param params Parameters
200 : : * of point in containing entity, unchanged if no containing entity \param num_evals If
201 : : * non-NULL, incremented each time reverse_eval is called \return Returns non-success only if
202 : : * evaulation failed for some reason (point not in element is NOT a reason for failure)
203 : : */
204 : : ErrorCode find_containing_entity( Range& entities, const double* point, const double iter_tol,
205 : : const double inside_tol, EntityHandle& containing_ent, double* params,
206 : : unsigned int* num_evals = NULL );
207 : :
208 : : /** \brief Given an entity set, return the contained entity the point is in, or none
209 : : * This function reverse-evaluates the entities, returning the first entity containing the
210 : : * point. If no entity contains the point, containing_ent is returned as 0 and params are
211 : : * unchanged. This function returns something other than MB_SUCCESS only when queries go wrong
212 : : * for some reason. num_evals, if non-NULL, is always incremented for each call to reverse_eval.
213 : : * This function calls set_ent_handle for each entity before calling reverse_eval, so the
214 : : * ElemEvaluator object is changed. \param ent_set Entity set containing the entities to be
215 : : * tested \param point Point tested, must have 3 dimensions, even for edge and face entities
216 : : * \param iter_tol Tolerance for non-linear reverse evaluation
217 : : * \param inside_tol Tolerance for is_inside test
218 : : * \param containing_ent Entity containing the point, returned 0 if no entity
219 : : * \param params Parameters of point in containing entity, unchanged if no containing entity
220 : : * \param num_evals If non-NULL, incremented each time reverse_eval is called
221 : : * \return Returns non-success only if evaulation failed for some reason (point not in element
222 : : * is NOT a reason for failure)
223 : : */
224 : : ErrorCode find_containing_entity( EntityHandle ent_set, const double* point, const double iter_tol,
225 : : const double inside_tol, EntityHandle& containing_ent, double* params,
226 : : unsigned int* num_evals = NULL );
227 : :
228 : : /** \brief Set the eval set for a given type entity
229 : : * \param tp Entity type for which to set the eval set
230 : : * \param eval_set Eval set object to use
231 : : */
232 : : ErrorCode set_eval_set( EntityType tp, const EvalSet& eval_set );
233 : :
234 : : /** \brief Set the eval set using a given entity handle
235 : : * This function queries the entity type and number of vertices on the entity to decide which
236 : : * type of shape function to use. NOTE: this function should only be called after setting a
237 : : * valid MOAB instance on the evaluator \param eh Entity handle whose type and #vertices are
238 : : * queried
239 : : */
240 : : ErrorCode set_eval_set( const EntityHandle eh );
241 : :
242 : : /** \brief Get the eval set for a given type entity */
243 : : inline EvalSet get_eval_set( EntityType tp )
244 : : {
245 : : return evalSets[tp];
246 : : }
247 : :
248 : : /** \brief Set entity handle & cache connectivty & vertex positions */
249 : : inline ErrorCode set_ent_handle( EntityHandle ent );
250 : :
251 : : /** \brief Get entity handle for this ElemEval */
252 : 12 : inline EntityHandle get_ent_handle() const
253 : : {
254 : 12 : return entHandle;
255 : : }
256 : :
257 : : /* \brief Get vertex positions cached on this EE
258 : : */
259 : 6 : inline double* get_vert_pos()
260 : : {
261 : 6 : return vertPos[0].array();
262 : : }
263 : :
264 : : /* \brief Get the vertex handles cached here */
265 : 3 : inline const EntityHandle* get_vert_handles() const
266 : : {
267 : 3 : return vertHandles;
268 : : }
269 : :
270 : : /* \brief Get the number of vertices for the cached entity */
271 : 9 : inline int get_num_verts() const
272 : : {
273 : 9 : return numVerts;
274 : : }
275 : :
276 : : /* \brief Get the tag handle cached on this entity */
277 : : inline Tag get_tag_handle() const
278 : : {
279 : : return tagHandle;
280 : : };
281 : :
282 : : /* \brief Set the tag handle to cache on this evaluator
283 : : * To designate that vertex coordinates are the desired tag, pass in a tag handle of 0
284 : : * and a tag dimension of 0.
285 : : * \param tag Tag handle to cache, or 0 to cache vertex positions
286 : : * \param tagged_ent_dim Dimension of entities tagged with this tag
287 : : */
288 : : inline ErrorCode set_tag_handle( Tag tag, int tagged_ent_dim = -1 );
289 : :
290 : : /* \brief Set the name of the tag to cache on this evaluator
291 : : * To designate that vertex coordinates are the desired tag, pass in "COORDS" as the name
292 : : * and a tag dimension of 0.
293 : : * \param tag_name Tag handle to cache, or 0 to cache vertex positions
294 : : * \param tagged_ent_dim Dimension of entities tagged with this tag
295 : : */
296 : : inline ErrorCode set_tag( const char* tag_name, int tagged_ent_dim = -1 );
297 : :
298 : : /* \brief Get the dimension of the entities on which tag is set */
299 : : inline int get_tagged_ent_dim() const
300 : : {
301 : : return taggedEntDim;
302 : : };
303 : :
304 : : /* \brief Set the dimension of entities having the tag */
305 : : inline ErrorCode set_tagged_ent_dim( int dim );
306 : :
307 : : /* \brief Get work space, sometimes this is useful for evaluating data you don't want to set as
308 : : * tags on entities Can't be const because most of the functions (evaluate, integrate, etc.)
309 : : * take non-const work space *'s
310 : : */
311 : 3 : inline double* get_work_space()
312 : : {
313 : 3 : return workSpace;
314 : : }
315 : :
316 : : /* \brief MOAB interface cached on this evaluator */
317 : 27 : inline Interface* get_moab()
318 : : {
319 : 27 : return mbImpl;
320 : : }
321 : :
322 : : private:
323 : : /** \brief Interface */
324 : : Interface* mbImpl;
325 : :
326 : : /** \brief Entity handle being evaluated */
327 : : EntityHandle entHandle;
328 : :
329 : : /** \brief Entity type */
330 : : EntityType entType;
331 : :
332 : : /** \brief Entity dimension */
333 : : int entDim;
334 : :
335 : : /** \brief Number of vertices cached here */
336 : : int numVerts;
337 : :
338 : : /** \brief Cached copy of vertex handle ptr */
339 : : const EntityHandle* vertHandles;
340 : :
341 : : /** \brief Cached copy of vertex positions */
342 : : CartVect vertPos[CN::MAX_NODES_PER_ELEMENT];
343 : :
344 : : /** \brief Tag being evaluated */
345 : : Tag tagHandle;
346 : :
347 : : /** \brief Whether tag is coordinates or something else */
348 : : bool tagCoords;
349 : :
350 : : /** \brief Number of values in this tag */
351 : : int numTuples;
352 : :
353 : : /** \brief Dimension of entities from which to grab tag */
354 : : int taggedEntDim;
355 : :
356 : : /** \brief Tag space */
357 : : std::vector< unsigned char > tagSpace;
358 : :
359 : : /** \brief Evaluation methods for elements of various topologies */
360 : : EvalSet evalSets[MBMAXTYPE];
361 : :
362 : : /** \brief Work space for element-specific data */
363 : : double* workSpace;
364 : :
365 : : }; // class ElemEvaluator
366 : :
367 : 11 : inline ElemEvaluator::ElemEvaluator( Interface* impl, EntityHandle ent, Tag tag, int tagged_ent_dim )
368 : : : mbImpl( impl ), entHandle( 0 ), entType( MBMAXTYPE ), entDim( -1 ), numVerts( 0 ), vertHandles( NULL ),
369 [ + + ][ + - ]: 440 : tagHandle( 0 ), tagCoords( false ), numTuples( 0 ), taggedEntDim( 0 ), workSpace( NULL )
[ + + ]
370 : : {
371 [ + + ][ + - ]: 11 : if( ent ) set_ent_handle( ent );
372 [ - + ][ # # ]: 11 : if( tag ) set_tag_handle( tag, tagged_ent_dim );
373 : 11 : }
374 : :
375 : 22 : inline ElemEvaluator::~ElemEvaluator()
376 : : {
377 [ + + ][ + - ]: 11 : if( workSpace ) delete[] workSpace;
378 : 11 : }
379 : :
380 : 6947 : inline ErrorCode ElemEvaluator::set_ent_handle( EntityHandle ent )
381 : : {
382 : 6947 : entHandle = ent;
383 [ + + ]: 6947 : if( workSpace )
384 : : {
385 [ + - ]: 15 : delete[] workSpace;
386 : 15 : workSpace = NULL;
387 : : }
388 : :
389 [ + - ]: 6947 : entType = mbImpl->type_from_handle( ent );
390 [ + - ]: 6947 : entDim = mbImpl->dimension_from_handle( ent );
391 : :
392 [ + - ]: 6947 : std::vector< EntityHandle > dum_vec;
393 [ + - ]: 6947 : ErrorCode rval = mbImpl->get_connectivity( ent, vertHandles, numVerts, false, &dum_vec );
394 [ - + ]: 6947 : if( MB_SUCCESS != rval ) return rval;
395 : :
396 [ + + ][ + - ]: 6947 : if( !evalSets[entType].evalFcn ) EvalSet::get_eval_set( entType, numVerts, evalSets[entType] );
397 : :
398 [ + - ][ + - ]: 6947 : rval = mbImpl->get_coords( vertHandles, numVerts, vertPos[0].array() );
399 [ - + ]: 6947 : if( MB_SUCCESS != rval ) return rval;
400 : :
401 [ + + ]: 6947 : if( tagHandle )
402 : : {
403 [ + - ]: 5 : rval = set_tag_handle( tagHandle );
404 [ - + ]: 5 : if( MB_SUCCESS != rval ) return rval;
405 : : }
406 [ + + ][ + - ]: 6947 : if( evalSets[entType].initFcn ) return ( *evalSets[entType].initFcn )( vertPos[0].array(), numVerts, workSpace );
[ + - ]
407 : 6947 : return MB_SUCCESS;
408 : : }
409 : :
410 : 17 : inline ErrorCode ElemEvaluator::set_tag_handle( Tag tag, int tagged_ent_dim )
411 : : {
412 : 17 : ErrorCode rval = MB_SUCCESS;
413 [ + + ][ + - ]: 17 : if( !tag && !tagged_ent_dim )
414 : : {
415 : 8 : tagCoords = true;
416 : 8 : numTuples = 3;
417 : 8 : taggedEntDim = 0;
418 : 8 : tagHandle = 0;
419 : 8 : return rval;
420 : : }
421 [ + + ]: 9 : else if( tagHandle != tag )
422 : : {
423 : 4 : tagHandle = tag;
424 [ + - ]: 4 : rval = mbImpl->tag_get_length( tagHandle, numTuples );
425 [ - + ]: 4 : if( MB_SUCCESS != rval ) return rval;
426 : : int sz;
427 [ + - ]: 4 : rval = mbImpl->tag_get_bytes( tag, sz );
428 [ - + ]: 4 : if( MB_SUCCESS != rval ) return rval;
429 [ + - ]: 4 : tagSpace.resize( CN::MAX_NODES_PER_ELEMENT * sz );
430 : 4 : tagCoords = false;
431 : : }
432 : :
433 [ + + ]: 9 : taggedEntDim = ( -1 == tagged_ent_dim ? 0 : tagged_ent_dim );
434 : :
435 [ + - ]: 9 : if( entHandle )
436 : : {
437 [ + - ]: 9 : if( 0 == taggedEntDim )
438 : : {
439 : 9 : rval = mbImpl->tag_get_data( tagHandle, vertHandles, numVerts, &tagSpace[0] );
440 [ - + ]: 9 : if( MB_SUCCESS != rval ) return rval;
441 : : }
442 [ # # ]: 0 : else if( taggedEntDim == entDim )
443 : : {
444 : 0 : rval = mbImpl->tag_get_data( tagHandle, &entHandle, 1, &tagSpace[0] );
445 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
446 : : }
447 : : }
448 : :
449 : 17 : return rval;
450 : : }
451 : :
452 : : inline ErrorCode ElemEvaluator::set_tag( const char* tag_name, int tagged_ent_dim )
453 : : {
454 : : ErrorCode rval = MB_SUCCESS;
455 : : if( !tag_name || strlen( tag_name ) == 0 ) return MB_FAILURE;
456 : : Tag tag;
457 : : if( !strcmp( tag_name, "COORDS" ) )
458 : : {
459 : : tagCoords = true;
460 : : taggedEntDim = 0;
461 : : numTuples = 3;
462 : : tagHandle = 0;
463 : : // can return here, because vertex coords already cached when entity handle set
464 : : return rval;
465 : : }
466 : : else
467 : : {
468 : : rval = mbImpl->tag_get_handle( tag_name, tag );
469 : : if( MB_SUCCESS != rval ) return rval;
470 : :
471 : : if( tagHandle != tag )
472 : : {
473 : : tagHandle = tag;
474 : : rval = mbImpl->tag_get_length( tagHandle, numTuples );
475 : : if( MB_SUCCESS != rval ) return rval;
476 : : int sz;
477 : : rval = mbImpl->tag_get_bytes( tag, sz );
478 : : if( MB_SUCCESS != rval ) return rval;
479 : : tagSpace.reserve( CN::MAX_NODES_PER_ELEMENT * sz );
480 : : tagCoords = false;
481 : : }
482 : :
483 : : taggedEntDim = ( -1 == tagged_ent_dim ? entDim : tagged_ent_dim );
484 : : }
485 : :
486 : : if( entHandle )
487 : : {
488 : : if( 0 == taggedEntDim )
489 : : {
490 : : rval = mbImpl->tag_get_data( tagHandle, vertHandles, numVerts, &tagSpace[0] );
491 : : if( MB_SUCCESS != rval ) return rval;
492 : : }
493 : : else if( taggedEntDim == entDim )
494 : : {
495 : : rval = mbImpl->tag_get_data( tagHandle, &entHandle, 1, &tagSpace[0] );
496 : : if( MB_SUCCESS != rval ) return rval;
497 : : }
498 : : }
499 : :
500 : : return rval;
501 : : }
502 : :
503 : 9 : inline ErrorCode ElemEvaluator::set_eval_set( EntityType tp, const EvalSet& eval_set )
504 : : {
505 : 9 : evalSets[tp] = eval_set;
506 [ + + ][ + + ]: 9 : if( entHandle && evalSets[entType].initFcn )
507 : : {
508 : 1 : ErrorCode rval = ( *evalSets[entType].initFcn )( vertPos[0].array(), numVerts, workSpace );
509 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
510 : : }
511 : 9 : return MB_SUCCESS;
512 : : }
513 : :
514 : 10714 : inline ErrorCode ElemEvaluator::eval( const double* params, double* result, int num_tuples ) const
515 : : {
516 [ + - ][ - + ]: 10714 : assert( entHandle && MBMAXTYPE != entType );
517 : : return ( *evalSets[entType].evalFcn )(
518 : 10714 : params, ( tagCoords ? (const double*)vertPos[0].array() : (const double*)&tagSpace[0] ), entDim,
519 [ + - ][ + - ]: 10714 : ( -1 == num_tuples ? numTuples : num_tuples ), workSpace, result );
520 : : }
521 : :
522 : 17631 : inline ErrorCode ElemEvaluator::reverse_eval( const double* posn, const double iter_tol, const double inside_tol,
523 : : double* params, int* ins ) const
524 : : {
525 [ + - ][ - + ]: 17631 : assert( entHandle && MBMAXTYPE != entType );
526 : : return ( *evalSets[entType].reverseEvalFcn )( evalSets[entType].evalFcn, evalSets[entType].jacobianFcn,
527 : : evalSets[entType].insideFcn, posn, vertPos[0].array(), numVerts,
528 : 17631 : entDim, iter_tol, inside_tol, workSpace, params, ins );
529 : : }
530 : :
531 : : /** \brief Evaluate the normal of the cached entity at a given facet */
532 : 32 : inline ErrorCode ElemEvaluator::get_normal( const int ientDim, const int facet, double normal[] ) const
533 : : {
534 [ + - ][ - + ]: 32 : assert( entHandle && MBMAXTYPE != entType );
535 : 32 : return ( *evalSets[entType].normalFcn )( ientDim, facet, numVerts, vertPos[0].array(), normal );
536 : : }
537 : :
538 : : /** \brief Evaluate the jacobian of the cached entity at a given parametric location */
539 : 10714 : inline ErrorCode ElemEvaluator::jacobian( const double* params, double* result ) const
540 : : {
541 [ + - ][ - + ]: 10714 : assert( entHandle && MBMAXTYPE != entType );
542 : 10714 : return ( *evalSets[entType].jacobianFcn )( params, vertPos[0].array(), numVerts, entDim, workSpace, result );
543 : : }
544 : :
545 : : /** \brief Integrate the cached tag over the cached entity */
546 : 8 : inline ErrorCode ElemEvaluator::integrate( double* result ) const
547 : : {
548 [ + - ][ + - ]: 8 : assert( entHandle && MBMAXTYPE != entType && ( tagCoords || tagHandle ) );
[ + - ][ - + ]
549 : 8 : ErrorCode rval = MB_SUCCESS;
550 [ + - ]: 8 : if( !tagCoords )
551 : : {
552 [ + - ]: 8 : if( 0 == taggedEntDim )
553 : 8 : rval = mbImpl->tag_get_data( tagHandle, vertHandles, numVerts, (void*)&tagSpace[0] );
554 : : else
555 : 0 : rval = mbImpl->tag_get_data( tagHandle, &entHandle, 1, (void*)&tagSpace[0] );
556 [ - + ]: 8 : if( MB_SUCCESS != rval ) return rval;
557 : : }
558 : 8 : return ( *evalSets[entType].integrateFcn )( ( tagCoords ? vertPos[0].array() : (const double*)&tagSpace[0] ),
559 [ - + ]: 8 : vertPos[0].array(), numVerts, entDim, numTuples, workSpace, result );
560 : : }
561 : :
562 : 2142 : inline ErrorCode ElemEvaluator::find_containing_entity( EntityHandle ent_set, const double* point,
563 : : const double iter_tol, const double inside_tol,
564 : : EntityHandle& containing_ent, double* params,
565 : : unsigned int* num_evals )
566 : : {
567 [ + - ][ - + ]: 2142 : assert( mbImpl->type_from_handle( ent_set ) == MBENTITYSET );
568 [ + - ]: 2142 : Range entities;
569 [ + - ]: 2142 : ErrorCode rval = mbImpl->get_entities_by_handle( ent_set, entities );
570 [ - + ]: 2142 : if( MB_SUCCESS != rval )
571 : 0 : return rval;
572 : : else
573 [ + - ]: 2142 : return find_containing_entity( entities, point, iter_tol, inside_tol, containing_ent, params, num_evals );
574 : : }
575 : :
576 : 11979 : inline int ElemEvaluator::inside( const double* params, const double tol ) const
577 : : {
578 : 11979 : return ( *evalSets[entType].insideFcn )( params, entDim, tol );
579 : : }
580 : :
581 : : inline ErrorCode ElemEvaluator::set_eval_set( const EntityHandle eh )
582 : : {
583 : : EvalSet eset;
584 : : ErrorCode rval = EvalSet::get_eval_set( mbImpl, eh, evalSets[mbImpl->type_from_handle( eh )] );
585 : : return rval;
586 : : }
587 : :
588 : : } // namespace moab
589 : :
590 : : #endif /*MOAB_ELEM_EVALUATOR_HPP*/
|
