// ============================================================================ // == == // == 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 Operator.h */ #ifndef AMDIS_OPERATOR_H #define AMDIS_OPERATOR_H #include #include "AMDiS_fwd.h" #include "FixVec.h" #include "Flag.h" #include "MatrixVector.h" #include "ElInfo.h" #include "AbstractFunction.h" #include "OpenMP.h" #include "SubAssembler.h" namespace AMDiS { /** * \ingroup Assembler * * \brief * Base class for ZeroOrderTerm, FirstOrderTerm and SecondOrderTerm. * OperatorTerms are the building blocks of an Operator. Each OperatorTerm * has its properties which are regarded, when constructing * an Assembler for the corresponding Operator. */ class OperatorTerm { public: /** \brief * Constructs an OperatorTerm with initially no properties. * degree_ is used to determine the degree of the needed quadrature * for the assemblage. */ OperatorTerm(int deg) : properties(0), degree(deg), auxFESpaces(0) {} /// Destructor. virtual ~OperatorTerm() {} /** \brief * Virtual method. It's called by SubAssembler::initElement() for * each OperatorTerm belonging to this SubAssembler. E.g., vectors * and coordinates at quadrature points can be calculated here. */ virtual void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL) {} virtual void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo, SubAssembler*, Quadrature *quad = NULL) {} /// Returs \auxFESpaces, the list of all aux fe spaces the operator makes use off. std::vector& getAuxFESpaces() { return auxFESpaces; } /// Specifies whether the matrix of the term is symmetric void setSymmetric(bool symm); /// Returns true, if the term is piecewise constant, returns false otherwise inline bool isPWConst() { return (degree == 0); } /// Returns true, if the term has a symmetric matrix, returns false otherwise. bool isSymmetric(); /// Returns \ref degree. inline int getDegree() { return degree; } /// Evaluation of the OperatorTerm at all quadrature points. virtual void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const = 0; /** \brief * Determines the value of a dof vector at the quadrature points of a given * element. It is used by all VecAtQP like operator terms. */ double *getVectorAtQPs(DOFVectorBase* vec, const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad); /** \brief * Determines the value of a dof vector at the quadrature points of a given * element. This function is used, if an operator is assembled on two different * meshes using the dual traverse. The element, i.e. the small or the large one, * is choosen, which corresponds to the mesh the dof vector is defined on. */ double *getVectorAtQPs(DOFVectorBase* vec, const ElInfo* smallElInfo, const ElInfo* largeElInfo, SubAssembler* subAssembler, Quadrature *quad); /// WorldVector* getGradientsAtQPs(DOFVectorBase* vec, const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad); /// WorldVector* getGradientsAtQPs(DOFVectorBase* vec, const ElInfo* smallElInfo, const ElInfo* largeElInfo, SubAssembler* subAssembler, Quadrature *quad); protected: /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$. void lalt(const DimVec >& Lambda, const WorldMatrix& matrix, DimMat& LALt, bool symm, double factor) const; /** \brief * Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for \f$ A \f$ * the matrix having a ONE in the position \f$ (K,L) \f$ * and ZEROS in all other positions. */ static void lalt_kl(const DimVec >& Lambda, int k, int l, DimMat& LALt, double factor); /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity. void l1lt(const DimVec >& Lambda, DimMat& LALt, double factor) const; /// Evaluation of \f$ \Lambda \cdot b\f$. static inline void lb(const DimVec >& Lambda, const WorldVector& b, DimVec& Lb, double factor) { const int dimOfWorld = Global::getGeo(WORLD); const int dim = Lambda.size() - 1; for (int i = 0; i <= dim; i++) { double val = 0.0; for (int j = 0; j < dimOfWorld; j++) val += Lambda[i][j] * b[j]; Lb[i] += val * factor; } } /** \brief * Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in * each component. */ static inline void l1(const DimVec >& Lambda, DimVec& Lb, double factor) { const int dimOfWorld = Global::getGeo(WORLD); const int dim = Lambda.size() - 1; for (int i = 0; i <= dim; i++) { double val = 0.0; for (int j = 0; j < dimOfWorld; j++) val += Lambda[i][j]; Lb[i] += val * factor; } } protected: /// Stores the properties of this OperatorTerm Flag properties; /// Polynomial degree of the term. Used to detemine the degree of the quadrature. int degree; /// List off all fe spaces, the operator term makes use off. std::vector auxFESpaces; /// Pointer to the Operator this OperatorTerm belongs to. Operator* operat; /// Flag for piecewise constant terms static const Flag PW_CONST; /// Flag for symmetric terms static const Flag SYMMETRIC; friend class SubAssembler; friend class ZeroOrderAssembler; friend class FirstOrderAssembler; friend class SecondOrderAssembler; friend class Operator; }; /** * \ingroup Assembler * * \brief * Describes the second order terms: \f$ \nabla \cdot A \nabla u(\vec{x}) \f$ */ class SecondOrderTerm : public OperatorTerm { public: /// Constructor. SecondOrderTerm(int deg) : OperatorTerm(deg) {} /// Destructor. virtual ~SecondOrderTerm() {} /// Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points. virtual void getLALt(const ElInfo *elInfo, int nPoints, DimMat **result) const = 0; /// Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points. virtual void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const = 0; }; /** * \ingroup Assembler * * \brief * Implements the laplace operator: \f$ \Delta u(\vec{x}) \f$ */ class Laplace_SOT : public SecondOrderTerm { public: /// Constructor. Laplace_SOT() : SecondOrderTerm(0) { setSymmetric(true); } /// Implenetation of SecondOrderTerm::getLALt(). inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const { const DimVec > &Lambda = elInfo->getGrdLambda(); for (int iq = 0; iq < nPoints; iq++) l1lt(Lambda, *(LALt[iq]), 1.0); } /// Implementation of SecondOrderTerm::eval(). inline void eval(int nPoints, const double *, const WorldVector *, const WorldMatrix *D2UhAtQP, double *result, double factor) const { int dow = Global::getGeo(WORLD); if (D2UhAtQP) { for (int iq = 0; iq < nPoints; iq++) { double resultQP = 0.0; for (int i = 0; i < dow; i++) resultQP += D2UhAtQP[iq][i][i]; result[iq] += factor * resultQP; } } } /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const { if (grdUhAtQP) for (int iq = 0; iq < nPoints; iq++) result[iq] += grdUhAtQP[iq]; } }; /** * \ingroup Assembler * * \brief * Implements the laplace operator multiplied with a scalar factor: * \f$ f \cdot \Delta u(\vec{x}) \f$ */ class FactorLaplace_SOT : public SecondOrderTerm { public: /// Constructor. FactorLaplace_SOT(double f) : SecondOrderTerm(0) { factor = new double; *factor = f; setSymmetric(true); } /// Constructor. FactorLaplace_SOT(double *fptr) : SecondOrderTerm(0), factor(fptr) { setSymmetric(true); } /// Implements SecondOrderTerm::getLALt(). inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const { const DimVec > &Lambda = elInfo->getGrdLambda(); for (int iq = 0; iq < nPoints; iq++) l1lt(Lambda, *(LALt[iq]), (*factor)); } /// Implenetation of SecondOrderTerm::eval(). void eval(int nPoints, const double *, const WorldVector *, const WorldMatrix *D2UhAtQP, double *result, double f) const { int dow = Global::getGeo(WORLD); if (D2UhAtQP) { for (int iq = 0; iq < nPoints; iq++) { double resultQP = 0.0; for (int i = 0; i < dow; i++) { resultQP += D2UhAtQP[iq][i][i]; } result[iq] += resultQP * f * (*factor); } } } /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const { if (grdUhAtQP) for (int iq = 0; iq < nPoints; iq++) axpy(*factor, grdUhAtQP[iq], result[iq]); } private: /// Pointer to the factor. double *factor; }; /** * \ingroup Assembler * * \brief * SecondOrderTerm where A is a function which maps a DOFVector evaluated at * a given quadrature point to a WolrdMatrix: * \f$ \nabla \cdot A(v(\vec{x})) \nabla u(\vec{x}) \f$ */ class MatrixFct_SOT : public SecondOrderTerm { public: /// Constructor. MatrixFct_SOT(DOFVectorBase *dv, AbstractFunction, double> *fct, AbstractFunction, WorldMatrix > *div, bool sym = false); /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /// Implements SecondOrderTerm::getLALt(). void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /// Implenetation of SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; /// Pointer to the values of the DOFVector at quadrature points. double* vecAtQPs; /// Function for A. AbstractFunction, double>* matrixFct; /// AbstractFunction, WorldMatrix >* divFct; /// True, if \ref matrixFct produces always symmetric matrices. bool symmetric; }; /** * \ingroup Assembler * * \brief * SecondOrderTerm where A is a given fixed WorldMatrix: * \f$ \nabla \cdot A \nabla u(\vec{x}) \f$ */ class Matrix_SOT : public SecondOrderTerm { public: /// Constructor Matrix_SOT(WorldMatrix mat) : SecondOrderTerm(0), matrix(mat) { symmetric = matrix.isSymmetric(); setSymmetric(symmetric); } /// Implements SecondOrderTerm::getLALt(). inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const { const DimVec >& Lambda = elInfo->getGrdLambda(); for (int iq = 0; iq < nPoints; iq++) lalt(Lambda, matrix, *(LALt[iq]), symmetric, 1.0); } /// Implenetation of SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: /// Matrix stroring A. WorldMatrix matrix; /// True, if \ref matrix is symmetric. bool symmetric; }; /** * \ingroup Assembler * * \brief * SecondOrderTerm where A is a WorldMatrix having a ONE in position IJ * and ZERO in all other positions * \f$ \nabla \cdot A \nabla u(\vec{x}) \f$ */ class FactorIJ_SOT : public SecondOrderTerm { public: /// Constructor. FactorIJ_SOT(int x_i, int x_j, double f) : SecondOrderTerm(0), xi(x_i), xj(x_j) { factor = new double; *factor = f; setSymmetric(xi == xj); } /// Constructor. FactorIJ_SOT(int x_i, int x_j, double *fptr) : SecondOrderTerm(0), xi(x_i), xj(x_j), factor(fptr) { setSymmetric(xi == xj); } /// Implements SecondOrderTerm::getLALt(). inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const { const DimVec > &Lambda = elInfo->getGrdLambda(); for (int iq = 0; iq < nPoints; iq++) lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*factor)); } /** \brief * Implenetation of SecondOrderTerm::eval(). */ void eval(int nPoints, const double *, const WorldVector *, const WorldMatrix *D2UhAtQP, double *result, double fac) const { if (D2UhAtQP) for (int iq = 0; iq < nPoints; iq++) result[iq] += (*factor) * D2UhAtQP[iq][xi][xj] * fac; } /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const { if (grdUhAtQP) for (int iq = 0; iq < nPoints; iq++) result[iq][xi] += (*factor) * grdUhAtQP[iq][xj]; } private: /** \brief * Directions for the derivatives. */ int xi, xj; /** \brief * Pointer to the factor. */ double *factor; }; /** * \ingroup Assembler * * \brief * Laplace operator multiplied with a function evaluated at the quadrature * points of a given DOFVector: * \f$ f(v(\vec{x})) \Delta u(\vec{x}) \f$ */ class VecAtQP_SOT : public SecondOrderTerm { public: /// Constructor. VecAtQP_SOT(DOFVectorBase *dv, AbstractFunction *af); /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes. void initElement(const ElInfo* smallElInfo, const ElInfo* largeElInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /// Implements SecondOrderTerm::getLALt(). void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /// Implenetation of SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; /// Pointer to an array containing the DOFVector evaluated at quadrature points. double* vecAtQPs; /// Function evaluated at \ref vecAtQPs. AbstractFunction *f; }; /** * \ingroup Assembler * * \brief * Laplace operator multiplied with a function evaluated at the quadrature * points of a given DOFVector: * \f$ f(v(\vec{x}), w(\vec{x})) \Delta u(\vec{x}) \f$ */ class Vec2AtQP_SOT : public SecondOrderTerm { public: /// Constructor. Vec2AtQP_SOT(DOFVectorBase *dv1, DOFVectorBase *dv2, BinaryAbstractFunction *af); /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /// Implements SecondOrderTerm::getLALt(). void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /// Implenetation of SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec1; DOFVectorBase* vec2; /// Pointer to an array containing the DOFVector evaluated at quadrature points. double* vecAtQPs1; double* vecAtQPs2; /// Function evaluated at \ref vecAtQPs. BinaryAbstractFunction *f; }; /** * \ingroup Assembler * * \brief * Laplace multiplied with a function evaluated at each quadrature point: * \f$ f(\vec{x}) \Delta u(\vec{x}) \f$ */ class CoordsAtQP_SOT : public SecondOrderTerm { public: /// Constructor. CoordsAtQP_SOT(AbstractFunction > *af) : SecondOrderTerm(af->getDegree()), g(af) { setSymmetric(true); } /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements SecondOrderTerm::getLALt(). */ void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /** \brief * Implenetation of SecondOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: /** \brief * Stores coordinates at quadrature points. Set in \ref initElement(). */ WorldVector* coordsAtQPs; /** \brief * Function evaluated at quadrature points. */ AbstractFunction > *g; }; /** * \ingroup Assembler * * \brief * SecondOrderTerm where A is a function which maps the gradient of a * DOFVector at each quadrature point to WorldMatrix: * \f$ \nabla \cdot A(\nabla v(\vec{x})) \nabla u(\vec{x})\f$ */ class MatrixGradient_SOT : public SecondOrderTerm { public: /// Constructor. MatrixGradient_SOT(DOFVectorBase *dv, AbstractFunction, WorldVector > *af, AbstractFunction, WorldMatrix > *divAf, bool symm = false); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements SecondOrderTerm::getLALt(). */ void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /** \brief * Implenetation of SecondOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: DOFVectorBase* vec; /** \brief * Function which maps each entry in \ref gradAtQPs to a WorldMatrix. */ AbstractFunction, WorldVector > *f; AbstractFunction, WorldMatrix > *divFct; /** \brief * Pointer to a WorldVector array containing the gradients of the DOFVector * at each quadrature point. */ WorldVector* gradAtQPs; /** \brief * True, if \ref f provides always a symmetric WorldMatrix. */ bool symmetric; }; /** * \ingroup Assembler * * \brief * Laplace multiplied with a function which maps the gradient of a DOFVector * at each quadrature point to a double: * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$ */ class FctGradient_SOT : public SecondOrderTerm { public: /// Constructor. FctGradient_SOT(DOFVectorBase *dv, AbstractFunction > *af); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements SecondOrderTerm::getLALt(). */ void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /** \brief * Implenetation of SecondOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: DOFVectorBase* vec; /** \brief * Function wich maps \ref gradAtQPs to a double. */ AbstractFunction > *f; /** \brief * Pointer to a WorldVector array containing the gradients of the DOFVector * at each quadrature point. */ WorldVector* gradAtQPs; }; /** * \ingroup Assembler * * \brief * Laplace multiplied with a function which maps the gradient of a DOFVector * at each quadrature point to a double: * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$ */ class VecAndGradientAtQP_SOT : public SecondOrderTerm { public: /// Constructor. VecAndGradientAtQP_SOT(DOFVectorBase *dv, BinaryAbstractFunction > *af); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements SecondOrderTerm::getLALt(). */ void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /** \brief * Implenetation of SecondOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: DOFVectorBase* vec; /** \brief * Function wich maps \ref gradAtQPs to a double. */ BinaryAbstractFunction > *f; /** \brief * Pointer to a WorldVector array containing the gradients of the DOFVector * at each quadrature point. */ double* vecAtQPs; WorldVector* gradAtQPs; }; /** * \ingroup Assembler * * \brief * Laplace multiplied with a function which maps the gradient of a DOFVector * at each quadrature point to a double: * \f$ \nabla \cdot A(v(\vec{x}), \nabla v(\vec{x})) \nabla u(\vec{x}) \f$ */ class VecMatrixGradientAtQP_SOT : public SecondOrderTerm { public: /// Constructor. VecMatrixGradientAtQP_SOT(DOFVectorBase *dv, BinaryAbstractFunction, double, WorldVector > *af, AbstractFunction, WorldMatrix > *divAf, bool symm = false); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements SecondOrderTerm::getLALt(). */ void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /** \brief * Implenetation of SecondOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: DOFVectorBase* vec; /** \brief * Function wich maps \ref gradAtQPs to a double. */ BinaryAbstractFunction, double, WorldVector > *f; AbstractFunction, WorldMatrix > *divFct; /** \brief * Pointer to a WorldVector array containing the gradients of the DOFVector * at each quadrature point. */ double* vecAtQPs; WorldVector* gradAtQPs; bool symmetric; }; /** * \ingroup Assembler * * \brief * Laplace multiplied with a function which maps the gradient of a DOFVector * at each quadrature point to a double: * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$ */ class VecGradCoordsAtQP_SOT : public SecondOrderTerm { public: /// Constructor. VecGradCoordsAtQP_SOT(DOFVectorBase *dv, TertiaryAbstractFunction, WorldVector > *af); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements SecondOrderTerm::getLALt(). */ void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /** \brief * Implenetation of SecondOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: DOFVectorBase* vec; /** \brief * Function wich maps \ref gradAtQPs to a double. */ TertiaryAbstractFunction, WorldVector > *f; /** \brief * Pointer to a WorldVector array containing the gradients of the DOFVector * at each quadrature point. */ double* vecAtQPs; WorldVector* gradAtQPs; WorldVector* coordsAtQPs; }; /** * \ingroup Assembler * * \brief * Laplace operator multiplied with a function evaluated at the quadrature * points of a given DOFVector: * \f$ -f(v(\vec{x})) \Delta u(\vec{x}) \f$ */ class VecAndCoordsAtQP_SOT : public SecondOrderTerm { public: /// Constructor. VecAndCoordsAtQP_SOT(DOFVectorBase *dv, BinaryAbstractFunction > *af); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements SecondOrderTerm::eval(). */ void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: /** \brief * DOFVector to be evaluated at quadrature points. */ DOFVectorBase* vec; /** \brief * Pointer to an array containing the DOFVector evaluated at quadrature * points. */ double* vecAtQPs; WorldVector* coordsAtQPs; /** \brief * Function evaluated at \ref vecAtQPs. */ BinaryAbstractFunction > *f; }; /** * \ingroup Assembler * * \brief * SecondOrderTerm where A is a function which maps the gradient of a * DOFVector at each quadrature point to WorldMatrix: * \f$ \nabla \cdot A(\nabla v(\vec{x})) \nabla u(\vec{x})\f$ */ class MatrixGradientAndCoords_SOT : public SecondOrderTerm { public: /// Constructor. MatrixGradientAndCoords_SOT(DOFVectorBase *dv, BinaryAbstractFunction, WorldVector , WorldVector > *af, AbstractFunction, WorldMatrix > *divAf, bool symm = false); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements SecondOrderTerm::getLALt(). */ void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /** \brief * Implenetation of SecondOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: DOFVectorBase* vec; /** \brief * Function which maps each entry in \ref gradAtQPs to a WorldMatrix. */ BinaryAbstractFunction, WorldVector, WorldVector > *f; AbstractFunction, WorldMatrix > *divFct; /** \brief * Pointer to a WorldVector array containing the gradients of the DOFVector * at each quadrature point. */ WorldVector* gradAtQPs; WorldVector* coordsAtQPs; /** \brief * True, if \ref f provides always a symmetric WorldMatrix. */ bool symmetric; }; class MatrixVec2_SOT : public SecondOrderTerm { public: /// Constructor. MatrixVec2_SOT(DOFVectorBase *dv1, DOFVectorBase *dv2, BinaryAbstractFunction *f, WorldMatrix Af, bool sym = false); /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /// Implements SecondOrderTerm::getLALt(). void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /// Implenetation of SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec1; DOFVectorBase* vec2; /// Pointer to the values of the DOFVector at quadrature points. double* vec1AtQPs; double* vec2AtQPs; /// Function for A. BinaryAbstractFunction * fct; /// WorldMatrix A; /// True, if \ref matrixFct produces always symmetric matrices. bool symmetric; }; class General_SOT : public SecondOrderTerm { public: /// Constructor. General_SOT(std::vector*> vecs, std::vector*> grads, TertiaryAbstractFunction, WorldVector, std::vector, std::vector > > *f, AbstractFunction, WorldMatrix > *divFct, bool symmetric); /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo*, SubAssembler*, Quadrature *quad= NULL); /// Implements SecondOrderTerm::getLALt(). void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /// Implenetation of SecondOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /// Implenetation of SecondOrderTerm::weakEval(). void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: std::vector*> vecs_; std::vector*> grads_; TertiaryAbstractFunction, WorldVector, std::vector, std::vector > > *f_; AbstractFunction, WorldMatrix > *divFct_; WorldVector *coordsAtQPs_; std::vector vecsAtQPs_; std::vector*> gradsAtQPs_; bool symmetric_; }; class GeneralParametric_SOT : public SecondOrderTerm { public: /// Constructor. GeneralParametric_SOT(std::vector*> vecs, std::vector*> grads, QuartAbstractFunction, WorldVector, WorldVector, std::vector, std::vector > > *f, AbstractFunction, WorldMatrix > *divFct, bool symmetric); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo*, SubAssembler* , Quadrature *quad= NULL); /** \brief * Implements SecondOrderTerm::getLALt(). */ void getLALt(const ElInfo *elInfo, int nPoints, DimMat **LALt) const; /** \brief * Implenetation of SecondOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; /** \brief * Implenetation of SecondOrderTerm::weakEval(). */ void weakEval(int nPoints, const WorldVector *grdUhAtQP, WorldVector *result) const; protected: std::vector*> vecs_; std::vector*> grads_; QuartAbstractFunction, WorldVector, WorldVector, std::vector, std::vector > > *f_; AbstractFunction, WorldMatrix > *divFct_; WorldVector *coordsAtQPs_; WorldVector elementNormal_; std::vector vecsAtQPs_; std::vector*> gradsAtQPs_; bool symmetric_; }; // ============================================================================ // ===== class FirstOrderTerm ================================================= // ============================================================================ /** * \ingroup Assembler * * \brief * Describes the first order terms: \f$ b \cdot \nabla u(\vec{x}) \f$ */ class FirstOrderTerm : public OperatorTerm { public: /// Constructor. FirstOrderTerm(int deg) : OperatorTerm(deg) {} /// Destructor. virtual ~FirstOrderTerm() {} /// Evaluation of \f$ \Lambda b \f$. virtual void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& result) const = 0; /// Implenetation of FirstOrderTerm::eval(). void eval(int nPoints, const double *, const WorldVector *grdUhAtQP, const WorldMatrix *, double *result, double factor) const { int dow = Global::getGeo(WORLD); if (grdUhAtQP) { for (int iq = 0; iq < nPoints; iq++) { double resultQP = 0.0; for (int i = 0; i < dow; i++) resultQP += grdUhAtQP[iq][i]; result[iq] += resultQP * factor; } } } }; /** * \ingroup Assembler * * \brief * First order term: \f$ b \cdot \nabla u(\vec{x}) \f$ with a vector b which * only consits of entries equal to one. */ class Simple_FOT : public FirstOrderTerm { public: /// Constructor. Simple_FOT() : FirstOrderTerm(0) {} /// Implements FirstOrderTerm::getLb(). inline void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& Lb) const { const DimVec > &Lambda = elInfo->getGrdLambda(); for (int iq = 0; iq < nPoints; iq++) l1(Lambda, Lb[iq], 1.0); } }; /** * \ingroup Assembler * * \brief * First order term: \f$ b \cdot \nabla u(\vec{x}) \f$ with a vector b which * only consits of entries equal to one. */ class FactorSimple_FOT : public FirstOrderTerm { public: /// Constructor. FactorSimple_FOT(double f) : FirstOrderTerm(0) { factor = new double; *factor = f; } /// Constructor. FactorSimple_FOT(double *fptr) : FirstOrderTerm(0), factor(fptr) {} /** \brief * Implements FirstOrderTerm::getLb(). */ inline void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& Lb) const { const DimVec > &Lambda = elInfo->getGrdLambda(); for (int iq = 0; iq < nPoints; iq++) l1(Lambda, Lb[iq], (*factor)); } private: /** \brief * Pointer to the factor. */ double *factor; }; /** * \ingroup Assembler * * \brief * First order term: \f$ b \cdot \nabla u(\vec{x}) \f$ with a given vector b. */ class Vector_FOT : public FirstOrderTerm { public: /// Constructor. Vector_FOT(WorldVector wv) : FirstOrderTerm(0), b(wv) {} /** \brief * Implements FirstOrderTerm::getLb(). */ inline void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >&Lb) const { const DimVec > &Lambda = elInfo->getGrdLambda(); for (int iq = 0; iq < nPoints; iq++) lb(Lambda, b, Lb[iq], 1.0); } /// Implements FirstOrderTerm::eval(). void eval(int nPoints, const double *, const WorldVector *grdUhAtQP, const WorldMatrix *, double *result, double factor) const { if (grdUhAtQP) for (int iq = 0; iq < nPoints; iq++) result[iq] += b * grdUhAtQP[iq] * factor; } protected: /// Vector which is multiplied with \f$ \nabla u(\vec{x}) \f$ WorldVector b; }; /** * \ingroup Assembler * * \brief * Simple_FOT multiplied with \f$ f(v(\vec{x})) \f$. */ class VecAtQP_FOT : public FirstOrderTerm { public: /// Constructor. VecAtQP_FOT(DOFVectorBase *dv, AbstractFunction *af, WorldVector *wv); /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /// Implements FirstOrderTerm::getLb(). void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& Lb) const; /// Implements FirstOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec; /// v at quadrature points. double *vecAtQPs; /// Function f. AbstractFunction *f; /// WorldVector *b; }; /** * \ingroup Assembler * * \brief * Simple_FOT multiplied with \f$ f(\vec{x}) \f$. */ class CoordsAtQP_FOT : public FirstOrderTerm { public: /// Constructor. CoordsAtQP_FOT(AbstractFunction > *af) : FirstOrderTerm(af->getDegree()), g(af) {} /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements FistOrderTerm::getLb(). */ void getLb(const ElInfo *elInfo, int nPoints,VectorOfFixVecs >&Lb) const; /** \brief * Implements FirstOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: /** \brief * Stores coordinates at quadrature points. Set in \ref initElement(). */ WorldVector* coordsAtQPs; /** \brief * Function avaluated at world coordinates. */ AbstractFunction > *g; }; /** * \ingroup Assembler * * \brief * Vector_FOT multiplied with \f$ f(\vec{x}) \f$. */ class VecCoordsAtQP_FOT : public FirstOrderTerm { public: /// Constructor. VecCoordsAtQP_FOT(AbstractFunction > *af, WorldVector wv) : FirstOrderTerm(af->getDegree()), g(af), b(wv) {} /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements FistOrderTerm::getLb(). */ void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& Lb) const; /** \brief * Implements FirstOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: /** \brief * Stores coordinates at quadrature points. Set in \ref initElement(). */ WorldVector* coordsAtQPs; /** \brief * Function evaluated at world coordinates. */ AbstractFunction > *g; /** \brief * Coefficient vector. */ WorldVector b; }; /** * \ingroup Assembler * * \brief * First order term: \f$ b(\nabla v(\vec{x})) \cdot \nabla u(\vec{x})\f$. */ class VectorGradient_FOT : public FirstOrderTerm { public: /// Constructor. VectorGradient_FOT(DOFVectorBase *dv, AbstractFunction, WorldVector > *af); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements FirstOrderTerm::getLb(). */ void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& Lb) const; /** \brief * Implements FirstOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: DOFVectorBase* vec; /** \brief * Function for b. */ AbstractFunction, WorldVector > *f; /** \brief * Gradient of v at quadrature points. */ WorldVector *gradAtQPs; }; /** * \ingroup Assembler * * \brief * First order term: \f$ b(v(\vec{x})) \cdot \nabla u(\vec{x})\f$. */ class VectorFct_FOT : public FirstOrderTerm { public: /// Constructor. VectorFct_FOT(DOFVectorBase *dv, AbstractFunction, double> *fct); /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements FirstOrderTerm::getLb(). */ void getLb(const ElInfo *elInfo, int qPoint, VectorOfFixVecs >& Lb) const; /** \brief * Implements FirstOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: /** \brief * DOFVector to be evaluated at quadrature points. */ DOFVectorBase* vec; /** \brief * Vector v at quadrature points. */ double *vecAtQPs; /** \brief * Function for b. */ AbstractFunction, double> *vecFct; }; /** * \ingroup Assembler * * \brief * */ class VecFctAtQP_FOT : public FirstOrderTerm { public: /// Constructor. VecFctAtQP_FOT(AbstractFunction, WorldVector > *af) : FirstOrderTerm(af->getDegree()), g(af) {} /** \brief * Implementation of \ref OperatorTerm::initElement(). */ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /** \brief * Implements FistOrderTerm::getLb(). */ void getLb(const ElInfo *elInfo, int nPoints,VectorOfFixVecs >&Lb) const; /** \brief * Implements FirstOrderTerm::eval(). */ void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: /** \brief * Stores coordinates at quadrature points. Set in \ref initElement(). */ WorldVector* coordsAtQPs; /** \brief * Function avaluated at world coordinates. */ AbstractFunction, WorldVector > *g; }; /** * \ingroup Assembler * * \brief * First order term: \f$ b(v(\vec{x}), \nabla w(\vec{x})) \cdot \nabla u(\vec{x})\f$. */ class VecGrad_FOT : public FirstOrderTerm { public: /// Constructor. VecGrad_FOT(DOFVectorBase *dv1, DOFVectorBase *dv2, BinaryAbstractFunction, double, WorldVector > *fct); /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes. void initElement(const ElInfo* smallElInfo, const ElInfo* largeElInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /// Implements FirstOrderTerm::getLb(). void getLb(const ElInfo *elInfo, int qPoint, VectorOfFixVecs >& Lb) const; /// Implements FirstOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: /** \brief * DOFVector to be evaluated at quadrature points. */ DOFVectorBase* vec1; DOFVectorBase* vec2; /** \brief * Vector v at quadrature points. */ double *vecAtQPs; WorldVector *gradAtQPs; /** \brief * Function for b. */ BinaryAbstractFunction, double, WorldVector > *vecFct; }; /** * \ingroup Assembler * * \brief * First order term: \f$ b(\nabla v(\vec{x}), \nabla w(\vec{x})) \cdot \nabla u(\vec{x})\f$. */ class FctGrad2_FOT : public FirstOrderTerm { public: /// Constructor FctGrad2_FOT(DOFVectorBase *dv1, DOFVectorBase *dv2, BinaryAbstractFunction, WorldVector, WorldVector > *vecFct_) : FirstOrderTerm(vecFct_->getDegree()), vec1(dv1), vec2(dv2), vecFct(vecFct_) {} /// Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); /// Implements FirstOrderTerm::getLb(). void getLb(const ElInfo *elInfo, int qPoint, VectorOfFixVecs >& Lb) const; /// Implements FirstOrderTerm::eval(). void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec1; /// DOFVector to be evaluated at quadrature points. DOFVectorBase* vec2; /// Gradient v at quadrature points. WorldVector *grad1AtQPs; /// Gradient v at quadrature points. WorldVector *grad2AtQPs; /// Function for b. BinaryAbstractFunction, WorldVector, WorldVector > *vecFct; }; class Vec2AndGradAtQP_FOT : public FirstOrderTerm { public: Vec2AndGradAtQP_FOT(DOFVectorBase *dv1, DOFVectorBase *dv2, TertiaryAbstractFunction > *f_, WorldVector *b_); void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& Lb) const; void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: DOFVectorBase* vec1; DOFVectorBase* vec2; double *vec1AtQPs; double *vec2AtQPs; WorldVector *gradAtQPs1; TertiaryAbstractFunction > *f; WorldVector *b; }; class FctVecAtQP_FOT : public FirstOrderTerm { public: FctVecAtQP_FOT(DOFVectorBase *dv, AbstractFunction > *f_, WorldVector *b_); void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& Lb) const; void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: DOFVectorBase* vec; double *vecAtQPs; WorldVector *coordsAtQPs; AbstractFunction > *f; WorldVector *b; }; class Vec2AtQP_FOT : public FirstOrderTerm { public: Vec2AtQP_FOT(DOFVectorBase *dv1, DOFVectorBase *dv2, WorldVector *b_); void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& Lb) const; void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: DOFVectorBase* vec1; DOFVectorBase* vec2; double *vec1AtQPs; double *vec2AtQPs; WorldVector *b; }; class Vec3FctAtQP_FOT : public FirstOrderTerm { public: Vec3FctAtQP_FOT(DOFVectorBase *dv1, DOFVectorBase *dv2,DOFVectorBase *dv3, BinaryAbstractFunction *f_, WorldVector *b_); void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, Quadrature *quad = NULL); void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs >& Lb) const; void eval(int nPoints, const double *uhAtQP, const WorldVector *grdUhAtQP, const WorldMatrix *D2UhAtQP, double *result, double factor) const; protected: DOFVectorBase* vec1; DOFVectorBase* vec2; DOFVectorBase* vec3; BinaryAbstractFunction* f; double *vec1AtQPs; double *vec2AtQPs; double *vec3AtQPs; WorldVector *b; }; class General_FOT : public FirstOrderTerm { public: /// Constructor General_FOT(std::vector*> vecs, std::vector*> grads, TertiaryAbstractFunction, WorldVector