DOFVector.h

// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file DOFVector.h */

#ifndef AMDIS_DOFVECTOR_H 
#define AMDIS_DOFVECTOR_H 
 
// ===========================================================================
// ===== includes ============================================================
// ===========================================================================

#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 "MemoryManager.h"
#include "Boundary.h"
#include "CreatorInterface.h"
#include "Serializable.h"
#include "DOFMatrix.h" 
#include "BasisFunction.h"

#ifdef HAVE_PARALLEL_AMDIS
#include "petscao.h"
#endif
 
namespace AMDiS {

  // ===========================================================================
  // ===== forward declarations ================================================
  // ===========================================================================

  template<typename T> 
  class DOFVectorBase : public DOFIndexed<T>
  {
  public:

    DOFVectorBase() 
      : feSpace(NULL),
	elementVector(NULL),
        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<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,
			  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 const std::string& getName() const { 
      return name; 
    } 

    inline BoundaryManager* getBoundaryManager() const { 
      return boundaryManager; 
    }

    inline void setBoundaryManager(BoundaryManager *bm) {
      boundaryManager = bm;
    }

  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;

    /// Are used to store temporary local dofs of an element. Thread safe.
    std::vector<DegreeOfFreedom*> localIndices;

    /// Are used to store temporary local values of an element. Thread safe.
    std::vector<T*> localVectors;

    /// Temporarly used in \ref getGrdAtQPs. Thread safe.
    std::vector<DimVec<double>*> grdTmp;

    /// Temporarly used in \ref getGrdAtQPs. Thread safe.
    std::vector<DimVec<double>*> grdPhis;

    /// Temporarly used in \ref getD2AtQPs. Thread safe.
    std::vector<DimMat<double>*> D2Phis;

    /// Dimension of the mesh this DOFVectorBase belongs to
    int dim;
  };


  // ===========================================================================
  // ===== defs ================================================================
  // ===========================================================================

  /** \brief
   * Specifies which operation should be done after coarsening
   */
  typedef enum{
    NO_OPERATION = 0,   
    COARSE_RESTRICT = 1,
    COARSE_INTERPOL = 2 
  } CoarsenOperation;
 
  // ===========================================================================
  // ===== class DOFVector =====================================================
  // ===========================================================================

  /** \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:
    MEMORY_MANAGED(DOFVector<T>);

    /** \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:
      MEMORY_MANAGED(Creator);

      Creator(FiniteElemSpace *feSpace_) 
        : feSpace(feSpace_) 
      {}

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

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

    private:
      FiniteElemSpace *feSpace;
    };

  public:
    /** \brief
     * Empty constructor. No initialization!
     */
    DOFVector() 
      : DOFVectorBase<T>(), 
	//elementVector(NULL), 
	feSpace(NULL),
	coarsenOperation(NO_OPERATION)
    {}

    /** \brief
     * Constructs a DOFVector with name n belonging to FiniteElemSpace f
     */
    DOFVector(const FiniteElemSpace* f, std::string n); 

    /** \brief
     * Initialization.
     */
    void init(const FiniteElemSpace* f, std::string n);

    /** \brief
     * Copy Constructor
     */
    DOFVector(const DOFVector& rhs) {
      *this = rhs;   
      name = rhs.name + "copy";
      if (feSpace && feSpace->getAdmin()) {
	(dynamic_cast<DOFAdmin*>(feSpace->getAdmin()))->addDOFIndexed(this);
      }
    }

    /** \brief
     * Destructor
     */
    virtual ~DOFVector();

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

    /** \brief
     * 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);
    /** \brief
     * Sets \ref coarsenOperation to op
     */
    inline void setCoarsenOperation(CoarsenOperation op) { 
      coarsenOperation = op; 
    }

    /** \brief
     * Returns \ref coarsenOperation
     */
    inline CoarsenOperation getCoarsenOperation() { 
      return coarsenOperation; 
    }

    /** \brief
     * Restriction after coarsening. Implemented for DOFVector<double>
     */
    inline void coarseRestrict(RCNeighbourList&, int) {}

    /** \brief
     * Interpolation after refinement.
     */
    inline void refineInterpol(RCNeighbourList&, int) {}

    /** \brief
     * Returns \ref vec
     */
    std::vector<T>& getVector() { 
      return vec;
    }

    /** \brief
     * Returns size of \ref vec
     */
    inline int getSize() const { 
      return vec.size();
    } 

    /** \brief
     * Returns used size of the vector.
     */
    inline int getUsedSize() const { 
      return feSpace->getAdmin()->getUsedSize(); 
    }

