cleanup and workaround of dune-typetree problems

dune/amdis/DirichletBC.hpp

@@ -7,6 +7,7 @@
 
 #include <dune/amdis/Output.hpp>
 #include <dune/amdis/common/Concepts.hpp>
+#include <dune/amdis/common/ValueCategory.hpp>
 #include <dune/amdis/utility/RangeType.hpp>
 #include <dune/amdis/utility/TreeData.hpp>
 
@@ -66,7 +67,10 @@ namespace AMDiS
                 VectorX& rhs,
                 VectorB& solution,
                 RowBasis const& rowBasis,
-                ColBasis const& colBasis);
+                ColBasis const& colBasis)
+    {
+      finishImpl(matrix, rhs, solution, rowBasis, colBasis, ValueCategory_t<Range>{});
+    }
 
   protected:
     template <class RowBasis>
@@ -78,6 +82,16 @@ namespace AMDiS
     template <class RB, class RowNodeTag>
     void initImpl(RB const&, RowNodeTag) {}
 
+    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis, class ValueCat>
+    void finishImpl(Matrix& matrix, VectorX& rhs, VectorB& solution,
+                    RowBasis const& rowBasis, ColBasis const& colBasis,
+                    ValueCat);
+
+    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
+    void finishImpl(Matrix& matrix, VectorX& rhs, VectorB& solution,
+                    RowBasis const& rowBasis, ColBasis const& colBasis,
+                    tag::unknown) {}
+
   private:
     std::function<bool(Domain)> predicate_;
     std::function<Range(Domain)> values_;
@@ -93,14 +107,7 @@ namespace AMDiS
     using WorldVector = typename GlobalBasis::GridView::template Codim<0>::Geometry::GlobalCoordinate;
 
     template <class RowNode, class ColNode>
-    using type = std::list< std::shared_ptr<DirichletBC<WorldVector, RangeType<RowNode>>> >;
-    // {
-    //   using Range = RangeType<RowNode>;
-    //   std::list<std::shared_ptr<DirichletBC<WorldVector, Range>>> bc;
-
-    //   void push_back(std::shared_ptr<DirichletBC<WorldVector, double>> const& bc) { scalar.push_back(bc); }
-    //   void push_back(std::shared_ptr<DirichletBC<WorldVector, WorldVector>> const& bc) { vector.push_back(bc); }
-    // };
+    using type = std::list< std::shared_ptr<DirichletBC<WorldVector, RangeType_t<RowNode>>> >;
   };
 
   template <class GlobalBasis>

dune/amdis/DirichletBC.inc.hpp

@@ -39,10 +39,10 @@ namespace AMDiS
 
 
   template <class WorldVector, class Range>
-    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis>
-  void DirichletBC<WorldVector, Range>::finish(
+    template <class Matrix, class VectorX, class VectorB, class RowBasis, class ColBasis, class ValueCat>
+  void DirichletBC<WorldVector, Range>::finishImpl(
       Matrix& matrix, VectorX& solution, VectorB& rhs,
-      RowBasis const& rowBasis, ColBasis const& colBasis)
+      RowBasis const& rowBasis, ColBasis const& colBasis, ValueCat)
       // Dune::TypeTree::LeafNodeTag, Dune::TypeTree::LeafNodeTag)
   {
     using Dune::Functions::interpolate;

dune/amdis/FileWriter.hpp

@@ -65,6 +65,12 @@ namespace AMDiS
                Vector const& dofvector)
       : FileWriterInterface(baseName)
       , dofvector_(dofvector)
+    {
+      init(baseName, ValueCategory_t<Range>{});
+    }
+
+    template <class ValueCat>
+    void init(std::string const& baseName, ValueCat)
     {
       int m = Parameters::get<int>(baseName + "->ParaView mode").value_or(0);
       mode_ = (m == 0 ? Dune::VTK::ascii : Dune::VTK::appendedraw);
@@ -80,11 +86,13 @@ namespace AMDiS
       vtkWriter_->addVertexData(dofvector_, Dune::VTK::FieldInfo(name_, VTKFieldType<Range>, VTKFieldSize<Range>));
     }
 
+    void init(std::string const&, tag::unknown) {}
 
     /// Implements \ref FileWriterInterface::writeFiles
     virtual void writeFiles(AdaptInfo& adaptInfo, bool force) override
     {
-      vtkSeqWriter_->write(adaptInfo.getTime(), mode_);
+      if (vtkSeqWriter_)
+        vtkSeqWriter_->write(adaptInfo.getTime(), mode_);
     }
 
   protected:

dune/amdis/GridFunctionOperator.hpp

@@ -5,7 +5,7 @@
 
 #include <dune/amdis/GridFunctions.hpp>
 #include <dune/amdis/common/Utility.hpp>
-#include <dune/amdis/utility/GetDegree.hpp>
+#include <dune/amdis/utility/FiniteElementType.hpp>
 
 namespace AMDiS
 {
@@ -110,7 +110,7 @@ namespace AMDiS
     {
       assert( bound_ );
 
-      int psiDegree = getPolynomialDegree(node);
+      int psiDegree = polynomialDegree(node);
 
       int degree = psiDegree + coeffDegree;
       if (isSimplex_)
@@ -133,8 +133,8 @@ namespace AMDiS
     {
       assert( bound_ );
 
-      int psiDegree = getPolynomialDegree(rowNode);
-      int phiDegree = getPolynomialDegree(colNode);
+      int psiDegree = polynomialDegree(rowNode);
+      int phiDegree = polynomialDegree(colNode);
 
       int degree = psiDegree + phiDegree + coeffDegree;
       if (isSimplex_)

dune/amdis/Initfile.hpp

@@ -7,6 +7,8 @@
 #include <list>
 #include <vector>
 
+#define BOOST_BIND_NO_PLACEHOLDERS
+
 #include <boost/lexical_cast.hpp>
 #include <boost/numeric/conversion/cast.hpp>
 #include <boost/tokenizer.hpp>

dune/amdis/LocalAssemblerBase.hpp

@@ -4,21 +4,41 @@
 #include <memory>
 #include <type_traits>
 
+#include <dune/amdis/common/ConceptsBase.hpp>
 #include <dune/amdis/common/TypeDefs.hpp>
-#include <dune/amdis/utility/GetEntity.hpp>
 #include <dune/amdis/utility/TreeData.hpp>
 
 namespace AMDiS
 {
+  namespace Impl
+  {
+    template <class E, class = Void_t<>>
+    struct ContextImpl
+    {
+      using Entity = E;
+      using Geometry = typename E::Geometry;
+    };
+
+    // specialization for intersections
+    template <class I>
+    struct ContextImpl<I, Void_t<typename I::Entity>>
+    {
+      using Entity = typename I::Entity;
+      using Geometry = typename I::LocalGeometry;
+    };
+
+  } // end namespace Impl
+
+
   /// Abstract base-class of a \ref LocalAssembler
   template <class LocalContext, class... Nodes>
   class LocalAssemblerBase
   {
   public:
 
-    using Element = typename Impl::Get<LocalContext>::Entity;
+    using Element = typename Impl::ContextImpl<LocalContext>::Entity;
     using Geometry = typename Element::Geometry;
-    using LocalGeometry = typename Impl::Get<LocalContext>::Geometry;
+    using LocalGeometry = typename Impl::ContextImpl<LocalContext>::Geometry;
 
     static constexpr int numNodes = sizeof...(Nodes);
     static_assert( numNodes == 1 || numNodes == 2,

dune/amdis/ProblemStat.hpp

@@ -10,7 +10,6 @@
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
-#include <dune/istl/matrix.hh>
 #include <dune/grid/common/grid.hh>
 
 #include <dune/amdis/AdaptInfo.hpp>
@@ -25,8 +24,6 @@
 #include <dune/amdis/ProblemStatTraits.hpp>
 #include <dune/amdis/StandardProblemIteration.hpp>
 
-#include <dune/amdis/common/OptionalDelete.hpp>
-#include <dune/amdis/common/Timer.hpp>
 #include <dune/amdis/common/TupleUtility.hpp>
 #include <dune/amdis/common/TypeDefs.hpp>
 #include <dune/amdis/common/Utility.hpp>
@@ -36,8 +33,6 @@
 
 #include <dune/amdis/io/FileWriterInterface.hpp>
 
-#include <dune/amdis/utility/DiscreteFunction.hpp>
-#include <dune/amdis/utility/RangeType.hpp>
 #include <dune/amdis/utility/TreePath.hpp>
 
 namespace AMDiS

dune/amdis/ProblemStat.inc.hpp

@@ -4,6 +4,7 @@
 #include <string>
 #include <utility>
 
+#include <dune/common/timer.hh>
 #include <dune/typetree/childextraction.hh>
 
 #include <dune/amdis/AdaptInfo.hpp>
@@ -12,7 +13,6 @@
 #include <dune/amdis/LocalAssembler.hpp>
 #include <dune/amdis/GridFunctionOperator.hpp>
 #include <dune/amdis/common/Loops.hpp>
-#include <dune/amdis/common/Timer.hpp>
 
 namespace AMDiS {
 
@@ -234,7 +234,7 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val
 
   auto valueGridFct = makeGridFunction(values, globalBasis->gridView());
 
-  using Range = typename decltype(valueGridFct)::Range;
+  using Range = RangeType_t<decltype(i)>;
   using BcType = DirichletBC<WorldVector,Range>;
   auto bc = std::make_shared<BcType>(predicate, valueGridFct);
   constraints[i][j].push_back(bc);
@@ -246,7 +246,7 @@ template <class Traits>
 void ProblemStat<Traits>::
 solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
 {
-  Timer t;
+  Dune::Timer t;
 
   SolverInfo solverInfo(name + "->solver");
   solverInfo.setCreateMatrixData(createMatrixData);
@@ -287,7 +287,7 @@ template <class Traits>
 void ProblemStat<Traits>::
 buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
 {
-  Timer t;
+  Dune::Timer t;
 
   Assembler<Traits> assembler(*globalBasis, matrixOperators, rhsOperators, constraints);
 
@@ -303,7 +303,7 @@ template <class Traits>
 void ProblemStat<Traits>::
 writeFiles(AdaptInfo& adaptInfo, bool force)
 {
-  Timer t;
+  Dune::Timer t;
   for (auto writer : filewriter)
     writer->writeFiles(adaptInfo, force);
   msg("writeFiles needed ", t.elapsed(), " seconds");

dune/amdis/assembler/FirstOrderPartialTestTrial.hpp

@@ -15,7 +15,7 @@ namespace AMDiS
   {
     struct partialtest_trial
     {
-      int comp;
+      std::size_t comp;
     };
   }
 

dune/amdis/assembler/FirstOrderTestPartialTrial.hpp

@@ -17,7 +17,7 @@ namespace AMDiS
   {
     struct test_partialtrial
     {
-      int comp;
+      std::size_t comp;
     };
   }
 
@@ -75,7 +75,7 @@ namespace AMDiS
 
         for (std::size_t i = 0; i < colPartial.size(); ++i) {
           colPartial[i] = jacobian[comp_][0] * referenceGradients[i][0][0];
-          for (std::size_t j = 1; j < jacobian.cols(); ++j)
+          for (std::size_t j = 1; j < jacobian.M(); ++j)
             colPartial[i] += jacobian[comp_][j] * referenceGradients[i][0][j];
         }
 
@@ -92,7 +92,7 @@ namespace AMDiS
     }
 
   private:
-    int comp_;
+    std::size_t comp_;
     std::vector<Dune::FieldVector<double,1>> rowShapeValues_;
   };
 

dune/amdis/assembler/SecondOrderPartialTestPartialTrial.hpp

@@ -17,8 +17,8 @@ namespace AMDiS
   {
     struct partialtest_partialtrial
     {
-      int comp_test; // i
-      int comp_trial; // j
+      std::size_t comp_test; // i
+      std::size_t comp_trial; // j
     };
   }
 
@@ -103,8 +103,8 @@ namespace AMDiS
     }
 
   private:
-    int compTest_;
-    int compTrial_;
+    std::size_t compTest_;
+    std::size_t compTrial_;
   };
 
   /** @} **/

dune/amdis/assembler/StokesOperator.hpp

+#pragma once
+
+#include <type_traits>
+#include <vector>
+
+#include <dune/amdis/GridFunctionOperator.hpp>
+#include <dune/amdis/common/ValueCategory.hpp>
+
+namespace AMDiS
+{
+  /**
+   * \addtogroup operators
+   * @{
+   **/
+
+  namespace tag
+  {
+    struct stokes {};
+  }
+
+
+  /// Stokes operator \f$ \langle\nabla\mathbf{v}, \nu \nabla\mathbf{u}\rangle + \langle\nabla\cdot\mathbf{v}, p\rangle + \langle q, \langle\nabla\cdot\mathbf{u}\rangle \f$
+  template <class ViscosityExpr, class QuadCreator>
+  class GridFunctionOperator<tag::stokes, ViscosityExpr, QuadCreator>
+      : public GridFunctionOperatorBase<ViscosityExpr, QuadCreator>
+  {
+    using Super = GridFunctionOperatorBase<ViscosityExpr, QuadCreator>;
+
+    static_assert( Category::Scalar<typename ViscosityExpr::Range>, "Viscosity must be of scalar type." );
+
+  public:
+    GridFunctionOperator(tag::stokes, ViscosityExpr const& expr, QuadCreator const& quadCreator)
+      : Super(expr, quadCreator, 1)
+    {}
+
+    template <class Context, class QuadratureRule,
+              class ElementMatrix, class Node, class Flag1, class Flag2>
+    void calculateElementMatrix(Context const& context,
+                                QuadratureRule const& quad,
+                                ElementMatrix& elementMatrix,
+                                Node const& tree,
+                                Node const& /*colNode*/,
+                                Flag1, Flag2)
+    {
+      using VelocityNode = typename Node::template Child<0>::type;
+      using PressureNode = typename Node::template Child<1>::type;
+      static_assert(VelocityNode::isPower && PressureNode::isLeaf, "");
+
+      using namespace Dune::Indices;
+
+      auto const& velocityLocalFE = tree.child(_0,0).finiteElement();
+      auto const& pressureLocalFE = tree.child(_1).finiteElement();
+
+      std::vector<Dune::FieldVector<double,1>> pressureValues;
+      std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
+
+      for (std::size_t iq = 0; iq < quad.size(); ++iq) {
+        // Position of the current quadrature point in the reference element
+        decltype(auto) local = context.position(quad[iq].position());
+
+        // The transposed inverse Jacobian of the map from the reference element to the element
+        const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+
+        // The multiplicative factor in the integral transformation formula
+        const double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
+        const double vel_factor = Super::coefficient(local) * factor;
+
+        // The gradients of the shape functions on the reference element
+        velocityLocalFE.localBasis().evaluateJacobian(local, referenceGradients);
+
+        // Compute the shape function gradients on the real element
+        std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
+
+        for (std::size_t i = 0; i < gradients.size(); ++i)
+          jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+        pressureLocalFE.localBasis().evaluateFunction(local, pressureValues);
+
+        // <viscosity * grad(u), grad(v)>
+        for (std::size_t i = 0; i < velocityLocalFE.size(); ++i) {
+          auto value_ii = vel_factor * (gradients[i] * gradients[i]);
+          for (std::size_t k = 0; k < Context::dow; ++k) {
+            std::size_t local_ki = tree.child(_0,k).localIndex(i);
+            elementMatrix[local_ki][local_ki] += value_ii;
+          }
+
+          for (std::size_t j = i+1; j < velocityLocalFE.size(); ++j) {
+            auto value = vel_factor * (gradients[i] * gradients[j]);
+            for (std::size_t k = 0; k < Context::dow; ++k) {
+              std::size_t local_ki = tree.child(_0,k).localIndex(i);
+              std::size_t local_kj = tree.child(_0,k).localIndex(j);
+              elementMatrix[local_ki][local_kj] += value;
+              elementMatrix[local_kj][local_ki] += value;
+            }
+          }
+        }
+
+        // <p, div(v)> + <div(u), q>
+        for (std::size_t i = 0; i < velocityLocalFE.size(); ++i) {
+          for (std::size_t j = 0; j < pressureLocalFE.size(); ++j) {
+            auto value = factor * pressureValues[j];
+            for (std::size_t k = 0; k < Context::dow; ++k) {
+              std::size_t local_v = tree.child(_0,k).localIndex(i);
+              std::size_t local_p = tree.child(_1).localIndex(j);
+
+              elementMatrix[local_v][local_p] += gradients[i][k] * value;
+              elementMatrix[local_p][local_v] += gradients[i][k] * value;
+            }
+          }
+        }
+      }
+
+    } // end calculateElementMatrix
+  };
+
+  /** @} **/
+
+} // end namespace AMDiS

dune/amdis/common/CMakeLists.txt

@@ -11,9 +11,9 @@ install(FILES
     Loops.hpp
     Math.hpp
     Mpl.hpp
-    OptionalDelete.hpp
+    MultiTypeFieldVector.hpp
     ScalarTypes.hpp
-    Timer.hpp
+    Size.hpp
     TupleUtility.hpp
     TypeDefs.hpp
     Utility.hpp

dune/amdis/common/Concepts.hpp

@@ -5,6 +5,7 @@
 #include <dune/functions/common/functionconcepts.hh>
 
 #include <dune/amdis/common/ConceptsBase.hpp>
+#include <dune/amdis/common/Mpl.hpp>
 #include <dune/amdis/common/Utility.hpp>
 
 namespace AMDiS

dune/amdis/common/FieldMatVec.hpp

@@ -97,6 +97,15 @@ namespace AMDiS
       return result;
     }
 
+    template <class T, int N, class Operation>
+    T accumulate(FieldMatrix<T, 1, N> const& x, Operation op)
+    {
+      T result = 0;
+      for (int i = 0; i < N; ++i)
+        result = op(result, x[0][i]);
+      return result;
+    }
+
   } // end namespace Impl
 
   /// Sum of vector entires.
@@ -106,6 +115,12 @@ namespace AMDiS
     return Impl::accumulate(x, Operation::Plus{});
   }
 
+  template <class T, int N>
+  T sum(FieldMatrix<T, 1, N> const& x)
+  {
+    return Impl::accumulate(x, Operation::Plus{});
+  }
+
 
   /// Dot-product with the vector itself
   template <class T, int N>
@@ -115,6 +130,13 @@ namespace AMDiS
     return Impl::accumulate(x, op);
   }
 
+  template <class T, int N>
+  auto unary_dot(FieldMatrix<T, 1, N> const& x)
+  {
+    auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); };
+    return Impl::accumulate(x, op);
+  }
+
   /// Maximum over all vector entries
   template <class T, int N>
   auto max(FieldVector<T, N> const& x)
@@ -122,6 +144,12 @@ namespace AMDiS
     return Impl::accumulate(x, Operation::Max{});
   }
 
+  template <class T, int N>
+  auto max(FieldMatrix<T, 1, N> const& x)
+  {
+    return Impl::accumulate(x, Operation::Max{});
+  }
+
   /// Minimum over all vector entries
   template <class T, int N>
   auto min(FieldVector<T, N> const& x)
@@ -129,6 +157,12 @@ namespace AMDiS
     return Impl::accumulate(x, Operation::Min{});
   }
 
+  template <class T, int N>
+  auto min(FieldMatrix<T, 1, N> const& x)
+  {
+    return Impl::accumulate(x, Operation::Min{});
+  }
+
   /// Maximum of the absolute values of vector entries
   template <class T, int N>
   auto abs_max(FieldVector<T, N> const& x)
@@ -136,6 +170,12 @@ namespace AMDiS
     return Impl::accumulate(x, Operation::AbsMax{});
   }
 
+  template <class T, int N>
+  auto abs_max(FieldMatrix<T, 1, N> const& x)
+  {
+    return Impl::accumulate(x, Operation::AbsMax{});
+  }
+
   /// Minimum of the absolute values of vector entries
   template <class T, int N>
   auto abs_min(FieldVector<T, N> const& x)
@@ -143,6 +183,12 @@ namespace AMDiS
     return Impl::accumulate(x, Operation::AbsMin{});
   }
 
+  template <class T, int N>
+  auto abs_min(FieldMatrix<T, 1, N> const& x)
+  {
+    return Impl::accumulate(x, Operation::AbsMin{});
+  }
+
   // ----------------------------------------------------------------------------
 
   /** \ingroup vector_norms
@@ -155,6 +201,13 @@ namespace AMDiS
     return Impl::accumulate(x, op);
   }
 
+  template <class T, int N>
+  auto one_norm(FieldMatrix<T, 1, N> const& x)
+  {
+    auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
+    return Impl::accumulate(x, op);
+  }
+
   /** \ingroup vector_norms
    *  \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 )
    **/
@@ -164,6 +217,12 @@ namespace AMDiS
     return std::sqrt(unary_dot(x));
   }
 
+  template <class T, int N>
+  auto two_norm(FieldMatrix<T, 1, N> const& x)
+  {
+    return std::sqrt(unary_dot(x));
+  }
+
   /** \ingroup vector_norms
    *  \brief The p-norm of a vector = ( sum_i |x_i|^p )^(1/p)
    **/
@@ -174,6 +233,13 @@ namespace AMDiS
     return std::pow( Impl::accumulate(x, op), 1.0/p );
   }
 
+  template <int p, class T, int N>
+  auto p_norm(FieldMatrix<T, 1, N> const& x)
+  {
+    auto op = [](auto const& a, auto const& b) { return a + Math::pow<p>(std::abs(b)); };
+    return std::pow( Impl::accumulate(x, op), 1.0/p );
+  }
+
   /** \ingroup vector_norms
    *  \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max
    **/
@@ -183,6 +249,12 @@ namespace AMDiS
     return abs_max(x);
   }
 
+  template <class T, int N>
+  auto infty_norm(FieldMatrix<T, 1, N> const& x)
+  {
+    return abs_max(x);
+  }
+
   // ----------------------------------------------------------------------------
 
   /// The euklidean distance between two vectors = |lhs-rhs|_2
@@ -342,7 +414,7 @@ namespace AMDiS
   }
 
   template <class T, int M, int N>
-  Dune::FieldMatrix<T,M,N>& multiplies_ABt(Dune::FieldMatrix<T, M, N> const& A,  Dune::DiagonalMatrix<T, N> const& B, Dune::FieldMatrix<T,M,N>& C)
+  FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A,  Dune::DiagonalMatrix<T, N> const& B, FieldMatrix<T,M,N>& C)
   {