DOFVector.h

// ============================================================================
// ==                                                                        ==
// == 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 <vector> 
#include <memory> 
#include <list> 

#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<typename T> 
  class DOFVectorBase : public DOFIndexed<T>
  {
  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<T>& 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<T>& vecAtQPs) const;

    void getVecAtQPs(const ElInfo *smallElInfo, 
		     const ElInfo *largeElInfo,
		     const Quadrature *quad,
		     const FastQuadrature *quadFast,
		     mtl::dense_vector<T>& vecAtQPs) const;

    const WorldVector<T> *getGrdAtQPs(const ElInfo *elInfo, 
				      const Quadrature *quad,
				      const FastQuadrature *quadFast,
				      WorldVector<T> *grdAtQPs) const;

    const WorldVector<T> *getGrdAtQPs(const ElInfo *smallElInfo, 
				      const ElInfo *largeElInfo,
				      const Quadrature *quad,
				      const FastQuadrature *quadFast,
				      WorldVector<T> *grdAtQPs) const;

    const WorldMatrix<T> *getD2AtQPs(const ElInfo *elInfo, 
				     const Quadrature *quad,
				     const FastQuadrature *quadFast,
				     WorldMatrix<T> *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<double*>::iterator getOperatorFactorBegin() 
    {
      return operatorFactor.begin();
    }

    inline std::vector<double*>::iterator getOperatorFactorEnd() 
    {
      return operatorFactor.end();
    }

    inline std::vector<double*>::iterator getOperatorEstFactorBegin() 
    {
      return operatorEstFactor.begin();
    }

    inline std::vector<double*>::iterator getOperatorEstFactorEnd() 
    {
      return operatorEstFactor.end();
    }


    inline std::vector<Operator*>::iterator getOperatorsBegin() 
    {
      return operators.begin();
    }

    inline std::vector<Operator*>::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<double>& lambda, DegreeOfFreedom* ind);

    inline std::vector<Operator*>& getOperators() 
    { 
      return operators; 
    }

    inline std::vector<double*>& getOperatorFactor() 
    { 
      return operatorFactor; 
    }

    inline std::vector<double*>& 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<DegreeOfFreedom, bool> &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<Operator*> operators;

    ///
    std::vector<double*> operatorFactor;

    ///
    std::vector<double*> 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<DegreeOfFreedom, bool> *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<typename T> 
  class DOFVector : public DOFVectorBase<T>, public Serializable
  {  
  public:
    /** \ingroup DOFAdministration
     * \brief
     * Enables the access of DOFVector<T>::Iterator. Alias for DOFIterator<T>
     */
    class Iterator : public DOFIterator<T> {
    public:
      Iterator(DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(c, type)
      {}

      Iterator(DOFAdmin *admin, DOFIndexed<T> *c, DOFIteratorType type)
	: DOFIterator<T>(admin, c, type)
      {}
    };

    class Creator : public CreatorInterface<DOFVector<T> > {
    public:
      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}

      DOFVector<T> *create() 
      { 
	return new DOFVector<T>(feSpace, ""); 
      }

      void free(DOFVector<T> *vec) 
      { 
	delete vec; 
      }

    private:
      FiniteElemSpace *feSpace;
    };

  public:
    /// Empty constructor. No initialization!
    DOFVector() 
      : DOFVectorBase<T>(), 
	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<DOFAdmin*>(this->feSpace->getAdmin()))->addDOFIndexed(this);
    }

    /// Destructor
    virtual ~DOFVector();

    /// Returns iterator to the begin of \ref vec
    typename std::vector<T>::iterator begin() 
    { 
      return vec.begin(); 
    }

    /// Returns iterator to the end of \ref vec
    typename std::vector<T>::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<DegreeOfFreedom> &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<double>
    inline void coarseRestrict(RCNeighbourList&, int) {}

    /// Interpolation after refinement. Implemented for DOFVector<double>
    inline void refineInterpol(RCNeighbourList&, int) {}

    /// Returns \ref vec
    std::vector<T>& 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<T>::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<T>::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<T>::operator[]");
      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size()))
	("Illegal vector index %d.\n", i);
      return vec[i];
    } 

    /// Returns \ref vec[i]
    inline T& operator[](DegreeOfFreedom i) 
    {
      FUNCNAME("DOFVector<T>::operator[]");

      TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(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<T>& operator=(T d) 
    {
      set(d); 
      return *this;
    } 

    /// Assignment operator between two vectors
    DOFVector<T>& operator=(const DOFVector<T>&); 

    /// vec[i] = v.vec[i]
    void copy(const DOFVector<T>& 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<T, WorldVector<double> > *fct);

    void interpol(DOFVector<T> *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<double> &p, ElInfo *oldElInfo = NULL, T* value = NULL) const 
    {
      FUNCNAME("DOFVector::evalAtPoint())");
      TEST_EXIT(false)("Please implement your evaluation\n");
      return *value;
    }

    /// Writes the data of the DOFVector to an output stream.
    void serialize(std::ostream &out) 
    {
      unsigned int size = vec.size();
      out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
      out.write(reinterpret_cast<const char*>(&(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<char*>(&size), sizeof(unsigned int));
      vec.resize(size);
      in.read(reinterpret_cast<char*>(&(vec[0])), size * sizeof(T));
    }

    DOFVector<WorldVector<T> > *getGradient(DOFVector<WorldVector<T> >*) const;

    WorldVector<DOFVector<T>*> *getGradient(WorldVector<DOFVector<T>*> *grad) const;

    DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;

  protected: 
    /// Data container
    std::vector<T> vec; 
 
    /// Specifies what operation should be performed after coarsening
    CoarsenOperation coarsenOperation;

    /// Used by \ref interpol
    AbstractFunction<T, WorldVector<double> > *interFct;

    /// Used for mesh traversal
    static DOFVector<T> *traverseVector;
  }; 


  template<>
  const double& DOFVector<double>::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, double* value) const;

  template<>
  const WorldVector<double>& DOFVector<WorldVector<double> >::evalAtPoint(WorldVector<double> &p, ElInfo *oldElInfo, WorldVector<double>* value) const;

  template<>
  void DOFVector<double>::refineInterpol(RCNeighbourList&, int);

  template<>
  void DOFVector<double>::coarseRestrict(RCNeighbourList&, int); 

  inline double min(const DOFVector<double>& v) 
  {
    return v.min();
  } 

  inline double max(const DOFVector<double>& v) 
  {
    return v.max();
  }


  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
    /** \brief
     * Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
     * as DOFContainer at DOFAdmin
     */
    DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
      : DOFVector<DegreeOfFreedom>(feSpace_, name_)
    {
      feSpace->getAdmin()->addDOFContainer(this);
    }
  
    /// Deregisters itself at DOFAdmin.
    ~DOFVectorDOF() 
    {
      feSpace->getAdmin()->removeDOFContainer(this);
    }

    /** \brief
     * Implements DOFContainer::operator[]() by calling 
     * DOFVector<DegreeOfFreedom>::operator[]()
     */
    DegreeOfFreedom& operator[](DegreeOfFreedom i) 
    {
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

    const DegreeOfFreedom& operator[](DegreeOfFreedom i) const 
    {
      return DOFVector<DegreeOfFreedom>::operator[](i);
    }

    /// Implements DOFIndexedBase::getSize()
    int getSize() const 
    {
      return DOFVector<DegreeOfFreedom>::getSize();
    }

    /// Implements DOFIndexedBase::resize()
    void resize(int size) 
    {
      DOFVector<DegreeOfFreedom>::resize(size);
    }

    void freeDOFContent(DegreeOfFreedom dof);

  protected:
    DOFVectorDOF();
  };


  template<typename T>
  double norm(DOFVector<T> *vec) 
  {
    return vec->nrm2();
  }

  template<typename T>
  double L2Norm(DOFVector<T> *vec) 
  {
    return vec->L2Norm();
  }

  template<typename T>
  double H1Norm(DOFVector<T> *vec) 
  {
    return vec->H1Norm();
  }

  template<typename T>
  void print(DOFVector<T> *vec) 
  {
    vec->print();
  }

  // point wise multiplication
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y);

  // multiplication with scalar
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal);

  // scalar product
  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y);

  // addition
  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y);

  // subtraction
  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y);

  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d);

  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v);

  template<typename T>
  const DOFVector<T>& operator+(const DOFVector<T>&v1 , const DOFVector<T>& v2);

  // y = a*x + y
  template<typename T>
  void axpy(double a,const DOFVector<T>& x, DOFVector<T>& y);

  // matrix vector product
  template<typename T>
  void mv(MatrixTranspose transpose, 
	  const DOFMatrix &a, 
	  const DOFVector<T> &x,
	  DOFVector<T> &result,
	  bool add = false); 

  template<typename T>
  void xpay(double a,const DOFVector<T>& x,DOFVector<T>& y);

  template<typename T>
  inline void scal(T a, DOFVector<T>& y) 
  {
    y *= a;
  }

  template<typename T>
  inline const DOFVector<T>& mult(double scal,
				  const DOFVector<T>& v,
				  DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v1,
				 const DOFVector<T>& v2,
				 DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& mod(const DOFVector<T>& v,
				 DOFVector<T>& result);

  template<typename T>
  inline const DOFVector<T>& Tanh(const DOFVector<T>& v,
				  DOFVector<T>& result);


  template<typename T>
  inline void set(DOFVector<T>& vec, T d) 
  {
    vec.set(d);
  }

  template<typename T>
  inline void setValue(DOFVector<T>& vec, T d) 
  {
    vec.set(d);
  }

  template<typename T>
  inline int size(DOFVector<T> *vec) 
  {
    return vec->getUsedSize();
  }

  template<typename T>
  inline void checkFeSpace(const FiniteElemSpace* feSpace, const std::vector<T>& vec)
  {
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       feSpace->getAdmin()->getUsedSize());
  }

  WorldVector<DOFVector<double>*> *transform(DOFVector<WorldVector<double> > *vec,
					     WorldVector<DOFVector<double>*> *result);

  /** \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<double> &vec1,
		   const DOFVector<double> &vec2,
		   BinaryAbstractFunction<double, double, double> *fct);

  /// Computes the integral of a function with two DOFVectors defined on the same mesh.
  double intSingleMesh(const DOFVector<double> &vec1,
		       const DOFVector<double> &vec2,
		       BinaryAbstractFunction<double, double, double> *fct);

  /// Computes the integral of a function with two DOFVectors defined on different meshes.
  double intDualMesh(const DOFVector<double> &vec1,
		     const DOFVector<double> &vec2,
		     BinaryAbstractFunction<double, double, double> *fct);

  double integrate(const DOFVector<double> &vec,
		   AbstractFunction<double, WorldVector<double> > *fct);

  double integrate(const FiniteElemSpace* feSpace,
		   AbstractFunction<double, WorldVector<double> > *fct);


}

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_