Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ElemUtil.hpp
Go to the documentation of this file.
1 #ifndef MOAB_ELEM_UTIL_HPP
2 #define MOAB_ELEM_UTIL_HPP
3 
4 #include "moab/CartVect.hpp"
5 #include <vector>
6 #include "moab/Matrix3.hpp"
7 
8 // to access data structures for spectral elements
9 
10 extern "C" {
11 #include "moab/FindPtFuncs.h"
12 }
13 
14 namespace moab
15 {
16 namespace ElemUtil
17 {
18 
19  bool nat_coords_trilinear_hex( const CartVect* hex_corners, const CartVect& x, CartVect& xi, double tol );
20  bool point_in_trilinear_hex( const CartVect* hex_corners, const CartVect& xyz, double etol );
21 
22  bool point_in_trilinear_hex( const CartVect* hex_corners,
23  const CartVect& xyz,
24  const CartVect& box_min,
25  const CartVect& box_max,
26  double etol );
27 
28  // wrapper to hex_findpt
29  void nat_coords_trilinear_hex2( const CartVect* hex_corners, const CartVect& x, CartVect& xi, double til );
30 
31  void hex_findpt( double* xm[3], int n, CartVect xyz, CartVect& rst, double& dist );
32 
33  void hex_eval( double* field, int n, CartVect rst, double& value );
34 
35  bool integrate_trilinear_hex( const CartVect* hex_corners, double* corner_fields, double& field_val, int num_pts );
36 
37 } // namespace ElemUtil
38 
39 namespace Element
40 {
41  /**\brief Class representing a map (diffeomorphism) F parameterizing a 3D element by its
42  * canonical preimage.*/
43  /*
44  Shape functions on the element can obtained by a pushforward (pullback by the inverse map)
45  of the shape functions on the canonical element. This is done by extending this class.
46 
47  We further assume that the parameterizing map is defined by the location of n vertices,
48  which can be set and retrieved on a Map instance. The number of vertices is fixed at
49  compile time.
50  */
51  class Map
52  {
53  public:
54  /**\brief Construct a Map defined by the given std::vector of vertices. */
55  Map( const std::vector< CartVect >& v )
56  {
57  this->vertex.resize( v.size() );
58  this->set_vertices( v );
59  };
60  /**\brief Construct a Map defined by n vertices. */
61  Map( const unsigned int n )
62  {
63  this->vertex = std::vector< CartVect >( n );
64  };
65  virtual ~Map();
66  /**\brief Evaluate the map on \f$x_i\f$ (calculate \f$\vec x = F($\vec \xi)\f$ )*/
67  virtual CartVect evaluate( const CartVect& xi ) const = 0;
68  /**\brief Evaluate the inverse map (calculate \f$\vec \xi = F^-1($\vec x)\f$ to given
69  * tolerance)*/
70  virtual CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
71  /**\brief decide if within the natural param space, with a tolerance*/
72  virtual bool inside_nat_space( const CartVect& xi, double& tol ) const = 0;
73  /* FIX: should evaluate and ievaluate return both the value and the Jacobian (first jet)? */
74  /**\brief Evaluate the map's Jacobi matrix. */
75  virtual Matrix3 jacobian( const CartVect& xi ) const = 0;
76  /* FIX: should this be evaluated in real coordinates and be obtained as part of a Newton
77  * solve? */
78  /**\brief Evaluate the inverse of the Jacobi matrix. */
79  virtual Matrix3 ijacobian( const CartVect& xi ) const
80  {
81  return this->jacobian( xi ).inverse();
82  };
83  /* det_jacobian and det_ijacobian should be overridden for efficiency. */
84  /**\brief Evaluate the determinate of the Jacobi matrix. */
85  virtual double det_jacobian( const CartVect& xi ) const
86  {
87  return this->jacobian( xi ).determinant();
88  };
89  /* FIX: should this be evaluated in real coordinates and be obtained as part of a Newton
90  * solve? */
91  /**\brief Evaluate the determinate of the inverse Jacobi matrix. */
92  virtual double det_ijacobian( const CartVect& xi ) const
93  {
94  return this->jacobian( xi ).inverse().determinant();
95  };
96 
97  /**\brief Evaluate a scalar field at a point given field values at the vertices. */
98  virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const = 0;
99  /**\brief Integrate a scalar field over the element given field values at the vertices. */
100  virtual double integrate_scalar_field( const double* field_vertex_values ) const = 0;
101 
102  /**\brief Size of the vertices vector. */
103  unsigned int size()
104  {
105  return this->vertex.size();
106  }
107  /**\brief Retrieve vertices. */
108  const std::vector< CartVect >& get_vertices();
109  /**\brief Set vertices. */
110  virtual void set_vertices( const std::vector< CartVect >& v );
111 
112  // will look at the box formed by vertex coordinates, and before doing any NR, bail out if
113  // necessary
114  virtual bool inside_box( const CartVect& xi, double& tol ) const;
115 
116  /* Exception thrown when an evaluation fails (e.g., ievaluate fails to converge). */
118  {
119  public:
120  EvaluationError( const CartVect& x, const std::vector< CartVect >& verts ) : p( x ), vertices( verts )
121  {
122 #ifndef NDEBUG
123  std::cout << "p:" << p << "\n vertices.size() " << vertices.size() << "\n";
124  for( size_t i = 0; i < vertices.size(); i++ )
125  std::cout << vertices[i] << "\n";
126 #endif
127  };
128 
129  private:
131  std::vector< CartVect > vertices;
132  }; // class EvaluationError
133 
134  /* Exception thrown when a bad argument is encountered. */
135  class ArgError
136  {
137  public:
138  ArgError(){};
139  }; // class ArgError
140  protected:
141  std::vector< CartVect > vertex;
142  }; // class Map
143 
144  /**\brief Shape function space for trilinear hexahedron, obtained by a pushforward of the
145  * canonical linear (affine) functions. */
146  class LinearHex : public Map
147  {
148  public:
149  LinearHex( const std::vector< CartVect >& vertices ) : Map( vertices ){};
150  LinearHex();
151  virtual ~LinearHex();
152 
153  virtual CartVect evaluate( const CartVect& xi ) const;
154  // virtual CartVect ievaluate(const CartVect& x, double tol) const ;
155  virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
156 
157  virtual Matrix3 jacobian( const CartVect& xi ) const;
158  virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
159  virtual double integrate_scalar_field( const double* field_vertex_values ) const;
160 
161  protected:
162  /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */
163  static const double corner[8][3];
164  static const double gauss[2][2];
165  static const unsigned int corner_count = 8;
166  static const unsigned int gauss_count = 2;
167 
168  }; // class LinearHex
169 
170  /**\brief Shape function space for trilinear hexahedron, obtained by a pushforward of the
171  * canonical linear (affine) functions. */
172  class QuadraticHex : public Map
173  {
174  public:
175  QuadraticHex( const std::vector< CartVect >& vertices ) : Map( vertices ){};
176  QuadraticHex();
177  virtual ~QuadraticHex();
178  virtual CartVect evaluate( const CartVect& xi ) const;
179  // virtual CartVect ievaluate(const CartVect& x, double tol) const ;
180  virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
181 
182  virtual Matrix3 jacobian( const CartVect& xi ) const;
183  virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
184  virtual double integrate_scalar_field( const double* field_vertex_values ) const;
185 
186  protected:
187  /* Preimages of the vertices -- "canonical vertices" -- are known as "corners".
188  * there are 27 vertices for a tri-quadratic xes*/
189  static const int corner[27][3];
190  static const double gauss[8][2]; // TODO fix me
191  static const unsigned int corner_count = 27;
192  static const unsigned int gauss_count = 2; // TODO fix me
193 
194  }; // class QuadraticHex
195  /**\brief Shape function space for a linear tetrahedron, obtained by a pushforward of the
196  * canonical affine shape functions. */
197  class LinearTet : public Map
198  {
199  public:
200  LinearTet( const std::vector< CartVect >& vertices ) : Map( vertices )
201  {
202  set_vertices( vertex );
203  };
204  LinearTet();
205  virtual ~LinearTet();
206  /* Override the evaluation routines to take advantage of the properties of P1. */
207  virtual CartVect evaluate( const CartVect& xi ) const
208  {
209  return this->vertex[0] + this->T * xi;
210  };
211  virtual CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
212  virtual Matrix3 jacobian( const CartVect& ) const
213  {
214  return this->T;
215  };
216  virtual Matrix3 ijacobian( const CartVect& ) const
217  {
218  return this->T_inverse;
219  };
220  virtual double det_jacobian( const CartVect& ) const
221  {
222  return this->det_T;
223  };
224  virtual double det_ijacobian( const CartVect& ) const
225  {
226  return this->det_T_inverse;
227  };
228  //
229  virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
230  virtual double integrate_scalar_field( const double* field_vertex_values ) const;
231  //
232  /* Override set_vertices so we can precompute the matrices effecting the mapping to and from
233  * the canonical simplex. */
234  virtual void set_vertices( const std::vector< CartVect >& v );
235  virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
236 
237  protected:
238  static const double corner[4][3];
241  }; // class LinearTet
242 
243  class SpectralHex : public Map
244  {
245  public:
246  SpectralHex( const std::vector< CartVect >& vertices ) : Map( vertices )
247  {
248  _xyz[0] = _xyz[1] = _xyz[2] = NULL;
249  };
250  SpectralHex( int order, double* x, double* y, double* z );
251  SpectralHex( int order );
252  SpectralHex();
253  virtual ~SpectralHex();
254  void set_gl_points( double* x, double* y, double* z );
255  virtual CartVect evaluate( const CartVect& xi ) const;
256  virtual CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
257  virtual Matrix3 jacobian( const CartVect& xi ) const;
258  double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
259  double integrate_scalar_field( const double* field_vertex_values ) const;
260  bool inside_nat_space( const CartVect& xi, double& tol ) const;
261 
262  // to compute the values that need to be cached for each element of order n
263  void Init( int order );
264  void freedata();
265 
266  protected:
267  /* values that depend only on the order of the element , cached */
268  /* the order in 3 directions */
269  static int _n;
270  static realType* _z[3];
271  static lagrange_data _ld[3];
273  static realType* _odwork; // work area
274 
275  // flag for initialization of data
276  static bool _init;
277 
279 
280  }; // class SpectralHex
281 
282  /**\brief Shape function space for bilinear quadrilateral, obtained from the canonical linear
283  * (affine) functions. */
284  class LinearQuad : public Map
285  {
286  public:
287  LinearQuad( const std::vector< CartVect >& vertices ) : Map( vertices ){};
288  LinearQuad();
289  virtual ~LinearQuad();
290  virtual CartVect evaluate( const CartVect& xi ) const;
291  // virtual CartVect ievaluate(const CartVect& x, double tol) const ;
292  virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
293 
294  virtual Matrix3 jacobian( const CartVect& xi ) const;
295  virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
296  virtual double integrate_scalar_field( const double* field_vertex_values ) const;
297 
298  protected:
299  /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */
300  static const double corner[4][3];
301  static const double gauss[1][2];
302  static const unsigned int corner_count = 4;
303  static const unsigned int gauss_count = 1;
304 
305  }; // class LinearQuad
306 
307  /**\brief Shape function space for bilinear quadrilateral on sphere, obtained from the
308  * canonical linear (affine) functions.
309  * It is mapped using gnomonic projection to a plane tangent at the first vertex
310  * It works well for edges that are great circle arcs; RLL meshes do not have this property,
311  * but HOMME or MPAS meshes do have it */
312  class SphericalQuad : public LinearQuad
313  {
314  public:
315  SphericalQuad( const std::vector< CartVect >& vertices );
316  virtual ~SphericalQuad(){};
317  virtual bool inside_box( const CartVect& pos, double& tol ) const;
318  CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
319 
320  protected:
322  Matrix3 transf; // so will have a lot of stuff, including the transf to a coordinate system
323  // double tangent_plane; // at first vertex; normal to the plane is first vertex
324 
325  }; // class SphericalQuad
326 
327  /**\brief Shape function space for linear triangle, similar to linear tet. */
328  class LinearTri : public Map
329  {
330  public:
331  LinearTri( const std::vector< CartVect >& vertices ) : Map( vertices )
332  {
333  set_vertices( vertex );
334  };
335  LinearTri();
336  virtual ~LinearTri();
337  /* Override the evaluation routines to take advantage of the properties of P1. */
338  /* similar to tets */
339  virtual CartVect evaluate( const CartVect& xi ) const
340  {
341  return this->vertex[0] + this->T * xi;
342  };
343  virtual CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
344  virtual Matrix3 jacobian( const CartVect& ) const
345  {
346  return this->T;
347  };
348  virtual Matrix3 ijacobian( const CartVect& ) const
349  {
350  return this->T_inverse;
351  };
352  virtual double det_jacobian( const CartVect& ) const
353  {
354  return this->det_T;
355  };
356  virtual double det_ijacobian( const CartVect& ) const
357  {
358  return this->det_T_inverse;
359  };
360 
361  /* Override set_vertices so we can precompute the matrices effecting the mapping to and from
362  * the canonical simplex. */
363  virtual void set_vertices( const std::vector< CartVect >& v );
364  virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
365 
366  virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
367  virtual double integrate_scalar_field( const double* field_vertex_values ) const;
368 
369  protected:
370  /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */
371  static const double corner[3][3];
374 
375  }; // class LinearTri
376 
377  /**\brief Shape function space for linear triangle on sphere, obtained from the
378  * canonical linear (affine) functions.
379  * It is mapped using gnomonic projection to a plane tangent at the first vertex
380  * It works well for edges that are great circle arcs; RLL meshes do not have this property,
381  * but HOMME or MPAS meshes do have it */
382  class SphericalTri : public LinearTri
383  {
384  public:
385  SphericalTri( const std::vector< CartVect >& vertices );
386  virtual ~SphericalTri(){};
387  virtual bool inside_box( const CartVect& pos, double& tol ) const;
388  CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
389 
390  protected:
392  Matrix3 transf; // so will have a lot of stuff, including the transf to a coordinate system
393  // double tangent_plane; // at first vertex; normal to the plane is first vertex
394 
395  }; // class SphericalTri
396 
397  /**\brief Shape function space for bilinear quadrilateral, obtained from the canonical linear
398  * (affine) functions. */
399  class LinearEdge : public Map
400  {
401  public:
402  LinearEdge( const std::vector< CartVect >& vertices ) : Map( vertices ){};
403  LinearEdge();
404  virtual CartVect evaluate( const CartVect& xi ) const;
405  // virtual CartVect ievaluate(const CartVect& x, double tol) const ;
406  virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
407 
408  virtual Matrix3 jacobian( const CartVect& xi ) const;
409  virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
410  virtual double integrate_scalar_field( const double* field_vertex_values ) const;
411 
412  protected:
413  /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */
414  static const double corner[2][3];
415  static const double gauss[1][2];
416  static const unsigned int corner_count = 2;
417  static const unsigned int gauss_count = 1;
418 
419  }; // class LinearEdge
420 
421  class SpectralQuad : public Map
422  {
423  public:
424  SpectralQuad( const std::vector< CartVect >& vertices ) : Map( vertices )
425  {
426  _xyz[0] = _xyz[1] = _xyz[2] = NULL;
427  };
428  SpectralQuad( int order, double* x, double* y, double* z );
429  SpectralQuad( int order );
430  SpectralQuad();
431  virtual ~SpectralQuad();
432  void set_gl_points( double* x, double* y, double* z );
433  virtual CartVect evaluate( const CartVect& xi ) const; // a 2d, so 3rd component is 0, always
434  virtual CartVect ievaluate(
435  const CartVect& x,
436  double tol = 1e-6,
437  const CartVect& x0 = CartVect( 0.0 ) ) const; // a 2d, so 3rd component is 0, always
438  virtual Matrix3 jacobian( const CartVect& xi ) const;
439  double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
440  double integrate_scalar_field( const double* field_vertex_values ) const;
441  bool inside_nat_space( const CartVect& xi, double& tol ) const;
442 
443  // to compute the values that need to be cached for each element of order n
444  void Init( int order );
445  void freedata();
446  // this will take node, vertex positions and compute the gl points
447  void compute_gl_positions();
448  void get_gl_points( double*& x, double*& y, double*& z, int& size );
449 
450  protected:
451  /* values that depend only on the order of the element , cached */
452  /* the order in all 3 directions ; it is also np in HOMME lingo*/
453  static int _n;
454  static realType* _z[2];
455  static lagrange_data _ld[2];
456  static opt_data_2 _data; // we should use only 2nd component
457  static realType* _odwork; // work area
458 
459  // flag for initialization of data
460  static bool _init;
461  static realType* _glpoints; // it is a space we can use to store gl positions for elements
462  // on the fly; we do not have a tag yet for them, as in Nek5000 application
463  // also, these positions might need to be moved on the sphere, for HOMME grids
464  // do we project them or how do we move them on the sphere?
465 
466  realType* _xyz[3]; // these are gl points; position?
467 
468  }; // class SpectralQuad
469 
470 } // namespace Element
471 
472 } // namespace moab
473 
474 #endif /*MOAB_ELEM_UTIL_HPP*/