redesign polynomialDegree and order of gridfunction and basis-nodes

amdis/Integrate.hpp

@@ -32,15 +32,6 @@ namespace AMDiS
       return gridView.comm().sum(result);
     }
 
-    template <class GF, class GV, class QP>
-    auto integrateImpl(GF&& gf, GV const& gv, QP makeQuad, std::true_type)
-    {
-      return integrateImpl(FWD(gf), gv, makeQuad);
-    }
-
-    template <class GF, class GV, class QP>
-    auto integrateImpl(GF&&, GV&&, QP&&, std::false_type) { return 0.0; }
-
   } // end namespace Impl
 
 
@@ -58,16 +49,9 @@ namespace AMDiS
   {
     auto&& gridFct = makeGridFunction(FWD(expr), gridView);
 
-    // test whether the gridFct model `Concepts::HasLocalFunctionOrder`
-    using GF = TYPEOF(gridFct);
-    static const bool expr_has_order = Concepts::HasLocalFunctionOrder<GF>;
-    static_assert(expr_has_order,
-      "Polynomial degree of expression can not be deduced. You need to provide an explicit value for the quadrature degree or a quadrature rule in `integrate()`.");
-
     using Rules = Dune::QuadratureRules<typename GridView::ctype, GridView::dimension>;
-    auto makeQuad = [](auto&& t, auto&& lf) { return Rules::rule(t, order(lf)); };
-
-    return Impl::integrateImpl(FWD(gridFct), gridView, makeQuad, bool_<expr_has_order>);
+    return Impl::integrateImpl(FWD(gridFct), gridView,
+      [](auto&& t, auto&& lf) { return Rules::rule(t, polynomialDegree(lf).value()); });
   }
 
 
@@ -84,7 +68,8 @@ namespace AMDiS
   auto integrate(Expr&& expr, GridView const& gridView, QuadratureRule const& quad)
   {
     auto&& gridFct = makeGridFunction(FWD(expr), gridView);
-    return Impl::integrateImpl(FWD(gridFct), gridView, [&](auto&&, auto&&) { return quad; });
+    return Impl::integrateImpl(FWD(gridFct), gridView,
+      [&](auto&&, auto&&) { return quad; });
   }
 
 

amdis/LocalOperator.hpp

@@ -132,8 +132,8 @@ namespace AMDiS
       test_warning(coeffDegree >= 0,
         "polynomial order of coefficient function not determined. Use order 4 by default.");
 
-      int psiDegree = polynomialDegree(rowNode);
-      int phiDegree = polynomialDegree(colNode);
+      int psiDegree = polynomialDegree(rowNode).value_or(1);
+      int phiDegree = polynomialDegree(colNode).value_or(1);
 
       int degree = psiDegree + phiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
       if (isSimplex_)
@@ -157,7 +157,7 @@ namespace AMDiS
       test_warning(coeffDegree >= 0,
         "polynomial order of coefficient function not determined. Use order 4 by default.");
 
-      int psiDegree = polynomialDegree(node);
+      int psiDegree = polynomialDegree(node).value_or(1);
 
       int degree = psiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
       if (isSimplex_)

amdis/common/CMakeLists.txt

