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