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