 |
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
|
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). */
117 class EvaluationError
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:
130 CartVect p;
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];
239 Matrix3 T, T_inverse;
240 double det_T, det_T_inverse;
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];
272 static opt_data_3 _data;
273 static realType* _odwork; // work area
274
275 // flag for initialization of data
276 static bool _init;
277
278 realType* _xyz[3];
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:
321 CartVect v1;
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];
372 Matrix3 T, T_inverse;
373 double det_T, det_T_inverse;
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:
391 CartVect v1;
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*/
1.9.1.