// ============================================================================
// ==                                                                        ==
// == 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 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"

namespace AMDiS {

  /** 
   * \ingroup Assembler
   * 
   * \brief
   * Base class for ZeroOrderTerm, FirstOrderTerm and SecondOrderTerm. 
   * OperatorTerms are the building blocks of an Operator. Each OperatorTerm
   * has its properties which are regarded, when constructing 
   * an Assembler for the corresponding Operator.
   */
  class OperatorTerm
  {
  public:
    /** \brief
     * Constructs an OperatorTerm with initially no properties.
     * degree_ is used to determine the degree of the needed quadrature
     * for the assemblage.  
     */
    OperatorTerm(int deg) 
      : properties(0), 
	degree(deg),
	dimOfWorld(Global::getGeo(WORLD)),
	bOne(-1)
    {}

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

    /** \brief
     * Virtual method. It's called by SubAssembler::initElement() for
     * each OperatorTerm belonging to this SubAssembler. E.g., vectors
     * and coordinates at quadrature points can be calculated here.
     */
    virtual void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL) 
    {}

    virtual void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
			     SubAssembler*, Quadrature *quad = NULL)
    {}


    /// Returs \auxFeSpaces, the list of all aux fe spaces the operator makes use off.
    inline std::set<const FiniteElemSpace*>& getAuxFeSpaces() 
    {
      return auxFeSpaces;
    }

    /// Specifies whether the matrix of the term is symmetric
    void setSymmetric(bool symm);

    /// Returns true, if the term is piecewise constant, returns false otherwise
    inline bool isPWConst() 
    { 
      return (degree == 0);
    }

    /// Returns true, if the term has a symmetric matrix, returns false otherwise.
    bool isSymmetric();

    /// Returns \ref degree.
    inline int getDegree() 
    { 
      return degree; 
    }

    /// Sets one component of the b vector to be one. See \ref bOne.
    void setB(int b)
    {
      bOne = b;
    }

    /// Evaluation of the OperatorTerm at all quadrature points.
    virtual void eval(int nPoints,
		      const double *uhAtQP,
		      const WorldVector<double> *grdUhAtQP,
		      const WorldMatrix<double> *D2UhAtQP,
		      double *result,
		      double factor) = 0;

    /** \brief
     * Determines the value of a dof vector at the quadrature points of a given 
     * element. It is used by all VecAtQP like operator terms.
     */
    double *getVectorAtQPs(DOFVectorBase<double>* vec,
			   const ElInfo* elInfo, 
			   SubAssembler* subAssembler,
			   Quadrature *quad);		   

    /** \brief
     * Determines the value of a dof vector at the quadrature points of a given 
     * element. This function is used, if an operator is assembled on two different
     * meshes using the dual traverse. The element, i.e. the small or the large one,
     * is choosen, which corresponds to the mesh the dof vector is defined on.
     */
    double *getVectorAtQPs(DOFVectorBase<double>* vec,
			   const ElInfo* smallElInfo, 
			   const ElInfo* largeElInfo, 
			   SubAssembler* subAssembler,
			   Quadrature *quad);	

    ///
    WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec,
					   const ElInfo* elInfo,
					   SubAssembler* subAssembler,
					   Quadrature *quad);

    ///
    WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec,
					   const ElInfo* smallElInfo, 
					   const ElInfo* largeElInfo,
					   SubAssembler* subAssembler,
					   Quadrature *quad);

  protected:
    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$.
    void lalt(const DimVec<WorldVector<double> >& Lambda,
	      const WorldMatrix<double>& matrix,
	      DimMat<double>& LALt,
	      bool symm,
	      double factor) const;

    /** \brief
     * Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for \f$ A \f$
     * the matrix having a ONE in the position \f$ (K,L) \f$
     * and ZEROS in all other positions.
     */
    static void lalt_kl(const DimVec<WorldVector<double> >& Lambda,
			int k, int l,
			DimMat<double>& LALt,
			double factor);

    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity.
    inline void l1lt(const DimVec<WorldVector<double> >& Lambda, 
		     DimMat<double>& LALt,
		     double factor) const
    {
      const int dim = LALt.getNumRows();
      
      for (int i = 0; i < dim; i++) {
	double val = 0.0;
	for (int k = 0; k < dimOfWorld; k++)
	  val += Lambda[i][k] * Lambda[i][k];
	val *= factor;
	LALt[i][i] += val;
	for (int j = i + 1; j < dim; j++) {
	  val = 0.0;
	  for (int k = 0; k < dimOfWorld; k++)
	    val += Lambda[i][k] * Lambda[j][k];
	  val *= factor;
	  LALt[i][j] += val;
	  LALt[j][i] += val;
	}    
      }      
    }

    /// Evaluation of \f$ \Lambda \cdot b\f$.
    inline void lb(const DimVec<WorldVector<double> >& Lambda,
		   const WorldVector<double>& b,
		   DimVec<double>& Lb,
		   double factor) const
    {
      const int dim = Lambda.getSize();

      for (int i = 0; i < dim; i++) {
	double val = 0.0;
	
	for (int j = 0; j < dimOfWorld; j++)
	  val += Lambda[i][j] * b[j];
	
	Lb[i] += val * factor;
      }    
    }

    /// Evaluation of \f$ \Lambda \cdot b\f$.
    inline void lb_one(const DimVec<WorldVector<double> >& Lambda,
		       DimVec<double>& Lb,
		       double factor) const
    {
      const int dim = Lambda.getSize();

      for (int i = 0; i < dim; i++)
	Lb[i] += Lambda[i][bOne] * factor;
    }

    /** \brief
     * Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in 
     * each component.
     */
    inline void l1(const DimVec<WorldVector<double> >& Lambda,
		   DimVec<double>& Lb,
		   double factor) const
    {
      const int dim = Lambda.getSize();

      for (int i = 0; i < dim; i++) {
	double val = 0.0;
      
	for (int j = 0; j < dimOfWorld; j++)
	  val += Lambda[i][j];

	Lb[i] += val * factor;
      }    
    }

  protected:
    /// Stores the properties of this OperatorTerm
    Flag properties;

    /// Polynomial degree of the term. Used to detemine the degree of the quadrature.
    int degree;

    /// Stores the dimension of the world.
    int dimOfWorld;

    /// List off all fe spaces, the operator term makes use off.
    std::set<const FiniteElemSpace*> auxFeSpaces;

    /// Pointer to the Operator this OperatorTerm belongs to.
    Operator* operat;

    /** \brief
     * In many cases, the vector b in the evaluation \f$ \Lambda \cdot b\f$ has zeros
     * in all components expect one that is set to one. Using the function \ref lb is
     * then unnecessary time consuming. Instead, this variable defines the component
     * of the vector b to be one. The function \ref lb_one is used if this variable is
     * not -1.
     */
    int bOne;

    /// Flag for piecewise constant terms
    static const Flag PW_CONST;

    /// Flag for symmetric terms
    static const Flag SYMMETRIC;

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

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Describes the second order terms: \f$ \nabla \cdot A \nabla u(\vec{x}) \f$
   */
  class SecondOrderTerm : public OperatorTerm
  {
  public:
    /// Constructor.
    SecondOrderTerm(int deg) 
      : OperatorTerm(deg) 
    {}

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

    /// Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points.
    virtual void getLALt(const ElInfo *elInfo, 
			 int nPoints, 
			 DimMat<double> **result) const = 0;

    /// Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
    virtual void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
			  std::vector<WorldVector<double> > &result) = 0;

  };

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Implements the laplace operator: \f$ \Delta u(\vec{x}) \f$
   */
  class Laplace_SOT : public SecondOrderTerm 
  {
  public:
    /// Constructor.
    Laplace_SOT() 
      : SecondOrderTerm(0) 
    {
      setSymmetric(true);
    }

    /// Implenetation of SecondOrderTerm::getLALt().
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
    {
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();

      for (int iq = 0; iq < nPoints; iq++)
	l1lt(Lambda, *(LALt[iq]), 1.0);
    }

    /// Implementation of SecondOrderTerm::eval().
    inline void eval(int nPoints,
		     const double *,
		     const WorldVector<double> *,
		     const WorldMatrix<double> *D2UhAtQP,
		     double *result,
		     double factor) 
    {
      int dow = Global::getGeo(WORLD);
    
      if (D2UhAtQP) {
	for (int iq = 0; iq < nPoints; iq++) {
	  double resultQP = 0.0;
	  for (int i = 0; i < dow; i++) 
	    resultQP += D2UhAtQP[iq][i][i];
	  result[iq] += factor * resultQP;
	}
      }
    }

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result) 
    {
      int nPoints = grdUhAtQP.size();
      for (int iq = 0; iq < nPoints; iq++)
	result[iq] += grdUhAtQP[iq];
    }
  };

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Implements the laplace operator multiplied with a scalar factor:
   * \f$ f \cdot \Delta u(\vec{x}) \f$
   */
  class FactorLaplace_SOT : public SecondOrderTerm 
  {
  public:
    /// Constructor.
    FactorLaplace_SOT(double f) 
      : SecondOrderTerm(0)  
    {
      factor = new double;
      *factor = f;

      setSymmetric(true);
    }

    /// Constructor.
    FactorLaplace_SOT(double *fptr) 
      : SecondOrderTerm(0), 
	factor(fptr)
    {
      setSymmetric(true);
    }

    /// Implements SecondOrderTerm::getLALt().
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
    {
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
      for (int iq = 0; iq < nPoints; iq++) 
	l1lt(Lambda, *(LALt[iq]), (*factor));
    }

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *,
	      const WorldVector<double> *,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double f) 
    {
      int dow = Global::getGeo(WORLD);

      if (D2UhAtQP) {
	for (int iq = 0; iq < nPoints; iq++) {
	  double resultQP = 0.0;
	  for (int i = 0; i < dow; i++) {
	    resultQP += D2UhAtQP[iq][i][i];
	  }
	  result[iq] += resultQP * f * (*factor);
	}
      }
    }

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result) 
    {
      int nPoints = grdUhAtQP.size();
      for (int iq = 0; iq < nPoints; iq++)
	axpy(*factor, grdUhAtQP[iq], result[iq]);
    }

  private:
    /// Pointer to the factor.
    double *factor;
  };

  /**
   * \ingroup Assembler
   * 
   * \brief
   * SecondOrderTerm where A is a function which maps a DOFVector evaluated at
   * a given quadrature point to a WolrdMatrix:
   * \f$ \nabla \cdot A(v(\vec{x})) \nabla u(\vec{x}) \f$
   */
  class MatrixFct_SOT : public SecondOrderTerm 
  {
  public:
    /// Constructor.
    MatrixFct_SOT(DOFVectorBase<double> *dv, 
		  AbstractFunction<WorldMatrix<double>, double> *fct,
		  AbstractFunction<WorldVector<double>, WorldMatrix<double> > *div,
		  bool sym = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
  
    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Pointer to the values of the DOFVector at quadrature points.
    double* vecAtQPs;

    /// Function for A.
    AbstractFunction<WorldMatrix<double>, double>* matrixFct;

    ///
    AbstractFunction<WorldVector<double>, WorldMatrix<double> >* divFct;

    /// True, if \ref matrixFct produces always symmetric matrices.
    bool symmetric;
  };

  /** 
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A is a given fixed WorldMatrix<double>:
   * \f$ \nabla \cdot A \nabla u(\vec{x}) \f$
   */
  class Matrix_SOT : public SecondOrderTerm {
  public:
    /// Constructor
    Matrix_SOT(WorldMatrix<double> mat) 
      : SecondOrderTerm(0), matrix(mat)
    {
      symmetric = matrix.isSymmetric();
      setSymmetric(symmetric);
    }

    /// Implements SecondOrderTerm::getLALt().
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
    {
      const DimVec<WorldVector<double> >& Lambda = elInfo->getGrdLambda();
      for (int iq = 0; iq < nPoints; iq++) 
	lalt(Lambda, matrix, *(LALt[iq]), symmetric, 1.0);
    }
  
    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// Matrix stroring A.
    WorldMatrix<double> matrix;

    /// True, if \ref matrix is symmetric.
    bool symmetric;
  };

  /** 
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A is a WorldMatrix<double> having a ONE in position IJ
   * and ZERO in all other positions
   * \f$ \nabla \cdot A \nabla u(\vec{x}) \f$
   */
  class FactorIJ_SOT : public SecondOrderTerm 
  {
  public:
    /// Constructor.
    FactorIJ_SOT(int x_i, int x_j, double f) 
      : SecondOrderTerm(0), xi(x_i), xj(x_j)
    {
      factor = new double;
      *factor = f;

      setSymmetric(xi == xj);
    }

    /// Constructor.
    FactorIJ_SOT(int x_i, int x_j, double *fptr) 
      : SecondOrderTerm(0), xi(x_i), xj(x_j), factor(fptr)
    {
      setSymmetric(xi == xj);
    }

    /// Implements SecondOrderTerm::getLALt().
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
    {
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
      for (int iq = 0; iq < nPoints; iq++)
	lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*factor));
    }

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *,
	      const WorldVector<double> *,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac)
    {
      if (D2UhAtQP)
	for (int iq = 0; iq < nPoints; iq++)
	  result[iq] += (*factor) * D2UhAtQP[iq][xi][xj] * fac;
    }

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result) 
    {
      int nPoints = grdUhAtQP.size();
      for (int iq = 0; iq < nPoints; iq++)
	result[iq][xi] += (*factor) * grdUhAtQP[iq][xj];
    }

  private:
    /// Directions for the derivatives.
    int xi, xj;

    /// Pointer to the factor.
    double *factor;
  };

  /**
   * \ingroup Assembler
   *  
   * \brief
   * Laplace operator multiplied with a function evaluated at the quadrature
   * points of a given DOFVector:
   * \f$ f(v(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class VecAtQP_SOT : public SecondOrderTerm {
  public:
    /// Constructor.
    VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo,
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
    double* vecAtQPs;

    /// Function evaluated at \ref vecAtQPs.
    AbstractFunction<double, double> *f;
  };

  /**
   * \ingroup Assembler
   *  
   * \brief
   * Laplace operator multiplied with a function evaluated at the quadrature
   * points of a given DOFVector:
   * \f$ f(v(\vec{x}), w(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class Vec2AtQP_SOT : public SecondOrderTerm {
  public:
    /// Constructor.
    Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, 
		 BinaryAbstractFunction<double, double, double> *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    
    
    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);
    
    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
    
  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;
    
    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
    double* vecAtQPs1;
    double* vecAtQPs2;
    
    /// Function evaluated at \ref vecAtQPs.
    BinaryAbstractFunction<double, double, double> *f;
  };

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Laplace multiplied with a function evaluated at each quadrature point:
   * \f$ f(\vec{x}) \Delta u(\vec{x}) \f$
   */
  class CoordsAtQP_SOT : public SecondOrderTerm {
  public:
    /// Constructor.
    CoordsAtQP_SOT(AbstractFunction<double, WorldVector<double> > *af)
      : SecondOrderTerm(af->getDegree()), g(af)
    {
      setSymmetric(true);
    }

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);


  protected:
    /// Stores coordinates at quadrature points. Set in \ref initElement().
    WorldVector<double>* coordsAtQPs;

    /// Function evaluated at quadrature points.
    AbstractFunction<double, WorldVector<double> > *g;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A is a function which maps the gradient of a 
   * DOFVector at each quadrature point to WorldMatrix<double>:
   * \f$ \nabla \cdot A(\nabla v(\vec{x})) \nabla u(\vec{x})\f$
   */
  class MatrixGradient_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    MatrixGradient_SOT(DOFVectorBase<double> *dv,
		       AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
		       AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
		       bool symm = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo, 
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    DOFVectorBase<double>* vec;

    /// Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
    AbstractFunction<WorldMatrix<double>, WorldVector<double> > *f;

    AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divFct;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    WorldVector<double>* gradAtQPs;

    /// True, if \ref f provides always a symmetric WorldMatrix<double>.
    bool symmetric;

  private:
    /// Temporary matrix used in \ref MatrixGradient_SOT::weakEval.
    WorldMatrix<double> tmpMat;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Laplace multiplied with a function which maps the gradient of a DOFVector
   * at each quadrature point to a double:
   * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class FctGradient_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    FctGradient_SOT(DOFVectorBase<double> *dv,  
		    AbstractFunction<double, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    DOFVectorBase<double>* vec;

    /// Function wich maps \ref gradAtQPs to a double.
    AbstractFunction<double, WorldVector<double> > *f;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    WorldVector<double>* gradAtQPs;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * Laplace multiplied with a function which maps the gradient of a DOFVector
   * at each quadrature point to a double:
   * \f$ \nabla \cdot A(v(\vec{x}), \nabla v(\vec{x})) \nabla u(\vec{x}) \f$
   */
  class VecMatrixGradientAtQP_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
			      BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
			      AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
			      bool symm = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
    
  protected:
    DOFVectorBase<double>* vec;
    
    /// Function wich maps \ref gradAtQPs to a double.
    BinaryAbstractFunction<WorldMatrix<double>, 
			   double, WorldVector<double> > *f;
    
    AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divFct;
    
    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    double* vecAtQPs;
    WorldVector<double>* gradAtQPs;
    
    bool symmetric;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Laplace multiplied with a function which maps the gradient of a DOFVector
   * at each quadrature point to a double:
   * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class VecGradCoordsAtQP_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    VecGradCoordsAtQP_SOT(DOFVectorBase<double> *dv,
			  TertiaryAbstractFunction<double, double,
			  WorldVector<double>, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    DOFVectorBase<double>* vec;

    /// Function wich maps \ref gradAtQPs to a double.
    TertiaryAbstractFunction<double, double, WorldVector<double>,
			     WorldVector<double> > *f;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    double* vecAtQPs;
    WorldVector<double>* gradAtQPs;
    WorldVector<double>* coordsAtQPs;
  
  };


  /**
   * \ingroup Assembler
   *  
   * \brief
   * Laplace operator multiplied with a function evaluated at the quadrature
   * points of a given DOFVector:
   * \f$ -f(v(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class VecAndCoordsAtQP_SOT : public SecondOrderTerm {
  public:
    /// Constructor.
    VecAndCoordsAtQP_SOT(DOFVectorBase<double> *dv, 
			 BinaryAbstractFunction<double, double, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::eval().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    

    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
    double* vecAtQPs;
  
    WorldVector<double>* coordsAtQPs;

    /// Function evaluated at \ref vecAtQPs.
    BinaryAbstractFunction<double, double,  WorldVector<double> > *f;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A is a function which maps the gradient of a 
   * DOFVector at each quadrature point to WorldMatrix<double>:
   * \f$ \nabla \cdot A(\nabla v(\vec{x})) \nabla u(\vec{x})\f$
   */
  class MatrixGradientAndCoords_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    MatrixGradientAndCoords_SOT(DOFVectorBase<double> *dv,
				BinaryAbstractFunction<WorldMatrix<double>,
				WorldVector<double> , 
				WorldVector<double> > *af,
				AbstractFunction<WorldVector<double>,
				WorldMatrix<double> > *divAf,
				bool symm = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    DOFVectorBase<double>* vec;

    /// Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
    BinaryAbstractFunction<WorldMatrix<double>,
			   WorldVector<double>, WorldVector<double> > *f;

    AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divFct;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    WorldVector<double>* gradAtQPs;

    WorldVector<double>* coordsAtQPs;

    /// True, if \ref f provides always a symmetric WorldMatrix<double>.
    bool symmetric;
  };

  class MatrixVec2_SOT : public SecondOrderTerm 
  {
  public:
    /// Constructor.
    MatrixVec2_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,  
		  BinaryAbstractFunction<double, double, double> *f,
		  WorldMatrix<double> Af,
		  bool sym = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
  
    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;

    /// Pointer to the values of the DOFVector at quadrature points.
    double* vec1AtQPs;
    double* vec2AtQPs;

    /// Function for A.
    BinaryAbstractFunction<double, double, double> * fct;

    ///
    WorldMatrix<double> A;

    /// True, if \ref matrixFct produces always symmetric matrices.
    bool symmetric;
  };


  class General_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    General_SOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
		TertiaryAbstractFunction<WorldMatrix<double>, 
		  WorldVector<double>,
		  std::vector<double>, 
		  std::vector<WorldVector<double> > > *f,
		AbstractFunction<WorldVector<double>, 
		  WorldMatrix<double> > *divFct,
		bool symmetric);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo*, SubAssembler*, Quadrature *quad= NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    std::vector<DOFVectorBase<double>*> vecs_; 

    std::vector<DOFVectorBase<double>*> grads_;

    TertiaryAbstractFunction<WorldMatrix<double>, 
			     WorldVector<double>,
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f_;

    AbstractFunction<WorldVector<double>, 
		     WorldMatrix<double> > *divFct_;

    WorldVector<double> *coordsAtQPs_;

    std::vector<double*> vecsAtQPs_;

    std::vector<WorldVector<double>*> gradsAtQPs_;

    bool symmetric_;
  };

  class GeneralParametric_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    GeneralParametric_SOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
			  QuartAbstractFunction<WorldMatrix<double>, 
			  WorldVector<double>,
			  WorldVector<double>,
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f,
			  AbstractFunction<WorldVector<double>, 
			  WorldMatrix<double> > *divFct,
			  bool symmetric);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo*, 
		     SubAssembler* ,
		     Quadrature *quad= NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    std::vector<DOFVectorBase<double>*> vecs_; 

    std::vector<DOFVectorBase<double>*> grads_;

    QuartAbstractFunction<WorldMatrix<double>, 
			  WorldVector<double>,
			  WorldVector<double>,
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f_;

    AbstractFunction<WorldVector<double>, 
		     WorldMatrix<double> > *divFct_;

    WorldVector<double> *coordsAtQPs_;
    WorldVector<double> elementNormal_;

    std::vector<double*> vecsAtQPs_;

    std::vector<WorldVector<double>*> gradsAtQPs_;

    bool symmetric_;
  };


  // ============================================================================
  // ===== class FirstOrderTerm =================================================
  // ============================================================================

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Describes the first order terms: \f$ b \cdot \nabla u(\vec{x}) \f$
   */
  class FirstOrderTerm : public OperatorTerm
  {
  public:
    /// Constructor.
    FirstOrderTerm(int deg) 
      : OperatorTerm(deg) 
    {}

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

    /// Evaluation of \f$ \Lambda b \f$.
    virtual void getLb(const ElInfo *elInfo,
		       int nPoints, 
		       VectorOfFixVecs<DimVec<double> >& result) const = 0;

    /// Implenetation of FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *,
	      double *result,
	      double factor)  
    {
      int dow = Global::getGeo(WORLD);

      if (grdUhAtQP) {
	for (int iq = 0; iq < nPoints; iq++) {
	  double resultQP = 0.0;
	  for (int i = 0; i < dow; i++)
	    resultQP += grdUhAtQP[iq][i];
	  
	  result[iq] += resultQP * factor;
	}
      }
    }

  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * First order term: \f$ b \cdot \nabla u(\vec{x}) \f$ with a vector b which 
   * only consits of entries equal to one.
   */
  class Simple_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    Simple_FOT() 
      : FirstOrderTerm(0) 
    {}

    /// Implements FirstOrderTerm::getLb().
    inline void getLb(const ElInfo *elInfo, 
		      int nPoints, 
		      VectorOfFixVecs<DimVec<double> >& Lb) const 
    {
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();

      for (int iq = 0; iq < nPoints; iq++)
	l1(Lambda, Lb[iq], 1.0);
    }
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * First order term: \f$ b \cdot \nabla u(\vec{x}) \f$ with a vector b which 
   * only consits of entries equal to one.
   */
  class FactorSimple_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    FactorSimple_FOT(double f) 
      : FirstOrderTerm(0)
    {
      factor = new double;
      *factor = f;
    }

    /// Constructor.
    FactorSimple_FOT(double *fptr) 
      : FirstOrderTerm(0), factor(fptr) 
    {}

    /// Implements FirstOrderTerm::getLb().
    inline void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
    {
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();

      for (int iq = 0; iq < nPoints; iq++)
	l1(Lambda, Lb[iq], (*factor));
    }

  private:
    /// Pointer to the factor.
    double *factor;
  };

  /**
   * \ingroup Assembler
   *
   * \brief 
   * First order term: \f$ b \cdot \nabla u(\vec{x}) \f$ with a given vector b.
   */
  class Vector_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    Vector_FOT(WorldVector<double> wv) 
      : FirstOrderTerm(0), b(wv) 
    {}

    /// Constructor.
    Vector_FOT(int bIdx) 
      : FirstOrderTerm(0) 
    {
      bOne = bIdx;
    }

    /// Implements FirstOrderTerm::getLb().
    inline void getLb(const ElInfo *elInfo, 
		      int nPoints, 
		      VectorOfFixVecs<DimVec<double> >&Lb) const 
    {
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();

      if (bOne > -1) {
	for (int iq = 0; iq < nPoints; iq++)
	  lb_one(Lambda, Lb[iq], 1.0);
      } else {
	for (int iq = 0; iq < nPoints; iq++)
	  lb(Lambda, b, Lb[iq], 1.0);
      }
    }

    /// Implements FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *,
	      double *result,
	      double factor)
    {
      if (grdUhAtQP)
	for (int iq = 0; iq < nPoints; iq++)
	  result[iq] += b * grdUhAtQP[iq] * factor; 
    }

  protected:
    /// Vector which is multiplied with \f$ \nabla u(\vec{x}) \f$
    WorldVector<double> b;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Simple_FOT multiplied with \f$ f(v(\vec{x})) \f$.
   */
  class VecAtQP_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    VecAtQP_FOT(DOFVectorBase<double> *dv, 
		AbstractFunction<double, double> *af,
		WorldVector<double> *wv);

    /// Constructor.
    VecAtQP_FOT(DOFVectorBase<double> *dv,
		AbstractFunction<double, double> *af,
		int bIdx);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements FirstOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int nPoints, 
	       VectorOfFixVecs<DimVec<double> >& Lb) const;

    /// Implements FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// v at quadrature points.
    double *vecAtQPs;

    /// Function f.
    AbstractFunction<double, double> *f;

    ///
    WorldVector<double> *b;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Simple_FOT multiplied with \f$ f(\vec{x}) \f$.
   */
  class CoordsAtQP_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    CoordsAtQP_FOT(AbstractFunction<double, WorldVector<double> > *af)
      : FirstOrderTerm(af->getDegree()), g(af) 
    {}

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements FistOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int nPoints,VectorOfFixVecs<DimVec<double> >&Lb) const;

    /// Implements FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

  protected:
    /// Stores coordinates at quadrature points. Set in \ref initElement().
    WorldVector<double>* coordsAtQPs;

    /// Function avaluated at world coordinates.
    AbstractFunction<double, WorldVector<double> > *g;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Vector_FOT multiplied with \f$ f(\vec{x}) \f$.
   */
  class VecCoordsAtQP_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    VecCoordsAtQP_FOT(AbstractFunction<double, WorldVector<double> > *af,
		      WorldVector<double> wv)
      : FirstOrderTerm(af->getDegree()), g(af), b(wv)
    {}

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements FistOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int nPoints,
	       VectorOfFixVecs<DimVec<double> >& Lb) const;

    /// Implements FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

  protected:
    /// Stores coordinates at quadrature points. Set in \ref initElement().
    WorldVector<double>* coordsAtQPs;

    /// Function evaluated at world coordinates.
    AbstractFunction<double, WorldVector<double> > *g;

    /// Coefficient vector.
    WorldVector<double> b;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * First order term: \f$ b(\nabla v(\vec{x})) \cdot \nabla u(\vec{x})\f$.
   */
  class VectorGradient_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    VectorGradient_FOT(DOFVectorBase<double> *dv,
		       AbstractFunction<WorldVector<double>, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements FirstOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int nPoints, 
	       VectorOfFixVecs<DimVec<double> >& Lb) const;

    /// Implements FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

  protected:
    DOFVectorBase<double>* vec;

    /// Function for b.
    AbstractFunction<WorldVector<double>, WorldVector<double> > *f;

    /// Gradient of v at quadrature points.
    WorldVector<double> *gradAtQPs;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * First order term: \f$ b(v(\vec{x})) \cdot \nabla u(\vec{x})\f$.
   */
  class VectorFct_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    VectorFct_FOT(DOFVectorBase<double> *dv,
		  AbstractFunction<WorldVector<double>, double> *fct);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements FirstOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int qPoint,
	       VectorOfFixVecs<DimVec<double> >& Lb) const; 

    /// Implements FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Vector v at quadrature points.
    double *vecAtQPs;

    /// Function for b.
    AbstractFunction<WorldVector<double>, double> *vecFct;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * 
   */
  class VecFctAtQP_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    VecFctAtQP_FOT(AbstractFunction<WorldVector<double>, WorldVector<double> > *af)
      : FirstOrderTerm(af->getDegree()), g(af) 
    {}

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements FistOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int nPoints,
	       VectorOfFixVecs<DimVec<double> >&Lb) const;

    /// Implements FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

  protected:
    /// Stores coordinates at quadrature points. Set in \ref initElement().
    WorldVector<double>* coordsAtQPs;

    /// Function avaluated at world coordinates.
    AbstractFunction<WorldVector<double>, WorldVector<double> > *g;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * First order term: \f$ b(v(\vec{x}), \nabla w(\vec{x})) \cdot \nabla u(\vec{x})\f$.
   */
  class VecGrad_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor.
    VecGrad_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
		BinaryAbstractFunction<WorldVector<double>, double, WorldVector<double> > *fct);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo,
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
    
    /// Implements FirstOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int qPoint, 
	       VectorOfFixVecs<DimVec<double> >& Lb) const; 
    
    /// Implements FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);
    
  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;
    
    /// Vector v at quadrature points.
    double *vecAtQPs;

    WorldVector<double> *gradAtQPs;
    
    /// Function for b.
    BinaryAbstractFunction<WorldVector<double>, double, WorldVector<double> > *vecFct;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * First order term: \f$ b(\nabla v(\vec{x}), \nabla w(\vec{x})) \cdot \nabla u(\vec{x})\f$.
   */
  class FctGrad2_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor
    FctGrad2_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
		BinaryAbstractFunction<WorldVector<double>, WorldVector<double>, WorldVector<double> > *vecFct_)
      : FirstOrderTerm(vecFct_->getDegree()), vec1(dv1), vec2(dv2), vecFct(vecFct_) 
    {}

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
    
    /// Implements FirstOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int qPoint, 
	       VectorOfFixVecs<DimVec<double> >& Lb) const; 
    
    /// Implements FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);
    
  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;

    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec2;
    
    /// Gradient v at quadrature points.
    WorldVector<double> *grad1AtQPs;

    /// Gradient v at quadrature points.
    WorldVector<double> *grad2AtQPs;

    /// Function for b.
    BinaryAbstractFunction<WorldVector<double>, WorldVector<double>, WorldVector<double> > *vecFct;
  };


  class Vec2AndGradAtQP_FOT : public FirstOrderTerm
  {
  public:
    Vec2AndGradAtQP_FOT(DOFVectorBase<double> *dv1, 
			DOFVectorBase<double> *dv2,
			TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f_, WorldVector<double> *b_);

    void initElement(const ElInfo* elInfo, 
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
    
    void getLb(const ElInfo *elInfo, 
	       int nPoints, 
	       VectorOfFixVecs<DimVec<double> >& Lb) const;
    
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

  protected:
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;

    double *vec1AtQPs;
    double *vec2AtQPs;  
    WorldVector<double> *gradAtQPs1;

    TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f;
    WorldVector<double> *b;
  };


  class FctVecAtQP_FOT : public FirstOrderTerm
  {
  public:
    FctVecAtQP_FOT(DOFVectorBase<double> *dv,
		   BinaryAbstractFunction<double, WorldVector<double>, double> *f_,
		   WorldVector<double> *b_);
    
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
    
    void getLb(const ElInfo *elInfo, int nPoints, 
	       VectorOfFixVecs<DimVec<double> >& Lb) const;
    
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);
    
  protected:
    DOFVectorBase<double>* vec;
    double *vecAtQPs;
    WorldVector<double> *coordsAtQPs;
    BinaryAbstractFunction<double, WorldVector<double>, double> *f;
    WorldVector<double> *b;
  };


  class Vec2AtQP_FOT : public FirstOrderTerm
  {
  public:
    Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
		 BinaryAbstractFunction<double, double, double> *af,
		 WorldVector<double> *b_);

    Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
		 BinaryAbstractFunction<double, double, double> *af,
		 int bIdx);
    
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
    
    void getLb(const ElInfo *elInfo, int nPoints, 
	       VectorOfFixVecs<DimVec<double> >& Lb) const;
    
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);
    
  protected:
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;
    double *vec1AtQPs;
    double *vec2AtQPs;  
    /// Function f.
    BinaryAbstractFunction<double, double, double> *f;
    WorldVector<double> *b;
  };


  class Vec3FctAtQP_FOT : public FirstOrderTerm
  {
  public:
    Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
		    TertiaryAbstractFunction<double, double, double, double> *f,
		    WorldVector<double> *b);
    
    Vec3FctAtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, DOFVectorBase<double> *dv3,
		    TertiaryAbstractFunction<double, double, double, double> *f,
		    int b);
    
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
    
    void getLb(const ElInfo *elInfo, int nPoints, 
	       VectorOfFixVecs<DimVec<double> >& Lb) const;
    
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);
    
  protected:
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;
    DOFVectorBase<double>* vec3;
    TertiaryAbstractFunction<double, double, double, double>* f;
    double *vec1AtQPs;
    double *vec2AtQPs;  
    double *vec3AtQPs;  
    WorldVector<double> *b;
  };


  class General_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor
    General_FOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
		TertiaryAbstractFunction<WorldVector<double>, 
		WorldVector<double>,
		std::vector<double>, 
		std::vector<WorldVector<double> > > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo*, 
		     SubAssembler* ,
		     Quadrature *quad= NULL);

    /// Implements FirstOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int nPoints, 
	       VectorOfFixVecs<DimVec<double> >& result) const;

    /// Implenetation of FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

  protected:
    std::vector<DOFVectorBase<double>*> vecs_; 

    std::vector<DOFVectorBase<double>*> grads_;

    TertiaryAbstractFunction<WorldVector<double>, 
			     WorldVector<double>,
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f_;

    WorldVector<double> *coordsAtQPs_;

    std::vector<double*> vecsAtQPs_;

    std::vector<WorldVector<double>*> gradsAtQPs_;
  };

  class GeneralParametric_FOT : public FirstOrderTerm
  {
  public:
    /// Constructor
    GeneralParametric_FOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
			  QuartAbstractFunction<WorldVector<double>, 
			  WorldVector<double>,
			  WorldVector<double>,
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL);

    /// Implements FirstOrderTerm::getLb().
    void getLb(const ElInfo *elInfo, int nPoints, 
	       VectorOfFixVecs<DimVec<double> >& result) const;

    /// Implenetation of FirstOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

  protected:
    std::vector<DOFVectorBase<double>*> vecs_; 

    std::vector<DOFVectorBase<double>*> grads_;

    QuartAbstractFunction<WorldVector<double>, 
			  WorldVector<double>,
			  WorldVector<double>,
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f_;

    WorldVector<double> *coordsAtQPs_;
    WorldVector<double> elementNormal_;

    std::vector<double*> vecsAtQPs_;

    std::vector<WorldVector<double>*> gradsAtQPs_;
  };


  // ============================================================================
  // ===== class ZeroOrderTerm ==================================================
  // ============================================================================

  /**
   * \ingroup Assembler
   *
   * \brief
   * Describes zero order terms: \f$ cu(\vec{x}) \f$.
   */
  class ZeroOrderTerm : public OperatorTerm
  {
  public:
    /// Constructor.
    ZeroOrderTerm(int deg) : OperatorTerm(deg) {}

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

    /// Evaluates \f$ c \f$
    virtual void getC(const ElInfo *elInfo, int nPoints, 
		      std::vector<double> &C) = 0;

  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Simple zero order term
   */
  class Simple_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    Simple_ZOT(double f = 1.0) : ZeroOrderTerm(0), factor(f) {}

    /// Implements ZeroOrderTerm::getC().
    inline void getC(const ElInfo *, int nPoints, std::vector<double> &C)  
    {
      for (int iq = 0; iq < nPoints; iq++)
	C[iq] += factor; 
    }

    /// Implements ZeroOrderTerm::eval().
    inline void eval(int nPoints,
		     const double *uhAtQP,
		     const WorldVector<double> *,
		     const WorldMatrix<double> *,
		     double *result,
		     double fac)  
    {
      for (int iq = 0; iq < nPoints; iq++)
	result[iq] += fac * factor * uhAtQP[iq];
    }
 
  protected:
    /// Constant factor of zero order term.
    double factor;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Zero order term: \f$ f(v(\vec{x})) u(\vec{x})\f$
   */
  class VecAtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    VecAtQP_ZOT(DOFVectorBase<double> *dv,
		AbstractFunction<double, double> *ab);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, 
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo,
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Vector v at quadrature points.
    double *vecAtQPs;

    /// Function for c.
    AbstractFunction<double, double> *f;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Zero order term: \f$ f(v(\vec{x})) g(w(\vec{x})) u(\vec{x})\f$
   */
  class MultVecAtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    MultVecAtQP_ZOT(DOFVectorBase<double> *dv1,
		    DOFVectorBase<double> *dv2,
		    AbstractFunction<double, double> *f1,
		    AbstractFunction<double, double> *f2);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// DOFVectorBase to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;

    /// Vector v at quadrature points.
    double *vecAtQPs1;
    double *vecAtQPs2;

    /// Function for c.
    AbstractFunction<double, double> *f1;
    AbstractFunction<double, double> *f2;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Zero order term: \f$ f(v(\vec{x}), w(\vec{x})) u(\vec{x})\f$
   */
  class Vec2AtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    Vec2AtQP_ZOT(DOFVectorBase<double> *dv1,
		 DOFVectorBase<double> *dv2,
		 BinaryAbstractFunction<double, double, double> *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo,
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// First DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;
    /// Second DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec2;

    /// Values of the first DOFVector at the quadrature points.
    double *vecAtQPs1;
    /// Values of the second DOFVector at the quadrature points.
    double *vecAtQPs2;

    /// Function for c.
    BinaryAbstractFunction<double, double, double> *f;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * Zero order term: \f$ f(v(\vec{x}), w(\vec{x})) u(\vec{x})\f$
   */
  class Vec3AtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    Vec3AtQP_ZOT(DOFVectorBase<double> *dv1,
		 DOFVectorBase<double> *dv2,
		 DOFVectorBase<double> *dv3,
		 TertiaryAbstractFunction<double, double, double, double> *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// DOFVectors to be evaluated at quadrature points.
    DOFVectorBase<double> *vec1, *vec2, *vec3;

    /// Vectors at quadrature points.
    double *vecAtQPs1, *vecAtQPs2, *vecAtQPs3;

    /// Function for c.
    TertiaryAbstractFunction<double, double, double, double> *f;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * 
   */
  class FctGradientCoords_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    FctGradientCoords_ZOT(DOFVectorBase<double> *dv,
			  BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getC().
    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    DOFVectorBase<double>* vec;

    /// Function wich maps \ref gradAtQPs to a double.
    BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> > *f;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    WorldVector<double>* gradAtQPs;

    WorldVector<double>* coordsAtQPs;

  };

  /**
   * \ingroup Assembler
   *
   * \brief
   */
  class VecGradCoordsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    VecGradCoordsAtQP_ZOT(DOFVectorBase<double> *dv,
			  TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Vector v at quadrature points.
    double *vecAtQPs;

    /// Gradient at quadrature points.
    WorldVector<double>* gradAtQPs;

    WorldVector<double>* coordsAtQPs;

    /// Function for c.
    TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *f;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   */
  class VecAndCoordsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    VecAndCoordsAtQP_ZOT(DOFVectorBase<double> *dv,
			 BinaryAbstractFunction<double, double, WorldVector<double> > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Vector v at quadrature points.
    double *vecAtQPs;

    /// Gradient at quadrature points.
    WorldVector<double>* coordsAtQPs;

    /// Function for c.
    BinaryAbstractFunction<double, double, WorldVector<double> > *f;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   */
  class Vec2AndGradAtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    Vec2AndGradAtQP_ZOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
			TertiaryAbstractFunction<double, double, WorldVector<double>, double > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;

    /// Vector v at quadrature points.
    double *vecAtQPs1;
    double *vecAtQPs2;

    /// Gradient at quadrature points.
    WorldVector<double> *gradAtQPs;

    /// Function for c.
    TertiaryAbstractFunction<double, double, WorldVector<double>, double > *f;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * 
   */
  class FctGradient_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    FctGradient_ZOT(DOFVectorBase<double> *dv,
		    AbstractFunction<double, WorldVector<double> > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getC().
    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    DOFVectorBase<double>* vec;

    /// Function wich maps \ref gradAtQPs to a double.
    AbstractFunction<double, WorldVector<double> > *f;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    WorldVector<double>* gradAtQPs;

  };

  /**
   * \ingroup Assembler
   *
   * \brief
   */
  class VecAndGradAtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    VecAndGradAtQP_ZOT(DOFVectorBase<double> *dv,
		       BinaryAbstractFunction<double, double, WorldVector<double> > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Vector v at quadrature points.
    double *vecAtQPs;

    /// Gradient at quadrature points.
    WorldVector<double> *gradAtQPs;

    /// Function for c.
    BinaryAbstractFunction<double, double, WorldVector<double> > *f;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * Laplace multiplied with a function which maps the gradient of a DOFVector
   * at each quadrature point to a double:
   * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class VecAndGradAtQP_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    VecAndGradAtQP_SOT(DOFVectorBase<double> *dv,
		       BinaryAbstractFunction<double, double, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);


    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implements SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implements SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
  
  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Vector v at quadrature points.
    double *vecAtQPs;

    /// Gradient at quadrature points.
    WorldVector<double> *gradAtQPs;

    /// Function for c.
    BinaryAbstractFunction<double, double, WorldVector<double> > *f;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * Second order term: \f$ f(v(\vec{x}),\nabla w(\vec{x})) \Delta u(\vec{x})\f$
   */
  class VecGrad_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    VecGrad_SOT( DOFVectorBase<double> *dv1,
			  DOFVectorBase<double> *dv2,
		       BinaryAbstractFunction<double, double, WorldVector<double> > *f_)
      : SecondOrderTerm(f_->getDegree()), vec1(dv1), vec2(dv2), f(f_) 
    {}

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implements SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implements SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
  
  protected:
    /// DOFVectors to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;

    /// Vector v at quadrature points.
    double *vecAtQPs;

    /// Gradient at quadrature points.
    WorldVector<double> *gradAtQPs;

    /// Function for c.
    BinaryAbstractFunction<double, double, WorldVector<double> > *f;
  };


  /* 
   * \ingroup Assembler
   *
   * \brief
   * Zero order term: \f$ f(\vec{x}) u(\vec{x})\f$
   */
  class CoordsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    CoordsAtQP_ZOT(AbstractFunction<double, WorldVector<double> > *af)
      : ZeroOrderTerm(af->getDegree()), g(af) 
    {}

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, 
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// Stores coordinates at quadrature points. Set in \ref initElement().
    WorldVector<double>* coordsAtQPs;

    /// Function for c.
    AbstractFunction<double, WorldVector<double> > *g;
  };

  /** 
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A(x) is a WorldMatrix having a in all positions
   * except possibly the position IJ
   * \f$ \nabla \cdot A(x) \nabla u(\vec{x}) \f$
   */
  class CoordsAtQP_IJ_SOT : public SecondOrderTerm 
  {
  public:
    /// Constructor.
    CoordsAtQP_IJ_SOT(AbstractFunction<double, WorldVector<double> > *af,
		      int x_i, int x_j) 
      : SecondOrderTerm(af->getDegree()), g(af), xi(x_i), xj(x_j)
    {
      setSymmetric(xi == xj);
    }

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;

    /// Implements SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implements SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
  
  private:
    /// Stores coordinates at quadrature points. Set in \ref initElement().
    WorldVector<double>* coordsAtQPs;

    /// Function evaluated at quadrature points.
    AbstractFunction<double, WorldVector<double> > *g;

    /// Directions for the derivatives.
    int xi, xj;
  };

  /**
   * \ingroup Assembler
   *  
   * \brief
   * SecondOrderTerm where A is a WorldMatrix having a in all positions
   * except possibly the position IJ, multiplied with a function
   * evaluated at the quadrature points of a given DOFVector:
   * \f$ \nabla \cdot f(v(\vec{x})) A \nabla u(\vec{x}) \f$
   */
  class VecAtQP_IJ_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    VecAtQP_IJ_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af,
		   int x_i, int x_j);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    

    /// Implements SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor);

    /// Implements SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result); 
  
  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
    double* vecAtQPs;

    /// Function evaluated at \ref vecAtQPs.
    AbstractFunction<double, double> *f;

  private:
    /// Directions for the derivatives.
    int xi, xj;
  };

  class VecAndGradVecAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
    /// Constructor. 
    VecAndGradVecAtQP_ZOT(DOFVectorBase<double> *dv, 
			  DOFVectorBase<double> *dGrd,
			  BinaryAbstractFunction<double, double, WorldVector<double> > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected: 
    /// DOFVector to be evaluated at quadrature points. 
    DOFVectorBase<double>* vec; 

    /// Vector v at quadrature points. 
    double *vecAtQPs; 

    /// First DOFVector whose gradient is evaluated at quadrature points. 
    DOFVectorBase<double>* vecGrd; 

    /// Gradient of first vector at quadrature points. 
    WorldVector<double> *gradAtQPs; 

    /// Function for c. 
    BinaryAbstractFunction<double, double, WorldVector<double> > *f; 
  }; 

  class VecAndGradVec2AtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
    /// Constructor. 
    VecAndGradVec2AtQP_ZOT(DOFVectorBase<double> *dv, 
			   DOFVectorBase<double> *dGrd1, 
			   DOFVectorBase<double> *dGrd2, 
			   TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected: 
    /// DOFVector to be evaluated at quadrature points. 
    DOFVectorBase<double>* vec; 

    /// Vector v at quadrature points. 
    double *vecAtQPs; 

    /// First DOFVector whose gradient is evaluated at quadrature points. 
    DOFVectorBase<double>* vecGrd1; 

    /// Gradient of first vector at quadrature points. 
    WorldVector<double> *grad1AtQPs; 

    /// Second DOFVector whose gradient is evaluated at quadrature points. 
    DOFVectorBase<double>* vecGrd2; 

    /// Gradient of second vector at quadrature points. 
    WorldVector<double> *grad2AtQPs; 

    /// Function for c. 
    TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *f; 
  }; 


  class VecOfDOFVecsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
    /// Constructor. 
    VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
			 AbstractFunction<double, std::vector<double> > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected: 
    /// Vector of DOFVectors to be evaluated at quadrature points. 
    std::vector<DOFVectorBase<double>*> vecs; 

    /// Vectors at quadrature points. 
    std::vector<double*> vecsAtQPs; 

    /// Function for c. 
    AbstractFunction<double, std::vector<double> > *f; 
  }; 

  class VecOfGradientsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
    /// Constructor. 
    VecOfGradientsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
			   AbstractFunction<double, std::vector<WorldVector<double>*> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected: 
    /// Vector of DOFVectors to be evaluated at quadrature points. 
    std::vector<DOFVectorBase<double>*> vecs; 

    /// Vectors at quadrature points. 
    std::vector<WorldVector<double>*> gradsAtQPs; 

    /// Function for c. 
    AbstractFunction<double, std::vector<WorldVector<double>*> > *f; 
  };


  class VecDivergence_ZOT : public ZeroOrderTerm
  {
  public: 
    /// Constructor. 
    VecDivergence_ZOT(int nComponents,
		      DOFVectorBase<double> *vec0,
		      DOFVectorBase<double> *vec1 = NULL,
		      DOFVectorBase<double> *vec2 = NULL);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected: 
    /// Vector of DOFVectors to be evaluated at quadrature points. 
    std::vector<DOFVectorBase<double>*> vecs; 

    /// Vectors at quadrature points. 
    std::vector<WorldVector<double>*> gradsAtQPs; 
  };



  class VecAndVecOfGradientsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    VecAndVecOfGradientsAtQP_ZOT(DOFVector<double> *v,
				 const std::vector<DOFVector<double>*>& dv,
				 BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVector<double>* vec;

    /// Vector v at quadrature points.
    double *vecAtQPs;

    /// Vector of DOFVectors to be evaluated at quadrature points.
    std::vector<DOFVector<double>*> vecs;

    /// Vectors at quadrature points.
    std::vector<WorldVector<double>*> gradsAtQPs;

    /// Function for c.
    BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *f;
  };


  class Vec2AndGrad2AtQP_ZOT : public ZeroOrderTerm
  {
  public:
    Vec2AndGrad2AtQP_ZOT(DOFVectorBase<double> *dv1, 
			 DOFVectorBase<double> *dv2,
			QuartAbstractFunction<double, double, double, WorldVector<double>, WorldVector<double> > *af);

    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    DOFVectorBase<double> *vec1, *vec2;

    double *vecAtQPs1, *vecAtQPs2;

    WorldVector<double> *gradAtQPs1, *gradAtQPs2;

    QuartAbstractFunction<double, double, double, WorldVector<double>,WorldVector<double> > *f;
  };


  class Vec2AndGradVecAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
    Vec2AndGradVecAtQP_ZOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, 
			  DOFVectorBase<double> *dGrd,
			  TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f);

    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected: 
    DOFVectorBase<double>* vec1; 
    DOFVectorBase<double>* vec2; 
    
    double *vec1AtQPs; 
    double *vec2AtQPs; 

    DOFVectorBase<double>* vecGrd; 
    WorldVector<double> *gradAtQPs; 

    TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f; 
  }; 



  class General_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
		TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    std::vector<DOFVectorBase<double>*> vecs_; 

    std::vector<DOFVectorBase<double>*> grads_;

    TertiaryAbstractFunction<double, 
			     WorldVector<double>,
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f_;

    WorldVector<double> *coordsAtQPs_;

    std::vector<double*> vecsAtQPs_;

    std::vector<WorldVector<double>*> gradsAtQPs_;

    std::vector<double> vecsArg;
    
    std::vector<WorldVector<double> > gradsArg;
  };

  class GeneralParametric_ZOT : public ZeroOrderTerm
  {
  public:
    /// Constructor.
    GeneralParametric_ZOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
			  QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, std::vector<double> &C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac);

  protected:
    std::vector<DOFVectorBase<double>*> vecs_; 

    std::vector<DOFVectorBase<double>*> grads_;

    QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f_;

    WorldVector<double> *coordsAtQPs_;

    WorldVector<double> elementNormal_;

    std::vector<double*> vecsAtQPs_;

    std::vector<WorldVector<double>*> gradsAtQPs_;
  };


  /** \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