@@ -23,6 +23,7 @@ install(FILES
     Literals.hpp
     Logical.hpp
     Math.hpp
+    PolynomialDegree.hpp
     Range.hpp
     QuadMath.hpp
     SharedPtr.hpp

amdis/common/PolynomialDegree.hpp

+#pragma once
+
+#include <optional>
+#include <type_traits>
+
+#include <dune/common/std/type_traits.hh>
+
+namespace AMDiS
+{
+  namespace Traits
+  {
+    template <class T>
+    using OrderFreeFunction = decltype(order(std::declval<T>()));
+
+    template <class T>
+    using OrderMemberFunction = decltype(std::declval<T>().order());
+
+  } // end namespace Traits
+
+
+  /// Polynomial degree of a local function or a finite-element space or anything
+  /// else that defines an order function.
+  /**
+   * Return the derived or defined polynomial degree of a function `f`
+   * If no order is defined, return an empty optional.
+   **/
+  template <class F>
+  std::optional<int> polynomialDegree(F const& f)
+  {
+    if constexpr (Dune::Std::is_detected_v<Traits::OrderFreeFunction,F>)
+      return order(f);
+    else if constexpr (Dune::Std::is_detected_v<Traits::OrderMemberFunction,F>)
+      return f.order();
+    else
+      return std::nullopt;
+  }
+
+} // end namespace AMDiS
\ No newline at end of file

amdis/functions/CMakeLists.txt

@@ -5,5 +5,6 @@ install(FILES
     Interpolate.hpp
     NodeIndices.hpp
     Nodes.hpp
+    Order.hpp
     ParallelGlobalBasis.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions)

amdis/functions/Order.hpp

+#pragma once
+
+#include <cmath>
+#include <type_traits>
+
+#include <dune/typetree/nodetags.hh>
+
+#include <amdis/common/ConceptsBase.hpp>
+#include <amdis/common/ForEach.hpp>
+
+namespace AMDiS
+{
+  template <class Node>
+  auto order(Node const& node)
+    -> decltype(node.finiteElement().localBasis().order())
+  {
+    return node.finiteElement().localBasis().order();
+  }
+
+  template <class Node,
+    REQUIRES(std::is_same_v<typename Node::NodeTag, Dune::TypeTree::PowerNodeTag>)>
+  int order(Node const& node)
+  {
+    return order(node.child(0u));
+  }
+
+  template <class Node,
+    REQUIRES(std::is_same_v<typename Node::NodeTag, Dune::TypeTree::CompositeNodeTag>)>
+  int order(Node const& node)
+  {
+    int degree = 0;
+    Tools::for_range<0,Node::CHILDREN>([&](auto const _i) {
+      degree = std::max(degree, int(order(node.child(_i))));
+    });
+    return degree;
+  }
+
+} // end namespace AMDiS

amdis/gridfunctions/DiscreteLocalFunction.inc.hpp

@@ -3,7 +3,9 @@
 #include <amdis/Output.hpp>
 #include <amdis/common/DerivativeTraits.hpp>
 #include <amdis/common/FieldMatVec.hpp>
+#include <amdis/common/Math.hpp>
 #include <amdis/functions/NodeIndices.hpp>
+#include <amdis/functions/Order.hpp>
 #include <amdis/utility/LocalBasisCache.hpp>
 #include <amdis/utility/LocalToGlobalAdapter.hpp>
 
@@ -130,11 +132,13 @@ public:
     return PartialLocalFunction{globalFunction_, type};
   }
 
-  /// \brief The \ref polynomialDegree() of the LocalFunctions
-  int order() const
+  /// \brief The polynomial degree of the LocalFunctions basis-functions
+  auto order() const
+    -> decltype(order(std::declval<SubTree>()))
   {
     assert( bound_ );
-    return polynomialDegree(*subTree_);
+    using AMDiS::order;
+    return order(*subTree_);
   }
 
   /// \brief Return the bound element
@@ -207,10 +211,12 @@ public:
     bound_ = false;
   }
 
-  int order() const
+  auto order() const
+    -> decltype(order(std::declval<SubTree>()))
   {
     assert( bound_ );
-    return std::max(0, polynomialDegree(*subTree_)-1);
+    using AMDiS::order;
+    return Math::max(0, order(*subTree_)-1);
   }
 
   /// Return the bound element

amdis/gridfunctions/FunctorGridFunction.hpp

