/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/



/** \file GenericOperatorTerm.h */

#ifndef AMDIS_GENERIC_OPERATOR_TERM_H
#define AMDIS_GENERIC_OPERATOR_TERM_H

#include "AMDiS_fwd.h"
#include "OperatorTerm.h"
#include "Functors.h"

#include <boost/static_assert.hpp>
#include <boost/type_traits.hpp>

#include "expressions/LazyOperatorTerm.h"
#include "expressions/expressions.h"


/** generic operator-terms provide an easy way of automated generation of
 * 'arbitrary' operator-terms out of some elementary operations, by using a 
 * recursive definition of the term. All necessary data will be initialized 
 * when an expression as part of the term uses this data.
 * Since no virtual functions, like in the AbstractFunction classes, are used
 * the overhead of a vtable is removed.
 * 
 * usage:
 * addZOT(Operator, Term) 
 *  ... add a zeroOrderTerm to Operator, i.e. (Term(x) * u, v)
 * 
 * addFOT(Operator, Term, FirstOrderType) 
 *  ... add a firstOrderTerm to Operator, (if Term::value_type = double) 
 *      i.e. (Term(x) * 1 * grad(u), v), rsp. (Term(x) * 1 * u, grad(v))
 * addFOT(Operator, Term, FirstOrderType)
 *  ... add a firstOrderTerm to Operator, (if Term::value_type = WorldVector)
 *      i.e. (Term(x) * b * grad(u), v), rsp. (Term(x) * u, grad(v))
 * addFOT<I>(Operator, Term, FirstOrderType)
 *  ... add a firstOrderTerm to Operator, 
 *      i.e. (Term(x) * e_I * grad(u), v), rsp. (Term(x) * e_I * u, grad(v))
 * 
 * addSOT(Operator, Term)
 *  ... add a secondOrderTerm to Operator, i.e. (Term(x) * grad(u), grad(v))
 * addSOT<I,J>(Operator, Term)
 *  ... add a secondOrderTerm to Operator, i.e. (E_IJ * Term(x) * grad(u), grad(v))
 * 
 * where Operator is eather a pointer or reference, FirstOrderType in {GRD_PHI, GRD_PSI}
 * and Term a componation of elementary terms by + - * /
 *  - constant(value) / value  ... a constant value
 *  - valueOf(DOFVector)       ... values of a DOFVector at QP
 *  - gradientOf(DOFVector)    ... gradient of a DOFVector at QP
 *  - derivative<I>(DOFVector) ... I'th partial derivative
 *  - X()                      ... coordinate at quadrature points
 *  - pow<I>(Term)             ... I'th power of a term
 *  - sqrt(Term)               ... square root of a term
 *  - Exp(Term)                ... exponential function of a term
 *  - function_<F>(Term)       ... evaluates F()(Term(iq))
 *  - function_(F f, Term)     ... evaluates f(Term(iq))
 * 
 * 
 * with F a functor that implements
 *   typedef (...) value_type;
 *   int getDegree(int d0);
 *   value_type operator()(const T0& v0) const;
 * 
 * respective
 *   int getDegree(int d0, int d1);
 *   value_type operator()(const T0& v0, const T1& v1) const;
 * 
 * respective
 *   int getDegree(int d0, int d1, int d2);
 *   value_type operator()(const T0& v0, const T1& v1, const T2& v2) const;
 * 
 * where the d0, d1, d2 give the polynomial degrees of the v0, v1, v2 terms.
 * */

namespace AMDiS {

// _______ ZeroOrderTerms ______________________________________________________

template<typename Term>
struct GenericZeroOrderTerm : public ZeroOrderTerm
{
  Term term;
  GenericZeroOrderTerm(const Term& term_)
    : ZeroOrderTerm(term_.getDegree()), term(term_) 
  {
    term.insertFeSpaces(auxFeSpaces);
  }

  void initElement(const ElInfo* elInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, elInfo, subAssembler, quad, NULL);
  }


  void initElement(const ElInfo* smallElInfo,
	    const ElInfo* largeElInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad, NULL);
  }


  void getC(const ElInfo *elInfo, int nPoints, ElementVector& C)
  {
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += term(iq);
  }


  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 * term(iq) * uhAtQP[iq];
  }
};

// _______ FirstOrderTerms _____________________________________________________


template<typename Term>
struct GenericFirstOrderTerm_1 : public FirstOrderTerm
{
  Term term;
  GenericFirstOrderTerm_1(const Term& term_)
    : FirstOrderTerm(term_.getDegree()), term(term_) 
  {
    term.insertFeSpaces(auxFeSpaces);
  }