    /** \brief
     * 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);
    } 

    /** \brief
     * 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);
    } 

    /** \brief
     * 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];
    } 

    /** \brief
     * 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];
    }
 
    /** \brief
     * Calculates Integral of this DOFVector
     */
    double Int(Quadrature* q = NULL) const;

    /** \brief
     * Calculates L1 norm of this DOFVector
     */
    double L1Norm(Quadrature* q = NULL) const;
 
    /** \brief
     * Calculates L2 norm of this DOFVector
     */
    inline double L2Norm(Quadrature* q = NULL) const {
      return sqrt(L2NormSquare());
    }

    /** \brief
     * Calculates square of L2 norm of this DOFVector
     */
    double L2NormSquare(Quadrature* q = NULL) const;

    /** \brief
     * Calculates H1 norm of this DOFVector
     */
    inline double H1Norm(Quadrature* q = NULL) const {
      return sqrt(H1NormSquare());
    }

    /** \brief
     * Calculates square of H1 norm of this DOFVector
     */
    double H1NormSquare(Quadrature* q = NULL) const;  

    /** \brief
     * Calculates euclidian norm of this DOFVector
     */
    double nrm2() const; 

    /** \brief
     * Returns square of the euclidian norm.
     */
    double squareNrm2() const;

    /** \brief
     * Calculates l2 norm of this DOFVector
     */
    inline double l2norm() const { 
      return nrm2();
    }

    /** \brief
     * Calculates the absolute sum of this DOFVector
     */
    T asum() const; 

    /** \brief
     * Calculates the l1 norm of this DOFVector
     */
    inline double l1norm() const { 
      return asum();
    } 

    /** \brief
     * Calculates doublewell of this DOFVector
     */
    double DoubleWell(Quadrature* q = NULL) const;
 
    /** \brief
     * Calculates the sum of this DOFVector
     */
    T sum() const; 
 
    /** \brief
     * Sets \ref vec[i] = val, i=0 , ... , size
     */
    void set(T val); 

    /** \brief
     * Assignment operator for setting complete vector to a certain value d
     */
    inline DOFVector<T>& operator=(T d) {
      set(d); 
      return *this;
    } 

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

    /** \brief
     * vec[i] = v.vec[i]
     */
    void copy(const DOFVector<T>& v); 

    /** \brief
     * Returns minimum of DOFVector
     */
    T min() const; 

    /** \brief
     * 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<T, WorldVector<double> > *fct);

    void interpol(DOFVector<T> *v, double factor);

    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));
    }

    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;

#ifdef HAVE_PARALLEL_AMDIS
    /// Sets the petsc application ordering object to map dof indices.
    void useApplicationOrdering(AO *ao) {
      applicationOrdering = ao;
    }
#endif

  protected: 
    /** \brief
     * Used by Int while mesh traversal
     */
    static int Int_fct(ElInfo* elinfo);

    /** \brief
     * Used by L1Norm while mesh traversal
     */
    static int L1Norm_fct(ElInfo* elinfo);
 
    /** \brief
     * Used by L2Norm while mesh traversal
     */
    static int L2Norm_fct(ElInfo* elinfo);

    /** \brief
     * Used by H1Norm while mesh traversal
     */
    static int H1Norm_fct(ElInfo* elinfo);

    /** \brief
     * Used by DoubleWell while mesh traversal
     */
    static int DoubleWell_fct(ElInfo* elinfo);

  protected: 
    /// Name of this DOFVector
    std::string name; 

    /// FiniteElemSpace of the vector
    const FiniteElemSpace *feSpace; 

    /// 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;

#ifdef HAVE_PARALLEL_AMDIS
    /// Petsc application ordering to map dof indices.
    AO *applicationOrdering;
#endif

  protected:
    /// Used while calculating vector norms
    static FastQuadrature *quad_fast;

    /// Stores the last calculated vector norm
    static double norm;
  }; 

  // ===============================================================================

  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();
  }

  // ===========================================================================
  // ===== class DOFVectorDOF ==================================================
  // ===========================================================================

  /** \ingroup DOFAdministration
   * \brief
   * A DOFVector that stores DOF indices.
   */
  class DOFVectorDOF : public DOFVector<DegreeOfFreedom>,
		       public DOFContainer
  {
  public:  
    MEMORY_MANAGED(DOFVectorDOF);

    /** \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();
  }

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

#include "DOFVector.hh"

#endif  // !_DOFVECTOR_H_