@@ -152,11 +152,11 @@ namespace AMDiS
   /**
    * **Requirements:**
    * - The functor `F` must model `Concepts::HasFunctorOrder`
-   * - All localFunctions `LFs...` must model `Concepts::HasOrder`
+   * - All localFunctions `LFs...` must model `Concepts::Polynomial`
    **/
   template <class Sig, class F, class... LFs,
     REQUIRES(Concepts::HasFunctorOrder<F,sizeof...(LFs)>
-             && (Concepts::HasOrder<LFs> &&...))>
+             && (Concepts::Polynomial<LFs> &&...))>
   int order(FunctorLocalFunction<Sig,F,LFs...> const& lf)
   {
     return Tools::apply([&lf](auto const&... lgfs) {

amdis/gridfunctions/Order.hpp

@@ -2,7 +2,7 @@
 
 #include <type_traits>
 
-#include <amdis/common/Concepts.hpp>
+#include <amdis/common/ConceptsBase.hpp>
 
 namespace AMDiS
 {
@@ -14,7 +14,6 @@ namespace AMDiS
     return lf.order();
   }
 
-
   namespace Concepts
   {
     /** \addtogroup Concepts
@@ -23,15 +22,7 @@ namespace AMDiS
 
     namespace Definition
     {
-      struct HasLocalFunctionOrder
-      {
-        template <class F>
-        auto require(F&& f) -> decltype(
-          order(localFunction(f))
-        );
-      };
-
-      struct HasOrder
+      struct Polynomial
       {
         template <class F>
         auto require(F&& f) -> decltype(
@@ -41,24 +32,14 @@ namespace AMDiS
 
     } // end namespace Definition
 
-
-    /// \brief GridFunction GF has free function `order(localFunction(F))`
-    template <class GF>
-    constexpr bool HasLocalFunctionOrder = models<Definition::HasLocalFunctionOrder(GF)>;
-
-    template <class GF>
-    using HasLocalFunctionOrder_t = models_t<Definition::HasLocalFunctionOrder(GF)>;
-
-
     /// \brief LocalFuncion LF has free function `order(F)`
     template <class LF>
-    constexpr bool HasOrder = models<Definition::HasOrder(LF)>;
+    constexpr bool Polynomial = models<Definition::Polynomial(LF)>;
 
     template <class LF>
-    using HasOrder_t = models_t<Definition::HasOrder(LF)>;
+    using Polynomial_t = models_t<Definition::Polynomial(LF)>;
 
     /** @} **/
 
   } // end namespace Concepts
-
 } // end namespace AMDiS

amdis/localoperators/ConvectionDiffusionOperator.hpp

@@ -176,16 +176,10 @@ namespace AMDiS
 
   private:
 
-    template <class LF>
-    using HasLocalFunctionOrder = decltype( order(std::declval<LF>()) );
-
     template <class LocalFct>
     int coeffOrder(LocalFct const& localFct)
     {
-      if constexpr (Dune::Std::is_detected<HasLocalFunctionOrder, LocalFct>::value)
-        return order(localFct);
-      else
-        return 0;
+      return polynomialDegree(localFct).value_or(0);
     }
 
     template <class T, int N>

amdis/typetree/FiniteElementType.hpp

@@ -5,6 +5,7 @@
 #include <dune/common/typetraits.hh>
 #include <dune/typetree/nodetags.hh>
 
+#include <amdis/common/Concepts.hpp>
 #include <amdis/common/ForEach.hpp>
 #include <amdis/common/Index.hpp>
 #include <amdis/common/Tags.hpp>
@@ -16,10 +17,8 @@ namespace AMDiS
     template <class Node, class NodeTag>
     struct FiniteElementTypeImpl
     {
-      static_assert( Dune::AlwaysFalse<NodeTag>::value, "Unknown node-type for range definition" );
-
-      // polynomial degree of finite-element basis functions
-      static int order(Node const&) { return 0; }
+      static_assert(Dune::AlwaysFalse<NodeTag>::value,
+        "Unknown node-type for range definition");
     };
   }
 
@@ -30,12 +29,6 @@ namespace AMDiS
   template <class Node>
   using FiniteElementType_t = typename FiniteElementType<Node>::type;
 
-  template <class Node>
-  int polynomialDegree(Node const& node)
-  {
-    return FiniteElementType<Node>::order(node);
-  }
-
 
   namespace Impl
   {
@@ -44,11 +37,6 @@ namespace AMDiS
     struct FiniteElementTypeImpl<Node, Dune::TypeTree::LeafNodeTag>
     {
       using type = typename Node::FiniteElement;
-
-      static int order(Node const& node)
-      {
-        return node.finiteElement().localBasis().order();
-      }
     };
 
     // Power node
@@ -57,11 +45,6 @@ namespace AMDiS
     {
       using ChildNode = typename Node::template Child<0>::type;
       using type = FiniteElementType_t<ChildNode>;
-
-      static int order(Node const& node)
-      {
-        return FiniteElementType<ChildNode>::order(node.child(0));
-      }
     };
 
     // Composite node
@@ -69,17 +52,7 @@ namespace AMDiS
     struct FiniteElementTypeImpl<Node, Dune::TypeTree::CompositeNodeTag>
     {
       using type = tag::unknown; // no common FiniteElement type
-
-      static int order(Node const& node)
-      {
-        int degree = 0;
-        Tools::for_range<0,Node::CHILDREN>([&](auto const _i) {
-          degree = std::max(degree, polynomialDegree(node.child(_i)));
-        });
-        return degree;
-      }
     };
 
   } // end namespace Impl
-
 } // end namespace AMDiS

