#include "Operator.h"
#include "ElInfo.h"
#include "Assembler.h"
#include "FixVec.h"
#include "DOFVector.h"
#include "Quadrature.h"
#include "OpenMP.h"

namespace AMDiS {

  const Flag OperatorTerm::PW_CONST = 1;
  const Flag OperatorTerm::SYMMETRIC = 2;

  const Flag Operator::MATRIX_OPERATOR = 1;
  const Flag Operator::VECTOR_OPERATOR = 2;


  int Operator::getQuadratureDegree(int order, FirstOrderType firstOrderType) 
  {
    std::vector<OperatorTerm*>* terms = NULL;
    int myRank = omp_get_thread_num();

    switch(order) {
    case 0:
      terms = &zeroOrder[myRank];
      break;
    case 1:
      if (firstOrderType == GRD_PHI)
	terms = &firstOrderGrdPhi[myRank];
      else 
	terms = &firstOrderGrdPsi[myRank];
      break;
    case 2:
      terms = &secondOrder[myRank];
      break;
    }

    const BasisFunction *psi = rowFESpace->getBasisFcts();
    const BasisFunction *phi = colFESpace->getBasisFcts();

    int psiDegree = psi->getDegree();
    int phiDegree = phi->getDegree();
    int maxTermDegree = 0;

    for (int i = 0; i < static_cast<int>(terms->size()); i++) {
      OperatorTerm *term = (*terms)[i];
      maxTermDegree = max(maxTermDegree, term->degree);
    }
   
    return psiDegree + phiDegree - order + maxTermDegree;
  }


  void OperatorTerm::setSymmetric(bool symm)
  {
    if (symm) {
      properties.setFlag(SYMMETRIC);
    } else {
      properties.unsetFlag(SYMMETRIC);   
    }
  }

  bool OperatorTerm::isSymmetric()
  {
    return properties.isSet(SYMMETRIC);
  }

  double *OperatorTerm::getVectorAtQPs(DOFVectorBase<double>* vec,
				       const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad)
  {
    FUNCNAME("OperatorTerm::getVectorAtQPs()");

    TEST_EXIT(elInfo->getMesh() == vec->getFESpace()->getMesh())
      ("There is something wrong!\n");
      
    return subAssembler->getVectorAtQPs(vec, elInfo, quad);
  }

  double *OperatorTerm::getVectorAtQPs(DOFVectorBase<double>* vec,
				       const ElInfo* smallElInfo, 
				       const ElInfo* largeElInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad)
  {
    FUNCNAME("OperatorTerm::getVectorAtQPs()");

    TEST_EXIT(smallElInfo->getMesh() == vec->getFESpace()->getMesh() ||
	      largeElInfo->getMesh() == vec->getFESpace()->getMesh())
      ("There is something wrong!\n");

    if (smallElInfo->getLevel() == largeElInfo->getLevel()) {

      // Both elements have the same size, so we can use the simple procedure
      // to determine the vecAtQPs.
      
      if (vec->getFESpace()->getMesh() == smallElInfo->getMesh()) {
	return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
      } else {
	return subAssembler->getVectorAtQPs(vec, largeElInfo, quad);
      }

    } else {

      // The two elements are different. If the vector is defined on the mesh of the
      // small element, we can still use the simple procedure to determine the vecAtQPs.

      if (vec->getFESpace()->getMesh() == largeElInfo->getMesh()) {
	return subAssembler->getVectorAtQPs(vec, smallElInfo, largeElInfo, quad);
      } else {
	return subAssembler->getVectorAtQPs(vec, smallElInfo, quad);
      }

    }
  }

  WorldVector<double>* OperatorTerm::getGradientsAtQPs(DOFVectorBase<double>* vec,
						       const ElInfo* elInfo,
						       SubAssembler* subAssembler,
						       Quadrature *quad)
  {
    FUNCNAME("OperatorTerm::getGradientsAtQPs()");

    TEST_EXIT_DBG(elInfo->getMesh() == vec->getFESpace()->getMesh())
      ("There is something wrong!\n");

    return subAssembler->getGradientsAtQPs(vec, elInfo, quad);
  }

  WorldVector<double>* OperatorTerm::getGradientsAtQPs(DOFVectorBase<double>* vec,
						       const ElInfo* smallElInfo, 
						       const ElInfo* largeElInfo,
						       SubAssembler* subAssembler,
						       Quadrature *quad)
  {
    FUNCNAME("OperatorTerm::getGradientsAtQPs()");

    TEST_EXIT(smallElInfo->getMesh() == vec->getFESpace()->getMesh() ||
	      largeElInfo->getMesh() == vec->getFESpace()->getMesh())
      ("There is something wrong!\n");

    if (smallElInfo->getLevel() == largeElInfo->getLevel()) {

      // Both elements have the same size, so we can use the simple procedure
      // to determine the gradients.
      
      if (vec->getFESpace()->getMesh() == smallElInfo->getMesh()) {
	return subAssembler->getGradientsAtQPs(vec, smallElInfo, quad);
      } else {
	return subAssembler->getGradientsAtQPs(vec, largeElInfo, quad);
      }

    } else {

      // The two elements are different. If the vector is defined on the mesh of the
      // small element, we can still use the simple procedure to determine the gradients.

      if (vec->getFESpace()->getMesh() == largeElInfo->getMesh()) {
	return subAssembler->getGradientsAtQPs(vec, smallElInfo, largeElInfo, quad);
      } else {
	return subAssembler->getGradientsAtQPs(vec, smallElInfo, quad);
      }

    }
  }

  void OperatorTerm::lalt(const DimVec<WorldVector<double> >& Lambda,
			  const WorldMatrix<double>& matrix,
			  DimMat<double>& LALt,
			  bool symm,
			  double factor)
  {
    int j, k, l;
    const int dimOfWorld = Global::getGeo(WORLD);
    int dim = LALt.getNumRows() - 1;
  
    double val = 0.0;

    if (symm) {
      for (int i = 0; i <= dim; i++) {
	val = 0.0;
	for (k = 0; k < dimOfWorld; k++)
	  for (l = 0; l < dimOfWorld; l++)
	    val += Lambda[i][k] * matrix[k][l] * Lambda[i][l];
	val *= factor;
	LALt[i][i] += val;
	for (j = i + 1; j <= dim; j++) {
	  val = 0.0;
	  for (k = 0; k < dimOfWorld; k++)
	    for (l = 0; l < dimOfWorld; l++)
	      val += Lambda[i][k] * matrix[k][l] * Lambda[j][l];
	  val *= factor;
	  LALt[i][j] += val;
	  LALt[j][i] += val;
	}
      }    
    } else {	
      for (int i = 0; i <= dim; i++) {
	for (j = 0; j <= dim; j++) {
	  val = 0.0;
	  for (k = 0; k < dimOfWorld; k++)
	    for (l = 0; l < dimOfWorld; l++) 
	      val += Lambda[i][k] * matrix[k][l] * Lambda[j][l];
	  val *= factor;
	  LALt[i][j] += val;
	}
      }    
    }      
  }

  void OperatorTerm::lalt_kl(const DimVec<WorldVector<double> >& Lambda,
			     int k, int l,
			     DimMat<double>& LALt,
			     double factor)
  {
    int dim = LALt.getNumRows() - 1;

    for (int i = 0; i <= dim; i++)
      for (int j = 0; j <= dim; j++)
	LALt[i][j] += factor * Lambda[i][k] * Lambda[j][l];
  }

  void OperatorTerm::l1lt(const DimVec<WorldVector<double> >& Lambda,
			  DimMat<double>& LALt,
			  double factor)
  {
    const int dimOfWorld = Global::getGeo(WORLD);
    int dim = LALt.getNumRows() - 1;

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

  void OperatorTerm::lb(const DimVec<WorldVector<double> >& Lambda,
			const WorldVector<double>& b,
			DimVec<double>& Lb,
			double factor)
  {
    int dim = Lb.getSize() - 1;
    const int dimOfWorld = Global::getGeo(WORLD);

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

  Operator::Operator(Flag operatorType,
		     const FiniteElemSpace *row,
		     const FiniteElemSpace *col)
    : rowFESpace(row), 
      colFESpace(col ? col : row),
      type(operatorType), 
      fillFlag(Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
	       Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA),
      needDualTraverse(false),
      uhOld(NULL),
      optimized(true)
  {
    int maxThreads = omp_get_overall_max_threads();

    assembler.resize(maxThreads);
    secondOrder.resize(maxThreads);
    firstOrderGrdPsi.resize(maxThreads);
    firstOrderGrdPhi.resize(maxThreads);
    zeroOrder.resize(maxThreads);

    for (int i = 0; i < maxThreads; i++) {
      assembler[i] = NULL;
      secondOrder[i].resize(0);
      firstOrderGrdPsi[i].resize(0);
      firstOrderGrdPhi[i].resize(0);
      zeroOrder[i].resize(0);
    }

    nRow = rowFESpace->getBasisFcts()->getNumber();
    nCol = colFESpace->getBasisFcts()->getNumber();
  }

  void Operator::setUhOld(const DOFVectorBase<double> *uhOld_)
  {
    uhOld = uhOld_;
  }

  void Operator::getElementMatrix(const ElInfo *elInfo, 
				  ElementMatrix& userMat, 
				  double factor)
  {
    int myRank = omp_get_thread_num();

    if (!assembler[myRank])
      initAssembler(myRank, NULL, NULL, NULL, NULL);

    assembler[myRank]->calculateElementMatrix(elInfo, userMat, factor);
  }

  void Operator::getElementMatrix(const ElInfo *rowElInfo, const ElInfo *colElInfo,
				  const ElInfo *smallElInfo, const ElInfo *largeElInfo,
				  ElementMatrix& userMat, 
				  double factor)
  {
    int myRank = omp_get_thread_num();

    if (!assembler[myRank])
      initAssembler(myRank, NULL, NULL, NULL, NULL);

    assembler[myRank]->calculateElementMatrix(rowElInfo, colElInfo, 
					      smallElInfo, largeElInfo,
    					      userMat, factor);
  }

  void Operator::getElementVector(const ElInfo *elInfo, 
				  ElementVector& userVec, 
				  double factor)
  {
    int myRank = omp_get_thread_num();

    if (!assembler[myRank])
      initAssembler(myRank, NULL, NULL, NULL, NULL);

    assembler[myRank]->calculateElementVector(elInfo, userVec, factor);
  }

  void Operator::getElementVector(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo, const ElInfo *largeElInfo,
				  ElementVector& userVec,
				  double factor)
  {
    int myRank = omp_get_thread_num();

    if (!assembler[myRank])
      initAssembler(myRank, NULL, NULL, NULL, NULL);

    assembler[myRank]->calculateElementVector(mainElInfo, auxElInfo, 
					      smallElInfo, largeElInfo,
					      userVec, factor);
  }

  void Operator::initAssembler(int rank,
			       Quadrature *quad2,
			       Quadrature *quad1GrdPsi,
			       Quadrature *quad1GrdPhi,
			       Quadrature *quad0) 
  {    
#ifdef _OPENMP
#pragma omp critical (initAssembler)
#endif
      {
	if (optimized) {
	  assembler[rank] = new OptimizedAssembler(this,
						   quad2, quad1GrdPsi, quad1GrdPhi, quad0,
						   rowFESpace, colFESpace);
	} else {
	  assembler[rank] = new StandardAssembler(this,
						  quad2, quad1GrdPsi, quad1GrdPhi, quad0,
						  rowFESpace, colFESpace);
	}
      }
  }

  void Operator::finishAssembling()
  {
    assembler[omp_get_thread_num()]->finishAssembling();
  }



  VecAtQP_ZOT::VecAtQP_ZOT(DOFVectorBase<double> *dv,
			   AbstractFunction<double, double> *af)
    : ZeroOrderTerm(af ? af->getDegree() : 0), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());        
  }

  MultVecAtQP_ZOT::MultVecAtQP_ZOT(DOFVectorBase<double> *dv1,
				   DOFVectorBase<double> *dv2,
				   AbstractFunction<double, double> *af1,
				   AbstractFunction<double, double> *af2)
    : ZeroOrderTerm(af1->getDegree() + af2->getDegree()), 
      vec1(dv1), vec2(dv2), f1(af1), f2(af2)
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

