// ============================================================================ // == == // == 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 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" namespace AMDiS { template class DOFVectorBase : public DOFIndexed { public: DOFVectorBase() : feSpace(NULL), elementVector(3), boundaryManager(NULL), nBasFcts(0) {} DOFVectorBase(const FiniteElemSpace *f, std::string n); virtual ~DOFVectorBase(); /** \brief * For the given element, this function returns an array of all DOFs of this * DOFVector that are defined on this element. */ virtual const T *getLocalVector(const Element *el, T* localVec) const; /** \brief * Evaluates the DOF vector at a set of quadrature points defined on the * given element. */ const T *getVecAtQPs(const ElInfo *elInfo, const Quadrature *quad, const FastQuadrature *quadFast, T *vecAtQPs) const; const T *getVecAtQPs(const ElInfo *smallElInfo, const ElInfo *largeElInfo, const Quadrature *quad, const FastQuadrature *quadFast, T *vecAtQPs) const; const WorldVector *getGrdAtQPs(const ElInfo *elInfo, const Quadrature *quad, const FastQuadrature *quadFast, WorldVector *grdAtQPs) const; const WorldVector *getGrdAtQPs(const ElInfo *smallElInfo, const ElInfo *largeElInfo, const Quadrature *quad, const FastQuadrature *quadFast, WorldVector *grdAtQPs) const; const WorldMatrix *getD2AtQPs(const ElInfo *elInfo, const Quadrature *quad, const FastQuadrature *quadFast, WorldMatrix *d2AtQPs) const; /// Returns the FE space the DOF vector is defined on. inline const FiniteElemSpace* getFESpace() const { return feSpace; } /** \brief * 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); /* \brief * 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 std::vector::iterator getOperatorFactorBegin() { return operatorFactor.begin(); } inline std::vector::iterator getOperatorFactorEnd() { return operatorFactor.end(); } inline std::vector::iterator getOperatorEstFactorBegin() { return operatorEstFactor.begin(); } inline std::vector::iterator getOperatorEstFactorEnd() { return operatorEstFactor.end(); } inline std::vector::iterator getOperatorsBegin() { return operators.begin(); } inline std::vector::iterator getOperatorsEnd() { return operators.end(); } Flag getAssembleFlag(); /** \brief * 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 std::vector& getOperators() { return operators; } inline std::vector& getOperatorFactor() { return operatorFactor; } inline std::vector& getOperatorEstFactor() { return operatorEstFactor; } /// Returns \ref name inline std::string getName() const { return name; } inline BoundaryManager* getBoundaryManager() const { return boundaryManager; } inline void setBoundaryManager(BoundaryManager *bm) { boundaryManager = bm; } #ifdef HAVE_PARALLEL_DOMAIN_AMDIS inline void setRankDofs(std::map dofmap) { rankDofs.clear(); // rankDofs = dofmap; } inline bool isRankDof(DegreeOfFreedom dof) { return rankDofs[dof]; } #endif protected: /// const FiniteElemSpace *feSpace; /// std::string name; /// ElementVector elementVector; /// std::vector operators; /// std::vector operatorFactor; /// std::vector operatorEstFactor; /// BoundaryManager *boundaryManager; /// Number of basis functions of the used finite element space. int nBasFcts; /// Are used to store temporary local dofs of an element. Thread safe. std::vector localIndices; /// Are used to store temporary local values of an element. Thread safe. std::vector localVectors; /// Temporarly used in \ref getGrdAtQPs. Thread safe. std::vector*> grdTmp; /// Temporarly used in \ref getGrdAtQPs. Thread safe. std::vector*> grdPhis; /// Temporarly used in \ref getD2AtQPs. Thread safe. std::vector*> D2Phis; /// Dimension of the mesh this DOFVectorBase belongs to int dim; #ifdef HAVE_PARALLEL_DOMAIN_AMDIS public: std::map rankDofs; #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, std::string n); /// Initialization. void init(const FiniteElemSpace* f, std::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 std::vector::iterator begin() { return vec.begin(); } /// Returns iterator to the end of \ref vec typename std::vector::iterator end() { return vec.end(); } /** \brief * Used by DOFAdmin to compress this DOFVector. Implementation of * DOFIndexedBase::compress() */ virtual void compressDOFIndexed(int first, int last, std::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. inline void refineInterpol(RCNeighbourList&, int) {} /// Returns \ref vec std::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; /// 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; /// Used by interpol while mesh traversal static int interpolFct(ElInfo* elinfo); /// Prints \ref vec to stdout void print() const; /// int calcMemoryUsage() const; /** \brief * 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); /// Writes the data of the DOFVector to an output stream. void serialize(std::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(std::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 > *getGradient(DOFVector >*) const; WorldVector*> *getGradient(WorldVector*> *grad) const; DOFVector >* getRecoveryGradient(DOFVector >*) const; protected: /// Used by Int while mesh traversal static int Int_fct(ElInfo* elinfo); /// Used by L1Norm while mesh traversal static int L1Norm_fct(ElInfo* elinfo); /// Used by L2Norm while mesh traversal static int L2Norm_fct(ElInfo* elinfo); /// Used by H1Norm while mesh traversal static int H1Norm_fct(ElInfo* elinfo); /// Used by DoubleWell while mesh traversal static int DoubleWell_fct(ElInfo* elinfo); protected: /// Data container std::vector vec; /// Specifies what operation should be performed after coarsening CoarsenOperation coarsenOperation; /// Used by \ref interpol AbstractFunction > *interFct; /// Used for mesh traversal static DOFVector *traverseVector; protected: /// Used while calculating vector norms static FastQuadrature *quad_fast; /// Stores the last calculated vector norm static double norm; }; 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: /** \brief * Calls constructor of DOFVector and registers itself * as DOFContainer at DOFAdmin */ DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_) : DOFVector(feSpace_, name_) { feSpace->getAdmin()->addDOFContainer(this); } /// Deregisters itself at DOFAdmin. ~DOFVectorDOF() { feSpace->getAdmin()->removeDOFContainer(this); } /** \brief * 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 void checkFeSpace(const FiniteElemSpace* feSpace, const std::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); } #include "DOFVector.hh" #endif // !_DOFVECTOR_H_