amdis/utility/QuadratureFactory.hpp

@@ -5,6 +5,7 @@
 #include <dune/geometry/quadraturerules.hh>
 #include <dune/geometry/type.hh>
 
+#include <amdis/common/PolynomialDegree.hpp>
 #include <amdis/gridfunctions/Order.hpp>
 
 namespace AMDiS
@@ -39,27 +40,16 @@ namespace AMDiS
   };
 
 
-  template <class LF>
-  using HasLocalFunctionOrder = decltype( order(std::declval<LF>()) );
-
-
   /// \brief Factory for quadrature rule, that calculates the coefficient order from
   /// a localFunction passed to the bind method.
   template <class ctype, int dim, class LocalFunction>
   class QuadFactoryFromLocalFunction
       : public QuadratureFactory<ctype, dim, LocalFunction>
   {
-    using Concept = Dune::Std::is_detected<HasLocalFunctionOrder, LocalFunction>;
-    static_assert(Concept::value,
-      "Polynomial order of GridFunction can not be extracted. Provide an explicit order parameter instead.");
-
   public:
     void bind(LocalFunction const& localFct) final
     {
-      if constexpr (Concept::value)
-        order_ = AMDiS::order(localFct);
-      else
-        order_ = -1;
+      order_ = polynomialDegree(localFct).value();
     }
 
     int order() const final { return order_; }

test/FiniteElementTypeTest.cpp

@@ -4,6 +4,8 @@
 #include <dune/functions/functionspacebases/powerbasis.hh>
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 
+#include <amdis/common/PolynomialDegree.hpp>
+#include <amdis/functions/Order.hpp>
 #include <amdis/typetree/FiniteElementType.hpp>
 #include "Tests.hpp"
 
@@ -36,14 +38,14 @@ int main()
     localView.bind(e);
 
     auto node = localView.tree();
-    AMDIS_TEST_EQ( polynomialDegree(node), k+1 ); // maximum over all polynomial degrees
+    AMDIS_TEST_EQ( polynomialDegree(node).value(), k+1 ); // maximum over all polynomial degrees
     static_assert( std::is_same_v<tag::unknown, FiniteElementType_t<decltype(node)>>, "" );
 
     auto v_node = TypeTree::child(node, _0);
-    AMDIS_TEST_EQ( polynomialDegree(v_node), k+1 );
+    AMDIS_TEST_EQ( polynomialDegree(v_node).value(), k+1 );
 
     auto p_node = TypeTree::child(node, _1);
-    AMDIS_TEST_EQ( polynomialDegree(p_node), k );
+    AMDIS_TEST_EQ( polynomialDegree(p_node).value(), k );
   }
 
   return report_errors();