ZeroOrderTerm.h

// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ==  http://www.amdis-fem.org                                              ==
// ==                                                                        ==
// ============================================================================
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.



/** \file ZeroOrderTerm.h */

#ifndef AMDIS_ZERO_ORDER_TERM_H
#define AMDIS_ZERO_ORDER_TERM_H

#include "AMDiS_fwd.h"
#include "OperatorTerm.h"
#include "AbstractFunction.h"

namespace AMDiS {

  /**
   * \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, ElementVector& 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, ElementVector& C)  
    {
      for (int iq = 0; iq < nPoints; iq++)
	C[iq] += factor; 
    }

    /// Implements ZeroOrderTerm::eval().
    inline void eval( int nPoints,
		      const mtl::dense_vector<double>& uhAtQP,
		      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
		      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
		      mtl::dense_vector<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 = NULL);

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

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

    /// Implements ZeroOrderTerm::getC().
    void getC(const ElInfo *, int nPoints, ElementVector& C);

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

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

    /// Vector v at quadrature points.
    mtl::dense_vector<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, ElementVector& C);

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

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

    /// Vector v at quadrature points.
    mtl::dense_vector<double> vecAtQPs1;
    mtl::dense_vector<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 = NULL);

    /// 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, ElementVector& C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const mtl::dense_vector<double>& uhAtQP,
	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	      mtl::dense_vector<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.
    mtl::dense_vector<double> vecAtQPs1;
    /// Values of the second DOFVector at the quadrature points.
    mtl::dense_vector<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 = 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, ElementVector& C);

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

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

    /// Vectors at quadrature points.
    mtl::dense_vector<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 ZeroOrderTerm::getC().
    void getC(const ElInfo *elInfo, int nPoints, ElementVector& C);

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

  protected:
    DOFVectorBase<double>* vec;

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

    /// Pointer to a WorldVector<double> array containing the gradients of the
    /// DOFVector at each quadrature point.
    mtl::dense_vector<WorldVector<double> > gradAtQPs;

    mtl::dense_vector<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, ElementVector& C);

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

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

    /// Vector v at quadrature points.
    mtl::dense_vector<double> vecAtQPs;

    /// Gradient at quadrature points.
    mtl::dense_vector<WorldVector<double> > gradAtQPs;

    mtl::dense_vector<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, ElementVector& C);

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

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

    /// Vector v at quadrature points.
    mtl::dense_vector<double> vecAtQPs;

    mtl::dense_vector<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, ElementVector& C);

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

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

    /// Vector v at quadrature points.
    mtl::dense_vector<double> vecAtQPs1, vecAtQPs2;

    /// Gradient at quadrature points.
    mtl::dense_vector<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 ZeroOrderTerm::getC().
    void getC(const ElInfo *elInfo, int nPoints, ElementVector& C);

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

  protected:
    DOFVectorBase<double>* vec;

    /// Function which 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.
     */
    mtl::dense_vector<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, ElementVector& C);

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

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

    /// Vector v at quadrature points.
    mtl::dense_vector<double> vecAtQPs;

    /// Gradient at quadrature points.
    mtl::dense_vector<WorldVector<double> > gradAtQPs;

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


  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, ElementVector& C);

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

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

    /// Vector v at quadrature points. 
    mtl::dense_vector<double> vecAtQPs; 

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

    /// Gradient of first vector at quadrature points. 
    mtl::dense_vector<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, ElementVector& C);

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

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

    /// Vector v at quadrature points. 
    mtl::dense_vector<double> vecAtQPs; 

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

    /// Gradient of first vector at quadrature points. 
    mtl::dense_vector<WorldVector<double> > grad1AtQPs;

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

    /// Gradient of second vector at quadrature points. 
    mtl::dense_vector<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, ElementVector& C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const mtl::dense_vector<double>& uhAtQP,
	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	      mtl::dense_vector<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<mtl::dense_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, ElementVector& C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const mtl::dense_vector<double>& uhAtQP,
	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	      mtl::dense_vector<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<mtl::dense_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, ElementVector& C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const mtl::dense_vector<double>& uhAtQP,
	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	      mtl::dense_vector<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<mtl::dense_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, ElementVector& C);

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

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

    /// Vector v at quadrature points.
    mtl::dense_vector<double> vecAtQPs;

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

    /// Vectors at quadrature points.
    std::vector<mtl::dense_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, ElementVector& C);

    void eval(int nPoints,
	      const mtl::dense_vector<double>& uhAtQP,
	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	      mtl::dense_vector<double>& result,
	      double fac);

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

    mtl::dense_vector<double> vecAtQPs1, vecAtQPs2;

    mtl::dense_vector<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, ElementVector& C);

    void eval(int nPoints,
	      const mtl::dense_vector<double>& uhAtQP,
	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	      mtl::dense_vector<double>& result,
	      double fac);

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

    DOFVectorBase<double>* vecGrd; 

    mtl::dense_vector<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, ElementVector& C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const mtl::dense_vector<double>& uhAtQP,
	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	      mtl::dense_vector<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_;

    mtl::dense_vector<WorldVector<double> > coordsAtQPs;

    std::vector<mtl::dense_vector<double> > vecsAtQPs;

    std::vector<mtl::dense_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, ElementVector& C);

    /// Implements ZeroOrderTerm::eval().
    void eval(int nPoints,
	      const mtl::dense_vector<double>& uhAtQP,
	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	      mtl::dense_vector<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_;

    mtl::dense_vector<WorldVector<double> > coordsAtQPs;

    WorldVector<double> elementNormal_;

    std::vector<mtl::dense_vector<double> > vecsAtQPs;

    std::vector<mtl::dense_vector<WorldVector<double> > > gradsAtQPs;
  };


  /* 
   * \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, ElementVector& C);

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

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

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

#endif