Simons changes.

AMDiS/src/ElementFunction.h

@@ -63,7 +63,7 @@ namespace AMDiS {
     {}
 
     /// evaluation at given coordinates.
-    const T& operator()(const DimVec<double>& bary) const {
+    T operator()(const DimVec<double>& bary) const {
       WorldVector<double> worldCoords;
       this->elInfo_->coordToWorld(bary, &worldCoords);
       return (*fct_)(worldCoords);
@@ -87,7 +87,7 @@ namespace AMDiS {
     {}
 
     /// evaluation at given coordinates.
-    const T& operator()(const DimVec<double>& bary) const 
+    T operator()(const DimVec<double>& bary) const 
     {
       static T t;
       const T* localVec = 

AMDiS/src/Operator.cc

@@ -2982,4 +2982,82 @@ namespace AMDiS {
     }
   }
 
+  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;
+      }
+    }
+  }
+
 }

AMDiS/src/Operator.h

@@ -1870,6 +1870,56 @@ namespace AMDiS {
     BinaryAbstractFunction<WorldVector<double>, double, WorldVector<double> > *vecFct;
   };
 
+
+  /**
+   * \ingroup Assembler
+   *
+   * \brief
+   * First order term: \f$ b(\nabla v(\vec{x}), \nabla w(\vec{x})) \cdot \nabla u(\vec{x})\f$.
+   */
+  class FctGrad2_FOT : public FirstOrderTerm
+  {
+  public:
+    /// Constructor
+    FctGrad2_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+		BinaryAbstractFunction<WorldVector<double>, WorldVector<double>, WorldVector<double> > *vecFct_)
+      : FirstOrderTerm(vecFct_->getDegree()), vec1(dv1), vec2(dv2), vecFct(vecFct_) 
+    {}
+
+    /// Implementation of \ref OperatorTerm::initElement().
+    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+		     Quadrature *quad = NULL);
+    
+    /// Implements FirstOrderTerm::getLb().
+    void getLb(const ElInfo *elInfo, int qPoint, 
+	       VectorOfFixVecs<DimVec<double> >& Lb) const; 
+    
+    /// Implements FirstOrderTerm::eval().
+    void eval(int numPoints,
+	      const double *uhAtQP,
+	      const WorldVector<double> *grdUhAtQP,
+	      const WorldMatrix<double> *D2UhAtQP,
+	      double *result,
+	      double factor) const;
+    
+  protected:
+    /// DOFVector to be evaluated at quadrature points.
+    DOFVectorBase<double>* vec1;
+
+    /// DOFVector to be evaluated at quadrature points.
+    DOFVectorBase<double>* vec2;
+    
+    /// Gradient v at quadrature points.
+    WorldVector<double> *grad1AtQPs;
+
+    /// Gradient v at quadrature points.
+    WorldVector<double> *grad2AtQPs;
+
+    /// Function for b.
+    BinaryAbstractFunction<WorldVector<double>, WorldVector<double>, WorldVector<double> > *vecFct;
+  };
+
+
   class General_FOT : public FirstOrderTerm
   {
   public:
@@ -2657,6 +2707,59 @@ namespace AMDiS {
     BinaryAbstractFunction<double, double, WorldVector<double> > *f;
   };
 
+
+  /**
+   * \ingroup Assembler
+   *
+   * \brief
+   * Second order term: \f$ f(v(\vec{x}),\nabla w(\vec{x})) \Delta u(\vec{x})\f$
+   */
+  class VecGrad_SOT : public SecondOrderTerm
+  {
+  public:
+    /// Constructor.
+    VecGrad_SOT( DOFVectorBase<double> *dv1,
+			  DOFVectorBase<double> *dv2,
+		       BinaryAbstractFunction<double, double, WorldVector<double> > *f_)
+      : SecondOrderTerm(f_->getDegree()), vec1(dv1), vec2(dv2), f(f_) 
+    {}
+
+    /// Implementation of \ref OperatorTerm::initElement().
+    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+		     Quadrature *quad = NULL);
+
+    /// Implements SecondOrderTerm::getLALt().
+    inline void getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const;
+
+    /// Implements SecondOrderTerm::eval().
+    void eval(int numPoints,
+	      const double              *uhAtQP,
+	      const WorldVector<double> *grdUhAtQP,
+	      const WorldMatrix<double> *D2UhAtQP,
+	      double *result,
+	      double factor) const;
+
+    /// Implements SecondOrderTerm::weakEval().
+    void weakEval(int numPoints,
+		  const WorldVector<double> *grdUhAtQP,
+		  WorldVector<double> *result) const;
+  
+  protected:
+    /// DOFVectors to be evaluated at quadrature points.
+    DOFVectorBase<double>* vec1;
+    DOFVectorBase<double>* vec2;
+
+    /// Vector v at quadrature points.
+    double *vecAtQPs;
+
+    /// Gradient at quadrature points.
+    WorldVector<double> *gradAtQPs;
+
+    /// Function for c.
+    BinaryAbstractFunction<double, double, WorldVector<double> > *f;
+  };
+
+
   /* 
    * \ingroup Assembler
    *