    auxFESpaces.push_back(dv1->getFESpace());
    auxFESpaces.push_back(dv2->getFESpace());
  }

  Vec2AtQP_ZOT::Vec2AtQP_ZOT(DOFVectorBase<double> *dv1,
			     DOFVectorBase<double> *dv2,
			     BinaryAbstractFunction<double, double, double> *af)
    : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af)
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

    auxFESpaces.push_back(dv1->getFESpace());
    auxFESpaces.push_back(dv2->getFESpace());
  }

  Vec3AtQP_ZOT::Vec3AtQP_ZOT(DOFVectorBase<double> *dv1,
			     DOFVectorBase<double> *dv2,
			     DOFVectorBase<double> *dv3,
			     TertiaryAbstractFunction<double, double, double, double> *af)
    : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), vec3(dv3), f(af)
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");
    TEST_EXIT(dv3)("No thierd vector!\n");

    auxFESpaces.push_back(dv1->getFESpace());
    auxFESpaces.push_back(dv2->getFESpace());
    auxFESpaces.push_back(dv3->getFESpace());
  }

  FctGradientCoords_ZOT::FctGradientCoords_ZOT(DOFVectorBase<double> *dv,
					       BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecGradCoordsAtQP_ZOT::VecGradCoordsAtQP_ZOT(DOFVectorBase<double> *dv,
					       TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAndCoordsAtQP_ZOT::VecAndCoordsAtQP_ZOT(DOFVectorBase<double> *dv,
					     BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  Vec2AndGradAtQP_ZOT::Vec2AndGradAtQP_ZOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
					   TertiaryAbstractFunction<double, double, WorldVector<double>, double > *af)
    : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af) 
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

    auxFESpaces.push_back(dv1->getFESpace());
    auxFESpaces.push_back(dv2->getFESpace());
  }

  FctGradient_ZOT::FctGradient_ZOT(DOFVectorBase<double> *dv,
				   AbstractFunction<double, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAndGradAtQP_ZOT::VecAndGradAtQP_ZOT(DOFVectorBase<double> *dv,
					 BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAndGradVecAtQP_ZOT::VecAndGradVecAtQP_ZOT(DOFVectorBase<double> *dv, 
					       DOFVectorBase<double> *dGrd,
					       BinaryAbstractFunction<double, double, WorldVector<double> > *af) 
    : ZeroOrderTerm(af->getDegree()), vec(dv), vecGrd(dGrd), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");
    TEST_EXIT(dGrd)("No gradient vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
    auxFESpaces.push_back(dGrd->getFESpace());
  }

  VecAndGradVec2AtQP_ZOT::VecAndGradVec2AtQP_ZOT(DOFVectorBase<double> *dv, 
						 DOFVectorBase<double> *dGrd1, 
						 DOFVectorBase<double> *dGrd2, 
						 TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *af) 
    : ZeroOrderTerm(af->getDegree()), vec(dv), vecGrd1(dGrd1), vecGrd2(dGrd2), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");
    TEST_EXIT(dGrd1)("No first gradient vector!\n");
    TEST_EXIT(dGrd2)("No second gradient vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
    auxFESpaces.push_back(dGrd1->getFESpace());
    auxFESpaces.push_back(dGrd2->getFESpace());
  }

  VecOfDOFVecsAtQP_ZOT::VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
					     AbstractFunction<double, std::vector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vecs(dv), f(af) 
  {
    vecsAtQPs.resize(vecs.size());

    for (int i = 0; i < static_cast<int>(dv.size()); i++) {
      TEST_EXIT(dv[i])("One vector is NULL!\n");

      auxFESpaces.push_back(dv[i]->getFESpace());
    }
  } 

  VecOfGradientsAtQP_ZOT::VecOfGradientsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv,
						 AbstractFunction<double, std::vector<WorldVector<double>*> > *af) 
    : ZeroOrderTerm(af->getDegree()), vecs(dv), f(af) 
  {
    gradsAtQPs.resize(vecs.size());

    for (int i = 0; i < static_cast<int>(dv.size()); i++) {
      TEST_EXIT(dv[i])("One vector is NULL!\n");

      auxFESpaces.push_back(dv[i]->getFESpace());
    }
  } 

  VecDivergence_ZOT::VecDivergence_ZOT(int nComponents,
				       DOFVectorBase<double> *vec0,
				       DOFVectorBase<double> *vec1,
				       DOFVectorBase<double> *vec2)
    : ZeroOrderTerm(0)
  {
    vecs.resize(nComponents);
    gradsAtQPs.resize(nComponents);
    vecs[0] = vec0;
    vecs[1] = vec1;
    vecs[2] = vec2;

    auxFESpaces.push_back(vec0->getFESpace());
    if (vec1) 
      auxFESpaces.push_back(vec1->getFESpace());
    if (vec2) 
      auxFESpaces.push_back(vec2->getFESpace());
  }


  VecAndVecOfGradientsAtQP_ZOT::VecAndVecOfGradientsAtQP_ZOT(DOFVector<double> *v,
							     const std::vector<DOFVector<double>*>& dv,
							     BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(v), vecs(dv), f(af)
  {
    gradsAtQPs.resize(vecs.size());

    TEST_EXIT(v)("No vector!\n");

    auxFESpaces.push_back(v->getFESpace());
    for (int i = 0; i < static_cast<int>(dv.size()); i++) {
      TEST_EXIT(dv[i])("One gradient vector is NULL!\n");

      auxFESpaces.push_back(dv[i]->getFESpace());
    }
  }

  General_ZOT::General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
			   std::vector<DOFVectorBase<double>*> grads,
			   TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *af)
    : ZeroOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
  {
    vecsAtQPs_.resize(vecs_.size());
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

      auxFESpaces.push_back(vecs[i]->getFESpace());
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

      auxFESpaces.push_back(grads[i]->getFESpace());
    }   
  }

  GeneralParametric_ZOT::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> > > *af)
    : ZeroOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
  {
    vecsAtQPs_.resize(vecs_.size());
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

      auxFESpaces.push_back(vecs[i]->getFESpace());
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

      auxFESpaces.push_back(grads[i]->getFESpace());
    }   
  }

  VecAtQP_FOT::VecAtQP_FOT(DOFVectorBase<double> *dv,
			   AbstractFunction<double, double> *af,
			   WorldVector<double> *wv)
    : FirstOrderTerm(af->getDegree()), vec(dv), f(af), b(wv)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VectorGradient_FOT::VectorGradient_FOT(DOFVectorBase<double> *dv,
					 AbstractFunction<WorldVector<double>, WorldVector<double> > *af)
    : FirstOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VectorFct_FOT::VectorFct_FOT(DOFVectorBase<double> *dv,
			       AbstractFunction<WorldVector<double>, double> *fct)
    : FirstOrderTerm(fct->getDegree()), vec(dv), vecFct(fct) 
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecGrad_FOT::VecGrad_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
			   BinaryAbstractFunction<WorldVector<double>, double, WorldVector<double> > *fct)
    : FirstOrderTerm(fct->getDegree()), vec1(dv1), vec2(dv2), vecFct(fct) 
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

    auxFESpaces.push_back(dv1->getFESpace());
    auxFESpaces.push_back(dv2->getFESpace());    
  }

  General_FOT::General_FOT(std::vector<DOFVectorBase<double>*> vecs,
			   std::vector<DOFVectorBase<double>*> grads,
			   TertiaryAbstractFunction<WorldVector<double>, 
			   WorldVector<double>,
			   std::vector<double>, 
			   std::vector<WorldVector<double> > > *af)
    : FirstOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
  {
    vecsAtQPs_.resize(vecs_.size());
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

      auxFESpaces.push_back(vecs[i]->getFESpace());
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

      auxFESpaces.push_back(grads[i]->getFESpace());
    }   
  }

  GeneralParametric_FOT::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> > > *af)
    : FirstOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
  {
    vecsAtQPs_.resize(vecs_.size());
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

      auxFESpaces.push_back(vecs[i]->getFESpace());
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

      auxFESpaces.push_back(grads[i]->getFESpace());
    }   
  }

  MatrixFct_SOT::MatrixFct_SOT(DOFVectorBase<double> *dv, 
			       AbstractFunction<WorldMatrix<double>, double> *fct,
			       AbstractFunction<WorldVector<double>, WorldMatrix<double> > *div,
			       bool sym)
    : SecondOrderTerm(fct->getDegree()), 
      vec(dv), 
      matrixFct(fct), 
      divFct(div),
      symmetric(sym)
  {
    setSymmetric(symmetric);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAtQP_SOT::VecAtQP_SOT(DOFVectorBase<double> *dv, 
			   AbstractFunction<double, double> *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    setSymmetric(true);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  Vec2AtQP_SOT::Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, 
			     BinaryAbstractFunction<double, double, double> *af)
    : SecondOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af)
  {
    setSymmetric(true);

    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

    auxFESpaces.push_back(dv1->getFESpace());
    auxFESpaces.push_back(dv2->getFESpace());
  }

  MatrixGradient_SOT::MatrixGradient_SOT(DOFVectorBase<double> *dv,
					 AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
					 AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
					 bool symm) 
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af), divFct(divAf), symmetric(symm)
  {
    setSymmetric(symmetric);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  FctGradient_SOT::FctGradient_SOT(DOFVectorBase<double> *dv,
				   AbstractFunction<double, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAndGradientAtQP_SOT::VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv,
						 BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecMatrixGradientAtQP_SOT::VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
						       BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
						       AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
						       bool symm)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af), 
      divFct(divAf), symmetric(symm)
  {
    setSymmetric(symmetric);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecGradCoordsAtQP_SOT::VecGradCoordsAtQP_SOT(DOFVectorBase<double> *dv,
					       TertiaryAbstractFunction<double, double,
					       WorldVector<double>, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAndCoordsAtQP_SOT::VecAndCoordsAtQP_SOT(DOFVectorBase<double> *dv, 
					     BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    setSymmetric(true);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  MatrixGradientAndCoords_SOT::MatrixGradientAndCoords_SOT(DOFVectorBase<double> *dv,
							   BinaryAbstractFunction<WorldMatrix<double>,
							   WorldVector<double>, WorldVector<double> > *af,
							   AbstractFunction<WorldVector<double>,
							   WorldMatrix<double> > *divAf,
							   bool symm) 
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af), divFct(divAf), symmetric(symm)
  {
    setSymmetric(symmetric);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  General_SOT::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)
      : SecondOrderTerm(f->getDegree()),
	vecs_(vecs),
	grads_(grads),
	f_(f),
	divFct_(divFct),
	symmetric_(symmetric)
  {
    vecsAtQPs_.resize(vecs_.size());
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

      auxFESpaces.push_back(vecs[i]->getFESpace());
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

      auxFESpaces.push_back(grads[i]->getFESpace());
    }
  }

  GeneralParametric_SOT::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)
    : SecondOrderTerm(f->getDegree()),
      vecs_(vecs),
      grads_(grads),
      f_(f),
      divFct_(divFct),
      symmetric_(symmetric)
  {
    vecsAtQPs_.resize(vecs_.size());
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

      auxFESpaces.push_back(vecs[i]->getFESpace());
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

      auxFESpaces.push_back(grads[i]->getFESpace());
    }
  }

  VecAndGradAtQP_SOT::VecAndGradAtQP_SOT(DOFVectorBase<double> *dv,
					 BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAtQP_IJ_SOT::VecAtQP_IJ_SOT(DOFVectorBase<double> *dv, 
				 AbstractFunction<double, double> *af,
				 int x_i, int x_j)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af), xi(x_i), xj(x_j)
  {
    setSymmetric(xi == xj);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }



  void MatrixFct_SOT::initElement(const ElInfo* elInfo, 
				  SubAssembler* subAssembler,
				  Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
  }
 
  void VecAtQP_SOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature* quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
  }

  void VecAtQP_SOT::initElement(const ElInfo* smallElInfo, 
				const ElInfo* largeElInfo,
				SubAssembler* subAssembler,
				Quadrature* quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad);
  }

  void Vec2AtQP_SOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
  }
 
  void CoordsAtQP_SOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad) 
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void MatrixGradient_SOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }

  void FctGradient_SOT::initElement(const ElInfo* elInfo, 
				    SubAssembler* subAssembler,
				    Quadrature *quad) 
  {
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }

  void VecAndGradientAtQP_SOT::initElement(const ElInfo* elInfo, 
					   SubAssembler* subAssembler,
					   Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }

  void VecMatrixGradientAtQP_SOT::initElement(const ElInfo* elInfo, 
					      SubAssembler* subAssembler,
					      Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }
  
  void VecAndCoordsAtQP_SOT::initElement(const ElInfo* elInfo, 
					 SubAssembler* subAssembler,
					 Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }
 
  void MatrixGradientAndCoords_SOT::initElement(const ElInfo* elInfo, 
						SubAssembler* subAssembler,
						Quadrature *quad) 
  {
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);  
  }

  void VecGradCoordsAtQP_SOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);  
  }

  void VecAtQP_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
  }

  void CoordsAtQP_FOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad) 
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void VecCoordsAtQP_FOT::initElement(const ElInfo* elInfo, 
				      SubAssembler* subAssembler,
				      Quadrature *quad) 
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void VectorGradient_FOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }

  void VectorFct_FOT::initElement(const ElInfo* elInfo, 
				  SubAssembler* subAssembler,
				  Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
  }
 
  void VecFctAtQP_FOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad) 
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void VecAtQP_ZOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
  }

  void VecAtQP_ZOT::initElement(const ElInfo* smallElInfo,
				const ElInfo* largeElInfo,
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
    vecAtQPs = getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad);
  }


  void MultVecAtQP_ZOT::initElement(const ElInfo* elInfo, 
				    SubAssembler* subAssembler,
				    Quadrature *quad) 
  {
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
  }

  void Vec2AtQP_ZOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
    TEST_EXIT(vec1->getFESpace() == vec2->getFESpace())
      ("Not yet implemented!\n");

    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
  }

  void Vec2AtQP_ZOT::initElement(const ElInfo* smallElInfo, 
				 const ElInfo* largeElInfo,
				 SubAssembler* subAssembler,
				 Quadrature *quad)
  {
    TEST_EXIT(vec1->getFESpace() == vec2->getFESpace())
      ("Not yet implemented!\n");

    vecAtQPs1 = getVectorAtQPs(vec1, smallElInfo, largeElInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, smallElInfo, largeElInfo, subAssembler, quad);
  }

  void Vec3AtQP_ZOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
    vecAtQPs3 = getVectorAtQPs(vec3, elInfo, subAssembler, quad);
  }

  void VecAndCoordsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					 SubAssembler* subAssembler,
					 Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }
  
  void Vec2AndGradAtQP_ZOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
  {
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec1, elInfo, subAssembler, quad);
  }

  void FctGradientCoords_ZOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }
 
  void VecGradCoordsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }
 
  void VecAndGradAtQP_ZOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }
 
  void VecAndGradVecAtQP_ZOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vecGrd, elInfo, subAssembler, quad);
  }

  void VecAndGradVec2AtQP_ZOT::initElement(const ElInfo* elInfo, 
					   SubAssembler* subAssembler,
					   Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    grad1AtQPs = getGradientsAtQPs(vecGrd1, elInfo, subAssembler, quad);
    grad2AtQPs = getGradientsAtQPs(vecGrd2, elInfo, subAssembler, quad);
  }
 
  void FctGradient_ZOT::initElement(const ElInfo* elInfo, 
				    SubAssembler* subAssembler,
				    Quadrature *quad) 
  {
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }
 
  void CoordsAtQP_ZOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad) {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void MatrixFct_SOT::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, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0);
  }

  void VecAtQP_SOT::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]), (*f)(vecAtQPs[iq]));
  }

  void Vec2AtQP_SOT::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]), (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
  }

  void CoordsAtQP_SOT::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]), (*g)(coordsAtQPs[iq]));
  }

  void MatrixGradient_SOT::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, (*f)(gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
    }
  }

  void FctGradient_SOT::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]), (*f)(gradAtQPs[iq]));
  }
  
  void VecAndGradientAtQP_SOT::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]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
  }

  void VecMatrixGradientAtQP_SOT::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, (*f)(vecAtQPs[iq], gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
  }

  void VecAndCoordsAtQP_SOT::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]), (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
  }
  
  void MatrixGradientAndCoords_SOT::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, (*f)(gradAtQPs[iq], coordsAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
  }

  void VecGradCoordsAtQP_SOT::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]), (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]));
  }

  void VecAtQP_FOT::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++) {
      if (b)
	lb(Lambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
      else
	l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
    }
  }

  void CoordsAtQP_FOT::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], (*g)(coordsAtQPs[iq]));
    }
  }

  void VecCoordsAtQP_FOT::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++) {
      lb(Lambda, b, Lb[iq], (*g)(coordsAtQPs[iq]));
    }
  }

  void VectorGradient_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    if (f) {
      for (int iq = 0; iq < nPoints; iq++) {
	lb(Lambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
      }
    } else {
      for (int iq = 0; iq < nPoints; iq++) {
	lb(Lambda, gradAtQPs[iq], Lb[iq], 1.0);
      }
    }
  }

  void VectorFct_FOT::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++) {
      lb(Lambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
    }
  }

  void VecFctAtQP_FOT::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++) {
      lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
    }
  }

  void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
    if (f) {
      for (int iq = 0; iq < nPoints; iq++) {
	C[iq] += (*f)(vecAtQPs[iq]);
      }
    } else {
      for (int iq = 0; iq < nPoints; iq++) {
	C[iq] += vecAtQPs[iq];
      }
    }
  }

  void VecAtQP_ZOT::eval(int nPoints,
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double fac) const
  {
    if (f) {
      for (int iq = 0; iq < nPoints; iq++) {
	result[iq] += fac * (*f)(vecAtQPs[iq]) * uhAtQP[iq];
      }
    } else {
      for (int iq = 0; iq < nPoints; iq++) {
	result[iq] += fac * vecAtQPs[iq] * uhAtQP[iq];
      }      
    }
  }


  void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
    for (int iq = 0; iq < nPoints; iq++) {
      C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);
    }
  }

  void MultVecAtQP_ZOT::eval(int nPoints,
			     const double *uhAtQP,
			     const WorldVector<double> *grdUhAtQP,
			     const WorldMatrix<double> *D2UhAtQP,
			     double *result,
			     double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++) {
      result[iq] += fac * (*f1)(vecAtQPs1[iq]) * 
	(*f2)(vecAtQPs2[iq]) * uhAtQP[iq];
    }
  }

  void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
    for (int iq = 0; iq < nPoints; iq++) {
      C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
    }
  }

  void Vec2AtQP_ZOT::eval(int nPoints,
			  const double *uhAtQP,
			  const WorldVector<double> *grdUhAtQP,
			  const WorldMatrix<double> *D2UhAtQP,
			  double *result,
			  double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++) {
      result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * 
	uhAtQP[iq];
    }
  }

  void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
    for (int iq = 0; iq < nPoints; iq++) {
      C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
    }
  }

  void Vec3AtQP_ZOT::eval(int nPoints,
			  const double *uhAtQP,
			  const WorldVector<double> *grdUhAtQP,
			  const WorldMatrix<double> *D2UhAtQP,
			  double *result,
			  double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++) {
      result[iq] += 
	fac * 
	(*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) * 
	uhAtQP[iq];
    }
  }


  void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, 
				  std::vector<double> &C) const { 
    for (int iq = 0; iq < nPoints; iq++) {
      C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
    }
  }

  void VecAndCoordsAtQP_ZOT::eval(int nPoints,
				  const double *uhAtQP,
				  const WorldVector<double> *grdUhAtQP,
				  const WorldMatrix<double> *D2UhAtQP,
				  double *result,
				  double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++) {
      result[iq] += fac * (*f)(vecAtQPs[iq], coordsAtQPs[iq]) * 
	uhAtQP[iq];
    }
  }

 
  void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints, 
				   std::vector<double> &C) const { 
    for (int iq = 0; iq < nPoints; iq++) {
      C[iq] += (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
    }
  }
 
  void FctGradientCoords_ZOT::eval(int nPoints,
				   const double *uhAtQP,
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++) {
      result[iq] += fac * (*f)(gradAtQPs[iq], coordsAtQPs[iq]) * 
	uhAtQP[iq];
    }
  }

  void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints,
				   std::vector<double> &C) const { 
    for (int iq = 0; iq < nPoints; iq++) {
      C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
    }
  }
 
  void VecGradCoordsAtQP_ZOT::eval(int nPoints,
				   const double *uhAtQP,
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * 
	(*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]) * 
	uhAtQP[iq];
  }

  void Vec2AndGradAtQP_ZOT::eval(int nPoints,
				 const double *uhAtQP,
				 const WorldVector<double> *grdUhAtQP,
				 const WorldMatrix<double> *D2UhAtQP,
				 double *result,
				 double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * 
	(*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]) * 
	uhAtQP[iq];
  }

  void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int nPoints,
				 std::vector<double> &C) const 
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]);
  }

  void VecAndGradAtQP_ZOT::getC(const ElInfo *, int nPoints,
				std::vector<double> &C) const 
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
  }
 
  void VecAndGradAtQP_ZOT::eval(int nPoints,
				const double *uhAtQP,
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
				double *result,
				double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];   
  }

  void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints,
				   std::vector<double> &C) const 
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
  }
 
  void VecAndGradVecAtQP_ZOT::eval(int nPoints,
				   const double *uhAtQP,
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];
  }

  void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int nPoints,
				    std::vector<double> &C) const { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]);
  }
 
  void VecAndGradVec2AtQP_ZOT::eval(int nPoints,
				    const double *uhAtQP,
				    const WorldVector<double> *grdUhAtQP,
				    const WorldMatrix<double> *D2UhAtQP,
				    double *result,
				    double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]) * uhAtQP[iq];  
  }

  void FctGradient_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(gradAtQPs[iq]);
  }
 
  void FctGradient_ZOT::eval(int nPoints,
			     const double *uhAtQP,
			     const WorldVector<double> *grdUhAtQP,
			     const WorldMatrix<double> *D2UhAtQP,
			     double *result,
			     double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(gradAtQPs[iq]) * uhAtQP[iq];    
  }

  void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, 
			    std::vector<double> &C) const { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*g)(coordsAtQPs[iq]);    
  }

  void CoordsAtQP_ZOT::eval(int nPoints,
			    const double *uhAtQP,
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double fac) const
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*g)(coordsAtQPs[iq]) * uhAtQP[iq];    
  }

  void MatrixFct_SOT::eval(int nPoints,
			   const double *uhAtQP,
			   const WorldVector<double> *grdUhAtQP,
			   const WorldMatrix<double> *D2UhAtQP,
			   double *result,
			   double factor) const
  {
    int dow = Global::getGeo(WORLD);

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

      WorldMatrix<double> A = (*matrixFct)(vecAtQPs[iq]);

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

      if (grdUhAtQP) {
	resultQP += (*divFct)(A) * grdUhAtQP[iq];
      }

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

  void MatrixFct_SOT::weakEval(int nPoints,
			       const WorldVector<double> *grdUhAtQP,
			       WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      WorldMatrix<double> A;
      for (int iq = 0; iq < nPoints; iq++) {
	A = (*matrixFct)(vecAtQPs[iq]);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }


  void Matrix_SOT::eval(int nPoints,
			const double *uhAtQP,
			const WorldVector<double> *grdUhAtQP,
			const WorldMatrix<double> *D2UhAtQP,
			double *result,
			double factor) const
  {
    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++) {
	  for (int j = 0; j < dow; j++) {
	    resultQP += matrix[i][j] * D2UhAtQP[iq][j][i];
	  }
	}
	result[iq] += resultQP * factor;
      }
    }
  }

  void Matrix_SOT::weakEval(int nPoints,
			    const WorldVector<double> *grdUhAtQP,
			    WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	result[iq] += matrix * grdUhAtQP[iq];
      }
    }
  }



  void MatrixGradient_SOT::eval(int nPoints,
				const double *uhAtQP,
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
				double *result,
				double factor) const
  {
    int dow = Global::getGeo(WORLD);

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

      WorldMatrix<double> A = (*f)(gradAtQPs[iq]);

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

      if (grdUhAtQP) {
	resultQP += (*divFct)(A) * grdUhAtQP[iq];
      }

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

  void MatrixGradient_SOT::weakEval(int nPoints,
				    const WorldVector<double> *grdUhAtQP,
				    WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      WorldMatrix<double> A;
      for (int iq = 0; iq < nPoints; iq++) {
	A = (*f)(gradAtQPs[iq]);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }

  void VecAtQP_SOT::eval(int nPoints,
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double fac) const
  {
    int dow = Global::getGeo(WORLD);

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

  void VecAtQP_SOT::weakEval(int nPoints,
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

  void Vec2AtQP_SOT::eval(int nPoints,
			  const double *uhAtQP,
			  const WorldVector<double> *grdUhAtQP,
			  const WorldMatrix<double> *D2UhAtQP,
			  double *result,
			  double fac) const
  {
    int dow = Global::getGeo(WORLD);
    
    if (D2UhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
	double resultQP = 0.0;
	for (int i = 0; i < dow; i++) {
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += fac * factor * resultQP;
      }
    }
  }
  
  void Vec2AtQP_SOT::weakEval(int nPoints,
			      const WorldVector<double> *grdUhAtQP,
			      WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

  void CoordsAtQP_SOT::eval(int nPoints,
			    const double *uhAtQP,
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double f) const
  {
    int dow = Global::getGeo(WORLD);

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

  void CoordsAtQP_SOT::weakEval(int nPoints,
				const WorldVector<double> *grdUhAtQP,
				WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*g)(coordsAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

  void FctGradient_SOT::eval(int nPoints,
			     const double *uhAtQP,
			     const WorldVector<double> *grdUhAtQP,
			     const WorldMatrix<double> *D2UhAtQP,
			     double *result,
			     double fac) const
  {
    int dow = Global::getGeo(WORLD);

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

  void FctGradient_SOT::weakEval(int nPoints,
				 const WorldVector<double> *grdUhAtQP,
				 WorldVector<double> *result) const
  {
    if (grdUhAtQP) {  
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(gradAtQPs[iq]);
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

  void VecAndGradientAtQP_SOT::eval(int nPoints,
				    const double *uhAtQP,
				    const WorldVector<double> *grdUhAtQP,
				    const WorldMatrix<double> *D2UhAtQP,
				    double *result,
				    double fac) const
  {
    int dow = Global::getGeo(WORLD);

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

  void VecAndGradientAtQP_SOT::weakEval(int nPoints,
					const WorldVector<double> *grdUhAtQP,
					WorldVector<double> *result) const
  {
    if (grdUhAtQP) {  
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

  void VecMatrixGradientAtQP_SOT::eval(int nPoints,
				       const double *uhAtQP,
				       const WorldVector<double> *grdUhAtQP,
				       const WorldMatrix<double> *D2UhAtQP,
				       double *result,
				       double factor) const
  {
    int dow = Global::getGeo(WORLD);
    
    for (int iq = 0; iq < nPoints; iq++) {
      double resultQP = 0.0;
      
      WorldMatrix<double> A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
      
      if (D2UhAtQP) {
	for (int i = 0; i < dow; i++) {
	  for (int j = 0; j < dow; j++) {
	    resultQP += A[i][j] * D2UhAtQP[iq][j][i];
	  }
	}
      }
      
      if (grdUhAtQP) {
	resultQP += (*divFct)(A) * grdUhAtQP[iq];
      }
      
      result[iq] += resultQP * factor;
    }
  }

  void VecMatrixGradientAtQP_SOT::weakEval(int nPoints,
					   const WorldVector<double> *grdUhAtQP,
					   WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      WorldMatrix<double> A;
      for (int iq = 0; iq < nPoints; iq++) {
	A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }

  void VecAndCoordsAtQP_SOT::eval(int nPoints,
				  const double *uhAtQP,
				  const WorldVector<double> *grdUhAtQP,
				  const WorldMatrix<double> *D2UhAtQP,
				  double *result,
				  double fac) const
  {
    int dow = Global::getGeo(WORLD);

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

  void VecAndCoordsAtQP_SOT::weakEval(int nPoints,
				      const WorldVector<double> *grdUhAtQP,
				      WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

  void MatrixGradientAndCoords_SOT::eval(int nPoints,
					 const double *uhAtQP,
					 const WorldVector<double> *grdUhAtQP,
					 const WorldMatrix<double> *D2UhAtQP,
					 double *result,
					 double factor) const
  {
    int dow = Global::getGeo(WORLD);

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

      WorldMatrix<double> A = (*f)(gradAtQPs[iq], coordsAtQPs[iq]);

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

      if(grdUhAtQP) {
	resultQP += (*divFct)(A) * grdUhAtQP[iq];
      }

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

  void MatrixGradientAndCoords_SOT::weakEval(int nPoints,
					     const WorldVector<double> *grdUhAtQP,
					     WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      WorldMatrix<double> A;
      for (int iq = 0; iq < nPoints; iq++) {
	A = (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }

  void VecGradCoordsAtQP_SOT::eval(int nPoints,
				   const double *uhAtQP,
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
    int  dow = Global::getGeo(WORLD);

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

  void VecGradCoordsAtQP_SOT::weakEval(int nPoints,
				       const WorldVector<double> *grdUhAtQP,
				       WorldVector<double> *result) const
  {
    if (grdUhAtQP) {  
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

  void VecAtQP_FOT::eval(int nPoints,
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double fac) const
  {
    int dow = Global::getGeo(WORLD);
  
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs[iq]);
	double resultQP = 0.0;
	for (int i = 0; i < dow; i++) {
	  resultQP += grdUhAtQP[iq][i];
	}
	result[iq] += fac * factor * resultQP;
      } 
    }
  }

  void CoordsAtQP_FOT::eval(int nPoints,
			    const double *uhAtQP,
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double f) const
  {
    int dow = Global::getGeo(WORLD);
  
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*g)(coordsAtQPs[iq]);
	double resultQP = 0.0;
	for (int i = 0; i < dow; i++) {
	  resultQP += grdUhAtQP[iq][i];
	}
	result[iq] += f * factor * resultQP;
      }
    }
  }

  void VecCoordsAtQP_FOT::eval(int nPoints,
			       const double *uhAtQP,
			       const WorldVector<double> *grdUhAtQP,
			       const WorldMatrix<double> *D2UhAtQP,
			       double *result,
			       double f) const
  {
    int dow = Global::getGeo(WORLD);
  
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*g)(coordsAtQPs[iq]);
	double resultQP = 0.0;
	for (int i = 0; i < dow; i++) {
	  resultQP += grdUhAtQP[iq][i];
	}
	result[iq] += f * factor * resultQP;
      }
    }
  }

  void VectorGradient_FOT::eval(int nPoints,
				const double *uhAtQP,
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
				double *result,
				double factor) const
  { 
    if (grdUhAtQP) {
      if (f) {
	for (int iq = 0; iq < nPoints; iq++) {
	  WorldVector<double> b = (*f)(gradAtQPs[iq]);
	  result[iq] += b * grdUhAtQP[iq] * factor;
	}
      } else {
	for (int iq = 0; iq < nPoints; iq++) {
	  result[iq] += gradAtQPs[iq] * grdUhAtQP[iq] * factor;
	}
      }
    }
  }

  void VectorFct_FOT::eval(int nPoints,
			   const double *uhAtQP,
			   const WorldVector<double> *grdUhAtQP,
			   const WorldMatrix<double> *D2UhAtQP,
			   double *result,
			   double factor) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	WorldVector<double> b = (*vecFct)(vecAtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }

  void VecFctAtQP_FOT::eval(int nPoints,
			    const double *uhAtQP,
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double f) const
  {
    int dow = Global::getGeo(WORLD);
  
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double resultQP = 0.0;
	const WorldVector<double> &b = (*g)(coordsAtQPs[iq]);
	for (int i = 0; i < dow; i++) {
	  resultQP += b[i] * grdUhAtQP[iq][i];
	}
	result[iq] += f * resultQP;
      }
    }
  }

  Assembler *Operator::getAssembler(int rank) 
  {
    if (!assembler[rank]) {
      initAssembler(rank, NULL, NULL, NULL, NULL);
    }

    return assembler[rank];
  }

  void Operator::setAssembler(int rank, Assembler *a)
  {
    assembler[rank] = a; 
  }

  void CoordsAtQP_IJ_SOT::initElement(const ElInfo* elInfo, 
				      SubAssembler* subAssembler,
				      Quadrature *quad)
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void CoordsAtQP_IJ_SOT::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]), (*g)(coordsAtQPs[iq]));
    }
  }

  void CoordsAtQP_IJ_SOT::eval(int nPoints,
			       const double *uhAtQP,
			       const WorldVector<double> *grdUhAtQP,
			       const WorldMatrix<double> *D2UhAtQP,
			       double *result,
			       double f) const
  {
    if (D2UhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*g)(coordsAtQPs[iq]);
	result[iq] = D2UhAtQP[iq][xi][xj] * factor * f;
      }
    }
  }

  void CoordsAtQP_IJ_SOT::weakEval(int nPoints,
				   const WorldVector<double> *grdUhAtQP,
				   WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*g)(coordsAtQPs[iq]);  
	result[iq][xi] += grdUhAtQP[iq][xj] * factor;
      }
    }
  }

  void VecAtQP_IJ_SOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad)
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
  }
 
  void VecAtQP_IJ_SOT::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]), (*f)(vecAtQPs[iq]));
    }
  }  

  void VecAtQP_IJ_SOT::eval(int nPoints,
			    const double *uhAtQP,
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double fac) const
  {
    if (D2UhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs[iq]);
	result[iq] = D2UhAtQP[iq][xi][xj] * factor * fac;
      }
    }
  }

  void VecAtQP_IJ_SOT::weakEval(int nPoints,
				const WorldVector<double> *grdUhAtQP,
				WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs[iq]);  
	result[iq][xi] += grdUhAtQP[iq][xj] * factor;
      }
    }
  }

 
  void VecOfDOFVecsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					 SubAssembler* subAssembler,
					 Quadrature *quad) 
  {
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++) {
      vecsAtQPs[i] = getVectorAtQPs(vecs[i], elInfo, subAssembler, quad);
    }
  }
 
  void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints,
				  std::vector<double> &C) const 
  { 
    int size = static_cast<int>(vecs.size());
    std::vector<double> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++) {
	arg[i] = vecsAtQPs[i][iq];
      }
      C[iq] += (*f)(arg);
    }
  }

  void VecOfDOFVecsAtQP_ZOT::eval(int nPoints,
				  const double *uhAtQP,
				  const WorldVector<double> *grdUhAtQP,
				  const WorldMatrix<double> *D2UhAtQP,
				  double *result,
				  double fac) const
  {
    int size = static_cast<int>(vecs.size());
    std::vector<double> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++) {
	arg[i] = vecsAtQPs[i][iq];
      }
      result[iq] += fac * (*f)(arg) * uhAtQP[iq];
    }
  }

 
  void VecOfGradientsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					   SubAssembler* subAssembler,
					   Quadrature *quad) 
  {
    int i ,size = static_cast<int>(vecs.size());
    for(i = 0; i < size; i++) {
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
    }
  }
 
  void VecDivergence_ZOT::initElement(const ElInfo* elInfo, 
				      SubAssembler* subAssembler,
				      Quadrature *quad) 
  {
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++) {
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
    }
  }
 

  void VecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
				    std::vector<double> &C) const 
  { 
    int size = static_cast<int>(vecs.size());
    std::vector<WorldVector<double>*> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++) {
	arg[i] = &(gradsAtQPs[i][iq]);
      }
      C[iq] += (*f)(arg);
    }
  }

  void VecOfGradientsAtQP_ZOT::eval(int nPoints,
				    const double *uhAtQP,
				    const WorldVector<double> *grdUhAtQP,
				    const WorldMatrix<double> *D2UhAtQP,
				    double *result,
				    double fac) const
  {
    int size = static_cast<int>(vecs.size());
    std::vector<WorldVector<double>*> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++) {
	arg[i] = &(gradsAtQPs[i][iq]);
      }
      result[iq] += fac * (*f)(arg) * uhAtQP[iq];
    }
  }

  void VecDivergence_ZOT::getC(const ElInfo *, int nPoints,
			       std::vector<double> &C) const 
  { 
    int size = static_cast<int>(vecs.size());

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++) {
	C[iq] += gradsAtQPs[i][iq][i];
      }
    }
  }

  void VecDivergence_ZOT::eval(int nPoints,
			       const double *uhAtQP,
			       const WorldVector<double> *grdUhAtQP,
			       const WorldMatrix<double> *D2UhAtQP,
			       double *result,
			       double fac) const
  {
    int size = static_cast<int>(vecs.size());

    for (int iq = 0; iq < nPoints; iq++) {
      double d = 0.0;
      for (int i = 0; i < size; i++) {
	d += gradsAtQPs[i][iq][i];
      }
      result[iq] += d * uhAtQP[iq] * fac;
    }
  }

  void VecAndGradAtQP_SOT::eval(int nPoints,
				const double *uhAtQP,
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
				double *result,
				double fac) const
  {
    int dow = Global::getGeo(WORLD);

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

  void VecAndGradAtQP_SOT::weakEval(int nPoints,
				    const WorldVector<double> *grdUhAtQP,
				    WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

  void VecAndGradAtQP_SOT::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]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
    }
  }
  
  void VecAndGradAtQP_SOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }
 
  void General_SOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
    }
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
    }
  }

  void General_SOT::getLALt(const ElInfo *elInfo, 
			    int nPoints, 
			    DimMat<double> **LALt) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      lalt(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg), 
	   *(LALt[iq]), symmetric_, 1.0);
    }
  }

  void General_SOT::eval(int nPoints,
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
  {
    int dow = Global::getGeo(WORLD);
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

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

      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }

      WorldMatrix<double> A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);

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

      if (grdUhAtQP) {
	resultQP += (*divFct_)(A) * grdUhAtQP[iq];
      }

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

  void General_SOT::weakEval(int nPoints,
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    if (grdUhAtQP) {
      WorldMatrix<double> A;
      for (int iq = 0; iq < nPoints; iq++) {
	for (int i = 0; i < nVecs; i++) {
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
	for (int i = 0; i < nGrads; i++) {
	  gradsArg[i] = gradsAtQPs_[i][iq];
	}
	A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }

  void General_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
    }
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
    }
  }

  void General_FOT::getLb(const ElInfo *elInfo, 
			  int nPoints, 
			  VectorOfFixVecs<DimVec<double> >& result) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      lb(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg), 
	 result[iq], 1.0);
    }
  }

  void General_FOT::eval(int nPoints,
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
  {
    int dow = Global::getGeo(WORLD);  
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

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

	for (int i = 0; i < nVecs; i++) {
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
	for (int i = 0; i < nGrads; i++) {
	  gradsArg[i] = gradsAtQPs_[i][iq];
	}

	const WorldVector<double> &b = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);

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

  void General_ZOT::initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
				Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
    }
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
    }
  }

  void General_ZOT::getC(const ElInfo *, int nPoints,
			 std::vector<double> &C) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      C[iq] += (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
    }
  }

  void General_ZOT::eval(int nPoints,
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double fac) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      result[iq] += fac * (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg) * uhAtQP[iq];
    }
  }



  //************************************ parametric general Operatorterms *****************
  void GeneralParametric_SOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    elInfo->getElementNormal(elementNormal_);
    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
    }
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
    }
  }

  void GeneralParametric_SOT::getLALt(const ElInfo *elInfo, 
				      int nPoints, 
				      DimMat<double> **LALt) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      lalt(Lambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg), 
	   *(LALt[iq]), symmetric_, 1.0);
    }
  }

  void GeneralParametric_SOT::eval(int nPoints,
				   const double *uhAtQP,
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double factor) const
  {
    int dow = Global::getGeo(WORLD);
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

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

      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }

      WorldMatrix<double> A = (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg);

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

      if (grdUhAtQP) {
	resultQP += (*divFct_)(A) * grdUhAtQP[iq];
      }

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

  void GeneralParametric_SOT::weakEval(int nPoints,
				       const WorldVector<double> *grdUhAtQP,
				       WorldVector<double> *result) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    if (grdUhAtQP) {
      WorldMatrix<double> A;
      for (int iq = 0; iq < nPoints; iq++) {
	for (int i = 0; i < nVecs; i++) {
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
	for (int i = 0; i < nGrads; i++) {
	  gradsArg[i] = gradsAtQPs_[i][iq];
	}
	A = (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }

  void GeneralParametric_FOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    elInfo->getElementNormal(elementNormal_);
    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
    }
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
    }
  }

  void GeneralParametric_FOT::getLb(const ElInfo *elInfo, 
				    int nPoints, 
				    VectorOfFixVecs<DimVec<double> >& result) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      lb(Lambda, (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg), result[iq], 1.0);
    }
  }

  void GeneralParametric_FOT::eval(int nPoints,
				   const double *uhAtQP,
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double factor) const
  {
    int dow = Global::getGeo(WORLD);
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

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

	for (int i = 0; i < nVecs; i++) {
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
	for (int i = 0; i < nGrads; i++) {
	  gradsArg[i] = gradsAtQPs_[i][iq];
	}

	const WorldVector<double> &b = (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg);

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

  void GeneralParametric_ZOT::initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
					  Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    elInfo->getElementNormal(elementNormal_);
    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
    }
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
    }
  }

  void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints, 
				   std::vector<double> &C) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      C[iq] += (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg);
    }
  }

  void GeneralParametric_ZOT::eval(int nPoints,
				   const double *uhAtQP,
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
      for (int i = 0; i < nGrads; i++) {
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      result[iq] += fac * (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg) * uhAtQP[iq];
    }
  }

  void VecAndVecOfGradientsAtQP_ZOT::initElement(const ElInfo* elInfo,
						 SubAssembler* subAssembler,
						 Quadrature *quad)
  {
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++)
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);

    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
  }

  void VecAndVecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
					  std::vector<double> &C) const
  {
    int size = static_cast<int>(vecs.size());
    std::vector<WorldVector<double>*> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++)
	arg[i] = &(gradsAtQPs[i][iq]);

      C[iq] += (*f)(vecAtQPs[iq], arg);
    }
  }

  void VecAndVecOfGradientsAtQP_ZOT::eval(int nPoints,
					  const double *uhAtQP,
					  const WorldVector<double> *grdUhAtQP,
					  const WorldMatrix<double> *D2UhAtQP,
					  double *result,
					  double fac) const
  {
    int size = static_cast<int>(vecs.size());
    std::vector<WorldVector<double>*> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++)
	arg[i] = &(gradsAtQPs[i][iq]);

      result[iq] += fac * (*f)(vecAtQPs[iq], arg) * uhAtQP[iq];
    }
  }

  void VecGrad_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
    vecAtQPs = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec2, elInfo, subAssembler, quad);
  }

  void VecGrad_FOT::initElement(const ElInfo* smallElInfo,
				const ElInfo* largeElInfo,
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
    vecAtQPs = getVectorAtQPs(vec1, smallElInfo, largeElInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec2, smallElInfo, largeElInfo, subAssembler, quad);
  }

  void VecGrad_FOT::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++)
      lb(Lambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
  }
  
  void VecGrad_FOT::eval(int nPoints,
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < nPoints; iq++) {
	WorldVector<double> b = (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }

  void VecGrad_SOT::eval(int numPoints,
				const double              *uhAtQP,
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
				double *result,
				double fac) const
  {
    int dow = Global::getGeo(WORLD);

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

  void VecGrad_SOT::weakEval(int numPoints,
				    const WorldVector<double> *grdUhAtQP,
				    WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < numPoints; iq++) {
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

  void VecGrad_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    for (int iq = 0; iq < numPoints; iq++) {
      l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
    }
  }
  
  void VecGrad_SOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
    vecAtQPs = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
    gradAtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
  }

  void FctGrad2_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
    grad1AtQPs = subAssembler->getGradientsAtQPs(vec1, elInfo, quad);
    grad2AtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
  }

  void FctGrad2_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    for (int iq = 0; iq < numPoints; iq++) {
      lb(Lambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
    }
  }
  
  void FctGrad2_FOT::eval(int numPoints,
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
  {
    if (grdUhAtQP) {
      for (int iq = 0; iq < numPoints; iq++) {
	WorldVector<double> b = (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }

}