Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
ElemUtil.hpp
Go to the documentation of this file.
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 #include "moab/FindPtFuncs.h"
00012 }
00013 
00014 namespace moab
00015 {
00016 namespace ElemUtil
00017 {
00018 
00019     bool nat_coords_trilinear_hex( const CartVect* hex_corners, const CartVect& x, CartVect& xi, double tol );
00020     bool point_in_trilinear_hex( const CartVect* hex_corners, const CartVect& xyz, double etol );
00021 
00022     bool point_in_trilinear_hex( const CartVect* hex_corners,
00023                                  const CartVect& xyz,
00024                                  const CartVect& box_min,
00025                                  const CartVect& box_max,
00026                                  double etol );
00027 
00028     // wrapper to hex_findpt
00029     void nat_coords_trilinear_hex2( const CartVect* hex_corners, const CartVect& x, CartVect& xi, double til );
00030 
00031     void hex_findpt( double* xm[3], int n, CartVect xyz, CartVect& rst, double& dist );
00032 
00033     void hex_eval( double* field, int n, CartVect rst, double& value );
00034 
00035     bool integrate_trilinear_hex( const CartVect* hex_corners, double* corner_fields, double& field_val, int num_pts );
00036 
00037 }  // namespace ElemUtil
00038 
00039 namespace Element
00040 {
00041     /**\brief Class representing a map (diffeomorphism) F parameterizing a 3D element by its
00042      * canonical preimage.*/
00043     /*
00044          Shape functions on the element can obtained by a pushforward (pullback by the inverse map)
00045          of the shape functions on the canonical element. This is done by extending this class.
00046 
00047          We further assume that the parameterizing map is defined by the location of n vertices,
00048          which can be set and retrieved on a Map instance.  The number of vertices is fixed at
00049          compile time.
00050     */
00051     class Map
00052     {
00053       public:
00054         /**\brief Construct a Map defined by the given std::vector of vertices. */
00055         Map( const std::vector< CartVect >& v )
00056         {
00057             this->vertex.resize( v.size() );
00058             this->set_vertices( v );
00059         };
00060         /**\brief Construct a Map defined by n vertices. */
00061         Map( const unsigned int n )
00062         {
00063             this->vertex = std::vector< CartVect >( n );
00064         };
00065         virtual ~Map();
00066         /**\brief Evaluate the map on \f$x_i\f$ (calculate \f$\vec x = F($\vec \xi)\f$ )*/
00067         virtual CartVect evaluate( const CartVect& xi ) const = 0;
00068         /**\brief Evaluate the inverse map (calculate \f$\vec \xi = F^-1($\vec x)\f$ to given
00069          * tolerance)*/
00070         virtual CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
00071         /**\brief decide if within the natural param space, with a tolerance*/
00072         virtual bool inside_nat_space( const CartVect& xi, double& tol ) const = 0;
00073         /* FIX: should evaluate and ievaluate return both the value and the Jacobian (first jet)? */
00074         /**\brief Evaluate the map's Jacobi matrix. */
00075         virtual Matrix3 jacobian( const CartVect& xi ) const = 0;
00076         /* FIX: should this be evaluated in real coordinates and be obtained as part of a Newton
00077          * solve? */
00078         /**\brief Evaluate the inverse of the Jacobi matrix. */
00079         virtual Matrix3 ijacobian( const CartVect& xi ) const
00080         {
00081             return this->jacobian( xi ).inverse();
00082         };
00083         /* det_jacobian and det_ijacobian should be overridden for efficiency. */
00084         /**\brief Evaluate the determinate of the Jacobi matrix. */
00085         virtual double det_jacobian( const CartVect& xi ) const
00086         {
00087             return this->jacobian( xi ).determinant();
00088         };
00089         /* FIX: should this be evaluated in real coordinates and be obtained as part of a Newton
00090          * solve? */
00091         /**\brief Evaluate the determinate of the inverse Jacobi matrix. */
00092         virtual double det_ijacobian( const CartVect& xi ) const
00093         {
00094             return this->jacobian( xi ).inverse().determinant();
00095         };
00096 
00097         /**\brief Evaluate a scalar field at a point given field values at the vertices. */
00098         virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const = 0;
00099         /**\brief Integrate a scalar field over the element given field values at the vertices. */
00100         virtual double integrate_scalar_field( const double* field_vertex_values ) const = 0;
00101 
00102         /**\brief Size of the vertices vector. */
00103         unsigned int size()
00104         {
00105             return this->vertex.size();
00106         }
00107         /**\brief Retrieve vertices. */
00108         const std::vector< CartVect >& get_vertices();
00109         /**\brief Set vertices.      */
00110         virtual void set_vertices( const std::vector< CartVect >& v );
00111 
00112         // will look at the box formed by vertex coordinates, and before doing any NR, bail out if
00113         // necessary
00114         virtual bool inside_box( const CartVect& xi, double& tol ) const;
00115 
00116         /* Exception thrown when an evaluation fails (e.g., ievaluate fails to converge). */
00117         class EvaluationError
00118         {
00119           public:
00120             EvaluationError( const CartVect& x, const std::vector< CartVect >& verts ) : p( x ), vertices( verts )
00121             {
00122 #ifndef NDEBUG
00123                 std::cout << "p:" << p << "\n vertices.size() " << vertices.size() << "\n";
00124                 for( size_t i = 0; i < vertices.size(); i++ )
00125                     std::cout << vertices[i] << "\n";
00126 #endif
00127             };
00128 
00129           private:
00130             CartVect p;
00131             std::vector< CartVect > vertices;
00132         };  // class EvaluationError
00133 
00134         /* Exception thrown when a bad argument is encountered. */
00135         class ArgError
00136         {
00137           public:
00138             ArgError(){};
00139         };  // class ArgError
00140       protected:
00141         std::vector< CartVect > vertex;
00142     };  // class Map
00143 
00144     /**\brief Shape function space for trilinear hexahedron, obtained by a pushforward of the
00145      * canonical linear (affine) functions. */
00146     class LinearHex : public Map
00147     {
00148       public:
00149         LinearHex( const std::vector< CartVect >& vertices ) : Map( vertices ){};
00150         LinearHex();
00151         virtual ~LinearHex();
00152 
00153         virtual CartVect evaluate( const CartVect& xi ) const;
00154         // virtual CartVect ievaluate(const CartVect& x, double tol) const ;
00155         virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
00156 
00157         virtual Matrix3 jacobian( const CartVect& xi ) const;
00158         virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
00159         virtual double integrate_scalar_field( const double* field_vertex_values ) const;
00160 
00161       protected:
00162         /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */
00163         static const double corner[8][3];
00164         static const double gauss[2][2];
00165         static const unsigned int corner_count = 8;
00166         static const unsigned int gauss_count  = 2;
00167 
00168     };  // class LinearHex
00169 
00170     /**\brief Shape function space for trilinear hexahedron, obtained by a pushforward of the
00171      * canonical linear (affine) functions. */
00172     class QuadraticHex : public Map
00173     {
00174       public:
00175         QuadraticHex( const std::vector< CartVect >& vertices ) : Map( vertices ){};
00176         QuadraticHex();
00177         virtual ~QuadraticHex();
00178         virtual CartVect evaluate( const CartVect& xi ) const;
00179         // virtual CartVect ievaluate(const CartVect& x, double tol) const ;
00180         virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
00181 
00182         virtual Matrix3 jacobian( const CartVect& xi ) const;
00183         virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
00184         virtual double integrate_scalar_field( const double* field_vertex_values ) const;
00185 
00186       protected:
00187         /* Preimages of the vertices -- "canonical vertices" -- are known as "corners".
00188          * there are 27 vertices for a tri-quadratic xes*/
00189         static const int corner[27][3];
00190         static const double gauss[8][2];  // TODO fix me
00191         static const unsigned int corner_count = 27;
00192         static const unsigned int gauss_count  = 2;  // TODO fix me
00193 
00194     };  // class QuadraticHex
00195     /**\brief Shape function space for a linear tetrahedron, obtained by a pushforward of the
00196      * canonical affine shape functions. */
00197     class LinearTet : public Map
00198     {
00199       public:
00200         LinearTet( const std::vector< CartVect >& vertices ) : Map( vertices )
00201         {
00202             set_vertices( vertex );
00203         };
00204         LinearTet();
00205         virtual ~LinearTet();
00206         /* Override the evaluation routines to take advantage of the properties of P1. */
00207         virtual CartVect evaluate( const CartVect& xi ) const
00208         {
00209             return this->vertex[0] + this->T * xi;
00210         };
00211         virtual CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
00212         virtual Matrix3 jacobian( const CartVect& ) const
00213         {
00214             return this->T;
00215         };
00216         virtual Matrix3 ijacobian( const CartVect& ) const
00217         {
00218             return this->T_inverse;
00219         };
00220         virtual double det_jacobian( const CartVect& ) const
00221         {
00222             return this->det_T;
00223         };
00224         virtual double det_ijacobian( const CartVect& ) const
00225         {
00226             return this->det_T_inverse;
00227         };
00228         //
00229         virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
00230         virtual double integrate_scalar_field( const double* field_vertex_values ) const;
00231         //
00232         /* Override set_vertices so we can precompute the matrices effecting the mapping to and from
00233          * the canonical simplex. */
00234         virtual void set_vertices( const std::vector< CartVect >& v );
00235         virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
00236 
00237       protected:
00238         static const double corner[4][3];
00239         Matrix3 T, T_inverse;
00240         double det_T, det_T_inverse;
00241     };  // class LinearTet
00242 
00243     class SpectralHex : public Map
00244     {
00245       public:
00246         SpectralHex( const std::vector< CartVect >& vertices ) : Map( vertices )
00247         {
00248             _xyz[0] = _xyz[1] = _xyz[2] = NULL;
00249         };
00250         SpectralHex( int order, double* x, double* y, double* z );
00251         SpectralHex( int order );
00252         SpectralHex();
00253         virtual ~SpectralHex();
00254         void set_gl_points( double* x, double* y, double* z );
00255         virtual CartVect evaluate( const CartVect& xi ) const;
00256         virtual CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
00257         virtual Matrix3 jacobian( const CartVect& xi ) const;
00258         double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
00259         double integrate_scalar_field( const double* field_vertex_values ) const;
00260         bool inside_nat_space( const CartVect& xi, double& tol ) const;
00261 
00262         // to compute the values that need to be cached for each element of order n
00263         void Init( int order );
00264         void freedata();
00265 
00266       protected:
00267         /* values that depend only on the order of the element , cached */
00268         /*  the order in 3 directions */
00269         static int _n;
00270         static realType* _z[3];
00271         static lagrange_data _ld[3];
00272         static opt_data_3 _data;
00273         static realType* _odwork;  // work area
00274 
00275         // flag for initialization of data
00276         static bool _init;
00277 
00278         realType* _xyz[3];
00279 
00280     };  // class SpectralHex
00281 
00282     /**\brief Shape function space for bilinear quadrilateral, obtained from the canonical linear
00283      * (affine) functions. */
00284     class LinearQuad : public Map
00285     {
00286       public:
00287         LinearQuad( const std::vector< CartVect >& vertices ) : Map( vertices ){};
00288         LinearQuad();
00289         virtual ~LinearQuad();
00290         virtual CartVect evaluate( const CartVect& xi ) const;
00291         // virtual CartVect ievaluate(const CartVect& x, double tol) const ;
00292         virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
00293 
00294         virtual Matrix3 jacobian( const CartVect& xi ) const;
00295         virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
00296         virtual double integrate_scalar_field( const double* field_vertex_values ) const;
00297 
00298       protected:
00299         /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */
00300         static const double corner[4][3];
00301         static const double gauss[1][2];
00302         static const unsigned int corner_count = 4;
00303         static const unsigned int gauss_count  = 1;
00304 
00305     };  // class LinearQuad
00306 
00307     /**\brief Shape function space for bilinear quadrilateral on sphere, obtained from the
00308      *  canonical linear (affine) functions.
00309      *  It is mapped using gnomonic projection to a plane tangent at the first vertex
00310      *  It works well for edges that are great circle arcs; RLL meshes  do not have this property,
00311      * but HOMME or MPAS meshes do have it */
00312     class SphericalQuad : public LinearQuad
00313     {
00314       public:
00315         SphericalQuad( const std::vector< CartVect >& vertices );
00316         virtual ~SphericalQuad(){};
00317         virtual bool inside_box( const CartVect& pos, double& tol ) const;
00318         CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
00319 
00320       protected:
00321         CartVect v1;
00322         Matrix3 transf;  // so will have a lot of stuff, including the transf to a coordinate system
00323         // double tangent_plane; // at first vertex; normal to the plane is first vertex
00324 
00325     };  // class SphericalQuad
00326 
00327     /**\brief Shape function space for linear triangle, similar to linear tet. */
00328     class LinearTri : public Map
00329     {
00330       public:
00331         LinearTri( const std::vector< CartVect >& vertices ) : Map( vertices )
00332         {
00333             set_vertices( vertex );
00334         };
00335         LinearTri();
00336         virtual ~LinearTri();
00337         /* Override the evaluation routines to take advantage of the properties of P1. */
00338         /* similar to tets */
00339         virtual CartVect evaluate( const CartVect& xi ) const
00340         {
00341             return this->vertex[0] + this->T * xi;
00342         };
00343         virtual CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
00344         virtual Matrix3 jacobian( const CartVect& ) const
00345         {
00346             return this->T;
00347         };
00348         virtual Matrix3 ijacobian( const CartVect& ) const
00349         {
00350             return this->T_inverse;
00351         };
00352         virtual double det_jacobian( const CartVect& ) const
00353         {
00354             return this->det_T;
00355         };
00356         virtual double det_ijacobian( const CartVect& ) const
00357         {
00358             return this->det_T_inverse;
00359         };
00360 
00361         /* Override set_vertices so we can precompute the matrices effecting the mapping to and from
00362          * the canonical simplex. */
00363         virtual void set_vertices( const std::vector< CartVect >& v );
00364         virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
00365 
00366         virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
00367         virtual double integrate_scalar_field( const double* field_vertex_values ) const;
00368 
00369       protected:
00370         /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */
00371         static const double corner[3][3];
00372         Matrix3 T, T_inverse;
00373         double det_T, det_T_inverse;
00374 
00375     };  // class LinearTri
00376 
00377     /**\brief Shape function space for linear triangle on sphere, obtained from the
00378      *  canonical linear (affine) functions.
00379      *  It is mapped using gnomonic projection to a plane tangent at the first vertex
00380      *  It works well for edges that are great circle arcs; RLL meshes  do not have this property,
00381      * but HOMME or MPAS meshes do have it */
00382     class SphericalTri : public LinearTri
00383     {
00384       public:
00385         SphericalTri( const std::vector< CartVect >& vertices );
00386         virtual ~SphericalTri(){};
00387         virtual bool inside_box( const CartVect& pos, double& tol ) const;
00388         CartVect ievaluate( const CartVect& x, double tol = 1e-6, const CartVect& x0 = CartVect( 0.0 ) ) const;
00389 
00390       protected:
00391         CartVect v1;
00392         Matrix3 transf;  // so will have a lot of stuff, including the transf to a coordinate system
00393         // double tangent_plane; // at first vertex; normal to the plane is first vertex
00394 
00395     };  // class SphericalTri
00396 
00397     /**\brief Shape function space for bilinear quadrilateral, obtained from the canonical linear
00398      * (affine) functions. */
00399     class LinearEdge : public Map
00400     {
00401       public:
00402         LinearEdge( const std::vector< CartVect >& vertices ) : Map( vertices ){};
00403         LinearEdge();
00404         virtual CartVect evaluate( const CartVect& xi ) const;
00405         // virtual CartVect ievaluate(const CartVect& x, double tol) const ;
00406         virtual bool inside_nat_space( const CartVect& xi, double& tol ) const;
00407 
00408         virtual Matrix3 jacobian( const CartVect& xi ) const;
00409         virtual double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
00410         virtual double integrate_scalar_field( const double* field_vertex_values ) const;
00411 
00412       protected:
00413         /* Preimages of the vertices -- "canonical vertices" -- are known as "corners". */
00414         static const double corner[2][3];
00415         static const double gauss[1][2];
00416         static const unsigned int corner_count = 2;
00417         static const unsigned int gauss_count  = 1;
00418 
00419     };  // class LinearEdge
00420 
00421     class SpectralQuad : public Map
00422     {
00423       public:
00424         SpectralQuad( const std::vector< CartVect >& vertices ) : Map( vertices )
00425         {
00426             _xyz[0] = _xyz[1] = _xyz[2] = NULL;
00427         };
00428         SpectralQuad( int order, double* x, double* y, double* z );
00429         SpectralQuad( int order );
00430         SpectralQuad();
00431         virtual ~SpectralQuad();
00432         void set_gl_points( double* x, double* y, double* z );
00433         virtual CartVect evaluate( const CartVect& xi ) const;  // a 2d, so 3rd component is 0, always
00434         virtual CartVect ievaluate(
00435             const CartVect& x,
00436             double tol         = 1e-6,
00437             const CartVect& x0 = CartVect( 0.0 ) ) const;  // a 2d, so 3rd component is 0, always
00438         virtual Matrix3 jacobian( const CartVect& xi ) const;
00439         double evaluate_scalar_field( const CartVect& xi, const double* field_vertex_values ) const;
00440         double integrate_scalar_field( const double* field_vertex_values ) const;
00441         bool inside_nat_space( const CartVect& xi, double& tol ) const;
00442 
00443         // to compute the values that need to be cached for each element of order n
00444         void Init( int order );
00445         void freedata();
00446         // this will take node, vertex positions and compute the gl points
00447         void compute_gl_positions();
00448         void get_gl_points( double*& x, double*& y, double*& z, int& size );
00449 
00450       protected:
00451         /* values that depend only on the order of the element , cached */
00452         /*  the order in all 3 directions ; it is also np in HOMME lingo*/
00453         static int _n;
00454         static realType* _z[2];
00455         static lagrange_data _ld[2];
00456         static opt_data_2 _data;   // we should use only 2nd component
00457         static realType* _odwork;  // work area
00458 
00459         // flag for initialization of data
00460         static bool _init;
00461         static realType* _glpoints;  // it is a space we can use to store gl positions for elements
00462         // on the fly; we do not have a tag yet for them, as in Nek5000 application
00463         // also, these positions might need to be moved on the sphere, for HOMME grids
00464         // do we project them or how do we move them on the sphere?
00465 
00466         realType* _xyz[3];  // these are gl points; position?
00467 
00468     };  // class SpectralQuad
00469 
00470 }  // namespace Element
00471 
00472 }  // namespace moab
00473 
00474 #endif /*MOAB_ELEM_UTIL_HPP*/
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines