Commit 661a7e84 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'issue/polynomial_degree' into 'master'

redesign polynomialDegree and order of gridfunction and basis-nodes

See merge request amdis/amdis-core!196
parents 605d1fc1 69096682
...@@ -32,15 +32,6 @@ namespace AMDiS ...@@ -32,15 +32,6 @@ namespace AMDiS
return gridView.comm().sum(result); 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 } // end namespace Impl
...@@ -58,16 +49,9 @@ namespace AMDiS ...@@ -58,16 +49,9 @@ namespace AMDiS
{ {
auto&& gridFct = makeGridFunction(FWD(expr), gridView); 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>; 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,
[](auto&& t, auto&& lf) { return Rules::rule(t, polynomialDegree(lf).value()); });
return Impl::integrateImpl(FWD(gridFct), gridView, makeQuad, bool_<expr_has_order>);
} }
...@@ -84,7 +68,8 @@ namespace AMDiS ...@@ -84,7 +68,8 @@ namespace AMDiS
auto integrate(Expr&& expr, GridView const& gridView, QuadratureRule const& quad) auto integrate(Expr&& expr, GridView const& gridView, QuadratureRule const& quad)
{ {
auto&& gridFct = makeGridFunction(FWD(expr), gridView); 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; });
} }
......
...@@ -132,8 +132,8 @@ namespace AMDiS ...@@ -132,8 +132,8 @@ namespace AMDiS
test_warning(coeffDegree >= 0, test_warning(coeffDegree >= 0,
"polynomial order of coefficient function not determined. Use order 4 by default."); "polynomial order of coefficient function not determined. Use order 4 by default.");
int psiDegree = polynomialDegree(rowNode); int psiDegree = polynomialDegree(rowNode).value_or(1);
int phiDegree = polynomialDegree(colNode); int phiDegree = polynomialDegree(colNode).value_or(1);
int degree = psiDegree + phiDegree + (coeffDegree >= 0 ? coeffDegree : 4); int degree = psiDegree + phiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
if (isSimplex_) if (isSimplex_)
...@@ -157,7 +157,7 @@ namespace AMDiS ...@@ -157,7 +157,7 @@ namespace AMDiS
test_warning(coeffDegree >= 0, test_warning(coeffDegree >= 0,
"polynomial order of coefficient function not determined. Use order 4 by default."); "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); int degree = psiDegree + (coeffDegree >= 0 ? coeffDegree : 4);
if (isSimplex_) if (isSimplex_)
......
...@@ -23,6 +23,7 @@ install(FILES ...@@ -23,6 +23,7 @@ install(FILES
Literals.hpp Literals.hpp
Logical.hpp Logical.hpp
Math.hpp Math.hpp
PolynomialDegree.hpp
Range.hpp Range.hpp
QuadMath.hpp QuadMath.hpp
SharedPtr.hpp SharedPtr.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
...@@ -5,5 +5,6 @@ install(FILES ...@@ -5,5 +5,6 @@ install(FILES
Interpolate.hpp Interpolate.hpp
NodeIndices.hpp NodeIndices.hpp
Nodes.hpp Nodes.hpp
Order.hpp
ParallelGlobalBasis.hpp ParallelGlobalBasis.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/functions)
#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
...@@ -3,7 +3,9 @@ ...@@ -3,7 +3,9 @@
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/common/DerivativeTraits.hpp> #include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/FieldMatVec.hpp> #include <amdis/common/FieldMatVec.hpp>
#include <amdis/common/Math.hpp>
#include <amdis/functions/NodeIndices.hpp> #include <amdis/functions/NodeIndices.hpp>
#include <amdis/functions/Order.hpp>
#include <amdis/utility/LocalBasisCache.hpp> #include <amdis/utility/LocalBasisCache.hpp>
#include <amdis/utility/LocalToGlobalAdapter.hpp> #include <amdis/utility/LocalToGlobalAdapter.hpp>
...@@ -130,11 +132,13 @@ public: ...@@ -130,11 +132,13 @@ public:
return PartialLocalFunction{globalFunction_, type}; return PartialLocalFunction{globalFunction_, type};
} }
/// \brief The \ref polynomialDegree() of the LocalFunctions /// \brief The polynomial degree of the LocalFunctions basis-functions
int order() const auto order() const
-> decltype(order(std::declval<SubTree>()))
{ {
assert( bound_ ); assert( bound_ );
return polynomialDegree(*subTree_); using AMDiS::order;
return order(*subTree_);
} }
/// \brief Return the bound element /// \brief Return the bound element
...@@ -207,10 +211,12 @@ public: ...@@ -207,10 +211,12 @@ public:
bound_ = false; bound_ = false;
} }
int order() const auto order() const
-> decltype(order(std::declval<SubTree>()))
{ {
assert( bound_ ); assert( bound_ );
return std::max(0, polynomialDegree(*subTree_)-1); using AMDiS::order;
return std::max(0, int(order(*subTree_))-1);
} }
/// Return the bound element /// Return the bound element
......
...@@ -152,11 +152,11 @@ namespace AMDiS ...@@ -152,11 +152,11 @@ namespace AMDiS
/** /**
* **Requirements:** * **Requirements:**
* - The functor `F` must model `Concepts::HasFunctorOrder` * - 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, template <class Sig, class F, class... LFs,
REQUIRES(Concepts::HasFunctorOrder<F,sizeof...(LFs)> REQUIRES(Concepts::HasFunctorOrder<F,sizeof...(LFs)>
&& (Concepts::HasOrder<LFs> &&...))> && (Concepts::Polynomial<LFs> &&...))>
int order(FunctorLocalFunction<Sig,F,LFs...> const& lf) int order(FunctorLocalFunction<Sig,F,LFs...> const& lf)
{ {
return Tools::apply([&lf](auto const&... lgfs) { return Tools::apply([&lf](auto const&... lgfs) {
......
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
#include <type_traits> #include <type_traits>
#include <amdis/common/Concepts.hpp> #include <amdis/common/ConceptsBase.hpp>
namespace AMDiS namespace AMDiS
{ {
...@@ -14,7 +14,6 @@ namespace AMDiS ...@@ -14,7 +14,6 @@ namespace AMDiS
return lf.order(); return lf.order();
} }
namespace Concepts namespace Concepts
{ {
/** \addtogroup Concepts /** \addtogroup Concepts
...@@ -23,15 +22,7 @@ namespace AMDiS ...@@ -23,15 +22,7 @@ namespace AMDiS
namespace Definition namespace Definition
{ {
struct HasLocalFunctionOrder struct Polynomial
{
template <class F>
auto require(F&& f) -> decltype(
order(localFunction(f))
);
};
struct HasOrder
{ {
template <class F> template <class F>
auto require(F&& f) -> decltype( auto require(F&& f) -> decltype(
...@@ -41,24 +32,14 @@ namespace AMDiS ...@@ -41,24 +32,14 @@ namespace AMDiS
} // end namespace Definition } // 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)` /// \brief LocalFuncion LF has free function `order(F)`
template <class LF> template <class LF>
constexpr bool HasOrder = models<Definition::HasOrder(LF)>; constexpr bool Polynomial = models<Definition::Polynomial(LF)>;
template <class 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 Concepts
} // end namespace AMDiS } // end namespace AMDiS
...@@ -176,16 +176,10 @@ namespace AMDiS ...@@ -176,16 +176,10 @@ namespace AMDiS
private: private:
template <class LF>
using HasLocalFunctionOrder = decltype( order(std::declval<LF>()) );
template <class LocalFct> template <class LocalFct>
int coeffOrder(LocalFct const& localFct) int coeffOrder(LocalFct const& localFct)
{ {
if constexpr (Dune::Std::is_detected<HasLocalFunctionOrder, LocalFct>::value) return polynomialDegree(localFct).value_or(0);
return order(localFct);
else
return 0;
} }
template <class T, int N> template <class T, int N>
......
...@@ -5,6 +5,7 @@ ...@@ -5,6 +5,7 @@
#include <dune/common/typetraits.hh> #include <dune/common/typetraits.hh>
#include <dune/typetree/nodetags.hh> #include <dune/typetree/nodetags.hh>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ForEach.hpp> #include <amdis/common/ForEach.hpp>
#include <amdis/common/Index.hpp> #include <amdis/common/Index.hpp>
#include <amdis/common/Tags.hpp> #include <amdis/common/Tags.hpp>
...@@ -16,10 +17,8 @@ namespace AMDiS ...@@ -16,10 +17,8 @@ namespace AMDiS
template <class Node, class NodeTag> template <class Node, class NodeTag>
struct FiniteElementTypeImpl struct FiniteElementTypeImpl
{ {
static_assert( Dune::AlwaysFalse<NodeTag>::value, "Unknown node-type for range definition" ); 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; }
}; };
} }
...@@ -30,12 +29,6 @@ namespace AMDiS ...@@ -30,12 +29,6 @@ namespace AMDiS
template <class Node> template <class Node>
using FiniteElementType_t = typename FiniteElementType<Node>::type; using FiniteElementType_t = typename FiniteElementType<Node>::type;
template <class Node>
int polynomialDegree(Node const& node)
{
return FiniteElementType<Node>::order(node);
}
namespace Impl namespace Impl
{ {
...@@ -44,11 +37,6 @@ namespace AMDiS ...@@ -44,11 +37,6 @@ namespace AMDiS
struct FiniteElementTypeImpl<Node, Dune::TypeTree::LeafNodeTag> struct FiniteElementTypeImpl<Node, Dune::TypeTree::LeafNodeTag>
{ {
using type = typename Node::FiniteElement; using type = typename Node::FiniteElement;
static int order(Node const& node)
{
return node.finiteElement().localBasis().order();
}
}; };
// Power node // Power node
...@@ -57,11 +45,6 @@ namespace AMDiS ...@@ -57,11 +45,6 @@ namespace AMDiS
{ {
using ChildNode = typename Node::template Child<0>::type; using ChildNode = typename Node::template Child<0>::type;
using type = FiniteElementType_t<ChildNode>; using type = FiniteElementType_t<ChildNode>;
static int order(Node const& node)
{
return FiniteElementType<ChildNode>::order(node.child(0));
}
}; };
// Composite node // Composite node
...@@ -69,17 +52,7 @@ namespace AMDiS ...@@ -69,17 +52,7 @@ namespace AMDiS
struct FiniteElementTypeImpl<Node, Dune::TypeTree::CompositeNodeTag> struct FiniteElementTypeImpl<Node, Dune::TypeTree::CompositeNodeTag>
{ {
using type = tag::unknown; // no common FiniteElement type 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 Impl
} // end namespace AMDiS } // end namespace AMDiS
...@@ -5,6 +5,7 @@ ...@@ -5,6 +5,7 @@
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh> #include <dune/geometry/type.hh>
#include <amdis/common/PolynomialDegree.hpp>
#include <amdis/gridfunctions/Order.hpp> #include <amdis/gridfunctions/Order.hpp>
namespace AMDiS namespace AMDiS
...@@ -39,27 +40,16 @@ 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 /// \brief Factory for quadrature rule, that calculates the coefficient order from
/// a localFunction passed to the bind method. /// a localFunction passed to the bind method.
template <class ctype, int dim, class LocalFunction> template <class ctype, int dim, class LocalFunction>
class QuadFactoryFromLocalFunction class QuadFactoryFromLocalFunction
: public QuadratureFactory<ctype, dim, LocalFunction> : 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: public:
void bind(LocalFunction const& localFct) final void bind(LocalFunction const& localFct) final
{ {
if constexpr (Concept::value) order_ = polynomialDegree(localFct).value();
order_ = AMDiS::order(localFct);
else
order_ = -1;
} }
int order() const final { return order_; } int order() const final { return order_; }
......
...@@ -4,6 +4,8 @@ ...@@ -4,6 +4,8 @@
#include <dune/functions/functionspacebases/powerbasis.hh> #include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh> #include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <amdis/common/PolynomialDegree.hpp>
#include <amdis/functions/Order.hpp>
#include <amdis/typetree/FiniteElementType.hpp> #include <amdis/typetree/FiniteElementType.hpp>
#include "Tests.hpp" #include "Tests.hpp"
...@@ -36,14 +38,14 @@ int main() ...@@ -36,14 +38,14 @@ int main()
localView.bind(e); localView.bind(e);
auto node = localView.tree(); 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)>>, "" ); static_assert( std::is_same_v<tag::unknown, FiniteElementType_t<decltype(node)>>, "" );
auto v_node = TypeTree::child(node, _0); 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); 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(); return report_errors();
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment