GenericOperatorTerm added - provides a math-like syntax for adding terms

AMDiS/src/Functors.h

@@ -123,14 +123,14 @@ template<typename T>
 struct Diff : public BinaryAbstractFunction<T,T,T>
 {
   Diff(int degree = 1) : BinaryAbstractFunction<T,T,T>(degree) {}
-  T operator()(const T &v1, const T &v2) const { return abs(v1-v2); }
+  T operator()(const T &v1, const T &v2) const { return std::abs(v1-v2); }
 };
 
 template<typename T=double>
 struct Abs : public AbstractFunction<T,T>
 {
   Abs(int degree = 1) : AbstractFunction<T,T>(degree) {}
-  T operator()(const T &v) const { return abs(v); }
+  T operator()(const T &v) const { return std::abs(v); }
 };
 
 template<typename T=double>
@@ -151,35 +151,36 @@ template<typename T=double>
 struct Sqrt : public AbstractFunction<T,T>
 {
   Sqrt(int degree = 4) : AbstractFunction<T,T>(degree) {}
-  T operator()(const T &v) const { return sqrt(v); }
+  T operator()(const T &v) const { return std::sqrt(v); }
 };
 
 namespace detail {
-  template<int p>
+  template<int p, typename T>
   struct Pow
   {
-    template<typename T>
-    static T eval(const T& v) { return v*Pow<p-1>::eval(v); }
+    typedef typename AMDiS::ProductType<T, typename Pow<p-1,T>::result_type>::type result_type;
+    static result_type eval(const T& v) { return v*Pow<p-1,T>::eval(v); }
   };
   
-  template<>
-  struct Pow<1> {
-    template<typename T>
-    static T eval(const T& v) { return v; }
+  template<typename T>
+  struct Pow<1,T> {
+    typedef T result_type;
+    static result_type eval(const T& v) { return v; }
   };
   
-  template<>
-  struct Pow<0> {
-    template<typename T>
-    static T eval(const T& v) { return 1.0; }
+  template<typename T>
+  struct Pow<0, T> {
+    typedef double result_type;
+    static result_type eval(const T& v) { return 1.0; }
   };
 }
 
 template<int p, typename T=double>
-struct Pow : public AbstractFunction<T,T>
+struct Pow : public AbstractFunction<typename detail::Pow<p,T>::result_type, T>
 {
-  Pow(double factor_=1.0, int degree = p) : AbstractFunction<T,T>(degree), factor(factor_) {}
-  T operator()(const T &v) const { return factor * detail::Pow<p>::eval(v); }
+  typedef typename detail::Pow<p,T>::result_type result_type;
+  Pow(double factor_=1.0, int degree = p) : AbstractFunction<result_type,T>(degree), factor(factor_) {}
+  result_type operator()(const T &v) const { return factor * detail::Pow<p,T>::eval(v); }
 private:
   double factor;
 };
@@ -188,7 +189,7 @@ template<typename T1, typename T2 = ProductType<T1, T1> >
 struct Norm2 : public AbstractFunction<T1, T2>
 {
   Norm2(int degree = 4) : AbstractFunction<T1, T2>(degree) {}
-  T1 operator()(const T2 &v) const { return sqrt(v*v); }
+  T1 operator()(const T2 &v) const { return std::sqrt(v*v); }
 };
 
 template<typename T1, typename T2>
@@ -202,7 +203,7 @@ template<typename T>
 struct Norm2_comp2 : public BinaryAbstractFunction<T,T,T>
 {
   Norm2_comp2(int degree = 4) : BinaryAbstractFunction<T,T,T>(degree) {}
-  T operator()(const T &v1, const T &v2) const { return sqrt(sqr(v1)+sqr(v2)); }
+  T operator()(const T &v1, const T &v2) const { return std::sqrt(sqr(v1)+sqr(v2)); }
 };
 
 template<typename T>
@@ -216,7 +217,7 @@ template<typename T>
 struct Norm2_comp3 : public TertiaryAbstractFunction<T,T,T,T>
 {
   Norm2_comp3(int degree = 4) : TertiaryAbstractFunction<T,T,T,T>(degree) {}
-  T operator()(const T &v1, const T &v2, const T &v3) const { return sqrt(sqr(v1)+sqr(v2)+sqr(v3)); }
+  T operator()(const T &v1, const T &v2, const T &v3) const { return std::sqrt(sqr(v1)+sqr(v2)+sqr(v3)); }
 };
 
 template<typename T>
@@ -229,7 +230,7 @@ struct Norm2Sqr_comp3 : public TertiaryAbstractFunction<T,T,T,T>
 template<typename T>
 struct L1Diff : public BinaryAbstractFunction<T,T,T>
 {
-  T operator()(const T &v1, const T &v2) const { return abs(v1-v2); }
+  T operator()(const T &v1, const T &v2) const { return std::abs(v1-v2); }
 };
 
 template<typename TOut, typename T=TOut>

AMDiS/src/GenericOperatorTerm.h

+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+
+
+/** \file GenericOperatorTerm.h */
+
+#ifndef AMDIS_GENERIC_OPERATOR_TERM_H
+#define AMDIS_GENERIC_OPERATOR_TERM_H
+
+#include "AMDiS_fwd.h"
+#include "OperatorTerm.h"
+#include "ValueTypes.h"
+#include "Functors.h"
+
+#include <boost/static_assert.hpp>
+
+
+#include "expressions/LazyOperatorTerm.h"
+#include "expressions/AtomicExpression.h"
+#include "expressions/LazyExpression.h"
+
+/** generic operator-terms provide an easy way of automated generation of
+ * 'arbitrary' operator-terms out of some elementary operations, by using a 
+ * recursive definition of the term. All necessary data will be initialized 
+ * when an expression as part of the term uses this data.
+ * Since no virtual functions, like in the AbstractFunction classes, are used
+ * the overhead of a vtable is removed.
+ * 
+ * usage:
+ * addZOT(Operator, Term) 
+ *  ... add a zeroOrderTerm to Operator, i.e. (Term(x) * u, v)
+ * 
+ * addFOT(Operator, Term, FirstOrderType) 
+ *  ... add a firstOrderTerm to Operator, (if Term::value_type = double) 
+ *      i.e. (Term(x) * 1 * grad(u), v), rsp. (Term(x) * 1 * u, grad(v))
+ * addFOT(Operator, Term, FirstOrderType)
+ *  ... add a firstOrderTerm to Operator, (if Term::value_type = WorldVector)
+ *      i.e. (Term(x) * b * grad(u), v), rsp. (Term(x) * u, grad(v))
+ * addFOT<I>(Operator, Term, FirstOrderType)
+ *  ... add a firstOrderTerm to Operator, 
+ *      i.e. (Term(x) * e_I * grad(u), v), rsp. (Term(x) * e_I * u, grad(v))
+ * 
+ * addSOT(Operator, Term)
+ *  ... add a secondOrderTerm to Operator, i.e. (Term(x) * grad(u), grad(v))
+ * addSOT<I,J>(Operator, Term)
+ *  ... add a secondOrderTerm to Operator, i.e. (E_IJ * Term(x) * grad(u), grad(v))
+ * 
+ * where Operator is eather a pointer or reference, FirstOrderType in {GRD_PHI, GRD_PSI}
+ * and Term a componation of elementary terms by + - * /
+ *  - constant(value) / value  ... a constant value
+ *  - valueOf(DOFVector)       ... values of a DOFVector at QP
+ *  - gradientOf(DOFVector)    ... gradient of a DOFVector at QP
+ *  - derivative<I>(DOFVector) ... I'th partial derivative
+ *  - X()                      ... coordinate at quadrature points
+ *  - pow<I>(Term)             ... I'th power of a term
+ *  - sqrt(Term)               ... square root of a term
+ *  - Exp(Term)                ... exponential function of a term
+ *  - function_<F>(Term)       ... evaluates F()(Term(iq))
+ *  - function_(F f, Term)     ... evaluates f(Term(iq))
+ * 
+ * 
+ * with F a functor that implements
+ *   typedef (...) value_type;
+ *   int getDegree(int d0);
+ *   value_type operator()(const T0& v0) const;
+ * 
+ * respective
+ *   int getDegree(int d0, int d1);
+ *   value_type operator()(const T0& v0, const T1& v1) const;
+ * 
+ * respective
+ *   int getDegree(int d0, int d1, int d2);
+ *   value_type operator()(const T0& v0, const T1& v1, const T2& v2) const;
+ * 
+ * where the d0, d1, d2 give the polynomial degrees of the v0, v1, v2 terms.
+ * */
+
+namespace AMDiS {
+
+// _______ ZeroOrderTerms ______________________________________________________
+
+template<typename Term>
+struct GenericZeroOrderTerm : public ZeroOrderTerm
+{
+  Term term;
+  GenericZeroOrderTerm(const Term& term_)
+    : ZeroOrderTerm(term_.getDegree()), term(term_) 
+  {
+    term.insertFeSpaces(auxFeSpaces);
+  }
+
+  void initElement(const ElInfo* elInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, elInfo, subAssembler, quad);
+  }
+
+
+  void initElement(const ElInfo* smallElInfo,
+	    const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+  }
+
+
+  void getC(const ElInfo *elInfo, int nPoints, ElementVector& C)
+  {
+    for (int iq = 0; iq < nPoints; iq++)
+      C[iq] += term(iq);
+  }
+
+
+  void eval( int nPoints,
+	      const mtl::dense_vector<double>& uhAtQP,
+	      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+	      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+	      mtl::dense_vector<double>& result,
+	      double fac)
+  {
+    for (int iq = 0; iq < nPoints; iq++)
+      result[iq] += fac * term(iq) * uhAtQP[iq];
+  }
+};
+
+// _______ FirstOrderTerms _____________________________________________________
+
+
+template<typename Term>
+struct GenericFirstOrderTerm_1 : public FirstOrderTerm
+{
+  Term term;
+  GenericFirstOrderTerm_1(const Term& term_)
+    : FirstOrderTerm(term_.getDegree()), term(term_) 
+  {
+    term.insertFeSpaces(auxFeSpaces);
+  }
+
+  void initElement(const ElInfo* elInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, elInfo, subAssembler, quad);
+  }
+
+
+  void initElement(const ElInfo* smallElInfo,
+	    const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+  }
+
+
+  /// Implements FirstOrderTerm::getLb().
+  void getLb(const ElInfo *elInfo,
+	     std::vector<mtl::dense_vector<double> >& Lb) const
+  {
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+    const int nPoints = static_cast<int>(Lb.size());
+
+    for (int iq = 0; iq < nPoints; iq++)
+      l1(grdLambda, Lb[iq], term(iq));
+  }
+};
+
+
+template<int I, typename Term>
+struct GenericFirstOrderTerm_i : public FirstOrderTerm
+{
+  Term term;
+  GenericFirstOrderTerm_i(const Term& term_)
+    : FirstOrderTerm(term_.getDegree()), term(term_) 
+  {
+    term.insertFeSpaces(auxFeSpaces);
+    FirstOrderTerm::bOne = I;
+  }
+  
+  GenericFirstOrderTerm_i(const Term& term_, int I0)
+    : FirstOrderTerm(term_.getDegree()), term(term_) 
+  {
+    term.insertFeSpaces(auxFeSpaces);
+    FirstOrderTerm::bOne = I0;
+    TEST_EXIT_DBG( I < 0 && I0 >= 0 )("You yould specify eather template<int I>, or constructor(int I0)\n");
+  }
+
+  void initElement(const ElInfo* elInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, elInfo, subAssembler, quad);
+  }
+
+
+  void initElement(const ElInfo* smallElInfo,
+	    const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+  }
+
+
+  /// Implements FirstOrderTerm::getLb().
+  void getLb(const ElInfo *elInfo,
+	     std::vector<mtl::dense_vector<double> >& Lb) const
+  {
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+    const int nPoints = static_cast<int>(Lb.size());
+
+    for (int iq = 0; iq < nPoints; iq++)
+      lb_one(grdLambda, Lb[iq], term(iq));
+  }
+};
+
+
+template<typename Term>
+struct GenericFirstOrderTerm_b : public FirstOrderTerm
+{
+  Term term;
+  
+  GenericFirstOrderTerm_b(const Term& term_)
+    : FirstOrderTerm(term_.getDegree()), term(term_)
+  {
+    term.insertFeSpaces(auxFeSpaces);
+  }
+
+  void initElement(const ElInfo* elInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, elInfo, subAssembler, quad);
+  }
+
+
+  void initElement(const ElInfo* smallElInfo,
+	    const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+  }
+
+
+  /// Implements FirstOrderTerm::getLb().
+  void getLb(const ElInfo *elInfo,
+	     std::vector<mtl::dense_vector<double> >& Lb) const
+  {
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+    const int nPoints = static_cast<int>(Lb.size());
+
+    for (int iq = 0; iq < nPoints; iq++)
+      lb(grdLambda, term(iq), Lb[iq], fac);
+  }
+};
+
+// _______ SecondOrderTerms ____________________________________________________
+
+
+template<typename Term>
+struct GenericSecondOrderTerm_1 : public SecondOrderTerm
+{
+  Term term;
+  
+  GenericSecondOrderTerm_1(const Term& term_)
+    : SecondOrderTerm(term_.getDegree()), term(term_) 
+  {
+    term.insertFeSpaces(auxFeSpaces);
+    setSymmetric(true);
+  }
+
+  void initElement(const ElInfo* elInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, elInfo, subAssembler, quad);
+  }
+
+
+  void initElement(const ElInfo* smallElInfo,
+	    const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+  }
+  
+  /// Implements SecondOrderTerm::getLALt().
+  void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const
+  {
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+    const int nPoints = static_cast<int>(LALt.size());
+
+    for (int iq = 0; iq < nPoints; iq++) 
+      l1lt(grdLambda, LALt[iq], term(iq));
+  }
+
+  /// Implenetation of SecondOrderTerm::eval().
+  void eval(int nPoints,
+	    const mtl::dense_vector<double>& uhAtQP,
+	    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+	    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+	    mtl::dense_vector<double>& result,
+	    double f) 
+  {
+    int dow = Global::getGeo(WORLD);
+
+    if (num_rows(D2UhAtQP) > 0) {
+      for (int iq = 0; iq < nPoints; iq++) {
+	double resultQP = 0.0;
+	for (int i = 0; i < dow; i++) {
+	  resultQP += D2UhAtQP[iq][i][i];
+	}
+	result[iq] += resultQP * f * term(iq);
+      }
+    }
+  }
+
+  /// Implenetation of SecondOrderTerm::weakEval().
+  void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+		std::vector<WorldVector<double> > &result) 
+  {
+    int nPoints = grdUhAtQP.size();
+    for (int iq = 0; iq < nPoints; iq++)
+      axpy(term(iq), grdUhAtQP[iq], result[iq]);
+  }
+};
+
+
+
+template<typename Term, bool symmetric = true>
+struct GenericSecondOrderTerm_A : public SecondOrderTerm
+{
+  Term term;
+  
+  GenericSecondOrderTerm_A(const Term& term_)
+    : SecondOrderTerm(term_.getDegree()), term(term_) 
+  {
+    term.insertFeSpaces(auxFeSpaces);
+    setSymmetric(symmetric);
+  }
+
+  void initElement(const ElInfo* elInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, elInfo, subAssembler, quad);
+  }
+
+
+  void initElement(const ElInfo* smallElInfo,
+	    const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+  }
+  
+  void getLALt(const ElInfo *elInfo, 
+	       std::vector<mtl::dense2D<double> > &LALt) const
+  {
+    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+    const int nPoints = static_cast<int>(LALt.size());
+
+    for (int iq = 0; iq < nPoints; iq++)
+      lalt(grdLambda, term(iq), LALt[iq], symmetric, 1.0);
+  }
+
+  void eval(int nPoints,
+	    const mtl::dense_vector<double>& uhAtQP,
+	    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+	    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+	    mtl::dense_vector<double>& result,
+	    double factor) 
+  {
+    int dow = Global::getGeo(WORLD);
+
+    for (int iq = 0; iq < nPoints; iq++) {
+      double resultQP = 0.0;
+
+      WorldMatrix<double> A = term(iq);
+
+      if (num_rows(D2UhAtQP) > 0)
+	for (int i = 0; i < dow; i++)
+	  for (int j = 0; j < dow; j++)
+	    resultQP += A[i][j] * D2UhAtQP[iq][j][i];
+
+//       if (num_rows(grdUhAtQP) > 0)
+// 	resultQP += (*divFct)(A) * grdUhAtQP[iq];
+
+      result[iq] += resultQP * factor;
+    }
+  }
+
+  void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+		std::vector<WorldVector<double> > &result)
+  {
+    int nPoints = grdUhAtQP.size();
+    WorldMatrix<double> A;
+    for (int iq = 0; iq < nPoints; iq++) {
+      A = term(iq);
+      result[iq] += term(iq) * grdUhAtQP[iq];
+    }    
+  }
+};
+
+
+
+template<int I, int J, typename Term>
+struct GenericSecondOrderTerm_ij : public SecondOrderTerm
+{
+  Term term;
+  int row, col;
+  
+  GenericSecondOrderTerm_ij(const Term& term_)
+    : SecondOrderTerm(term_.getDegree()), term(term_), row(I), col(J)
+  {
+    term.insertFeSpaces(auxFeSpaces);
+    setSymmetric(row == col);
+  }
+  
+  GenericSecondOrderTerm_ij(const Term& term_, int I0, int J0)
+    : SecondOrderTerm(term_.getDegree()), term(term_), row(I0), col(J0)
+  {
+    term.insertFeSpaces(auxFeSpaces);
+    setSymmetric(row == col);
+    TEST_EXIT_DBG( I < 0 && I0 >= 0 && J < 0 && J0 >= 0 ) 
+      ("You yould specify eather template<int I, int J>, or constructor(int I0, int J0)\n");
+  }
+
+  void initElement(const ElInfo* elInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)
+  {
+    term.initElement(*this, elInfo, subAssembler, quad);
+  }
+
+
+  void initElement(const ElInfo* smallElInfo,
+	    const ElInfo* largeElInfo,
+	    SubAssembler* subAssembler,
+	    Quadrature *quad)