added virtual destructors, wrong extension string in mesh reader, forward...

added virtual destructors, wrong extension string in mesh reader, forward declaration in operation::composer

src/amdis/LocalAssemblerBase.hpp

@@ -51,6 +51,8 @@ namespace AMDiS
 
   public:
 
+    virtual ~LocalAssemblerBase() {}
+
     virtual void bind(Element const& element, Geometry const& geometry) = 0;
     virtual void unbind() = 0;
 
@@ -97,6 +99,21 @@ namespace AMDiS
       {
         return flag && (!assembled || changing);
       }
+
+      template <class Geo>
+      void bind(Element const& elem, Geo const& geo)
+      {
+        for (auto& scaled : element) scaled.op->bind(elem,geo);
+        for (auto& scaled : boundary) scaled.op->bind(elem,geo);
+        for (auto& scaled : intersection) scaled.op->bind(elem,geo);
+      }
+
+      void unbind()
+      {
+        for (auto& scaled : element) scaled.op->unbind();
+        for (auto& scaled : boundary) scaled.op->unbind();
+        for (auto& scaled : intersection) scaled.op->unbind();
+      }
     };
 
   public:

src/amdis/Mesh.hpp

@@ -113,14 +113,14 @@ namespace AMDiS
         auto ext = fn.extension();
 
 #if HAVE_ALBERTA
-        if (ext == "1d" || ext == "2d" || ext == "3d") {
+        if (ext == ".1d" || ext == ".2d" || ext == ".3d") {
           Dune::GridFactory<Grid> factory;
           Dune::AlbertaReader<Grid> reader;
           reader.readGrid(filename, factory);
           return std::unique_ptr<Grid>{factory.createGrid()};
         }
 #endif
-        if (ext == "msh") {
+        if (ext == ".msh") {
           Dune::GmshReader<Grid> reader;
           return std::unique_ptr<Grid>{reader.read(filename)};
         }

src/amdis/Operations.hpp

@@ -14,5 +14,6 @@
 #include <amdis/operations/Basic.hpp>
 #include <amdis/operations/CMath.hpp>
 #include <amdis/operations/Composer.hpp>
+#include <amdis/operations/Composer.impl.hpp>
 #include <amdis/operations/FieldMatVec.hpp>
 #include <amdis/operations/MaxMin.hpp>

src/amdis/ProblemInstat.hpp

@@ -79,6 +79,23 @@ namespace AMDiS
     std::unique_ptr<SystemVector> oldSolution;
   };
 
+
+
+#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
+  // Deduction rule
+  template <class Traits>
+  ProblemInstat(std::string name, ProblemStat<Traits>& prob)
+    -> ProblemInstat<Traits>;
+#endif
+
+  // Generator for ProblemInstat with given ProblemStat
+  template <class Traits>
+  ProblemInstat<Traits>
+  makeProblemInstat(std::string name, ProblemStat<Traits>& prob)
+  {
+    return {std::move(name), prob};
+  }
+
 } // end namespace AMDiS
 
 #include "ProblemInstat.inc.hpp"

src/amdis/assembler/SecondOrderGradTestGradTrial.hpp

@@ -65,11 +65,11 @@ namespace AMDiS
         const auto exprValue = Super::coefficient(local);
 
         // The gradients of the shape functions on the reference element
-        std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
+        thread_local std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
         rowLocalFE.localBasis().evaluateJacobian(local, referenceGradients);
 
         // Compute the shape function gradients on the real element
-        std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
+        thread_local 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]);
@@ -113,7 +113,7 @@ namespace AMDiS
                            tag::scalar)
     {
       auto const& localFE = node.finiteElement();
-      std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
+      thread_local 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
@@ -129,7 +129,7 @@ namespace AMDiS
         localFE.localBasis().evaluateJacobian(local, referenceGradients);
 
         // Compute the shape function gradients on the real element
-        std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
+        thread_local 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]);
@@ -160,7 +160,7 @@ namespace AMDiS
                            tag::matrix)
     {
       auto const& localFE = node.finiteElement();
-      std::vector<Dune::FieldMatrix<double,1,Context::dim> > referenceGradients;
+      thread_local 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
@@ -177,7 +177,7 @@ namespace AMDiS
         localFE.localBasis().evaluateJacobian(local, referenceGradients);
 
         // Compute the shape function gradients on the real element
-        std::vector<Dune::FieldVector<double,Context::dow> > gradients(referenceGradients.size());
+        thread_local 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]);

src/amdis/gridfunctions/FunctorGridFunction.hpp

@@ -139,7 +139,7 @@ namespace AMDiS
              && all_of_v<Concepts::HasOrder<LFs>...>)>
   int order(FunctorLocalFunction<Sig,F,LFs...> const& lf)
   {
-    return Dune::Std::apply([&lf](auto const&... gfs) { return order(lf.fct(), order(gfs)...); },
+    return Dune::Std::apply([&lf](auto const&... lgfs) { return order(lf.fct(), order(lgfs)...); },
       lf.localFcts());
   }
 

src/amdis/io/FileWriterInterface.hpp

@@ -24,6 +24,8 @@ namespace AMDiS
       Parameters::get(base + "->name", name_);
     }
 
+    virtual ~FileWriterInterface() {}
+
 
     // pure virtual method to be implemented by derived classes
     virtual void writeFiles(AdaptInfo& adaptInfo, bool force) = 0;

src/amdis/operations/Composer.hpp

@@ -67,34 +67,10 @@ namespace AMDiS
       return Dune::Std::apply([&](auto const&... gs) { return order(c.f_, deg(gs)...); }, c.gs_);
     }
 
-    // forward declaration
-    struct Multiplies;
-    struct Plus;
-
     /// Partial derivative of composed function:
     /// Implements: // sum_i [ d_i(f)[g...] * d_j(g_i) ]
     template <int J, class F, class... Gs>
-    auto partial(Composer<F,Gs...> const& c, index_t<J> _j)
-    {
-      auto index_seq = MakeSeq_t<sizeof...(Gs)>{};
-
-      // d_i(f)[g...] * d_j(g_i)
-      auto term_i = [&](auto const _i)
-      {
-        auto di_f = Dune::Std::apply([&](auto const&... gs) {
-          return compose(partial(c.f_, _i), gs...);
-        }, c.gs_);
-
-        auto const& g_i = std::get<_i>(c.gs_);
-        return compose(Multiplies{}, di_f, partial(g_i, _j));
-      };
-
-      // sum_i [ d_i(f)[g...] * d_j(g_i) ]
-      return Dune::Std::apply([&](auto const... _i)
-      {
-        return compose(Plus{}, term_i(_i)...);
-      }, index_seq);
-    }
+    auto partial(Composer<F,Gs...> const& c, index_t<J> _j);
 
 
 #ifndef DOXYGEN

src/amdis/utility/Traversal.hpp

@@ -4,219 +4,114 @@
 #include <dune/typetree/treepath.hh>
 #include <dune/typetree/visitor.hh>
 
