// ============================================================================
// ==                                                                        ==
// == 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 FirstOrderTerm.h */

#ifndef AMDIS_FIRST_ORDER_TERM_H
#define AMDIS_FIRST_ORDER_TERM_H

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

namespace AMDiS {

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

  protected:
    /** \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;
      }    
    }

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

  };

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

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

}

#endif