  void initElement(const ElInfo* elInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, elInfo, subAssembler, quad);
  }


  void initElement(const ElInfo* smallElInfo,
	    const ElInfo* largeElInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
  }


  /// Implements FirstOrderTerm::getLb().
  void getLb(const ElInfo *elInfo,
	     std::vector<mtl::dense_vector<double> >& Lb) const
  {
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(Lb.size());

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


template<int I, typename Term>
struct GenericFirstOrderTerm_i : public FirstOrderTerm
{
  Term term;
  GenericFirstOrderTerm_i(const Term& term_)
    : FirstOrderTerm(term_.getDegree()), term(term_) 
  {
    term.insertFeSpaces(auxFeSpaces);
    FirstOrderTerm::bOne = I;
  }
  
  GenericFirstOrderTerm_i(const Term& term_, int I0)
    : FirstOrderTerm(term_.getDegree()), term(term_) 
  {
    term.insertFeSpaces(auxFeSpaces);
    FirstOrderTerm::bOne = I0;
    TEST_EXIT_DBG( I < 0 && I0 >= 0 )("You yould specify eather template<int I>, or constructor(int I0)\n");
  }

  void initElement(const ElInfo* elInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, elInfo, subAssembler, quad);
  }


  void initElement(const ElInfo* smallElInfo,
	    const ElInfo* largeElInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
  }


  /// Implements FirstOrderTerm::getLb().
  void getLb(const ElInfo *elInfo,
	     std::vector<mtl::dense_vector<double> >& Lb) const
  {
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(Lb.size());

    for (int iq = 0; iq < nPoints; iq++)
      lb_one(grdLambda, Lb[iq], term(iq));
  }
};


template<typename Term>
struct GenericFirstOrderTerm_b : public FirstOrderTerm
{
  Term term;
  
  GenericFirstOrderTerm_b(const Term& term_)
    : FirstOrderTerm(term_.getDegree()), term(term_)
  {
    term.insertFeSpaces(auxFeSpaces);
  }

  void initElement(const ElInfo* elInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, elInfo, subAssembler, quad);
  }


  void initElement(const ElInfo* smallElInfo,
	    const ElInfo* largeElInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
  }


  /// Implements FirstOrderTerm::getLb().
  void getLb(const ElInfo *elInfo,
	     std::vector<mtl::dense_vector<double> >& Lb) const
  {
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(Lb.size());

    for (int iq = 0; iq < nPoints; iq++)
      lb(grdLambda, term(iq), Lb[iq], 1.0);
  }
};

// _______ SecondOrderTerms ____________________________________________________


template<typename Term>
struct GenericSecondOrderTerm_1 : public SecondOrderTerm
{
  Term term;
  
  GenericSecondOrderTerm_1(const Term& term_)
    : SecondOrderTerm(term_.getDegree()), term(term_) 
  {
    term.insertFeSpaces(auxFeSpaces);
    setSymmetric(true);
  }

  void initElement(const ElInfo* elInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, elInfo, subAssembler, quad);
  }


  void initElement(const ElInfo* smallElInfo,
	    const ElInfo* largeElInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
  }
  
  /// Implements SecondOrderTerm::getLALt().
  void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const
  {
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(LALt.size());

    for (int iq = 0; iq < nPoints; iq++) 
      l1lt(grdLambda, LALt[iq], term(iq));
  }

  /// Implenetation of SecondOrderTerm::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 f) 
  {
    int dow = Global::getGeo(WORLD);

    if (num_rows(D2UhAtQP) > 0) {
      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 * term(iq);
      }
    }
  }

  /// 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(term(iq), grdUhAtQP[iq], result[iq]);
  }
};



template<typename Term, bool symmetric = false>
struct GenericSecondOrderTerm_A : public SecondOrderTerm
{
  Term term;
  
  GenericSecondOrderTerm_A(const Term& term_)
    : SecondOrderTerm(term_.getDegree()), term(term_) 
  {
    term.insertFeSpaces(auxFeSpaces);
    setSymmetric(symmetric);
  }

