// ============================================================================ // == == // == AMDiS - Adaptive multidimensional simulations == // == == // ============================================================================ // == == // == TU Dresden == // == == // == Institut für Wissenschaftliches Rechnen == // == Zellescher Weg 12-14 == // == 01069 Dresden == // == germany == // == == // ============================================================================ // == == // == https://gforge.zih.tu-dresden.de/projects/amdis/ == // == == // ============================================================================ /** \file Quadrature.h */ #ifndef AMDIS_QUADRATURE_H #define AMDIS_QUADRATURE_H #include #include "AMDiS_fwd.h" #include "BasisFunction.h" #include "Flag.h" #include "FixVec.h" namespace AMDiS { /** * \ingroup Assembler * * \brief * For the assemblage of the system matrix and right hand side vector of the * linear system, we have to compute integrals, for example: * \f[ \int_{\Omega} f(x)\varphi_i(x) dx \f] * For general data A, b, c, and f, these integrals can not be calculated * exactly. Quadrature formulas have to be used in order to calculate the * integrals approximately. Numerical integration in finite element methods is * done by looping over all grid elements and using a quadrature formula on * each element. */ class Quadrature { protected: /// Avoids call of default constructor Quadrature(); /** \brief * Constructs a Quadrature with name name_ of degree degree_ for dim dim_. * The Quadrature has n_points_ quadrature points with barycentric * coordinates lambda_ and weights w_. The constructor is protected because * access to a Quadrature should be done via \ref provideQuadrature. */ Quadrature(const char* name_, int degree_, int dim_, int n_points_, VectorOfFixVecs > *lambda_, double* w_) : name(name_), degree(degree_), dim(dim_), n_points(n_points_), lambda(lambda_), w(w_) {} public: /// Copy constructor Quadrature(const Quadrature&); /// Destructor ~Quadrature(); /// Returns a Quadrature for dimension dim exact for degree degree. static Quadrature *provideQuadrature(int dim, int degree); /** \brief * Approximates an integral by the numerical quadrature described by quad; * f is a pointer to an AbstractFunction to be integrated, evaluated in * barycentric coordinates; the return value is * \f[ \sum_{k = 0}^{n_points-1} w[k] * (*f)(lambda[k]) \f] * For the approximation of \f$\int_S f\f$ we have to multiply this value * with d!|S| for a simplex S; for a parametric simplex f should be a pointer * to a function which calculates * \f$f(\lambda)|det DF_S(\hat{x}(\lambda))| \f$. */ double integrateStdSimplex(AbstractFunction > *f); /// Returns \ref name inline std::string getName() { return name; } /// Returns \ref n_points inline int getNumPoints() const { return n_points; } /// Returns \ref w[p] inline double getWeight(int p) const { return w[p]; } /// Returns \ref w. inline double* getWeight() const { return w; } /// Returns \ref dim inline int getDim() const { return dim; } /// Returns \ref degree inline int getDegree() const { return degree; } /** \brief * Returns a pointer to a vector storing the values of a doubled valued * function at all quadrature points; f is that AbstractFunction * , evaluated in barycentric coordinates; if vec is not NULL, the values are * stored in this vector, otherwise the values are stored in some static * local vector, which is overwritten on the next call */ const double *fAtQp(const AbstractFunction >& f, double *vec) const ; /** \brief * Returns a pointer to a vector storing the gradient (with respect to world * coordinates) of a double valued function at all quadrature points; * grdF is a pointer to a AbstractFunction, evaluated in barycentric * coordinates and returning a pointer to a WorldVector storing the gradient; * if vec is not NULL, the values are stored in this vector, otherwise the * values are stored in some local static vector, which is overwritten on the * next call */ const WorldVector *grdFAtQp(const AbstractFunction, DimVec >& grdF, WorldVector* vec) const; /** \brief * Returns \ref lambda[a][b] which is the b-th coordinate entry of the a-th * quadrature point */ inline double getLambda(int a, int b) const { return (lambda ? (*lambda)[a][b] : 0.0); } /** \brief * Returns \ref lambda[a] which is a DimVec containing the * coordiantes of the a-th quadrature point */ inline const DimVec& getLambda(int a) const { return (*lambda)[a]; } /// Returns \ref lambda which is a VectorOfFixvecs >. VectorOfFixVecs > *getLambda() const { return lambda; } public: /// Maximal number of quadrature points for the different dimensions static const int maxNQuadPoints[4]; protected: /// Name of this Quadrature std::string name; /// Quadrature is exact of this degree int degree; /// Quadrature for dimension dim int dim; /// Number of quadrature points int n_points; /// Vector of quadrature points given in barycentric coordinates VectorOfFixVecs > *lambda; /// Vector of quadrature weights double *w; protected: /** \brief * Initialisation of all static Quadrature objects which will be returned * by \ref provideQuadrature() */ static void initStaticQuadratures(); /** \name static quadratures, used weights, and barycentric coords * \{ */ static Quadrature **quad_nd[4]; static Quadrature *quad_0d[1]; static Quadrature *quad_1d[20]; static Quadrature *quad_2d[18]; static Quadrature *quad_3d[8]; static VectorOfFixVecs > *x_0d; static double *w_0d; static VectorOfFixVecs > *x0_1d; static VectorOfFixVecs > *x1_1d; static VectorOfFixVecs > *x2_1d; static VectorOfFixVecs > *x3_1d; static VectorOfFixVecs > *x4_1d; static VectorOfFixVecs > *x5_1d; static VectorOfFixVecs > *x6_1d; static VectorOfFixVecs > *x7_1d; static VectorOfFixVecs > *x8_1d; static VectorOfFixVecs > *x9_1d; static double *w0_1d; static double *w1_1d; static double *w2_1d; static double *w3_1d; static double *w4_1d; static double *w5_1d; static double *w6_1d; static double *w7_1d; static double *w8_1d; static double *w9_1d; static VectorOfFixVecs > *x1_2d; static VectorOfFixVecs > *x2_2d; static VectorOfFixVecs > *x3_2d; static VectorOfFixVecs > *x4_2d; static VectorOfFixVecs > *x5_2d; static VectorOfFixVecs > *x7_2d; static VectorOfFixVecs > *x8_2d; static VectorOfFixVecs > *x9_2d; static VectorOfFixVecs > *x10_2d; static VectorOfFixVecs > *x11_2d; static VectorOfFixVecs > *x12_2d; static VectorOfFixVecs > *x17_2d; static double *w1_2d; static double *w2_2d; static double *w3_2d; static double *w4_2d; static double *w5_2d; static double *w7_2d; static double *w8_2d; static double *w9_2d; static double *w10_2d; static double *w11_2d; static double *w12_2d; static double *w17_2d; static VectorOfFixVecs > *x1_3d; static VectorOfFixVecs > *x2_3d; static VectorOfFixVecs > *x3_3d; static VectorOfFixVecs > *x4_3d; static VectorOfFixVecs > *x5_3d; static VectorOfFixVecs > *x7_3d; static double *w1_3d; static double *w2_3d; static double *w3_3d; static double *w4_3d; static double *w5_3d; static double *w7_3d; /** \} */ }; /** \brief * Pre-compute the values of all basis functions at all quadrature nodes; */ const Flag INIT_PHI=1; /** \brief * Pre-compute the gradients (with respect to the barycentric coordinates) of * all basis functions at all quadrature nodes */ const Flag INIT_GRD_PHI=2; /** \brief * pre-compute all 2nd derivatives (with respect to the barycentric * coordinates) of all basis functions at all quadrature nodes; */ const Flag INIT_D2_PHI=4; /** * \ingroup Integration * *\brief * Often numerical integration involves basis functions, such as the assembling * of the system matrix and right hand side, or the integration of finite * element functions. Since numerical quadrature involves only the values at * the quadrature points and the values of basis functions and its derivatives * are the same at these points for all elements of the grid, such routines can * be much more efficient, if they can use pre-computed values of the basis * functions at the quadrature points. In this case the basis functions do not * have to be evaluated for each quadrature point on every element newly. * Information that should be pre-computed can be specified by the following * symbolic constants: * \ref INIT_PHI, \ref INIT_GRD_PHI, \ref INIT_D2_PHI */ class FastQuadrature { protected: /** \brief * Constructs a FastQuadrature for the given Quadrature, BasisFunction, and * flag. */ FastQuadrature(BasisFunction* basFcts, Quadrature* quad, Flag flag) : init_flag(flag), phi(NULL), grdPhi(NULL), D2Phi(NULL), quadrature(quad), basisFunctions(basFcts) {} /// Copy constructor FastQuadrature(const FastQuadrature&); /// Extended copy constructor FastQuadrature(const FastQuadrature&, const Flag); /// Destructor ~FastQuadrature(); public: /// Returns a FastQuadrature for the given BasisFunction, Quadrature, and flags. static FastQuadrature* provideFastQuadrature(const BasisFunction*, const Quadrature&, Flag); /// inits FastQuadrature like speciefied in flag void init(Flag init_flag); inline bool initialized(Flag flag) { if (flag == INIT_PHI) return (phi != NULL); if (flag == INIT_GRD_PHI) return (grdPhi != NULL); if (flag == INIT_D2_PHI) return (D2Phi != NULL); ERROR_EXIT("invalid flag\n"); return false; } /// Returns \ref quadrature inline const Quadrature* getQuadrature() const { return quadrature; } /// Returns \ref max_points inline int getMaxQuadPoints() { return max_points; } /// Returns (*\ref D2Phi)[q][i][j][m] const double getSecDer(int q, int i, int j, int m) const; /// Returns (*\ref D2Phi)[q] const VectorOfFixVecs > *getSecDer(int q) const; /// Returns (*\ref grdPhi)[q][i][j] inline const double getGradient(int q, int i ,int j) const { return (grdPhi) ? (*grdPhi)[q][i][j] : 0.0; } /// Returns (*\ref grdPhi)[q] inline VectorOfFixVecs >* getGradient(int q) const { return (grdPhi) ? &((*grdPhi)[q]) : NULL; } inline DimVec& getGradient(int q, int i) const { return (*grdPhi)[q][i]; } /// Returns \ref phi[q][i] inline const double getPhi(int q, int i) const { return (phi) ? phi[q][i] : 0.0; } /// Returns \ref phi[q] inline const double *getPhi(int q) const { return (phi) ? phi[q] : NULL; } /// Returns \ref quadrature ->integrateStdSimplex(f) inline double integrateStdSimplex(AbstractFunction > *f) { return quadrature->integrateStdSimplex(f); } /// Returns \ref quadrature ->getNumPoints() inline int getNumPoints() const { return quadrature->getNumPoints(); } /// Returns \ref quadrature ->getWeight(p) inline double getWeight(int p) const { return quadrature->getWeight(p); } /// Returns \ref quadrature ->getDim() inline int getDim() const { return quadrature->getDim(); } /// Returns \ref quadrature ->getDegree() inline int getDegree() const { return quadrature->getDegree(); } /// Returns \ref quadrature ->grdFAtQp(f, vec) inline const WorldVector *grdFAtQp(const AbstractFunction, DimVec >& f, WorldVector* vec) const { return quadrature->grdFAtQp(f, vec); } /// Returns \ref quadrature ->fAtQp(f, vec) inline const double *fAtQp(const AbstractFunction >& f, double *vec) const { return quadrature->fAtQp(f, vec); } /// Returns \ref quadrature ->getLambda(a,b) inline double getLambda(int a,int b) const { return quadrature->getLambda(a,b); } /// Returns \ref quadrature ->getLambda(a) inline const DimVec& getLambda(int a) const { return quadrature->getLambda(a); } /// Returns \ref basisFunctions inline BasisFunction* getBasisFunctions() const { return basisFunctions; } protected: /** \brief * Specifies which information should be pre-computed. Can be \ref INIT_PHI, * \ref INIT_GRD_PHI, or \ref INIT_D2_PHI */ Flag init_flag; /** \brief * Matrix storing function values if the flag \ref INIT_PHI is set; * phi[i][j] stores the value \ref basisFunctions->phi[j] * (quadrature->lambda[i]), 0 <= j < basisFunctions->getNumber() and * 0 <= i < n_points */ double **phi; /** \brief * Matrix storing all gradients (with respect to the barycentric coordinates) * if the flag \ref INIT_GRD_PHI is set; grdPhi[i][j][k] stores the value * basisFunctions->grdPhi[j](quadrature->lambda[i])[k] * for 0 <= j < basisFunctions->getNumber(), * 0 <= i < . . . , n_points, and 0 <= k < DIM */ MatrixOfFixVecs > *grdPhi; /** \brief * Matrix storing all second derivatives (with respect to the barycentric * coordinates) if the flag \ref INIT_D2_PHI is set; D2Phi[i][j][k][l] stores * the value basisFunctions->D2Phi[j](quadrature->lambda[i])[k][l] * for 0 <= j < basisFunctions->getNumber(), * 0 <= i < n_points, and 0 <= k,l < DIM */ MatrixOfFixVecs > *D2Phi; /** \brief * List of all used FastQuadratures */ static std::list fastQuadList; /** \brief * Maximal number of quadrature points for all yet initialised FastQuadrature * objects. This value may change after a new initialisation of a * FastQuadrature */ static int max_points; /// This FastQuadrature stores values for Quadrature quadrature Quadrature* quadrature; /// Values stored for basis functions basisFunctions BasisFunction* basisFunctions; }; } #endif // AMDIS_QUADRATURE_H