Operator.h

// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  TU Dresden                                                            ==
// ==                                                                        ==
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
// ==                                                                        ==
// ============================================================================

/** \file Operator.h */

#ifndef AMDIS_OPERATOR_H
#define AMDIS_OPERATOR_H

#include <vector>
#include "AMDiS_fwd.h"
#include "FixVec.h"
#include "Flag.h"
#include "MatrixVector.h"
#include "ElInfo.h"
#include "AbstractFunction.h"
#include "OpenMP.h"
#include "SubAssembler.h"
#include "OperatorTerm.h"
#include "ZeroOrderTerm.h"
#include "FirstOrderTerm.h"
#include "SecondOrderTerm.h"

namespace AMDiS {

  /** \brief
   * An Operator holds all information needed to assemble the system matrix
   * and the right hand side. It consists of four OperatorTerm lists each storing
   * Terms of a specific order and type. You can define your own Operator by 
   * creating an empty Operator and than adding OperatorTerms to it.
   * An Operator can by of type MATRIX_OPERATOR, if it's used to assemble the
   * system matrix on the left hand side, or it can be of type VECTOR_OPERATOR,
   * if it's used to assemble the right hand side vector. If an Operator gives
   * contributions to both sides of the system it is a MATRIX_OPERATOR and a
   * VECTOR_OPERATOR in one instance. This allows to efficiently reuse element 
   * matrices once calculated.
   * By calling \ref getElementMatrix() or \ref getElementVector() one can 
   * initiate the assembling procedure. Therefor each Operator has its own
   * Assembler, especially created for this Operator, by the first call of
   * \ref getElementMatrix() or \ref getElementVector(). 
   */
  class Operator
  {
  public:
    /** \brief
     * Constructs an empty Operator of type operatorType for the given 
     * FiniteElemSpace. 
     * The type is given by a Flag that can contain the values MATRIX_OPERATOR,
     * VECTOR_OPERATOR, or MATRIX_OPERATOR | VECTOR_OPERATOR. This type specifies 
     * whether the Operator is used on the left hand side, the right hand side,
     * or on both sides of the system. 
     */
    Operator(Flag operatorType,
	     const FiniteElemSpace *rowFESpace,
	     const FiniteElemSpace *colFESpace = NULL);

    /// Destructor.
    virtual ~Operator() {}

    /// Sets \ref optimized.
    inline void useOptimizedAssembler(bool opt) 
    {
      optimized = opt;
    }

    /// Returns \ref optimized.
    inline bool isOptimized() 
    {
      return optimized;
    }

    /// Adds a SecondOrderTerm to the Operator
    template <typename T>
    void addSecondOrderTerm(T *term);

    /// Adds a FirstOrderTerm to the Operator
    template <typename T>
    void addFirstOrderTerm(T *term, FirstOrderType type = GRD_PHI);

    /// Adds a ZeroOrderTerm to the Operator
    template <typename T>
    void addZeroOrderTerm(T *term);


    /** \brief
     * Calculates the element matrix for this ElInfo and adds it multiplied by
     * factor to userMat.
     */
    virtual void getElementMatrix(const ElInfo *elInfo, 
				  ElementMatrix& userMat, 
				  double factor = 1.0);

    virtual void getElementMatrix(const ElInfo *rowElInfo,
				  const ElInfo *colElInfo,
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
				  ElementMatrix& userMat,
				  double factor = 1.0);

    /** \brief
     * Calculates the element vector for this ElInfo and adds it multiplied by
     * factor to userVec.
     */
    virtual void getElementVector(const ElInfo *elInfo, 
				  ElementVector& userVec, 
				  double factor = 1.0);

    virtual void getElementVector(const ElInfo *mainElInfo, 
				  const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
				  ElementVector& userVec,
				  double factor = 1.0);

    /// That function must be called after one assembling cycle has been finished.
    void finishAssembling();

    /// Returns \ref rowFESpace.
    inline const FiniteElemSpace *getRowFESpace() 
    { 
      return rowFESpace; 
    }

    /// Returns \ref colFESpace.
    inline const FiniteElemSpace *getColFESpace() 
    { 
      return colFESpace; 
    }

    /// Returns \ref auxFeSpaces.
    inline std::set<const FiniteElemSpace*>& getAuxFeSpaces()
    {
      return auxFeSpaces;
    }

    /// Sets \ref uhOld.
    void setUhOld(const DOFVectorBase<double> *uhOld);

    /// Returns \ref uhOld.
    inline const DOFVectorBase<double> *getUhOld() 
    {
      return uhOld;
    }    

    /// Returns \ref assembler
    Assembler *getAssembler(int rank = 0);

    /// Sets \ref assembler
    void setAssembler(int rank, Assembler *ass);

    /// Returns whether this is a matrix operator.
    inline const bool isMatrixOperator() 
    {
      return type.isSet(MATRIX_OPERATOR);
    }

    /// Returns whether this is a vector operator
    inline const bool isVectorOperator() 
    {
      return type.isSet(VECTOR_OPERATOR);
    }

    /// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
    inline void setFillFlag(Flag f) 
    { 
      fillFlag = f; 
    }

    /// Sets \ref needDualTraverse.
    void setNeedDualTraverse(bool b) 
    {
      needDualTraverse = b;
    }

    /// Returns \ref fillFlag
    inline Flag getFillFlag() 
    { 
      return fillFlag; 
    }

    /// Returns \ref needDualTraverse
    bool getNeedDualTraverse() 
    {
      return needDualTraverse;
    }

    /// Initialization of the needed SubAssemblers using the given quadratures. 
    void initAssembler(int rank,
		       Quadrature *quad2,
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


    /// Calculates the needed quadrature degree for the given order. 
    int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);

    /// Evaluation of all terms in \ref zeroOrder. 
    void evalZeroOrder(int nPoints,
		       const double *uhAtQP,
		       const WorldVector<double> *grdUhAtQP,
		       const WorldMatrix<double> *D2UhAtQP,
		       double *result,
		       double factor) const
    {
      int myRank = omp_get_thread_num();

      std::vector<OperatorTerm*>::const_iterator termIt;
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
	   ++termIt)
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
    }

    /// Evaluation of all terms in \ref firstOrderGrdPsi. 
    void evalFirstOrderGrdPsi(int nPoints,
			      const double *uhAtQP,
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
      int myRank = omp_get_thread_num();

      std::vector<OperatorTerm*>::const_iterator termIt;
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
	   ++termIt)
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
    }

    /// Evaluation of all terms in \ref firstOrderGrdPhi. 
    void evalFirstOrderGrdPhi(int nPoints,
			      const double *uhAtQP,
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
      int myRank = omp_get_thread_num();

      std::vector<OperatorTerm*>::const_iterator termIt;
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
	   ++termIt)
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
    }

    /// Evaluation of all terms in \ref secondOrder. 
    void evalSecondOrder(int nPoints,
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
    {
      int myRank = omp_get_thread_num();

      std::vector<OperatorTerm*>::const_iterator termIt;
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
	   ++termIt)
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
    }

    /// Weak evaluation of all terms in \ref secondOrder. 
    void weakEvalSecondOrder(const std::vector<WorldVector<double> > &grdUhAtQP,
			     std::vector<WorldVector<double> > &result) const
    {
      int myRank = omp_get_thread_num();

      for (std::vector<OperatorTerm*>::const_iterator termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
	   ++termIt)
	static_cast<SecondOrderTerm*>(*termIt)->weakEval(grdUhAtQP, result);
    }
  
    /// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
    {
      int myRank = omp_get_thread_num();

      std::vector<OperatorTerm*>::const_iterator termIt;
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
	   ++termIt)
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
    }
  
    /// Calls getLb() for each term in \ref firstOrderGrdPsi and adds the results to Lb.
    void getLbGrdPsi(const ElInfo *elInfo, 
		     int nPoints, 
		     VectorOfFixVecs<DimVec<double> >& Lb) const
    {
      int myRank = omp_get_thread_num();

      std::vector<OperatorTerm*>::const_iterator termIt;
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
	   ++termIt)
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
    }

    /// Calls getLb() for each term in \ref firstOrderGrdPhi and adds the results to Lb.
    void getLbGrdPhi(const ElInfo *elInfo, 
		     int nPoints, 
		     VectorOfFixVecs<DimVec<double> >& Lb) const
    {
      int myRank = omp_get_thread_num();

      std::vector<OperatorTerm*>::const_iterator termIt;
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
	   ++termIt)
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
    }

    /// Calls getC() for each term in \ref zeroOrder and adds the results to c.
    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c)
    {
      int myRank = omp_get_thread_num();

      std::vector<OperatorTerm*>::const_iterator termIt;
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
	   ++termIt)
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
    }

    /// Returns true, if there are second order terms. Returns false otherwise.
    inline bool secondOrderTerms() 
    {
      return secondOrder[omp_get_thread_num()].size() != 0;
    }

    /// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
    inline bool firstOrderTermsGrdPsi() 
    {
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
    }

    /// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
    inline bool firstOrderTermsGrdPhi() 
    {
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
    }

    /// Returns true, if there are zero order terms. Returns false otherwise.
    inline bool zeroOrderTerms() 
    {
      return zeroOrder[omp_get_thread_num()].size() != 0;
    }

  public:
    /// Constant type flag for matrix operators
    static const Flag MATRIX_OPERATOR;

    /// Constant type flag for vector operators
    static const Flag VECTOR_OPERATOR;

  protected:
    /// FiniteElemSpace for matrix rows and element vector
    const FiniteElemSpace *rowFESpace;

    /// FiniteElemSpace for matrix columns. Can be equal to \rowFESpace.
    const FiniteElemSpace *colFESpace;

    /// List of aux fe spaces, e.g., if a term is multiplied with a DOF vector
    std::set<const FiniteElemSpace*> auxFeSpaces;

    /// Number of rows in the element matrix
    int nRow;

    /// Number of columns in the element matrix
    int nCol;

    /// Type of this Operator.
    Flag type;

    /// Flag for mesh traversal
    Flag fillFlag;

    /// If true, the operator needs a dual traverse over two meshes for assembling.
    bool needDualTraverse;

    /** \brief
     * Calculates the element matrix and/or the element vector. It is
     * created especially for this Operator, when \ref getElementMatrix()
     * or \ref getElementVector is called for the first time.
     */
    std::vector<Assembler*> assembler;

    /// List of all second order terms
    std::vector< std::vector<OperatorTerm*> > secondOrder;

    /// List of all first order terms derived to psi
    std::vector< std::vector<OperatorTerm*> > firstOrderGrdPsi;

    /// List of all first order terms derived to phi
    std::vector< std::vector<OperatorTerm*> > firstOrderGrdPhi;

    /// List of all zero order terms
    std::vector< std::vector<OperatorTerm*> > zeroOrder;

    /** \brief
     * Pointer to the solution of the last timestep. Can be used if the 
     * Operator is MATRIX_OPERATOR and VECTOR_OPERATOR for a more efficient
     * assemblage of the element vector when the element matrix was already
     * computed.
     */
    const DOFVectorBase<double> *uhOld;

    /// Spezifies whether optimized assemblers are used or not.
    bool optimized;

    friend class Assembler;
    friend class SubAssembler;
    friend class ZeroOrderAssembler;
    friend class FirstOrderAssembler;
    friend class SecondOrderAssembler;
  };

}

#include "Operator.hh"

#endif // AMDIS_OPERATOR_H