-namespace AMDiS {
-
-#ifndef DOXYGEN // these are all internals and not public API. Only access is using applyToTree().
-
+namespace AMDiS
+{
   // forward declaration of main engine struct
-  template <typename tag = Dune::TypeTree::StartTag, bool doApply = true>
+  template <typename NodeTag, bool visit = true>
   struct TraverseTree;
 
-  // This struct is the core of the algorithm. While this specialization simply serves as the starting point
-  // of the traversal and takes care of some setup work, the struct has to be specialized for each TreeType node type it
-  // should support.
-  // The first parameter is the tag of the node type
-  // and the second one must always be specialized as true, as a value of false means the node should in fact not be visited.
-  // That case is already handled by a specialization of the struct.
-  template <bool doApply>
-  struct TraverseTree<Dune::TypeTree::StartTag, doApply>
-  {
-
-    template <typename Node, typename Visitor>
-    static void apply(Node&& node, Visitor&& visitor)
-    {
-      TraverseTree<Dune::TypeTree::NodeTag<Node>>::apply(std::forward<Node>(node),
-                                          std::forward<Visitor>(visitor),
-                                          Dune::TypeTree::hybridTreePath());
-    }
-
-  };
-
 
   // Do not visit nodes the visitor is not interested in
   template <typename NodeTag>
   struct TraverseTree<NodeTag, false>
   {
-    // we won't do anything with the objects, so having them all const
-    // works fine.
     template <typename Node, typename Visitor, typename TreePath>
-    static void apply(const Node& node, const Visitor& visitor, TreePath treePath)
+    static void apply(const Node& node, const Visitor& visitor, TreePath const& tp)
     {}
   };
 
+#ifndef DOXYGEN
 
+  // some implementation details
 
-  // ********************************************************************************
-  // LeafNode
-  // ********************************************************************************
+  template <class Node, class Index>
+  struct HybridChildType
+      : HybridChildType<std::remove_const_t<Node>, std::remove_const_t<Index>> {};
 
-  // LeafNode - just call the leaf() callback
-  template <>
-  struct TraverseTree<Dune::TypeTree::LeafNodeTag, true>
+  template <class Node>
+  struct HybridChildType<Node, std::size_t>
   {
-    template <typename N, typename V, typename TreePath>
-    static void apply(N&& n, V&& v, TreePath tp)
-    {
-      v.leaf(std::forward<N>(n),tp);
-    }
+    using type = typename Node::template Child<0>::Type;
   };
 
+  template <class Node, std::size_t K>
+  struct HybridChildType<Node, Dune::index_constant<K>>
+  {
+    using type = typename Node::template Child<K>::Type;
+  };
 
-  // ********************************************************************************
-  // PwerNode
-  // ********************************************************************************
+  template <class NodeTag, class Node>
+  constexpr std::size_t hybridDegree(NodeTag, Node const& node)
+  {
+    return Dune::TypeTree::degree(node);
+  }
 
-  template <>
-  struct TraverseTree<Dune::TypeTree::PowerNodeTag, true>
+  template <class Node>
+  constexpr auto hybridDegree(Dune::TypeTree::CompositeNodeTag, Node const& node)
   {
-    template <typename N, typename V, typename TreePath>
-    static void apply(N&& n, V&& v, TreePath tp)
-    {
-      v.pre(std::forward<N>(n),tp);
+    return Dune::index_constant<Node::CHILDREN>{};
+  }
 
-      typedef typename std::remove_reference<N>::type Node;
-      typedef typename std::remove_reference<V>::type Visitor;
 
-      // get child type
-      typedef typename Node::template Child<0>::Type C;
+  template <std::size_t k, std::size_t n>
+  constexpr bool notLastChild(Dune::index_constant<k> const&, Dune::index_constant<n> const&)
+  {
+    return k < n-1;
+  }
+
+  constexpr bool notLastChild(std::size_t k, std::size_t n)
+  {
+    return k < n-1;
+  }
 
-      // Do we have to visit the children? As the TreePath is dynamic, it does not
-      // contain any information that could be evaluated at compile time, so we only
-      // have to query the visitor once.
-      const bool visit = Visitor::template VisitChild<Node,C,TreePath>::value;
+#endif
 
-      auto child_tp = push_back(tp, std::size_t(0));
-      static const std::size_t last = treePathSize(tp);
 
-      // iterate over children
-      for (std::size_t k = 0; k < degree(n); ++k)
+  template <class NodeTag>
+  struct TraverseTree<NodeTag, true>
+  {
+    template <typename N, typename V, typename TreePath>
+    static void apply(N&& n, V&& v, TreePath const& tp)
+    {
+      using Node = std::remove_reference_t<N>;
+      using Visitor = std::remove_reference_t<V>;
+
+      v.pre(std::forward<N>(n),tp);
+
+      auto const deg = hybridDegree(NodeTag{}, n);
+      forEach(Dune::range(deg), [&](auto const _k)
       {
         // always call beforeChild(), regardless of the value of visit
-        v.beforeChild(std::forward<N>(n),n.child(k),tp,k);
-
-        // update TreePath
-        std::get<last>(child_tp._data) = k;
+        v.beforeChild(std::forward<N>(n),n.child(_k),tp,_k);
 
         // descend to child
-        TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.child(k),std::forward<V>(v),child_tp);
+        using C = typename HybridChildType<Node, decltype(_k)>::type;
+        const bool visit = Visitor::template VisitChild<Node,C,TreePath>::value;
+        TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.child(_k),std::forward<V>(v),push_back(tp, _k));
 
         // always call afterChild(), regardless of the value of visit
-        v.afterChild(std::forward<N>(n),n.child(k),tp,k);
+        v.afterChild(std::forward<N>(n),n.child(_k),tp,_k);
 
         // if this is not the last child, call infix callback
-        if (k < degree(n) - 1)
+        if (notLastChild(_k, deg))
           v.in(std::forward<N>(n),tp);
-      }
+      });
 
       v.post(std::forward<N>(n),tp);
     }
   };
 
