diff --git a/AMDiS/src/Functors.h b/AMDiS/src/Functors.h
index 64a2e7b66ab82dfb8f5e1b52d7892ecd81de18d6..815c6b03b3751f113d53c8b29c74514a1650090e 100644
--- a/AMDiS/src/Functors.h
+++ b/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>
diff --git a/AMDiS/src/GenericOperatorTerm.h b/AMDiS/src/GenericOperatorTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..5f19af70c36a8ffa8a85538850b600f0f949ad4e
--- /dev/null
+++ b/AMDiS/src/GenericOperatorTerm.h
@@ -0,0 +1,666 @@
+/******************************************************************************
+ *
+ * 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)
+ {
+ term.initElement(*this, smallElInfo, largeElInfo, subAssembler, quad);
+ }
+
+ void getLALt(const ElInfo *elInfo,
+ std::vector<mtl::dense2D<double> > &LALt) const
+ {
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
+ for (int iq = 0; iq < nPoints; iq++)
+ lalt_kl(grdLambda, row, col, LALt[iq], term(iq));
+ }
+
+ void eval(int nPoints,
+ const mtl::dense_vector<double>& uhAtQP,
+ const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+ const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+ mtl::dense_vector<double>& result,
+ double fac)
+ {
+ if (num_rows(D2UhAtQP) > 0) {
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += D2UhAtQP[iq][row][col] * term(iq) * fac;
+ }
+ }
+
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result)
+ {
+ int nPoints = grdUhAtQP.size();
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq][row] += grdUhAtQP[iq][col] * term(iq);
+ }
+};
+
+
+// _____________________________________________________________________________
+
+template <typename Term>
+void addZOT(Operator* op, const Term& term)
+{
+ op->addZeroOrderTerm(new GenericZeroOrderTerm<Term>(term));
+}
+
+template <typename Term>
+void addZOT(Operator& op, const Term& term)
+{
+ op.addZeroOrderTerm(new GenericZeroOrderTerm<Term>(term));
+}
+
+// .............................................................................
+// Rosenbrock method
+
+// template <typename Term>
+// void addZOT_(Operator* op, const Term& term)
+// {
+// if (op->getUhOld()) {
+// addZOT(op, jacobian(term,0) * ...
+//
+// }
+
+// template <typename Term>
+// void addZOT(Operator& op, const Term& term)
+// {
+// op.addZeroOrderTerm(new GenericZeroOrderTerm<Term>(term));
+// }
+
+
+// _____________________________________________________________________________
+
+// first order term using FirstOrderTerm::l1
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
+addFOT(Operator* op, const Term& term, FirstOrderType type)
+{
+ op->addFirstOrderTerm(new GenericFirstOrderTerm_1<Term>(term), type);
+}
+
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
+addFOT(Operator& op, const Term& term, FirstOrderType type)
+{
+ op.addFirstOrderTerm(new GenericFirstOrderTerm_1<Term>(term), type);
+}
+
+// first order term using FirstOrderTerm::lb_one
+template <int I, typename Term>
+inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
+addFOT(Operator* op, const Term& term, FirstOrderType type)
+{
+ op->addFirstOrderTerm(new GenericFirstOrderTerm_i<I,Term>(term), type);
+}
+
+template <int I, typename Term>
+inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
+addFOT(Operator& op, const Term& term, FirstOrderType type)
+{
+ op.addFirstOrderTerm(new GenericFirstOrderTerm_i<I,Term>(term), type);
+}
+
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
+addFOT(Operator* op, const Term& term, int I, FirstOrderType type)
+{
+ op->addFirstOrderTerm(new GenericFirstOrderTerm_i<-1,Term>(term,I), type);
+}
+
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
+addFOT(Operator& op, const Term& term, int I, FirstOrderType type)
+{
+ op.addFirstOrderTerm(new GenericFirstOrderTerm_i<-1,Term>(term,I), type);
+}
+
+// first order term using FirstOrderTerm::lb
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<WorldVector<double>, typename Term::value_type>::type >::type
+addFOT(Operator* op, const Term& term, FirstOrderType type)
+{
+ op->addFirstOrderTerm(new GenericFirstOrderTerm_b<Term>(term), type);
+}
+
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<WorldVector<double>, typename Term::value_type>::type >::type
+addFOT(Operator& op, const Term& term, FirstOrderType type)
+{
+ op.addFirstOrderTerm(new GenericFirstOrderTerm_b<Term>(term), type);
+}
+
+// _____________________________________________________________________________
+
+// second order term using matrix functions
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<WorldMatrix<double>, typename Term::value_type>::type >::type
+addSOT(Operator* op, const Term& term)
+{
+ op->addSecondOrderTerm(new GenericSecondOrderTerm_A<Term>(term));
+}
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<WorldMatrix<double>, typename Term::value_type>::type >::type
+addSOT(Operator& op, const Term& term)
+{
+ op.addSecondOrderTerm(new GenericSecondOrderTerm_A<Term>(term));
+}
+
+template <bool Symmetric, typename Term>
+inline typename boost::enable_if< typename boost::is_same<WorldMatrix<double>, typename Term::value_type>::type >::type
+addSOT(Operator* op, const Term& term)
+{
+ op->addSecondOrderTerm(new GenericSecondOrderTerm_A<Term, Symmetric>(term));
+}
+
+template <bool Symmetric, typename Term>
+inline typename boost::enable_if< typename boost::is_same<WorldMatrix<double>, typename Term::value_type>::type >::type
+addSOT(Operator& op, const Term& term)
+{
+ op.addSecondOrderTerm(new GenericSecondOrderTerm_A<Term, Symmetric>(term));
+}
+
+// second order term using scalar functions with identity matrix
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
+addSOT(Operator* op, const Term& term)
+{
+ op->addSecondOrderTerm(new GenericSecondOrderTerm_1<Term>(term));
+}
+
+template <typename Term>
+inline typename boost::enable_if< typename boost::is_same<double, typename Term::value_type>::type >::type
+addSOT(Operator& op, const Term& term)
+{
+ op.addSecondOrderTerm(new GenericSecondOrderTerm_1<Term>(term));
+}
+
+// second order term using matrix=0 with matrix_ij = function value
+template <int I, int J, typename Term>
+void addSOT(Operator* op, const Term& term)
+{
+ op->addSecondOrderTerm(new GenericSecondOrderTerm_ij<I,J,Term>(term));
+}
+
+template <int I, int J, typename Term>
+void addSOT(Operator& op, const Term& term)
+{
+ op.addSecondOrderTerm(new GenericSecondOrderTerm_ij<I,J,Term>(term));
+}
+
+template <typename Term>
+void addSOT(Operator* op, const Term& term, int I, int J)
+{
+ op->addSecondOrderTerm(new GenericSecondOrderTerm_ij<-1,-1,Term>(term,I,J));
+}
+
+template <typename Term>
+void addSOT(Operator& op, const Term& term, int I, int J)
+{
+ op.addSecondOrderTerm(new GenericSecondOrderTerm_ij<-1,-1,Term>(term,I,J));
+}
+
+}
+
+#endif // AMDIS_GENERIC_OPERATOR_TERM_H
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index 8880653851446587980e1cefc480feb9c7a21a13..85ba5974b53cf81f47706e893fe749e6fe8b0282 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -71,7 +71,9 @@ namespace AMDiS {
deserialized(false),
computeExactError(false),
boundaryConditionSet(false),
- writeAsmInfo(false)
+ writeAsmInfo(false),
+ solutionTime(0.0),
+ buildTime(0.0)
{
Parameters::get(name + "->components", nComponents);
TEST_EXIT(nComponents > 0)("No value set for parameter \"%s->components\"!\n",
@@ -666,6 +668,7 @@ namespace AMDiS {
INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
t.elapsed());
+ solutionTime = t.elapsed();
adaptInfo->setSolverIterations(solver->getIterations());
adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
@@ -924,6 +927,7 @@ namespace AMDiS {
MPI::COMM_WORLD.Barrier();
#endif
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", t.elapsed());
+ buildTime = t.elapsed();
}
diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h
index 8ae9b1dddd5f55f4ac2deedea3a604dbeeb961aa..055bb53ff6635907d65d54c223f36370f4262b57 100644
--- a/AMDiS/src/ProblemStat.h
+++ b/AMDiS/src/ProblemStat.h
@@ -547,6 +547,16 @@ namespace AMDiS {
{
return fileWriters;
}
+
+ double getSolutionTime()
+ {
+ return solutionTime;
+ }
+
+ double getBuildTime()
+ {
+ return buildTime;
+ }
protected:
/// If the exact solution is known, the problem can compute the exact
@@ -664,6 +674,9 @@ namespace AMDiS {
bool writeAsmInfo;
std::map<Operator*, std::vector<OperatorPos> > operators;
+
+ double solutionTime;
+ double buildTime;
};
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
diff --git a/AMDiS/src/Recovery.h b/AMDiS/src/Recovery.h
index 81b1ee77ab76b8138b9ce307dc44e9b904d2e509..f8375a14d0e5e81963e7bdf611e559bde1c2f17e 100644
--- a/AMDiS/src/Recovery.h
+++ b/AMDiS/src/Recovery.h
@@ -51,9 +51,9 @@ namespace AMDiS {
double operator()(const WorldVector<double>& y,
const WorldVector<double>& z) const
{
- double result = pow(y[0] - z[0], exponent[0]);
+ double result = std::pow(y[0] - z[0], double(exponent[0]));
for (int i = 1; i < exponent.size(); i++)
- result *= pow(y[i] - z[i], exponent[i]);
+ result *= std::pow(y[i] - z[i], double(exponent[i]));
return result;
}
diff --git a/AMDiS/src/ZeroOrderTerm.h b/AMDiS/src/ZeroOrderTerm.h
index c7c8adec2e59ab4b89add07fa9bf54aef0b1f1cf..5c0dc44632e62ab1030972839218b6c835a84bbf 100644
--- a/AMDiS/src/ZeroOrderTerm.h
+++ b/AMDiS/src/ZeroOrderTerm.h
@@ -99,20 +99,20 @@ namespace AMDiS {
public:
/// Constructor.
VecAtQP_ZOT(DOFVectorBase<double> *dv,
- AbstractFunction<double, double> *ab = nullptr,
+ AbstractFunction<double, double> *ab = NULL,
double factor_ = 1.0
);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implementation of \ref OperatorTerm::initElement() for multiple meshes.
void initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -160,7 +160,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -203,19 +203,19 @@ namespace AMDiS {
/// Constructor.
Vec2AtQP_ZOT(DOFVectorBase<double> *dv1,
DOFVectorBase<double> *dv2,
- BinaryAbstractFunction<double, double, double> *f = nullptr,
+ BinaryAbstractFunction<double, double, double> *f = NULL,
double factor_ = 1.0
);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
void initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -260,13 +260,13 @@ namespace AMDiS {
Vec3AtQP_ZOT(DOFVectorBase<double> *dv1,
DOFVectorBase<double> *dv2,
DOFVectorBase<double> *dv3,
- TertiaryAbstractFunction<double, double, double, double> *f = nullptr,
+ TertiaryAbstractFunction<double, double, double, double> *f = NULL,
double factor_ = 1.0
);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -309,7 +309,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *elInfo, int nPoints, ElementVector& C);
@@ -349,7 +349,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -392,7 +392,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -433,7 +433,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -477,7 +477,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *elInfo, int nPoints, ElementVector& C);
@@ -519,7 +519,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -557,7 +557,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -599,7 +599,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -645,7 +645,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -679,7 +679,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -710,12 +710,12 @@ namespace AMDiS {
/// Constructor.
VecDivergence_ZOT(int nComponents,
DOFVectorBase<double> *vec0,
- DOFVectorBase<double> *vec1 = nullptr,
- DOFVectorBase<double> *vec2 = nullptr);
+ DOFVectorBase<double> *vec1 = NULL,
+ DOFVectorBase<double> *vec2 = NULL);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -747,7 +747,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -786,7 +786,7 @@ namespace AMDiS {
QuartAbstractFunction<double, double, double, WorldVector<double>, WorldVector<double> > *af);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -816,7 +816,7 @@ namespace AMDiS {
TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f);
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -850,7 +850,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -899,7 +899,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
@@ -950,7 +950,7 @@ namespace AMDiS {
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
- Quadrature *quad = nullptr);
+ Quadrature *quad = NULL);
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, ElementVector& C);
diff --git a/AMDiS/src/expressions/AtomicExpression.h b/AMDiS/src/expressions/AtomicExpression.h
new file mode 100644
index 0000000000000000000000000000000000000000..b1c2dec1e83241b3929ee89315eada335a119fd2
--- /dev/null
+++ b/AMDiS/src/expressions/AtomicExpression.h
@@ -0,0 +1,189 @@
+/******************************************************************************
+ *
+ * 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 AtomicExpression.h */
+
+#ifndef AMDIS_ATOMIC_EXPRESSION_H
+#define AMDIS_ATOMIC_EXPRESSION_H
+
+#include "AMDiS_fwd.h"
+#include "LazyOperatorTerm.h"
+#include "ValueTypes.h"
+
+#include <boost/static_assert.hpp>
+
+#include "ValueOf.h"
+#include "GradientOf.h"
+
+namespace AMDiS {
+
+namespace result_of {
+
+ struct Coords : public LazyOperatorTermBase
+ {
+ typedef WorldVector<double> value_type;
+ mutable mtl::dense_vector<WorldVector<double> > x;
+
+ Coords() {}
+
+ template<typename List>
+ void insertFeSpaces(List& feSpaces) const {}
+
+ int getDegree() const
+ {
+ return 1;
+ }
+
+ template<typename OT>
+ void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ subAssembler->getCoordsAtQPs(elInfo, quad, x);
+ }
+
+
+ template<typename OT>
+ void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ subAssembler->getCoordsAtQPs(smallElInfo, quad, x);
+ }
+
+ inline value_type operator()(const int& iq) const { return x[iq]; }
+
+ inline value_type eval(const int& iq) const { return x[iq]; }
+ inline value_type derivative(const int& iq, int identifier) const { WorldVector<double> vec0; vec0.set(0.0); return vec0; }
+ };
+
+
+ template<int I>
+ struct Coord : public LazyOperatorTermBase
+ {
+ typedef double value_type;
+ mutable mtl::dense_vector<WorldVector<double> > x;
+
+ Coord() {}
+
+ template<typename List>
+ void insertFeSpaces(List& feSpaces) const {}
+
+ int getDegree() const
+ {
+ return 1;
+ }
+
+ template<typename OT>
+ void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ subAssembler->getCoordsAtQPs(elInfo, quad, x);
+ }
+
+ template<typename OT>
+ void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ subAssembler->getCoordsAtQPs(smallElInfo, quad, x);
+ }
+
+ inline double operator()(const int& iq) const { return x[iq][I]; }
+
+ inline double eval(const int& iq) const { return x[iq][I]; }
+ inline double derivative(const int& iq, int identifier) const { return 0.0; }
+ };
+
+
+ template<typename T>
+ struct Const : public LazyOperatorTermBase
+ {
+ typedef T value_type;
+ T value;
+ Const(const T& value_) : value(value_) {}
+
+ template<typename List>
+ void insertFeSpaces(List& feSpaces) const {}
+
+ int getDegree() const
+ {
+ return 0;
+ }
+
+ template<typename OT>
+ void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad) {}
+
+ template<typename OT>
+ void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad) {}
+
+ inline value_type operator()(const int& iq) const { return value; }
+
+ inline value_type eval(const int& iq) const { return value; }
+ inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
+ };
+
+
+
+
+ ///
+ template<typename Term>
+ struct Jacobian : public LazyOperatorTerm1<Term>
+ {
+ typedef LazyOperatorTerm1<Term> super;
+ typedef typename Term::value_type value_type;
+
+ const int I;
+
+ Jacobian(const Term& term_, int I_) : super(term_), I(I_) {}
+
+ int getDegree() const
+ {
+ // TODO: eventuell noch Polynomgrad der Ableitung extra bestimmen
+ return super::term.getDegree();
+ }
+
+ inline value_type operator()(const int& iq) const { return super::term.derivative(iq, I); }
+ };
+
+}
+
+result_of::Coords X() { return result_of::Coords(); }
+
+template<int I>
+result_of::Coord<I> X() { return result_of::Coord<I>(); }
+
+template<typename T>
+result_of::Const<T> constant(const T& value) { return result_of::Const<T>(value); }
+
+namespace Private {
+
+ template<typename Term>
+ inline typename boost::enable_if<
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
+ result_of::Jacobian<Term> >::type
+ jacobian(const Term& t, int I) { return result_of::Jacobian<Term>(t, I); }
+
+}
+
+} // end namespace AMDiS
+
+#endif
\ No newline at end of file
diff --git a/AMDiS/src/expressions/GradientOf.h b/AMDiS/src/expressions/GradientOf.h
new file mode 100644
index 0000000000000000000000000000000000000000..d6aba4ff3e341e71e76c07b4019ccf1a54f4ef31
--- /dev/null
+++ b/AMDiS/src/expressions/GradientOf.h
@@ -0,0 +1,157 @@
+/******************************************************************************
+ *
+ * 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 GradientOf.h */
+
+#ifndef AMDIS_GRADIENT_OF_H
+#define AMDIS_GRADIENT_OF_H
+
+#include "AMDiS_fwd.h"
+#include "LazyOperatorTerm.h"
+#include "ValueTypes.h"
+#include "DOFVector.h"
+
+#include <boost/static_assert.hpp>
+
+namespace AMDiS {
+
+namespace result_of {
+
+ template<typename Vector>
+ struct GradientOf : public LazyOperatorTermBase
+ {
+ typedef typename ValueType<Vector>::type T;
+ typedef typename GradientType<T>::type value_type;
+
+ DOFVector<T>* vecDV;
+ mutable mtl::dense_vector<typename GradientType<T>::type> vec;
+
+ GradientOf(Vector& vector) : vecDV(&vector) {}
+ GradientOf(Vector* vector) : vecDV(vector) {}
+
+ template<typename List>
+ void insertFeSpaces(List& feSpaces) const
+ {
+ feSpaces.insert(vecDV->getFeSpace());
+ }
+
+ int getDegree() const
+ {
+ return vecDV->getFeSpace()->getBasisFcts()->getDegree() /* -1 */;
+ }
+
+ template<typename OT>
+ void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ ot.getGradientsAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ }
+
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ ot.getGradientsAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ }
+
+ value_type operator()(const int& iq) const { return vec[iq]; }
+ };
+
+
+ template<int I, typename Vector>
+ struct DerivativeOf : public LazyOperatorTermBase
+ {
+ typedef typename ValueType<Vector>::type T;
+ typedef T value_type;
+
+ DOFVector<T>* vecDV;
+ mutable mtl::dense_vector<typename GradientType<T>::type> vec;
+ int comp;
+
+ DerivativeOf(Vector& vector) : vecDV(&vector), comp(I) {}
+ DerivativeOf(Vector* vector) : vecDV(vector), comp(I) {}
+
+ DerivativeOf(Vector& vector, int I0) : vecDV(&vector), comp(I0)
+ {
+ TEST_EXIT_DBG( I < 0 && I0 >= 0 )
+ ("You yould specify eather template<int I>, or constructor(int I0)\n");
+ }
+ DerivativeOf(Vector* vector, int I0) : vecDV(vector), comp(I0)
+ {
+ TEST_EXIT_DBG( I < 0 && I0 >= 0 )
+ ("You yould specify eather template<int I>, or constructor(int I0)\n");
+ }
+
+ template<typename List>
+ void insertFeSpaces(List& feSpaces) const
+ {
+ feSpaces.insert(vecDV->getFeSpace());
+ }
+
+ int getDegree() const
+ {
+ return vecDV->getFeSpace()->getBasisFcts()->getDegree() /* -1 */;
+ }
+
+ template<typename OT>
+ void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ ot.getGradientsAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ }
+
+
+ template<typename OT>
+ void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ ot.getGradientsAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ }
+
+ value_type operator()(const int& iq) const { return vec[iq][comp]; }
+ };
+
+} // end namespace result_of
+
+template<typename T>
+result_of::GradientOf<DOFVector<T> > gradientOf(DOFVector<T>& vector) { return result_of::GradientOf<DOFVector<T> >(vector); }
+
+template<typename T>
+result_of::GradientOf<DOFVector<T> > gradientOf(DOFVector<T>* vector) { return result_of::GradientOf<DOFVector<T> >(vector); }
+
+template<int I, typename T>
+result_of::DerivativeOf<I, DOFVector<T> > derivativeOf(DOFVector<T>& vector) { return result_of::DerivativeOf<I, DOFVector<T> >(vector); }
+
+template<int I, typename T>
+result_of::DerivativeOf<I, DOFVector<T> > derivativeOf(DOFVector<T>* vector) { return result_of::DerivativeOf<I, DOFVector<T> >(vector); }
+
+template<typename T>
+result_of::DerivativeOf<-1, DOFVector<T> > derivativeOf(DOFVector<T>& vector, int I0) { return result_of::DerivativeOf<-1, DOFVector<T> >(vector, I0); }
+
+template<typename T>
+result_of::DerivativeOf<-1, DOFVector<T> > derivativeOf(DOFVector<T>* vector, int I0) { return result_of::DerivativeOf<-1, DOFVector<T> >(vector, I0); }
+
+} // end namespace AMDiS
+
+
+#endif
\ No newline at end of file
diff --git a/AMDiS/src/expressions/LazyExpression.h b/AMDiS/src/expressions/LazyExpression.h
new file mode 100644
index 0000000000000000000000000000000000000000..da406a1f67285c8767b4143a7bcab7413813ccd2
--- /dev/null
+++ b/AMDiS/src/expressions/LazyExpression.h
@@ -0,0 +1,648 @@
+/******************************************************************************
+ *
+ * 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 LazyExpression.h */
+
+#ifndef AMDIS_LAZY_EXPRESSION_H
+#define AMDIS_LAZY_EXPRESSION_H
+
+#include "AMDiS_fwd.h"
+#include "LazyOperatorTerm.h"
+#include "ValueTypes.h"
+#include "Functors.h"
+
+#include <boost/static_assert.hpp>
+
+
+namespace AMDiS {
+
+struct NoType {};
+
+struct FunctorBase
+{
+ typedef NoType value_type;
+protected:
+ virtual void helper() {};
+};
+
+
+namespace result_of {
+
+ ///
+ template<int I, typename Term>
+ struct Pow : public LazyOperatorTerm1<Term>
+ {
+ typedef LazyOperatorTerm1<Term> super;
+ typedef typename AMDiS::detail::Pow<I,typename Term::value_type>::result_type value_type;
+
+ Pow(const Term& term_) : super(term_) {}
+
+ int getDegree() const
+ {
+ return I * (super::term.getDegree());
+ }
+
+ inline value_type operator()(const int& iq) const { return AMDiS::detail::Pow<I,typename Term::value_type>::eval(super::term(iq)); }
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term.derivative(iq, identifier) * I*AMDiS::detail::Pow<I-1,typename Term::value_type>::eval(super::term(iq));
+ }
+ };
+
+
+ ///
+ template<typename Term>
+ struct Sqrt : public LazyOperatorTerm1<Term>
+ {
+ typedef LazyOperatorTerm1<Term> super;
+ typedef typename Term::value_type value_type;
+
+ Sqrt(const Term& term_) : super(term_) {}
+
+ int getDegree() const
+ {
+ return 2 * (super::term.getDegree()); // stimmt nicht ganz
+ }
+
+ inline value_type operator()(const int& iq) const { return sqrt(super::term(iq)); }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term.derivative(iq, identifier) / (2.0 * sqrt(super::term(iq)));
+ }
+ };
+
+
+ ///
+ template<typename Term>
+ struct Exp : public LazyOperatorTerm1<Term>
+ {
+ typedef LazyOperatorTerm1<Term> super;
+ typedef typename Term::value_type value_type;
+ int degree;
+
+ Exp(const Term& term_, int degree_ = 1) : super(term_), degree(degree_) {}
+
+ int getDegree() const
+ {
+ return degree * (super::term.getDegree());
+ }
+
+ inline value_type operator()(const int& iq) const { return exp(super::term(iq)); }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term.derivative(iq, identifier) * exp(super::term(iq));
+ }
+ };
+
+
+ ///
+ template<typename Term>
+ struct Abs : public LazyOperatorTerm1<Term>
+ {
+ typedef LazyOperatorTerm1<Term> super;
+ typedef typename Term::value_type value_type;
+
+ Abs(const Term& term_) : super(term_) {}
+
+ int getDegree() const
+ {
+ return super::term.getDegree();
+ }
+
+ inline value_type operator()(const int& iq) const { return std::abs(super::term(iq)); }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term(iq) > 0.0 ? super::term.derivative(iq, identifier) : -super::term.derivative(iq, identifier);
+ }
+ };
+
+
+ ///
+ template<typename Term>
+ struct Signum : public LazyOperatorTerm1<Term>
+ {
+ typedef LazyOperatorTerm1<Term> super;
+ typedef typename Term::value_type value_type;
+
+ Signum(const Term& term_) : super(term_) {}
+
+ int getDegree() const
+ {
+ return 0;
+ }
+
+ inline value_type operator()(const int& iq) const { return (super::term(iq) > 0.0 ? 1.0 : -1.0); }
+
+ inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
+ };
+
+
+ ///
+ template<typename Term1, typename Term2>
+ struct Add : public LazyOperatorTerm2<Term1, Term2>
+ {
+ typedef LazyOperatorTerm2<Term1, Term2> super;
+ typedef typename Term1::value_type value_type;
+
+ Add(const Term1& term1_, const Term2& term2_)
+ : super(term1_, term2_) {}
+
+ int getDegree() const
+ {
+ return std::max(super::term1.getDegree(), super::term2.getDegree());
+ }
+
+ inline value_type operator()(const int& iq) const { return super::term1(iq) + super::term2(iq); }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term1.derivative(iq, identifier) + super::term2.derivative(iq, identifier);
+ }
+ };
+
+
+ ///
+ template<typename Term1, typename Term2>
+ struct Mult : public LazyOperatorTerm2<Term1, Term2>
+ {
+ typedef LazyOperatorTerm2<Term1, Term2> super;
+ typedef typename ProductType<typename Term1::value_type,
+ typename Term2::value_type>::type value_type;
+
+ Mult(const Term1& term1_, const Term2& term2_)
+ : super(term1_, term2_) {}
+
+ int getDegree() const
+ {
+ return super::term1.getDegree() + super::term2.getDegree();
+ }
+
+ inline value_type operator()(const int& iq) const { return super::term1(iq) * super::term2(iq); }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term1.derivative(iq, identifier) * super::term2(iq) + super::term1(iq) * super::term2.derivative(iq, identifier);
+ }
+ };
+
+
+ ///
+ template<typename Term1, typename Term2>
+ struct Divide : public LazyOperatorTerm2<Term1, Term2>
+ {
+ typedef LazyOperatorTerm2<Term1, Term2> super;
+ typedef typename ProductType<typename Term1::value_type,
+ typename Term2::value_type>::type value_type;
+
+ Divide(const Term1& term1_, const Term2& term2_)
+ : super(term1_, term2_) {}
+
+ int getDegree() const
+ {
+ return super::term1.getDegree() + super::term2.getDegree(); // stimmt nicht ganz :)
+ }
+
+ inline value_type operator()(const int& iq) const { return super::term1(iq) / super::term2(iq); }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ typename Term2::value_type term2_eval = super::term2(iq);
+ return (super::term1.derivative(iq, identifier) * term2_eval - super::term1(iq) * super::term2.derivative(iq, identifier)) / (sqr(term2_eval));
+ }
+ };
+
+
+ ///
+ template<typename Term1, typename Term2>
+ struct Max : public LazyOperatorTerm2<Term1, Term2>
+ {
+ typedef LazyOperatorTerm2<Term1, Term2> super;
+ typedef typename Term1::value_type value_type;
+
+ Max(const Term1& term1_, const Term2& term2_)
+ : super(term1_, term2_) {}
+
+ int getDegree() const
+ {
+ return std::max(super::term1.getDegree(), super::term2.getDegree());
+ }
+
+ inline value_type operator()(const int& iq) const
+ {
+ return std::max(super::term1(iq), super::term2(iq));
+ }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term1(iq) > super::term2(iq) ? super::term1.derivative(iq, identifier) : super::term2.derivative(iq, identifier);
+ }
+ };
+
+
+ ///
+ template<typename Term1, typename Term2>
+ struct Min : public LazyOperatorTerm2<Term1, Term2>
+ {
+ typedef LazyOperatorTerm2<Term1, Term2> super;
+ typedef typename Term1::value_type value_type;
+
+ Min(const Term1& term1_, const Term2& term2_)
+ : super(term1_, term2_) {}
+
+ int getDegree() const
+ {
+ return std::max(super::term1.getDegree(), super::term2.getDegree());
+ }
+
+ inline value_type operator()(const int& iq) const
+ {
+ return std::min(super::term1(iq), super::term2(iq));
+ }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term1(iq) < super::term2(iq) ? super::term1.derivative(iq, identifier) : super::term2.derivative(iq, identifier);
+ }
+ };
+
+
+ ///
+ template<typename F, typename Term>
+ struct Function1 : public LazyOperatorTerm1<Term>
+ {
+ typedef LazyOperatorTerm1<Term> super;
+ BOOST_STATIC_ASSERT_MSG( (boost::is_base_of<FunctorBase, F>::value), "********** ERROR: Only functors with base FunctorBase allowed **********" );
+
+ typedef typename F::value_type value_type;
+ BOOST_STATIC_ASSERT_MSG( (!boost::is_same<value_type, NoType>::value), "********** ERROR: You have to define a value_type for your Functor **********" );
+
+ F f;
+ Function1(const F& f_, const Term& term_) : super(term_), f(f_) {}
+
+ int getDegree() const
+ {
+ return f.getDegree(super::term.getDegree());
+ }
+
+ inline value_type operator()(const int& iq) const { return f(super::term(iq)); }
+
+ inline value_type derivative(const int& iq, int identifier) const { return super::term.derivative(iq, identifier) * f.derivative<0>(super::term(iq)); }
+ };
+
+
+ ///
+ template<typename F, typename Term1, typename Term2>
+ struct Function2 : public LazyOperatorTerm2<Term1, Term2>
+ {
+ typedef LazyOperatorTerm2<Term1, Term2> super;
+ BOOST_STATIC_ASSERT_MSG( (boost::is_base_of<FunctorBase, F>::value), "********** ERROR: Only functors with base FunctorBase allowed **********" );
+
+ typedef typename F::value_type value_type;
+ BOOST_STATIC_ASSERT_MSG( (!boost::is_same<value_type, NoType>::value), "********** ERROR: You have to define a value_type for your Functor **********" );
+
+ F f;
+
+ Function2(const F& f_, const Term1& term1_, const Term2& term2_) : super(term1_, term2_), f(f_) {}
+
+ int getDegree() const
+ {
+ return f.getDegree(super::term1.getDegree(), super::term2.getDegree());
+ }
+
+ inline value_type operator()(const int& iq) const { return f(super::term1(iq), super::term2(iq)); }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term1.derivative(iq, identifier) * f.derivative<0>(super::term1(iq), super::term2(iq))
+ +super::term2.derivative(iq, identifier) * f.derivative<1>(super::term1(iq), super::term2(iq));
+ }
+ };
+
+
+ ///
+ template<typename F, typename Term1, typename Term2, typename Term3>
+ struct Function3 : public LazyOperatorTerm3<Term1, Term2, Term3>
+ {
+ typedef LazyOperatorTerm3<Term1, Term2, Term3> super;
+ BOOST_STATIC_ASSERT_MSG( (boost::is_base_of<FunctorBase, F>::value), "********** ERROR: Only functors with base FunctorBase allowed **********" );
+
+ typedef typename F::value_type value_type;
+ BOOST_STATIC_ASSERT_MSG( (!boost::is_same<value_type, NoType>::value), "********** ERROR: You have to define a value_type for your Functor **********" );
+
+ F f;
+
+ Function3(const F& f_, const Term1& term1_, const Term2& term2_, const Term3& term3_)
+ : super(term1_, term2_, term3_), f(f_) {}
+
+ int getDegree() const
+ {
+ return f.getDegree(super::term1.getDegree(), super::term2.getDegree(), super::term3.getDegree());
+ }
+
+ inline value_type operator()(const int& iq) const { return f(super::term1(iq), super::term2(iq), super::term3(iq)); }
+
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return super::term1.derivative(iq, identifier) * f.derivative<0>(super::term1(iq), super::term2(iq), super::term3(iq))
+ +super::term2.derivative(iq, identifier) * f.derivative<1>(super::term1(iq), super::term2(iq), super::term3(iq))
+ +super::term3.derivative(iq, identifier) * f.derivative<2>(super::term1(iq), super::term2(iq), super::term3(iq));
+ }
+ };
+}
+
+// add two terms
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Add<Term1, Term2> >::type
+add(const Term1& t1, const Term2& t2) { return result_of::Add<Term1, Term2>(t1, t2); }
+
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Add<Term1, Term2> >::type
+operator+(const Term1& t1, const Term2& t2) { return result_of::Add<Term1, Term2>(t1, t2); }
+
+template<typename T, typename Term>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Add<result_of::Const<T>, Term> >::type
+operator+(const T& t1, const Term& t2) { return result_of::Add<result_of::Const<T>, Term>(constant(t1), t2); }
+
+template<typename Term, typename T>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Add<Term, result_of::Const<T> > >::type
+operator+(const Term& t1, const T& t2) { return result_of::Add<Term, result_of::Const<T> >(t1, constant(t2)); }
+
+
+// subtract two terms
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Add<Term1, result_of::Mult<result_of::Const<double>, Term2> > >::type
+subtract(const Term1& t1, const Term2& t2)
+{
+ return result_of::Add<Term1, result_of::Mult<result_of::Const<double>, Term2> >
+ (t1, result_of::Mult<result_of::Const<double>, Term2>
+ (result_of::Const<double>(-1.0), t2) );
+}
+
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Add<Term1, result_of::Mult<result_of::Const<double>, Term2> > >::type
+operator-(const Term1& t1, const Term2& t2) { return subtract(t1, t2); }
+
+template<typename T, typename Term>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Add<result_of::Const<T>, result_of::Mult<result_of::Const<double>, Term> > >::type
+operator-(const T& t1, const Term& t2) { return subtract(constant(t1), t2); }
+
+template<typename Term, typename T>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Add<Term, result_of::Const<T> > >::type
+operator-(const Term& t1, const T& t2) { return result_of::Add<Term, result_of::Const<T> >(t1, constant(-t2)); }
+
+
+// multiply two terms
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Mult<Term1, Term2> >::type
+mult(const Term1& t1, const Term2& t2) { return result_of::Mult<Term1, Term2>(t1, t2); }
+
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Mult<Term1, Term2> >::type
+operator*(const Term1& t1, const Term2& t2) { return result_of::Mult<Term1, Term2>(t1, t2); }
+
+template<typename T, typename Term>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Mult<result_of::Const<T>, Term> >::type
+operator*(const T& t1, const Term& t2) { return result_of::Mult<result_of::Const<T>, Term>(constant(t1), t2); }
+
+template<typename Term, typename T>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Mult<Term, result_of::Const<T> > >::type
+operator*(const Term& t1, const T& t2) { return result_of::Mult<Term, result_of::Const<T> >(t1, constant(t2)); }
+
+
+// divide two terms
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Divide<Term1, Term2> >::type
+divide(const Term1& t1, const Term2& t2) { return result_of::Divide<Term1, Term2>(t1, t2); }
+
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Divide<Term1, Term2> >::type
+operator/(const Term1& t1, const Term2& t2) { return divide(t1, t2); }
+
+template<typename T, typename Term>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Divide<result_of::Const<T>, Term> >::type
+operator/(const T& t1, const Term& t2) { return divide(constant(t1), t2); }
+
+template<typename Term, typename T>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Divide<Term, result_of::Const<T> > >::type
+operator/(const Term& t1, const T& t2) { return divide(t1, constant(t2)); }
+
+// maximum of two terms
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Max<Term1, Term2> >::type
+max(const Term1& t1, const Term2& t2) { return result_of::Max<Term1, Term2>(t1, t2); }
+
+template<typename T, typename Term>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Max<result_of::Const<T>, Term> >::type
+max(const T& t1, const Term& t2) { return result_of::Max<result_of::Const<T>, Term>(constant(t1), t2); }
+
+template<typename Term, typename T>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Max<Term, result_of::Const<T> > >::type
+max(const Term& t1, const T& t2) { return result_of::Max<Term, result_of::Const<T> >(t1, constant(t2)); }
+
+
+
+// minimum of two terms
+template<typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Min<Term1, Term2> >::type
+min(const Term1& t1, const Term2& t2) { return result_of::Min<Term1, Term2>(t1, t2); }
+
+template<typename T, typename Term>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Min<result_of::Const<T>, Term> >::type
+min(const T& t1, const Term& t2) { return result_of::Max<result_of::Const<T>, Term>(constant(t1), t2); }
+
+template<typename Term, typename T>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_fundamental<T>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type>::type,
+ result_of::Min<Term, result_of::Const<T> > >::type
+min(const Term& t1, const T& t2) { return result_of::Min<Term, result_of::Const<T> >(t1, constant(t2)); }
+
+
+// I'th power of a term
+template<int I, typename Term>
+inline typename boost::enable_if<
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
+ result_of::Pow<I, Term> >::type
+pow(const Term& t) { return result_of::Pow<I, Term>(t); }
+
+
+// square root of a term
+template<typename Term>
+inline typename boost::enable_if<
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
+ result_of::Sqrt<Term> >::type
+sqrt(const Term& t) { return result_of::Sqrt<Term>(t); }
+
+
+// exponential function of a term
+template<typename Term>
+inline typename boost::enable_if<
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
+ result_of::Exp<Term> >::type
+exp(const Term& t) { return result_of::Exp<Term>(t); }
+
+
+// absolute value of a term
+template<typename Term>
+inline typename boost::enable_if<
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
+ result_of::Abs<Term> >::type
+abs(const Term& t) { return result_of::Abs<Term>(t); }
+
+
+// signum of a term
+template<typename Term>
+inline typename boost::enable_if<
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
+ result_of::Signum<Term> >::type
+signum(const Term& t) { return result_of::Signum<Term>(t); }
+
+
+
+// call a function with 1 argument
+template<typename F, typename Term>
+inline typename boost::enable_if<
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
+ result_of::Function1<F,Term> >::type
+function_(const F& f, const Term& t) { return result_of::Function1<F,Term>(f,t); }
+
+template<typename F, typename Term>
+inline typename boost::enable_if<
+ typename boost::is_base_of<LazyOperatorTermBase, Term>::type,
+ result_of::Function1<F,Term> >::type
+function_(const Term& t) { return result_of::Function1<F,Term>(F(),t); }
+
+
+// call a function with 2 arguments
+template<typename F, typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Function2<F, Term1, Term2> >::type
+function_(const F& f, const Term1& t1, const Term2& t2) { return result_of::Function2<F,Term1,Term2>(f,t1,t2); }
+
+template<typename F, typename Term1, typename Term2>
+inline typename boost::enable_if<
+ typename boost::mpl::and_<typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type>::type,
+ result_of::Function2<F, Term1, Term2> >::type
+function_(const Term1& t1, const Term2& t2) { return result_of::Function2<F,Term1,Term2>(F(),t1,t2); }
+
+
+// call a function with 3 arguments
+template<typename F, typename Term1, typename Term2, typename Term3>
+inline
+typename boost::enable_if<
+ typename boost::mpl::and_<
+ typename boost::mpl::and_<
+ typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type
+ >::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term3>::type
+ >::type,
+ result_of::Function3<F, Term1, Term2, Term3>
+ >::type
+function_(const F& f, const Term1& t1, const Term2& t2, const Term3& t3) { return result_of::Function3<F,Term1,Term2,Term3>(f,t1,t2,t3); }
+
+template<typename F, typename Term1, typename Term2, typename Term3>
+inline
+typename boost::enable_if<
+ typename boost::mpl::and_<
+ typename boost::mpl::and_<
+ typename boost::is_base_of<LazyOperatorTermBase, Term1>::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term2>::type
+ >::type,
+ typename boost::is_base_of<LazyOperatorTermBase, Term3>::type
+ >::type,
+ result_of::Function3<F, Term1, Term2, Term3>
+ >::type
+function_(const Term1& t1, const Term2& t2, const Term3& t3) { return result_of::Function3<F,Term1,Term2,Term3>(F(),t1,t2,t3); }
+
+} // end namespace AMDiS
+
+#endif
\ No newline at end of file
diff --git a/AMDiS/src/expressions/LazyOperatorTerm.h b/AMDiS/src/expressions/LazyOperatorTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..67b15eeaa29e43a2d787ec4188c23b9b4b324ed4
--- /dev/null
+++ b/AMDiS/src/expressions/LazyOperatorTerm.h
@@ -0,0 +1,141 @@
+/******************************************************************************
+ *
+ * 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 LazyOperatorTerm.h */
+
+#ifndef AMDIS_LAZY_OPERATOR_TERM_H
+#define AMDIS_LAZY_OPERATOR_TERM_H
+
+#include "AMDiS_fwd.h"
+#include <boost/static_assert.hpp>
+
+namespace AMDiS {
+
+struct LazyOperatorTermBase
+{
+protected:
+ virtual void helper() {}
+};
+
+template<typename Term>
+struct LazyOperatorTerm1 : public LazyOperatorTermBase
+{
+ Term term;
+ LazyOperatorTerm1(const Term& term_) : term(term_) {}
+
+ template<typename List>
+ inline void insertFeSpaces(List& feSpaces)
+ {
+ term.insertFeSpaces(feSpaces);
+ }
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ term.initElement(ot, elInfo, subAssembler, quad);
+ }
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ term.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
+ }
+
+ inline double operator()(const int& iq) const;
+};
+
+template<typename Term1, typename Term2>
+struct LazyOperatorTerm2 : public LazyOperatorTermBase
+{
+ Term1 term1;
+ Term2 term2;
+ LazyOperatorTerm2(const Term1& term1_, const Term2& term2_) : term1(term1_), term2(term2_) {}
+
+ template<typename List>
+ inline void insertFeSpaces(List& feSpaces)
+ {
+ term1.insertFeSpaces(feSpaces);
+ term2.insertFeSpaces(feSpaces);
+ }
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ term1.initElement(ot, elInfo, subAssembler, quad);
+ term2.initElement(ot, elInfo, subAssembler, quad);
+ }
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ term1.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
+ term2.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
+ }
+
+ double operator()(const int& iq) const;
+};
+
+template<typename Term1, typename Term2, typename Term3>
+struct LazyOperatorTerm3 : public LazyOperatorTermBase
+{
+ Term1 term1;
+ Term2 term2;
+ Term3 term3;
+ LazyOperatorTerm3(const Term1& term1_, const Term2& term2_, const Term3& term3_)
+ : term1(term1_), term2(term2_), term3(term3_) {}
+
+ template<typename List>
+ inline void insertFeSpaces(List& feSpaces)
+ {
+ term1.insertFeSpaces(feSpaces);
+ term2.insertFeSpaces(feSpaces);
+ term3.insertFeSpaces(feSpaces);
+ }
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ term1.initElement(ot, elInfo, subAssembler, quad);
+ term2.initElement(ot, elInfo, subAssembler, quad);
+ term3.initElement(ot, elInfo, subAssembler, quad);
+ }
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ term1.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
+ term2.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
+ term3.initElement(ot, smallElInfo, largeElInfo, subAssembler, quad);
+ }
+
+ inline double operator()(const int& iq) const;
+};
+
+} // end namespace AMDiS
+
+#endif
\ No newline at end of file
diff --git a/AMDiS/src/expressions/ValueOf.h b/AMDiS/src/expressions/ValueOf.h
new file mode 100644
index 0000000000000000000000000000000000000000..6914ecd18d7829b05af5f6b95aa825b2ee8b16b8
--- /dev/null
+++ b/AMDiS/src/expressions/ValueOf.h
@@ -0,0 +1,209 @@
+/******************************************************************************
+ *
+ * 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_VALUE_OF_H
+#define AMDIS_VALUE_OF_H
+
+#include "AMDiS_fwd.h"
+#include "LazyOperatorTerm.h"
+#include "ValueTypes.h"
+#include "DOFVector.h"
+
+#include <boost/static_assert.hpp>
+
+
+namespace AMDiS {
+
+namespace result_of {
+
+ template<typename Vector>
+ struct ValueOf : public LazyOperatorTermBase {};
+
+ template<typename T>
+ struct ValueOf<DOFVector<T> > : public LazyOperatorTermBase
+ {
+ typedef T value_type;
+
+ DOFVector<T>* vecDV;
+ mutable mtl::dense_vector<T> vec;
+
+ ValueOf(DOFVector<T>& vector) : vecDV(&vector) {}
+ ValueOf(DOFVector<T>* vector) : vecDV(vector) {}
+
+ template<typename List>
+ inline void insertFeSpaces(List& feSpaces) const
+ {
+ feSpaces.insert(vecDV->getFeSpace());
+ }
+
+ inline int getDegree() const
+ {
+ return vecDV->getFeSpace()->getBasisFcts()->getDegree();
+ }
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ ot.getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ }
+
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ ot.getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ }
+
+ inline value_type operator()(const int& iq) const { return vec[iq]; }
+ inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
+ };
+
+
+ template<typename T, template<class> class Vector>
+ struct ValueOf<Vector<DOFVector<T>*> > : public LazyOperatorTermBase
+ {
+ typedef Vector<T> value_type;
+
+ Vector<DOFVector<T>*> vecDV;
+ mutable mtl::dense_vector<value_type> vec;
+
+ ValueOf(Vector<DOFVector<T>*>& vector) : vecDV(vector) {}
+
+ template<typename List>
+ inline void insertFeSpaces(List& feSpaces) const
+ {
+ for (size_t i = 0; i < num_rows(vecDV); i++)
+ feSpaces.insert(vecDV[i]->getFeSpace());
+ }
+
+ inline int getDegree() const
+ {
+ return vecDV[0]->getFeSpace()->getBasisFcts()->getDegree();
+ }
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ Vector<mtl::dense_vector<T> > helper(num_rows(vecDV));
+ for (size_t i = 0; i < num_rows(vecDV); i++)
+ ot.getVectorAtQPs(vecDV[i], elInfo, subAssembler, quad, helper[i]);
+ vec.change_dim(num_rows(helper[0]));
+ for (size_t iq = 0; iq < num_rows(helper[0]); iq++) {
+ value_type tmp(num_rows(vecDV));
+ for (size_t i = 0; i < num_rows(vecDV); i++)
+ tmp[i] = helper[i][iq];
+ vec[iq] = tmp;
+ }
+ }
+
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ Vector<mtl::dense_vector<T> > helper(num_rows(vecDV));
+ for (size_t i = 0; i < num_rows(vecDV); i++)
+ ot.getVectorAtQPs(vecDV[i], smallElInfo, largeElInfo, subAssembler, quad, helper[i]);
+ vec.change_dim(num_rows(helper[0]));
+ for (size_t iq = 0; iq < num_rows(helper[0]); iq++) {
+ value_type tmp(num_rows(vecDV));
+ for (size_t i = 0; i < num_rows(vecDV); i++)
+ tmp[i] = helper[i][iq];
+ vec[iq] = tmp;
+ }
+ }
+
+ inline value_type operator()(const int& iq) const { return vec[iq]; }
+ inline value_type derivative(const int& iq, int identifier) const { return 0.0; }
+ };
+
+
+ struct ValueOf2 : public LazyOperatorTermBase
+ {
+ typedef double value_type;
+
+ DOFVector<value_type>* vecDV;
+
+ mutable mtl::dense_vector<value_type> vec;
+
+ const int I; // component of problem
+
+ ValueOf2(ProblemStat& prob, int I_) : vecDV(prob.getSolution(I_)), I(I_) {}
+ ValueOf2(ProblemStat* prob, int I_) : vecDV(prob->getSolution(I_)), I(I_) {}
+
+ // for Rosenbrock problems
+ ValueOf2(RosenbrockStationary& prob, int I_)
+ : vecDV(prob.getStageSolution(I_)), I(I_) {}
+ ValueOf2(RosenbrockStationary* prob, int I_)
+ : vecDV(prob->getStageSolution(I_)), I(I_) {}
+
+ template<typename List>
+ inline void insertFeSpaces(List& feSpaces) const
+ {
+ feSpaces.insert(vecDV->getFeSpace());
+ }
+
+ inline int getDegree() const
+ {
+ return vecDV->getFeSpace()->getBasisFcts()->getDegree();
+ }
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* elInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ ot.getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
+ }
+
+
+ template<typename OT>
+ inline void initElement(OT& ot, const ElInfo* smallElInfo, const ElInfo* largeElInfo,
+ SubAssembler* subAssembler, Quadrature *quad)
+ {
+ ot.getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
+ }
+
+ inline value_type operator()(const int& iq) const { return vec[iq]; }
+ inline value_type derivative(const int& iq, int identifier) const
+ {
+ return static_cast<value_type>(identifier == I);
+ }
+ };
+
+
+} // end namespace result_of
+
+
+template<typename T>
+result_of::ValueOf<DOFVector<T> > valueOf(DOFVector<T>& vector) { return result_of::ValueOf<DOFVector<T> >(vector); }
+
+template<typename T>
+result_of::ValueOf<DOFVector<T> > valueOf(DOFVector<T>* vector) { return result_of::ValueOf<DOFVector<T> >(vector); }
+
+} // end namespace AMDiS
+
+#endif // AMDIS_VALUE_OF_H
\ No newline at end of file
diff --git a/AMDiS/src/io/detail/VtkWriter.h b/AMDiS/src/io/detail/VtkWriter.h
index 5abb4280ca0fa88be6e5f24fbdfd7fbf105b774f..1d831b350fbe2400ba94ef4cbe206796a3c86061 100644
--- a/AMDiS/src/io/detail/VtkWriter.h
+++ b/AMDiS/src/io/detail/VtkWriter.h
@@ -103,6 +103,16 @@ namespace AMDiS { namespace io {
namespace detail
{
+ inline std::string timeToString(double time, char delimiter = '.')
+ {
+ std::stringstream ss; ss << time;
+ std::string s = ss.str();
+
+ size_t pos = s.find_first_of ('.', 0);
+ if (pos != std::string::npos)
+ s[pos] = delimiter;
+ return s;
+ }
/// Adds a new entry to a ParaView animation file.
int updateAnimationFile(AdaptInfo *adaptInfo,
diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc
index 2d2b2b9c5c0dd17a3e4984cecd73cc91217d8234..e6d00f378d9673d6cbb68c9130e94854e3443724 100644
--- a/AMDiS/src/parallel/MeshDistributor.cc
+++ b/AMDiS/src/parallel/MeshDistributor.cc
@@ -361,7 +361,7 @@ namespace AMDiS { namespace Parallel {
#endif
// === Global refinements. ===
-
+
int globalRefinement = 0;
Parameters::get(mesh->getName() + "->global refinements", globalRefinement);