// ============================================================================ // == == // == AMDiS - Adaptive multidimensional simulations == // == == // == http://www.amdis-fem.org == // == == // ============================================================================ // // Software License for AMDiS // // Copyright (c) 2010 Dresden University of Technology // All rights reserved. // Authors: Simon Vey, Thomas Witkowski et al. // // This file is part of AMDiS // // See also license.opensource.txt in the distribution. /** \file DOFVector.h */ #ifndef AMDIS_DOFVECTOR_H #define AMDIS_DOFVECTOR_H #include #include #include #include "AMDiS_fwd.h" #include "FixVec.h" #include "Global.h" #include "Flag.h" #include "RCNeighbourList.h" #include "DOFIterator.h" #include "DOFIndexed.h" #include "DOFContainer.h" #include "Boundary.h" #include "CreatorInterface.h" #include "Serializable.h" #include "DOFMatrix.h" #include "BasisFunction.h" #include "FiniteElemSpace.h" #include "SurfaceQuadrature.h" #ifdef HAVE_PARALLEL_DOMAIN_AMDIS #include #include "parallel/ParallelDofMapping.h" #endif namespace AMDiS { using namespace std; template class DOFVectorBase : public DOFIndexed { public: DOFVectorBase() : feSpace(NULL), elementVector(3), boundaryManager(NULL), nBasFcts(0) { #ifdef HAVE_PARALLEL_DOMAIN_AMDIS dofMap = NULL; #endif } DOFVectorBase(const FiniteElemSpace *f, string n); virtual ~DOFVectorBase(); /// For the given element, this function returns an array of all DOFs of /// this DOFVector that are defined on this element. virtual void getLocalVector(const Element *el, mtl::dense_vector& localVec) const; /// Evaluates the DOF vector at a set of quadrature points defined on the /// given element. void getVecAtQPs(const ElInfo *elInfo, const Quadrature *quad, const FastQuadrature *quadFast, mtl::dense_vector& vecAtQPs) const; void getVecAtQPs(const ElInfo *smallElInfo, const ElInfo *largeElInfo, const Quadrature *quad, const FastQuadrature *quadFast, mtl::dense_vector& vecAtQPs) const; /// Evaluates the gradient of a DOF vector at a set of quadrature points /// defined on the given element. void getGrdAtQPs( const ElInfo *elInfo, const Quadrature *quad, const FastQuadrature *quadFast, mtl::dense_vector::type> &grdAtQPs) const; void getGrdAtQPs( const ElInfo *smallElInfo, const ElInfo *largeElInfo, const Quadrature *quad, const FastQuadrature *quadFast, mtl::dense_vector::type> &grdAtQPs) const; /// Evaluates the comp'th component of the derivative of a DOF vector at a /// set of quadrature points defined on the given element. void getDerivativeAtQPs( const ElInfo *elInfo, const Quadrature *quad, const FastQuadrature *quadFast, int comp, mtl::dense_vector &derivativeAtQPs) const; void getDerivativeAtQPs( const ElInfo *smallElInfo, const ElInfo *largeElInfo, const Quadrature *quad, const FastQuadrature *quadFast, int comp, mtl::dense_vector &derivativeAtQPs) const; /// Evaluates the jacobian of a DOF vector at a set of quadrature points /// defined on the given element. void getD2AtQPs(const ElInfo *elInfo, const Quadrature *quad, const FastQuadrature *quadFast, mtl::dense_vector::type> &d2AtQPs) const; /// Returns the FE space the DOF vector is defined on. inline const FiniteElemSpace* getFeSpace() const { return feSpace; } /// Assembles the element vector for the given ellement and adds the /// element matrix to the current DOF vector. void assemble(T factor, ElInfo *elInfo, const BoundaryType *bound, Operator *op = NULL); void assemble2(T factor, ElInfo *mainElInfo, ElInfo *auxElInfo, ElInfo *smallElInfo, ElInfo *largeElInfo, const BoundaryType *bound, Operator *op = NULL); void addElementVector(T sign, const ElementVector& elVec, const BoundaryType *bound, ElInfo *elInfo, bool add = true); /// That function must be called after the matrix assembling has been /// finished. This makes it possible to start some cleanup or matrix /// data compressing procedures. void finishAssembling(); inline void addOperator(Operator* op, double *factor = NULL, double *estFactor = NULL) { operators.push_back(op); operatorFactor.push_back(factor); operatorEstFactor.push_back(estFactor); } inline vector::iterator getOperatorFactorBegin() { return operatorFactor.begin(); } inline vector::iterator getOperatorFactorEnd() { return operatorFactor.end(); } inline vector::iterator getOperatorEstFactorBegin() { return operatorEstFactor.begin(); } inline vector::iterator getOperatorEstFactorEnd() { return operatorEstFactor.end(); } inline vector::iterator getOperatorsBegin() { return operators.begin(); } inline vector::iterator getOperatorsEnd() { return operators.end(); } Flag getAssembleFlag(); /// Evaluates \f[ u_h(x(\lambda)) = \sum_{i=0}^{m-1} vec[ind[i]] * /// \varphi^i(\lambda) \f] where \f$\varphi^i \f$ is the i-th basis /// function, \f$x(\lambda) \f$ are the world coordinates of lambda /// and \f$m \f$ is the number of basis functions T evalUh(const DimVec& lambda, DegreeOfFreedom* ind); inline vector& getOperators() { return operators; } inline vector& getOperatorFactor() { return operatorFactor; } inline vector& getOperatorEstFactor() { return operatorEstFactor; } /// Returns \ref name inline string getName() const { return name; } inline void setName(string n) { name = n; } inline BoundaryManager* getBoundaryManager() const { return boundaryManager; } inline void setBoundaryManager(BoundaryManager *bm) { boundaryManager = bm; } #ifdef HAVE_PARALLEL_DOMAIN_AMDIS void setDofMap(FeSpaceDofMap& m) { dofMap = &m; } inline bool isRankDof(DegreeOfFreedom dof) { TEST_EXIT_DBG(dofMap)("No rank dofs set!\n"); return dofMap->isRankDof(dof); } inline void setDirichletDofValue(DegreeOfFreedom dof, T value) { dirichletDofValues[dof] = value; } map& getDirichletValues() { return dirichletDofValues; } #endif protected: /// const FiniteElemSpace *feSpace; /// string name; /// ElementVector elementVector; /// vector operators; /// vector operatorFactor; /// vector operatorEstFactor; /// BoundaryManager *boundaryManager; /// Number of basis functions of the used finite element space. int nBasFcts; /// Dimension of the mesh this DOFVectorBase belongs to int dim; #ifdef HAVE_PARALLEL_DOMAIN_AMDIS FeSpaceDofMap *dofMap; map dirichletDofValues; #endif }; /// Specifies which operation should be done after coarsening typedef enum{ NO_OPERATION = 0, COARSE_RESTRICT = 1, COARSE_INTERPOL = 2 } CoarsenOperation; /** \ingroup DOFAdministration * \brief * The DOFs described above are just integers that can be used as indices into * vectors and matrices. During refinement and coarsening of the mesh, the * number of used DOFs, the meaning of one integer index, and even the total * range of DOFs change. To be able to handle these changes automatically for * all vectors, which are indexed by the DOFs, special data structures are * used which contain such vector data. Lists of these structures are kept in * DOFAdmin, so that all vectors in the lists can be resized together with the * range of DOFs. During refinement and coarsening of elements, values can be * interpolated automatically to new DOFs, and restricted from old DOFs. */ template class DOFVector : public DOFVectorBase, public Serializable { public: /** \ingroup DOFAdministration * \brief * Enables the access of DOFVector::Iterator. Alias for DOFIterator */ class Iterator : public DOFIterator { public: Iterator(DOFIndexed *c, DOFIteratorType type) : DOFIterator(c, type) {} Iterator(DOFAdmin *admin, DOFIndexed *c, DOFIteratorType type) : DOFIterator(admin, c, type) {} }; class Creator : public CreatorInterface > { public: Creator(FiniteElemSpace *feSpace_) : feSpace(feSpace_) {} DOFVector *create() { return new DOFVector(feSpace, ""); } void free(DOFVector *vec) { delete vec; } private: FiniteElemSpace *feSpace; }; public: /// Empty constructor. No initialization! DOFVector() : DOFVectorBase(), coarsenOperation(NO_OPERATION) {} /// Constructs a DOFVector with name n belonging to FiniteElemSpace f DOFVector(const FiniteElemSpace* f, string n); /// Initialization. void init(const FiniteElemSpace* f, string n); /// Copy Constructor DOFVector(const DOFVector& rhs) { *this = rhs; this->name = rhs.name + "copy"; if (this->feSpace && this->feSpace->getAdmin()) (dynamic_cast(this->feSpace->getAdmin()))->addDOFIndexed(this); } /// Destructor virtual ~DOFVector(); /// Returns iterator to the begin of \ref vec typename vector::iterator begin() { return vec.begin(); } /// Returns iterator to the end of \ref vec typename vector::iterator end() { return vec.end(); } /// Returns const_iterator to the begin of \ref vec typename vector::const_iterator begin() const { return vec.begin(); } /// Returns const_iterator to the end of \ref vec typename vector::const_iterator end() const { return vec.end(); } /// Used by DOFAdmin to compress this DOFVector. Implementation of /// \ref DOFIndexedBase::compress() virtual void compressDOFIndexed(int first, int last, vector &newDof); /// Sets \ref coarsenOperation to op inline void setCoarsenOperation(CoarsenOperation op) { coarsenOperation = op; } /// Returns \ref coarsenOperation inline CoarsenOperation getCoarsenOperation() { return coarsenOperation; } /// Restriction after coarsening. Implemented for DOFVector inline void coarseRestrict(RCNeighbourList&, int) {} /// Interpolation after refinement. Implemented for DOFVector inline void refineInterpol(RCNeighbourList&, int) {} /// Returns \ref vec vector& getVector() { return vec; } /// Returns size of \ref vec inline int getSize() const { return vec.size(); } /// Returns used size of the vector. inline int getUsedSize() const { return this->feSpace->getAdmin()->getUsedSize(); } /// Resizes \ref vec to n inline void resize(int n) { FUNCNAME("DOFVector::resize()"); TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); vec.resize(n); } /// Resizes \ref vec to n and inits new values with init inline void resize(int n, T init) { FUNCNAME("DOFVector::resize()"); TEST_EXIT_DBG(n >= 0)("Can't resize DOFVector to negative size\n"); vec.resize(n, init); } /// Returns \ref vec[i] inline const T& operator[](DegreeOfFreedom i) const { FUNCNAME("DOFVector::operator[]"); TEST_EXIT_DBG(i >= 0 && i < static_cast(vec.size())) ("Illegal vector index %d.\n", i); return vec[i]; } /// Returns \ref vec[i] inline T& operator[](DegreeOfFreedom i) { FUNCNAME("DOFVector::operator[]"); TEST_EXIT_DBG(i >= 0 && i < static_cast(vec.size())) ("Illegal vector index %d.\n", i); return vec[i]; } /// Calculates Integral of this DOFVector double Int(Quadrature* q = NULL) const { return Int(-1, q); } /** \brief * Calculates Integral of this DOFVector. * * \param[in] meshlevel Then mesh level on which the integral should be * calculated. Usually, only -1 for the leaf level * makes sence, but other values can be used, e.g., * for debugging. * \param[in] q Quadrature object. If not specified, the function * creates a new quadrature object. */ double Int(int meshLevel, Quadrature* q = NULL) const; /// Calculates Integral of this DOFVector over parts of the domain /// boundary, indicated by boundaryType. Implemented for DOFVector T IntOnBoundary(BoundaryType boundary, Quadrature* q = NULL) const { FUNCNAME("DOFVector::IntOnBoundary())"); TEST_EXIT(false)("Please implement your integration\n"); return 0.0; } /// Calculates Integral of this DOFVector times normal vector over parts /// of the domain boundary, indicated by boundaryType. Implemented for /// DOFVector > double IntOnBoundaryNormal(BoundaryType boundary, Quadrature* q = NULL) const { FUNCNAME("DOFVector::IntOnBoundaryNormal())"); TEST_EXIT(false)("Please implement your integration\n"); return 0.0; } /// Calculates L1 norm of this DOFVector double L1Norm(Quadrature* q = NULL) const; /// Calculates L2 norm of this DOFVector inline double L2Norm(Quadrature* q = NULL) const { return sqrt(L2NormSquare()); } /// Calculates square of L2 norm of this DOFVector double L2NormSquare(Quadrature* q = NULL) const; /// Calculates H1 norm of this DOFVector inline double H1Norm(Quadrature* q = NULL) const { return sqrt(H1NormSquare()); } /// Calculates square of H1 norm of this DOFVector double H1NormSquare(Quadrature* q = NULL) const; /// Calculates euclidian norm of this DOFVector double nrm2() const; /// Returns square of the euclidian norm. double squareNrm2() const; /// Calculates l2 norm of this DOFVector inline double l2norm() const { return nrm2(); } /// Calculates the absolute sum of this DOFVector T asum() const; /// Calculates the l1 norm of this DOFVector inline double l1norm() const { return asum(); } /// Calculates doublewell of this DOFVector double DoubleWell(Quadrature* q = NULL) const; /// Calculates the sum of this DOFVector T sum() const; /// Sets \ref vec[i] = val, i=0 , ... , size void set(T val); /// Assignment operator for setting complete vector to a certain value d inline DOFVector& operator=(T d) { set(d); return *this; } /// Assignment operator between two vectors DOFVector& operator=(const DOFVector&); /// vec[i] = v.vec[i] void copy(const DOFVector& v); /// Returns minimum of DOFVector. T min() const; /// Returns maximum of DOFVector. T max() const; /// Returns absolute maximum of DOFVector. T absMax() const; /// Returns the average value of the DOFVector. T average() const; /// Prints \ref vec to stdout. void print() const; /// int calcMemoryUsage() const; /// Computes the coefficients of the interpolant of the function fct and /// stores these in the DOFVector void interpol(AbstractFunction > *fct); void interpol(DOFVector *v, double factor = 1.0); /// Eval DOFVector at given point p. If oldElInfo != NULL the search for /// the element, where p is inside, starts from oldElInfo. implemented for: /// double, WorldVector< double > T evalAtPoint(WorldVector &p, ElInfo *oldElInfo = NULL) const { FUNCNAME("DOFVector::evalAtPoint())"); TEST_EXIT(false)("Please implement your evaluation\n"); } /// Determine the DegreeOfFreedom that has coords with minimal euclidean /// distance to WorldVector p. return true if DOF is found, and false /// otherwise. const bool getDofIdxAtPoint(WorldVector &p, DegreeOfFreedom &idx, ElInfo *oldElInfo = NULL, bool useOldElInfo = false) const; /// Writes the data of the DOFVector to an output stream. void serialize(ostream &out) { unsigned int size = vec.size(); out.write(reinterpret_cast(&size), sizeof(unsigned int)); out.write(reinterpret_cast(&(vec[0])), size * sizeof(T)); } /// Reads data of an DOFVector from an input stream. void deserialize(istream &in) { unsigned int size; in.read(reinterpret_cast(&size), sizeof(unsigned int)); vec.resize(size); in.read(reinterpret_cast(&(vec[0])), size * sizeof(T)); } DOFVector::type>* getGradient(DOFVector::type> *grad) const; WorldVector*> *getGradient(WorldVector*> *grad) const; DOFVector::type>* getRecoveryGradient(DOFVector::type> *grad) const; protected: /// Data container vector vec; /// Specifies what operation should be performed after coarsening CoarsenOperation coarsenOperation; }; template<> double DOFVector::IntOnBoundary( BoundaryType boundaryType, Quadrature* q) const; template<> double DOFVector >::IntOnBoundaryNormal( BoundaryType boundaryType, Quadrature* q) const; template<> double DOFVector::evalAtPoint(WorldVector &p, ElInfo *oldElInfo) const; template<> WorldVector DOFVector >::evalAtPoint(WorldVector &p, ElInfo *oldElInfo) const; template<> void DOFVector::refineInterpol(RCNeighbourList&, int); template<> void DOFVector::coarseRestrict(RCNeighbourList&, int); inline double min(const DOFVector& v) { return v.min(); } inline double max(const DOFVector& v) { return v.max(); } /** \ingroup DOFAdministration * \brief * A DOFVector that stores DOF indices. */ class DOFVectorDOF : public DOFVector, public DOFContainer { public: /// Calls constructor of DOFVector and registers itself /// as DOFContainer at DOFAdmin DOFVectorDOF(const FiniteElemSpace* feSpace_, string name_) : DOFVector(feSpace_, name_) { feSpace->getAdmin()->addDOFContainer(this); } /// Deregisters itself at DOFAdmin. ~DOFVectorDOF() { feSpace->getAdmin()->removeDOFContainer(this); } /// Implements DOFContainer::operator[]() by calling /// DOFVector::operator[]() DegreeOfFreedom& operator[](DegreeOfFreedom i) { return DOFVector::operator[](i); } const DegreeOfFreedom& operator[](DegreeOfFreedom i) const { return DOFVector::operator[](i); } /// Implements DOFIndexedBase::getSize() int getSize() const { return DOFVector::getSize(); } /// Implements DOFIndexedBase::resize() void resize(int size) { DOFVector::resize(size); } void freeDOFContent(DegreeOfFreedom dof); protected: DOFVectorDOF(); }; template double norm(DOFVector *vec) { return vec->nrm2(); } template double L2Norm(DOFVector *vec) { return vec->L2Norm(); } template double H1Norm(DOFVector *vec) { return vec->H1Norm(); } template void print(DOFVector *vec) { vec->print(); } // point wise multiplication template const DOFVector& operator*=(DOFVector& x, const DOFVector& y); // multiplication with scalar template const DOFVector& operator*=(DOFVector& x, T scal); // scalar product template T operator*(DOFVector& x, DOFVector& y); // addition template const DOFVector& operator+=(DOFVector& x, const DOFVector& y); // subtraction template const DOFVector& operator-=(DOFVector& x, const DOFVector& y); template const DOFVector& operator*(const DOFVector& v, double d); template const DOFVector& operator*(double d, const DOFVector& v); template const DOFVector& operator+(const DOFVector&v1 , const DOFVector& v2); // y = a*x + y template void axpy(double a, const DOFVector& x, DOFVector& y); // matrix vector product template void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector &x, DOFVector &result, bool add = false); template void xpay(double a,const DOFVector& x,DOFVector& y); template inline void scal(T a, DOFVector& y) { y *= a; } template inline const DOFVector& mult(double scal, const DOFVector& v, DOFVector& result); template inline const DOFVector& add(const DOFVector& v, double scal, DOFVector& result); template inline const DOFVector& add(const DOFVector& v1, const DOFVector& v2, DOFVector& result); template inline const DOFVector& mod(const DOFVector& v, DOFVector& result); template inline const DOFVector& Tanh(const DOFVector& v, DOFVector& result); template inline void set(DOFVector& vec, T d) { vec.set(d); } template inline void setValue(DOFVector& vec, T d) { vec.set(d); } template inline int size(DOFVector *vec) { return vec->getUsedSize(); } template inline int size(const DOFVector& vec) { return vec.getUsedSize(); } template inline void checkFeSpace(const FiniteElemSpace* feSpace, const vector& vec) { TEST_EXIT_DBG(feSpace)("feSpace is NULL\n"); TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n"); TEST_EXIT_DBG(static_cast(vec.size()) >= feSpace->getAdmin()->getUsedSize()) ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), feSpace->getAdmin()->getUsedSize()); } WorldVector*> *transform(DOFVector > *vec, WorldVector*> *result); template vector*> *transform(DOFVector::type> *vec, vector*> *res); /// Computes the integral: \f$\int f(\vec{x})\,\text{d}\vec{x}\f$ template TOut integrate_Coords(const FiniteElemSpace* feSpace, AbstractFunction > *fct); template TOut integrate(const FiniteElemSpace* feSpace, AbstractFunction > *fct) { return integrate_Coords(feSpace, fct); } /// Computes the integral: \f$\int f(v(\vec{x}))\,\text{d}\vec{x}\f$ template TOut integrate_Vec(const DOFVector &vec, AbstractFunction *fct); template TOut integrate(const DOFVector &vec, AbstractFunction *fct = NULL) { return fct == NULL ? vec.Int() : integrate_Vec(vec, fct); } /** \brief * Computes the integral of a function that includes two different DOFVectors: * \f$\int f(v(\vec{x}), w(\vec{x}))\,\text{d}\vec{x}\f$ * This function works also for the case that the DOFVectors are defined on * two different meshes. */ template TOut integrate_Vec2(const DOFVector &vec1, const DOFVector &vec2, BinaryAbstractFunction *fct); template TOut integrate(const DOFVector &vec1, const DOFVector &vec2, BinaryAbstractFunction *fct) { return integrate_Vec2(vec1, vec2, fct); } /// Computes the integral of a function with two DOFVectors defined on the /// same mesh. template TOut int_Vec2_SingleMesh(const DOFVector &vec1, const DOFVector &vec2, BinaryAbstractFunction *fct); /// Computes the integral of a function with two DOFVectors defined on /// different meshes. template TOut int_Vec2_DualMesh(const DOFVector &vec1, const DOFVector &vec2, BinaryAbstractFunction *fct); /// Computes the integral: \f$\int v(\vec{x}) f(\vec{x})\,\text{d}\vec{x}\f$ template typename ProductType::type integrate_VecTimesCoords(const DOFVector &vec, AbstractFunction > *fct); /// Computes the integral: \f$\int f(v(\vec{x}), \vec{x})\,\text{d}\vec{x}\f$ template TOut integrate_VecAndCoords(const DOFVector &vec, BinaryAbstractFunction > *fct); /// Computes the integral: \f$\int f(\{v_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$ double integrateGeneral(const vector*> &vecs, AbstractFunction > *fct); /// Computes the integral: /// \f$\int f(\{v_i(\vec{x})\}_i,\{\nabla w_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$ double integrateGeneralGradient(const vector*> &vecs, const vector*> &grds, BinaryAbstractFunction, vector > > *fct); } #include "DOFVector.hh" #endif // !_DOFVECTOR_H_