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