// ============================================================================ // == == // == 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" namespace AMDiS { template class DOFVectorBase : public DOFIndexed { public: DOFVectorBase() : feSpace(NULL), elementVector(3), boundaryManager(NULL), nBasFcts(0) { #ifdef HAVE_PARALLEL_DOMAIN_AMDIS rankDofs = NULL; #endif } 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 void getLocalVector(const Element *el, mtl::dense_vector& localVec) const; /** \brief * 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; /** \brief * 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; /** \brief * 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; /** \brief * 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; } /** \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 void setName(std::string n) { name = n; } 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 = &dofmap; } inline bool isRankDof(DegreeOfFreedom dof) { TEST_EXIT_DBG(rankDofs)("No rank dofs set!\n"); 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; /// Dimension of the mesh this DOFVectorBase belongs to int dim; #ifdef HAVE_PARALLEL_DOMAIN_AMDIS 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. Implemented for DOFVector 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; /** \brief * 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; } /** \brief * 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; /// Used by interpol while mesh traversal. void 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 = 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 > */ inline const T& evalAtPoint( WorldVector &p, ElInfo *oldElInfo = NULL, T* value = NULL) const { FUNCNAME("DOFVector::evalAtPoint())"); TEST_EXIT(false)("Please implement your evaluation\n"); return *value; } /** 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(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; DOFVector::type>* getGradient(DOFVector::type> *grad) const; WorldVector*> *getGradient(WorldVector*> *grad) const; // DOFVector >* getRecoveryGradient(DOFVector >*) const; DOFVector::type>* getRecoveryGradient(DOFVector::type> *grad) const; 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; }; template<> double DOFVector::IntOnBoundary( BoundaryType boundaryType, Quadrature* q) const; template<> double DOFVector >::IntOnBoundaryNormal( BoundaryType boundaryType, Quadrature* q) const; template<> const double& DOFVector::evalAtPoint( WorldVector &p, ElInfo *oldElInfo, double* value) const; template<> const WorldVector& DOFVector >::evalAtPoint( WorldVector &p, ElInfo *oldElInfo, WorldVector* value) 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: /** \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); template std::vector*> *transform(DOFVector::type> *vec, std::vector*> *res); /** \brief * Computes the integral of a function that includes two different DOFVectors. This * function works also for the case that the DOFVectors are defined on two different * meshes. */ double integrate(const DOFVector &vec1, const DOFVector &vec2, BinaryAbstractFunction *fct); /// Computes the integral of a function with two DOFVectors defined on the same mesh. double intSingleMesh(const DOFVector &vec1, const DOFVector &vec2, BinaryAbstractFunction *fct); /// Computes the integral of a function with two DOFVectors defined on different meshes. double intDualMesh(const DOFVector &vec1, const DOFVector &vec2, BinaryAbstractFunction *fct); double integrate(const DOFVector &vec, AbstractFunction > *fct); double integrate(const FiniteElemSpace* feSpace, AbstractFunction > *fct); double integrateGeneral(const std::vector*> &vecs, AbstractFunction > *fct); double integrateGeneralGradient(const std::vector*> &vecs, const std::vector*> &grds, BinaryAbstractFunction, std::vector > > *fct); template T integrate_VecAndCoords(const DOFVector &vec, BinaryAbstractFunction > *fct); template T2 integrate_Vec(const DOFVector &vec, AbstractFunction *fct); } #include "DOFVector.hh" #endif // !_DOFVECTOR_H_