moab
|
00001 #ifndef MOAB_ELEM_UTIL_HPP 00002 #define MOAB_ELEM_UTIL_HPP 00003 00004 #include "moab/CartVect.hpp" 00005 #include <vector> 00006 #include "moab/Matrix3.hpp" 00007 00008 // to access data structures for spectral elements 00009 00010 extern "C" 00011 { 00012 #include "moab/FindPtFuncs.h" 00013 } 00014 00015 namespace moab { 00016 namespace ElemUtil { 00017 00018 bool nat_coords_trilinear_hex(const CartVect* hex_corners, 00019 const CartVect& x, 00020 CartVect& xi, 00021 double tol); 00022 bool point_in_trilinear_hex(const CartVect *hex_corners, 00023 const CartVect& xyz, 00024 double etol); 00025 00026 bool point_in_trilinear_hex(const CartVect *hex_corners, 00027 const CartVect& xyz, 00028 const CartVect& box_min, 00029 const CartVect& box_max, 00030 double etol); 00031 00032 //wrapper to hex_findpt 00033 void nat_coords_trilinear_hex2(const CartVect* hex_corners, 00034 const CartVect& x, 00035 CartVect& xi, 00036 double til); 00037 00038 00039 00040 void hex_findpt(double *xm[3], 00041 int n, 00042 CartVect xyz, 00043 CartVect& rst, 00044 double& dist); 00045 00046 void hex_eval(double *field, 00047 int n, 00048 CartVect rst, 00049 double &value); 00050 00051 bool integrate_trilinear_hex(const CartVect* hex_corners, 00052 double *corner_fields, 00053 double& field_val, 00054 int num_pts); 00055 00056 } // namespace ElemUtil 00057 00058 namespace Element { 00060 /* 00061 Shape functions on the element can obtained by a pushforward (pullback by the inverse map) 00062 of the shape functions on the canonical element. This is done by extending this class. 00063 00064 We further assume that the parameterizing map is defined by the location of n vertices, 00065 which can be set and retrived on a Map instance. The number of vertices is fixed at 00066 compile time. 00067 */ 00068 class Map { 00069 public: 00071 Map(const std::vector<CartVect>& v) {this->vertex.resize(v.size()); this->set_vertices(v);}; 00073 Map(const unsigned int n) {this->vertex = std::vector<CartVect>(n);}; 00074 virtual ~Map(); 00076 virtual CartVect evaluate( const CartVect& xi ) const = 0; 00078 virtual CartVect ievaluate( const CartVect& x, double tol, const CartVect& x0 = CartVect(0.0)) const ; 00080 virtual bool inside_nat_space(const CartVect & xi, double & tol) const = 0; 00081 /* FIX: should evaluate and ievaluate return both the value and the Jacobian (first jet)? */ 00083 virtual Matrix3 jacobian( const CartVect& xi ) const = 0; 00084 /* FIX: should this be evaluated in real coordinates and be obtained as part of a Newton solve? */ 00086 virtual Matrix3 ijacobian( const CartVect& xi ) const {return this->jacobian(xi).inverse();}; 00087 /* det_jacobian and det_ijacobian should be overriden for efficiency. */ 00089 virtual double det_jacobian(const CartVect& xi) const {return this->jacobian(xi).determinant();}; 00090 /* FIX: should this be evaluated in real coordinates and be obtained as part of a Newton solve? */ 00092 virtual double det_ijacobian(const CartVect& xi) const {return this->jacobian(xi).inverse().determinant();}; 00093 00095 virtual double evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const = 0; 00097 virtual double integrate_scalar_field(const double *field_vertex_values) const = 0; 00098 00100 unsigned int size() {return this->vertex.size();} 00102 const std::vector<CartVect>& get_vertices(); 00104 virtual void set_vertices(const std::vector<CartVect>& v); 00105 00106 /* Exception thrown when an evaluation fails (e.g., ievaluate fails to converge). */ 00107 class EvaluationError { 00108 public: 00109 EvaluationError(){}; 00110 };// class EvaluationError 00111 00112 /* Exception thrown when a bad argument is encountered. */ 00113 class ArgError { 00114 public: 00115 ArgError(){}; 00116 };// class ArgError 00117 protected: 00118 std::vector<CartVect> vertex; 00119 };// class Map 00120 00122 class LinearHex : public Map { 00123 public: 00124 LinearHex(const std::vector<CartVect>& vertices) : Map(vertices){}; 00125 LinearHex(); 00126 virtual ~LinearHex(); 00127 00128 virtual CartVect evaluate( const CartVect& xi ) const; 00129 //virtual CartVect ievaluate(const CartVect& x, double tol) const ; 00130 virtual bool inside_nat_space(const CartVect & xi, double & tol) const; 00131 00132 virtual Matrix3 jacobian(const CartVect& xi) const; 00133 virtual double evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const; 00134 virtual double integrate_scalar_field(const double *field_vertex_values) const; 00135 00136 protected: 00137 /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */ 00138 static const double corner[8][3]; 00139 static const double gauss[1][2]; 00140 static const unsigned int corner_count = 8; 00141 static const unsigned int gauss_count = 1; 00142 00143 };// class LinearHex 00144 00146 class QuadraticHex : public Map { 00147 public: 00148 QuadraticHex(const std::vector<CartVect>& vertices) : Map(vertices){}; 00149 QuadraticHex(); 00150 virtual ~QuadraticHex(); 00151 virtual CartVect evaluate( const CartVect& xi ) const; 00152 //virtual CartVect ievaluate(const CartVect& x, double tol) const ; 00153 virtual bool inside_nat_space(const CartVect & xi, double & tol) const; 00154 00155 virtual Matrix3 jacobian(const CartVect& xi) const; 00156 virtual double evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const; 00157 virtual double integrate_scalar_field(const double *field_vertex_values) const; 00158 00159 protected: 00160 /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". 00161 * there are 27 vertices for a tri-quadratic xes*/ 00162 static const int corner[27][3]; 00163 static const double gauss[8][2];// TODO fix me 00164 static const unsigned int corner_count = 27; 00165 static const unsigned int gauss_count = 8; // TODO fix me 00166 00167 };// class QuadraticHex 00169 class LinearTet : public Map { 00170 public: 00171 LinearTet(const std::vector<CartVect>& vertices) : Map(vertices){ LinearTet::set_vertices(vertex);}; 00172 LinearTet(); 00173 virtual ~LinearTet(); 00174 /* Override the evaluation routines to take advantage of the properties of P1. */ 00175 virtual CartVect evaluate(const CartVect& xi) const {return this->vertex[0] + this->T*xi;}; 00176 virtual CartVect ievaluate(const CartVect& x) const {return this->T_inverse*(x-this->vertex[0]);}; 00177 virtual Matrix3 jacobian(const CartVect& ) const {return this->T;}; 00178 virtual Matrix3 ijacobian(const CartVect& ) const {return this->T_inverse;}; 00179 virtual double det_jacobian(const CartVect& ) const {return this->det_T;}; 00180 virtual double det_ijacobian(const CartVect& ) const {return this->det_T_inverse;}; 00181 // 00182 virtual double evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const; 00183 virtual double integrate_scalar_field(const double *field_vertex_values) const; 00184 // 00185 /* Override set_vertices so we can precompute the matrices effecting the mapping to and from the canonical simplex. */ 00186 virtual void set_vertices(const std::vector<CartVect>& v); 00187 virtual bool inside_nat_space(const CartVect & xi, double & tol) const; 00188 protected: 00189 static const double corner[4][3]; 00190 Matrix3 T, T_inverse; 00191 double det_T, det_T_inverse; 00192 };// class LinearTet 00193 00194 class SpectralHex : public Map { 00195 public: 00196 SpectralHex(const std::vector<CartVect>& vertices) : Map(vertices){}; 00197 SpectralHex(int order, double * x, double * y, double *z) ; 00198 SpectralHex(int order); 00199 SpectralHex(); 00200 virtual ~SpectralHex(); 00201 void set_gl_points( double * x, double * y, double *z) ; 00202 virtual CartVect evaluate( const CartVect& xi ) const; 00203 virtual CartVect ievaluate(const CartVect& x) const; 00204 virtual Matrix3 jacobian(const CartVect& xi) const; 00205 double evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const; 00206 double integrate_scalar_field(const double *field_vertex_values) const; 00207 bool inside_nat_space(const CartVect & xi, double & tol) const; 00208 00209 // to compute the values that need to be cached for each element of order n 00210 void Init(int order); 00211 void freedata(); 00212 protected: 00213 /* values that depend only on the order of the element , cached */ 00214 /* the order in 3 directions */ 00215 static int _n; 00216 static real *_z[3]; 00217 static lagrange_data _ld[3]; 00218 static opt_data_3 _data; 00219 static real * _odwork;// work area 00220 00221 // flag for initialization of data 00222 static bool _init; 00223 00224 real * _xyz[3]; 00225 00226 };// class SpectralHex 00227 00229 class LinearQuad : public Map { 00230 public: 00231 LinearQuad(const std::vector<CartVect>& vertices) : Map(vertices){}; 00232 LinearQuad(); 00233 virtual ~LinearQuad(); 00234 virtual CartVect evaluate( const CartVect& xi ) const; 00235 //virtual CartVect ievaluate(const CartVect& x, double tol) const ; 00236 virtual bool inside_nat_space(const CartVect & xi, double & tol) const; 00237 00238 virtual Matrix3 jacobian(const CartVect& xi) const; 00239 virtual double evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const; 00240 virtual double integrate_scalar_field(const double *field_vertex_values) const; 00241 00242 protected: 00243 /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */ 00244 static const double corner[4][3]; 00245 static const double gauss[1][2]; 00246 static const unsigned int corner_count = 4; 00247 static const unsigned int gauss_count = 1; 00248 00249 };// class LinearQuad 00250 00252 class LinearEdge : public Map { 00253 public: 00254 LinearEdge(const std::vector<CartVect>& vertices) : Map(vertices){}; 00255 LinearEdge(); 00256 virtual CartVect evaluate( const CartVect& xi ) const; 00257 //virtual CartVect ievaluate(const CartVect& x, double tol) const ; 00258 virtual bool inside_nat_space(const CartVect & xi, double & tol) const; 00259 00260 virtual Matrix3 jacobian(const CartVect& xi) const; 00261 virtual double evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const; 00262 virtual double integrate_scalar_field(const double *field_vertex_values) const; 00263 00264 protected: 00265 /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */ 00266 static const double corner[2][3]; 00267 static const double gauss[1][2]; 00268 static const unsigned int corner_count = 2; 00269 static const unsigned int gauss_count = 1; 00270 00271 };// class LinearEdge 00272 00273 class SpectralQuad : public Map { 00274 public: 00275 SpectralQuad(const std::vector<CartVect>& vertices) : Map(vertices){}; 00276 SpectralQuad(int order, double * x, double * y, double *z) ; 00277 SpectralQuad(int order); 00278 SpectralQuad(); 00279 virtual ~SpectralQuad(); 00280 void set_gl_points( double * x, double * y, double *z) ; 00281 virtual CartVect evaluate( const CartVect& xi ) const;// a 2d, so 3rd component is 0, always 00282 virtual CartVect ievaluate(const CartVect& x) const; //a 2d, so 3rd component is 0, always 00283 virtual Matrix3 jacobian(const CartVect& xi) const; 00284 double evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const; 00285 double integrate_scalar_field(const double *field_vertex_values) const; 00286 bool inside_nat_space(const CartVect & xi, double & tol) const; 00287 00288 // to compute the values that need to be cached for each element of order n 00289 void Init(int order); 00290 void freedata(); 00291 // this will take node, vertex positions and compute the gl points 00292 void compute_gl_positions(); 00293 void get_gl_points( double *& x, double *& y, double *& z, int & size) ; 00294 protected: 00295 /* values that depend only on the order of the element , cached */ 00296 /* the order in all 3 directions ; it is also np in HOMME lingo*/ 00297 static int _n; 00298 static real *_z[2]; 00299 static lagrange_data _ld[2]; 00300 static opt_data_2 _data; // we should use only 2nd component 00301 static real * _odwork;// work area 00302 00303 // flag for initialization of data 00304 static bool _init; 00305 static real * _glpoints; // it is a space we can use to store gl positions for elements 00306 // on the fly; we do not have a tag yet for them, as in Nek5000 application 00307 // also, these positions might need to be moved on the sphere, for HOMME grids 00308 // do we project them or how do we move them on the sphere? 00309 00310 real * _xyz[3]; // these are gl points; position? 00311 00312 00313 };// class SpectralQuad 00314 00315 00316 }// namespace Element 00317 00318 } // namespace moab 00319 00320 #endif /*MOAB_ELEM_UTIL_HPP*/