  void initElement(const ElInfo* elInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, elInfo, subAssembler, quad);
  }


  void initElement(const ElInfo* smallElInfo,
	    const ElInfo* largeElInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
  }
  
  void getLALt(const ElInfo *elInfo, 
	       std::vector<mtl::dense2D<double> > &LALt) const
  {
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(LALt.size());

    for (int iq = 0; iq < nPoints; iq++)
      lalt(grdLambda, term(iq), LALt[iq], symmetric, 1.0);
  }

  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 factor) 
  {
    int dow = Global::getGeo(WORLD);

    for (int iq = 0; iq < nPoints; iq++) {
      double resultQP = 0.0;

      WorldMatrix<double> A = term(iq);

      if (num_rows(D2UhAtQP) > 0)
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++)
	    resultQP += A[i][j] * D2UhAtQP[iq][j][i];

//       if (num_rows(grdUhAtQP) > 0)
// 	resultQP += (*divFct)(A) * grdUhAtQP[iq];

      result[iq] += resultQP * factor;
    }
  }

  void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		std::vector<WorldVector<double> > &result)
  {
    int nPoints = grdUhAtQP.size();
    WorldMatrix<double> A;
    for (int iq = 0; iq < nPoints; iq++) {
      A = term(iq);
      result[iq] += term(iq) * grdUhAtQP[iq];
    }    
  }
};



template<int I, int J, typename Term>
struct GenericSecondOrderTerm_ij : public SecondOrderTerm
{
  Term term;
  int row, col;
  
  GenericSecondOrderTerm_ij(const Term& term_)
    : SecondOrderTerm(term_.getDegree()), term(term_), row(I), col(J)
  {
    term.insertFeSpaces(auxFeSpaces);
    setSymmetric(row == col);
  }
  
  GenericSecondOrderTerm_ij(const Term& term_, int I0, int J0)
    : SecondOrderTerm(term_.getDegree()), term(term_), row(I0), col(J0)
  {
    term.insertFeSpaces(auxFeSpaces);
    setSymmetric(row == col);
    TEST_EXIT_DBG( I < 0 && I0 >= 0 && J < 0 && J0 >= 0 ) 
      ("You yould specify eather template<int I, int J>, or constructor(int I0, int J0)\n");
  }

  void initElement(const ElInfo* elInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, elInfo, subAssembler, quad);
  }


  void initElement(const ElInfo* smallElInfo,
	    const ElInfo* largeElInfo,
	    SubAssembler* subAssembler,
	    Quadrature *quad)
  {
    term.initElement(this, smallElInfo, largeElInfo, subAssembler, quad);
  }

  void getLALt(const ElInfo *elInfo, 
	       std::vector<mtl::dense2D<double> > &LALt) const
  {
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(LALt.size());

    for (int iq = 0; iq < nPoints; iq++)
      lalt_kl(grdLambda, row, col, LALt[iq], term(iq));
  }  

  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) 
  {
    if (num_rows(D2UhAtQP) > 0) {
      for (int iq = 0; iq < nPoints; iq++)
	result[iq] += D2UhAtQP[iq][row][col] * term(iq) * fac;
    }
  }

  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][row] += grdUhAtQP[iq][col] * term(iq);
  }
};


// _____________________________________________________________________________

template <typename Term>
void addZOT(Operator* op, const Term& term)
{
  op->addZeroOrderTerm(new GenericZeroOrderTerm<Term>(term));
}

template <typename Term>
void addZOT(Operator& op, const Term& term)
{
  op.addZeroOrderTerm(new GenericZeroOrderTerm<Term>(term));
}

inline void addZOT(Operator* op, double term)
{
  addZOT(op, constant(term));
}

inline void addZOT(Operator& op, double term)
{
  addZOT(op, constant(term));
}


// .............................................................................
// Rosenbrock method

// template <typename Term>
// void addZOT_(Operator* op, const Term& term)
// {
//   if (op->getUhOld()) {
//     addZOT(op, jacobian(term,0) * ...
//     
// }

// template <typename Term>
// void addZOT(Operator& op, const Term& term)
// {
//   op.addZeroOrderTerm(new GenericZeroOrderTerm<Term>(term));
// }


// _____________________________________________________________________________

// first order term using FirstOrderTerm::l1
template <typename Term>
inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
addFOT(Operator* op, const Term& term, FirstOrderType type)
{
  op->addFirstOrderTerm(new GenericFirstOrderTerm_1<Term>(term), type);
}

template <typename Term>
inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
addFOT(Operator& op, const Term& term, FirstOrderType type)
{
  op.addFirstOrderTerm(new GenericFirstOrderTerm_1<Term>(term), type);
}

inline void addFOT(Operator* op, double term, FirstOrderType type)
{
  addFOT(op, constant(term), type);
}

inline void addFOT(Operator& op, double term, FirstOrderType type)
{
  addFOT(op, constant(term), type);
}

// first order term using FirstOrderTerm::lb_one
template <int I, typename Term>
inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
addFOT(Operator* op, const Term& term, FirstOrderType type)
{
  op->addFirstOrderTerm(new GenericFirstOrderTerm_i<I,Term>(term), type);
}

template <int I, typename Term>
inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
addFOT(Operator& op, const Term& term, FirstOrderType type)
{
  op.addFirstOrderTerm(new GenericFirstOrderTerm_i<I,Term>(term), type);
}

template <int I>
inline void addFOT(Operator* op, double term, FirstOrderType type)
{
  addFOT<I>(op, constant(term), type);
}

template <int I>
inline void addFOT(Operator& op, double term, FirstOrderType type)
{
  addFOT<I>(op, constant(term), type);
}

template <typename Term>
inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
addFOT(Operator* op, const Term& term, int I, FirstOrderType type)
{
  op->addFirstOrderTerm(new GenericFirstOrderTerm_i<-1,Term>(term,I), type);
}

template <typename Term>
inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
addFOT(Operator& op, const Term& term, int I, FirstOrderType type)
{
  op.addFirstOrderTerm(new GenericFirstOrderTerm_i<-1,Term>(term,I), type);
}

inline void addFOT(Operator* op, double term, int I, FirstOrderType type)
{
  addFOT(op, constant(term), I, type);
}

inline void addFOT(Operator& op, double term, int I, FirstOrderType type)
{
  addFOT(op, constant(term), I, type);
}

// first order term using FirstOrderTerm::lb
template <typename Term>
inline typename boost::enable_if< typename boost::is_same<WorldVector<double>, typename Term::value_type>::type >::type
addFOT(Operator* op, const Term& term, FirstOrderType type)
{
  op->addFirstOrderTerm(new GenericFirstOrderTerm_b<Term>(term), type);
}

template <typename Term>
inline typename boost::enable_if< typename boost::is_same<WorldVector<double>, typename Term::value_type>::type >::type
addFOT(Operator& op, const Term& term, FirstOrderType type)
{
  op.addFirstOrderTerm(new GenericFirstOrderTerm_b<Term>(term), type);
}

inline void addFOT(Operator* op, WorldVector<double> term, FirstOrderType type)
{
  addFOT(op, constant(term), type);
}

inline void addFOT(Operator& op, WorldVector<double> term, FirstOrderType type)
{
  addFOT(op, constant(term), type);
}

// _____________________________________________________________________________

// second order term using matrix functions
template <typename Term>
inline typename boost::enable_if< typename boost::is_same<WorldMatrix<double>, typename Term::value_type>::type >::type
addSOT(Operator* op, const Term& term)
{
  op->addSecondOrderTerm(new GenericSecondOrderTerm_A<Term>(term));
}
template <typename Term>
inline typename boost::enable_if< typename boost::is_same<WorldMatrix<double>, typename Term::value_type>::type >::type
addSOT(Operator& op, const Term& term)
{
  op.addSecondOrderTerm(new GenericSecondOrderTerm_A<Term>(term));
}

inline void addSOT(Operator* op, WorldMatrix<double> term)
{
  addSOT(op, constant(term));
}

inline void addSOT(Operator& op, WorldMatrix<double> term)
{
  addSOT(op, constant(term));
}

template <bool Symmetric, typename Term>
inline typename boost::enable_if< typename boost::is_same<WorldMatrix<double>, typename Term::value_type>::type >::type
addSOT(Operator* op, const Term& term)
{
  op->addSecondOrderTerm(new GenericSecondOrderTerm_A<Term, Symmetric>(term));
}

template <bool Symmetric, typename Term>
inline typename boost::enable_if< typename boost::is_same<WorldMatrix<double>, typename Term::value_type>::type >::type
addSOT(Operator& op, const Term& term)
{
  op.addSecondOrderTerm(new GenericSecondOrderTerm_A<Term, Symmetric>(term));
}

template <bool Symmetric>
inline void addSOT(Operator* op, WorldMatrix<double> term)
{
  addSOT<Symmetric>(op, constant(term));
}

template <bool Symmetric>
inline void addSOT(Operator& op, WorldMatrix<double> term)
{
  addSOT<Symmetric>(op, constant(term));
}

// second order term using scalar functions with identity matrix
template <typename Term>
inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
addSOT(Operator* op, const Term& term)
{
  op->addSecondOrderTerm(new GenericSecondOrderTerm_1<Term>(term));
}

template <typename Term>
inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
addSOT(Operator& op, const Term& term)
{
  op.addSecondOrderTerm(new GenericSecondOrderTerm_1<Term>(term));
}

inline void addSOT(Operator* op, double term)
{
  addSOT(op, constant(term));
}

inline void addSOT(Operator& op, double term)
{
  addSOT(op, constant(term));
}

// second order term using matrix=0 with matrix_ij = function value
template <int I, int J, typename Term>
void addSOT(Operator* op, const Term& term)
{
  op->addSecondOrderTerm(new GenericSecondOrderTerm_ij<I,J,Term>(term));
}

template <int I, int J, typename Term>
void addSOT(Operator& op, const Term& term)
{
  op.addSecondOrderTerm(new GenericSecondOrderTerm_ij<I,J,Term>(term));
}

template <int I, int J>
inline void addSOT(Operator* op, double term)
{
  addSOT<I,J>(op, constant(term));
}
template <int I, int J>
inline void addSOT(Operator& op, double term)
{
  addSOT<I,J>(op, constant(term));
}

template <typename Term>
void addSOT(Operator* op, const Term& term, int I, int J)
{
  op->addSecondOrderTerm(new GenericSecondOrderTerm_ij<-1,-1,Term>(term,I,J));
}

template <typename Term>
void addSOT(Operator& op, const Term& term, int I, int J)
{
  op.addSecondOrderTerm(new GenericSecondOrderTerm_ij<-1,-1,Term>(term,I,J));
}

inline void addSOT(Operator* op, double term, int I, int J)
{
  addSOT(op, constant(term), I, J);
}

inline void addSOT(Operator& op, double term, int I, int J)
{
  addSOT(op, constant(term), I, J);
}


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

/// Create an expression functor by wrapping an AbstractFunction and evaluate it a coordinates.
template<typename TOut>
inline expressions::Function1<expressions::Wrapper<TOut,WorldVector<double> >, expressions::Coords>
eval(AbstractFunction<TOut, WorldVector<double> >* fct) { return function_(wrap(fct), X()); }

/// Integrate an expression over a domain.
/** IF the expression does not contain any DOFVector, a mesh must be given as second argument */
template<typename Term>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
integrate(Term term, Mesh* mesh_ = NULL);

// -----------------------------------------------------------------------------

/// Accumulate the values of an expression at the Degrees of freedom
template<typename Term, typename Functor>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
accumulate(Term term, Functor f, typename Term::value_type value0);

/// Maximum of an expression at DOFs, using the \ref accumulate function.
template<typename Term>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
max(Term term)
{
  typename Term::value_type value0 = -1.e25;
  value0 = accumulate(term, functors::max<typename Term::value_type>(), value0);
  
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  Parallel::mpi::globalMax(value0);
#endif
  return value0;
}

/// Minimum of an expression at DOFs, using the \ref accumulate function.
template<typename Term>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
min(Term term)
{
  typename Term::value_type value0 = 1.e25;
  value0 = accumulate(term, functors::min<typename Term::value_type>(), value0);
  
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  Parallel::mpi::globalMin(value0);
#endif
  return value0;
}

/// Maximum of  absolute values of an expression at DOFs, using the \ref accumulate function.
template<typename Term>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
abs_max(Term term)
{
  typename Term::value_type value0 = 0.0;
  value0 = accumulate(term, functors::abs_max<typename Term::value_type>(), value0);
  
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  Parallel::mpi::globalMax(value0);
#endif
  return value0;
}

/// Minimum of  absolute values of an expression at DOFs, using the \ref accumulate function.
template<typename Term>
inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
abs_min(Term term)
{
  typename Term::value_type value0 = 1.e25;
  value0 = accumulate(term, functors::abs_min<typename Term::value_type>(), value0);
  
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  Parallel::mpi::globalMin(value0);
#endif
  return value0;
}

// -----------------------------------------------------------------------------

/// Assign an expression to a DOFVector
template<typename T, typename Term>
inline typename boost::enable_if<
  typename boost::mpl::and_<typename traits::is_expr<Term>::type, 
			    typename boost::is_convertible<typename Term::value_type, T>::type
			    >::type
  >::type
transformDOF(Term term, DOFVector<T>* result);

/// Assign an expression to a DOFVector
template<typename T, typename Term>
typename boost::enable_if<
  typename boost::mpl::and_<typename traits::is_expr<Term>::type, 
			    typename boost::is_convertible<typename Term::value_type, T>::type
			    >::type,
  DOFVector<T>& >::type
operator<<(DOFVector<T>& result, const Term& term);

// -----------------------------------------------------------------------------

/// Print an expression to an output stream
template<typename Term>
typename boost::enable_if<typename traits::is_expr<Term>::type, 
			  std::ostream& >::type
operator<<(std::ostream& result, const Term& term);

} // end namespace AMDiS

#include "GenericOperatorTerm.hh"

#endif // AMDIS_GENERIC_OPERATOR_TERM_H