-
-  namespace Impl
-  {
-    // ********************************************************************************
-    // Static Version
-    // ********************************************************************************
-
-    template <std::size_t inverse_k, std::size_t count>
-    struct TraverseTreeStatic
-    {
-      template <typename N, typename V, typename TreePath>
-      static void apply(N&& n, V&& v, TreePath tp)
-      {
-        // make sure we do not try to work with references to the actual types
-        typedef typename std::remove_reference<N>::type Node;
-        typedef typename std::remove_reference<V>::type Visitor;
-
-        // get child type
-        typedef typename Node::template Child<count-inverse_k>::Type C;
-
-        // extend TreePath by child index
-        auto k = std::integral_constant<std::size_t, count-inverse_k>{};
-        auto child_tp = push_back(tp, k);
-
-        // is the visitor interested in this child?
-        const bool visit = Visitor::template VisitChild<Node,C,decltype(child_tp)>::value;
-
-        // beforeChild() gets called regardless of the value of visit
-        v.beforeChild(std::forward<N>(n),n.template child<count-inverse_k>(),tp,k);
-
-        // traverse to child
-        TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.template child<count-inverse_k>(),
-                                              std::forward<V>(v),
-                                              child_tp);
-
-        // afterChild() gets called regardless of the value of visit
-        v.afterChild(std::forward<N>(n),n.template child<count-inverse_k>(),tp,k);
-
-        // we are not at the last child (that is specialized), so call infix visitor callback
-        v.in(std::forward<N>(n),tp);
-
-        // continue with next child
-        TraverseTreeStatic<inverse_k-1,count>::apply(std::forward<N>(n), std::forward<V>(v), tp);
-      }
-    };
-
-    // Specialization for last child. This specialization stops the recursion and
-    // does not call the infix visitor on the CompositeNode.
-    template <std::size_t count>
-    struct TraverseTreeStatic<1,count>
-    {
-      template<typename N, typename V, typename TreePath>
-      static void apply(N&& n, V&& v, TreePath tp)
-      {
-        typedef typename std::remove_reference<N>::type Node;
-        typedef typename std::remove_reference<V>::type Visitor;
-        typedef typename Node::template Child<count-1>::Type C;
-        auto k = std::integral_constant<std::size_t, count-1>{};
-        auto child_tp = push_back(tp, k);
-
-        const bool visit = Visitor::template VisitChild<Node,C,decltype(child_tp)>::value;
-        v.beforeChild(std::forward<N>(n),n.template child<count-1>(),tp,k);
-        TraverseTree<Dune::TypeTree::NodeTag<C>,visit>::apply(n.template child<count-1>(),
-                                              std::forward<V>(v),
-                                              child_tp);
-        v.afterChild(std::forward<N>(n),n.template child<count-1>(),tp,k);
-      }
-    };
-
-    // Specialization for CompositeNode without any children.
-    template <>
-    struct TraverseTreeStatic<0,0>
-    {
-      template <typename N, typename V, typename TreePath>
-      static void apply(N&& n, V&& v, TreePath tp) {}
-    };
-
-  } // end namespace Impl
-
-
-  // ********************************************************************************
-  // CompositeNode
-  // ********************************************************************************
-
-  // Traverse CompositeNode - just forward to the generic algorithm
+  // LeafNode - just call the leaf() callback
   template <>
-  struct TraverseTree<Dune::TypeTree::CompositeNodeTag,true>
+  struct TraverseTree<Dune::TypeTree::LeafNodeTag, true>
   {
     template <typename N, typename V, typename TreePath>
-    static void apply(N&& n, V&& v, TreePath tp)
+    static void apply(N&& n, V&& v, TreePath const& tp)
     {
-      v.pre(std::forward<N>(n),tp);
-      typedef typename std::remove_reference<N>::type Node;
-      static const std::size_t C = Dune::TypeTree::StaticDegree<Node>::value;
-      Impl::TraverseTreeStatic<C,C>::apply(std::forward<N>(n), std::forward<V>(v), tp);
-      v.post(std::forward<N>(n),tp);
+      v.leaf(std::forward<N>(n),tp);
     }
   };
 
-  #endif // DOXYGEN
-
-  // ********************************************************************************
-  // Public Interface
-  // ********************************************************************************
-
   //! Apply visitor to TypeTree.
   /**
     * This function applies the given visitor to the given tree. Both visitor and tree may be const
@@ -231,7 +126,10 @@ namespace AMDiS {
   template <typename Tree, typename Visitor>
   void traverseTree(Tree&& tree, Visitor&& visitor)
   {
-    TraverseTree<>::apply(std::forward<Tree>(tree), std::forward<Visitor>(visitor));
+    using Node = std::remove_reference_t<Tree>;
+    using NodeTag = Dune::TypeTree::NodeTag<Node>;
+    using TreePath = Dune::TypeTree::HybridTreePath<>;
+    TraverseTree<NodeTag>::apply(std::forward<Tree>(tree), std::forward<Visitor>(visitor), TreePath{});
   }
 
 } // end namespace AMDiS

src/amdis/utility/TreePath.hpp

@@ -145,6 +145,26 @@ namespace AMDiS
     return "";
   }
 
+  namespace Impl
+  {
+    template <class TreePath, std::size_t... I>
+    std::array<std::size_t, sizeof...(I)> toArrayImpl(TreePath const& tp, std::index_sequence<I...>)
+    {
+      return {{std::size_t(Dune::TypeTree::treePathEntry<I>(tp))...}};
+    }
+  }
+
+  template <class T0, class... T>
+  auto to_array(Dune::TypeTree::HybridTreePath<T0,T...> const& tp)
+  {
+    return Impl::toArrayImpl(tp, std::make_index_sequence<1+sizeof...(T)>{});
+  }
+
+  std::array<std::size_t,1> to_array(Dune::TypeTree::HybridTreePath<> const& tp)
+  {
+    return {{0u}};
+  }
+
   /// \brief Generate a TreePath from a sequence of integers and integral-constants
   /**
    * Converts a sequence of arguments to a \ref Dune::TypeTree::HybridTreePath.
@@ -158,4 +178,33 @@ namespace AMDiS
     return Dune::TypeTree::HybridTreePath<T...>(t...);
   }
 
+
+  namespace Impl
+  {
+    template <class TreePath, std::size_t... I>
+    auto popFrontImpl(TreePath const& tp, std::index_sequence<I...>)
+    {
+      return Dune::TypeTree::hybridTreePath(Dune::TypeTree::treePathEntry<I+1>(tp)...);
+    }
+
+    template <class TreePath, std::size_t... I>
+    auto popBackImpl(TreePath const& tp, std::index_sequence<I...>)
+    {
+      return Dune::TypeTree::hybridTreePath(Dune::TypeTree::treePathEntry<I>(tp)...);
+    }
+  }
+
+  template <class T0, class... T>
+  Dune::TypeTree::HybridTreePath<T...> pop_front(Dune::TypeTree::HybridTreePath<T0,T...> const& tp)
+  {
+    return Impl::popFrontImpl(tp, std::make_index_sequence<sizeof...(T)>{});
+  }
+
+  template <class... T, class TN>
+  Dune::TypeTree::HybridTreePath<T...> pop_front(Dune::TypeTree::HybridTreePath<T...,TN> const& tp)
+  {
+    return Impl::popBackImpl(tp, std::make_index_sequence<sizeof...(T)>{});
+  }
+
+
 } // end namespace AMDiS

src/amdis/utility/Visitor.hpp

@@ -9,22 +9,9 @@ namespace AMDiS
   // from dune-typetree merge-request !2
   namespace Impl
   {
-    struct CallInnerNode
-    {
-      template <class Node, class /*TreePath*/>
-      using Visit = bool_t<Node::isPower && Node::template Child<0>::type::isLeaf>;
-    };
-
-    struct CallAnyNode
-    {
-      template <class /*Node*/, class /*TreePath*/>
-      using Visit = std::true_type;
-    };
-
     template <class PreFunc, class LeafFunc, class PostFunc>
     class CallbackVisitor
         : public Dune::TypeTree::TreeVisitor
-        , public CallAnyNode
     {
     public:
       CallbackVisitor(PreFunc& preFunc, LeafFunc& leafFunc, PostFunc& postFunc)
@@ -36,9 +23,7 @@ namespace AMDiS
       template <typename Node, typename TreePath>
       void pre(Node&& node, TreePath treePath)
       {
-        static_assert(std::decay_t<Node>::isPower || std::decay_t<Node>::isComposite, "");
-        Dune::Hybrid::ifElse(Visit<std::decay_t<Node>, TreePath>{},
-          [&](auto id) { id(preFunc_)(node, treePath); });
+        preFunc_(node, treePath);
       }
 
       template <typename Node, typename TreePath>
@@ -50,9 +35,7 @@ namespace AMDiS
       template <typename Node, typename TreePath>
       void post(Node&& node, TreePath treePath)
       {
-        static_assert(std::decay_t<Node>::isPower || std::decay_t<Node>::isComposite, "");
-        Dune::Hybrid::ifElse(Visit<std::decay_t<Node>, TreePath>{},
-          [&](auto id) { id(postFunc_)(node, treePath); });
+        postFunc_(node, treePath);
       }
 
     private: