From e54ba68f61a094d79214d36c5aba6fa9d105f86e Mon Sep 17 00:00:00 2001 From: Simon Praetorius <simon.praetorius@tu-dresden.de> Date: Mon, 13 Nov 2017 17:32:44 +0100 Subject: [PATCH] Concepts redesigned and FieldVector/FieldMatrix operations added --- dune/amdis/Assembler.hpp | 13 +- dune/amdis/DirichletBC.hpp | 6 +- dune/amdis/FirstOrderAssembler.hpp | 8 +- dune/amdis/Operator.hpp | 23 ++- dune/amdis/OperatorTermBase.hpp | 89 +++++---- dune/amdis/SecondOrderAssembler.hpp | 8 +- dune/amdis/ZeroOrderAssembler.hpp | 8 +- dune/amdis/common/Concepts.hpp | 149 +++++++++++++++ dune/amdis/common/ConceptsBase.hpp | 188 ++++--------------- dune/amdis/common/FieldMatVec.hpp | 269 ++++++++++++++++++++++++++++ dune/amdis/common/Operations.hpp | 215 ++++++++++++++++++++++ dune/amdis/common/Traits.hpp | 121 ------------- dune/amdis/terms/CoordsTerm.hpp | 2 +- dune/amdis/terms/DOFVectorTerm.hpp | 1 - dune/amdis/terms/FunctorTerm.hpp | 4 +- dune/amdis/terms/TermGenerator.hpp | 2 +- dune/amdis/utility/GetEntity.hpp | 13 +- 17 files changed, 772 insertions(+), 347 deletions(-) create mode 100644 dune/amdis/common/Concepts.hpp create mode 100644 dune/amdis/common/FieldMatVec.hpp create mode 100644 dune/amdis/common/Operations.hpp delete mode 100644 dune/amdis/common/Traits.hpp diff --git a/dune/amdis/Assembler.hpp b/dune/amdis/Assembler.hpp index a812338c..227a967a 100644 --- a/dune/amdis/Assembler.hpp +++ b/dune/amdis/Assembler.hpp @@ -6,24 +6,25 @@ namespace AMDiS { - template <class GridView, class Element> + // the LocalContext is either an Codim=0-EntityType or an IntersectionType + template <class GridView, class LocalContext> class Assembler { static constexpr int dim = GridView::dimension; using ctype = typename GridView::ctype; - using QuadratureRule = QuadratureRuleFactory_t<Element, ctype, dim>; - using Geometry = typename Impl::Get<Element>::Geometry; + using QuadratureRule = QuadratureRuleFactory_t<LocalContext, ctype, dim>; + using Geometry = typename Impl::Get<LocalContext>::Geometry; public: /// Constructor, initializes the geometry of the element and a /// quadrature-rule wrapper template <class Operator> - Assembler(Operator& op, Element const& element, int degree, FirstOrderType type = GRD_PHI) + Assembler(Operator& op, LocalContext const& element, int degree, FirstOrderType type = GRD_PHI) : geometry(get_geometry(element)) { int order = op.getQuadratureDegree(geometry.type(), geometry, degree, type); - auto const& quad = Dune::QuadratureRules<ctype, Element::mydimension>::rule(geometry.type(), order); + auto const& quad = Dune::QuadratureRules<ctype, LocalContext::mydimension>::rule(geometry.type(), order); quadrature.reset(new QuadratureRule(element, quad)); } @@ -56,7 +57,7 @@ namespace AMDiS template <class Coordinate> decltype(auto) map(Coordinate const& p) const { - return get_position<Element>(geometry, p); + return get_position<LocalContext>(geometry, p); } private: diff --git a/dune/amdis/DirichletBC.hpp b/dune/amdis/DirichletBC.hpp index b9e3d910..ab301bcb 100644 --- a/dune/amdis/DirichletBC.hpp +++ b/dune/amdis/DirichletBC.hpp @@ -4,7 +4,7 @@ #include <type_traits> #include <vector> -#include <dune/amdis/common/ConceptsBase.hpp> +#include <dune/amdis/common/Concepts.hpp> #include <dune/amdis/Output.hpp> @@ -30,8 +30,8 @@ namespace AMDiS { public: template <class Predicate, class Values, - std::enable_if_t< Concepts::Functor<Predicate, bool(WorldVector)> && - Concepts::Functor<Values, double(WorldVector)>, int*> = nullptr> + REQUIRES(Concepts::Functor<Predicate, bool(WorldVector)> && + Concepts::Functor<Values, double(WorldVector)>) > DirichletBC(Predicate&& predicate, Values&& values) : predicate(std::forward<Predicate>(predicate)) , values(std::forward<Values>(values)) diff --git a/dune/amdis/FirstOrderAssembler.hpp b/dune/amdis/FirstOrderAssembler.hpp index b6621351..20467769 100644 --- a/dune/amdis/FirstOrderAssembler.hpp +++ b/dune/amdis/FirstOrderAssembler.hpp @@ -9,11 +9,11 @@ namespace AMDiS { - template <class GridView, class Element, FirstOrderType type> + template <class GridView, class LocalContext, FirstOrderType type> class FirstOrderAssembler - : public Assembler<GridView, Element> + : public Assembler<GridView, LocalContext> { - using Super = Assembler<GridView, Element>; + using Super = Assembler<GridView, LocalContext>; static constexpr int dim = GridView::dimension; static constexpr int dow = GridView::dimensionworld; @@ -23,7 +23,7 @@ namespace AMDiS public: template <class Operator> - FirstOrderAssembler(Operator& op, Element const& element) + FirstOrderAssembler(Operator& op, LocalContext const& element) : Super(op, element, 1, type) { if (type == GRD_PHI) diff --git a/dune/amdis/Operator.hpp b/dune/amdis/Operator.hpp index 32698ba6..af42a851 100644 --- a/dune/amdis/Operator.hpp +++ b/dune/amdis/Operator.hpp @@ -16,11 +16,13 @@ namespace AMDiS { - template <class MeshView, class Element> + + // the LocalContext is either an Codim=0-EntityType or an IntersectionType + template <class MeshView, class LocalContext> class Operator { using Self = Operator; - using OperatorTermType = OperatorTerm<Element>; + using OperatorTermType = OperatorTerm<LocalContext>; using ElementVector = Impl::ElementVector; using ElementMatrix = Impl::ElementMatrix; @@ -32,8 +34,13 @@ namespace AMDiS template <class, class> friend class SecondOrderAssembler; public: - /// Add coefficients for zero-order operator < c * u, v >. - /// The coefficient \p c must be a scalar expression. + /// \brief Add coefficients for zero-order operator < C(phi), psi >. + /** + * If `c` is a scalar expression, then C(phi) := c * phi, + * otherwise if `c` is derived from \ref ZeroOrderOperatorTerm, it is evaluated + * using its `evalZot` method at quadrature points and represents the local + * evaluation of the scalar product `c = C(phi(lambda))*v(lambda)` with lambda a local coordinate. + **/ template <class Term> Self& addZOT(Term const& c) { @@ -49,8 +56,8 @@ namespace AMDiS } - /// Add coefficients for first-order operator < psi, term * grad(phi) > or - /// < grad(psi), term * phi >, where the first corresponds to + /// Add coefficients for first-order operator < psi, b * grad(phi) > or + /// < grad(psi), b * phi >, where the first corresponds to /// \p type == GRD_PHI and the second one to \p type == GRD_PSI. /// The coefficient \p b must be a vector expression. template <class Term> @@ -122,14 +129,14 @@ namespace AMDiS template <class RowView, class ColView> - bool getElementMatrix(Element const& element, + bool getElementMatrix(LocalContext const& element, RowView const& rowView, ColView const& colView, ElementMatrix& elementMatrix, double* factor = NULL); template <class RowView> - bool getElementVector(Element const& element, + bool getElementVector(LocalContext const& element, RowView const& rowView, ElementVector& elementVector, double* factor = NULL); diff --git a/dune/amdis/OperatorTermBase.hpp b/dune/amdis/OperatorTermBase.hpp index cd1a6164..3989dc58 100644 --- a/dune/amdis/OperatorTermBase.hpp +++ b/dune/amdis/OperatorTermBase.hpp @@ -16,108 +16,137 @@ namespace AMDiS { /// Abstract base class of all operator terms - template <class Element> + template <class LocalContext> class OperatorTerm { protected: - static constexpr int dim = Element::dimension; - static constexpr int dow = Element::Geometry::coorddimension; + static constexpr int dim = LocalContext::dimension; + static constexpr int dow = LocalContext::Geometry::coorddimension; - using QuadratureRule = QuadratureRuleFactory_t<Element,double,dim>; + using QuadratureRule = QuadratureRuleFactory_t<LocalContext,double,dim>; public: - // operator-terms, when used in operators, must implement this method - virtual void init(Element const& element, QuadratureRule const& points) = 0; + // initialize operator-term on the current element + virtual void init(LocalContext const& element, QuadratureRule const& points) = 0; + // evaluate operator-term at quadrature point as zero-order term, + // with no gradients involved virtual double evalZot(std::size_t iq, Dune::FieldVector<double,1> const& test, Dune::FieldVector<double,1> const trial = 1.0) const = 0; + // evaluate operator-term at quadrature point as first-order term, + // with gradients on the trial-function virtual double evalFot1(std::size_t iq, Dune::FieldVector<double,1> const& test, Dune::FieldVector<double,dow> const& grad_trial) const = 0; + // evaluate operator-term at quadrature point as first-order term, + // with gradients on the test-function virtual double evalFot2(std::size_t iq, Dune::FieldVector<double,dow> const& grad_test, Dune::FieldVector<double,1> const trial = 1.0) const = 0; + // evaluate operator-term at quadrature point as second-order term, + // with gradients on the trial and test-function virtual double evalSot(std::size_t iq, Dune::FieldVector<double,dow> const& grad_test, Dune::FieldVector<double,dow> const& grad_trial) const = 0; + // return the polynomial degree, necessary to intergrate the operator-term + // on an element accurately virtual int getDegree(Dune::GeometryType const& t) const = 0; }; - - /// Base class for all operators based on expressions - template <class Element, class Term, class Traits = tag::none> + /// \brief Base class for all operator-terms based on expressions + /** + * Represents an operator-term with a coefficient the scales the classical operator, e.g. + * - zero order terms `<expr * u, v>` + * - first order terms (GRD_PHI) `<(vecexpr * grad(u)), v>`, `<expr * d_i(u), v>` (a) + * - first order terms (GRD_PSI) `<vecexpr * u, grad(v)>`, `<expr * u, d_i(v)>` (a) + * - second order terms `<expr * grad(u), grad(v)>`, `<(matexpr *grad(u)), grad(v)>`, `<expr * d_i(u), d_j(v)>` (b) + **/ + // LocalContext ... The type of the element where to evaluate, either entity (codim=0), or intersection (codim=1) + // Expression ..... The type of the expression to evaluate + // Traits ......... Specialization for the operator evaluation [tag:none | VectorComponent (a) | MatrixComponent (b)] + template <class LocalContext, class Expression, class Traits = tag::none> class GenericOperatorTerm - : public OperatorTerm<Element> + : public OperatorTerm<LocalContext> { - using Super = OperatorTerm<Element>; + using Super = OperatorTerm<LocalContext>; using QuadratureRule = typename Super::QuadratureRule; - static constexpr int dim = Super::dim; static constexpr int dow = Super::dow; public: - GenericOperatorTerm(Term const& term, Traits traits = {}) - : term(term) + /// Constructor, stores a copy of the expression `expr` and of the traits `traits`. + GenericOperatorTerm(Expression const& expr, Traits traits = {}) + : expr(expr) , traits(traits) {} - virtual void init(Element const& element, QuadratureRule const& points) override + /// \brief Initialize the expression on the current (inside-)element, in the quadrature points. + /** + * This initialization does not only bind the operator-term to the inside-element, + * but also calculates values in all quadrature points. These pre-calculated values + * can be accessed using the `operator[]` method of the expression. + * + * Note: For entities the expression is bound to the `element` directly, but + * for intersections, the expression is bound to inside-entity of the element. + **/ + virtual void init(LocalContext const& element, QuadratureRule const& points) override { - term.init(get_entity(element), points); - - // cache term evaluation - // values.resize(points.size()); - // for (std::size_t iq = 0; iq < points.size(); ++iq) - // values[iq] = term[iq]; + expr.init(get_entity(element), points); } + /// Calculates `expr[iq] * test * trial` virtual double evalZot(std::size_t iq, Dune::FieldVector<double,1> const& test, Dune::FieldVector<double,1> const trial = 1.0) const override { - return Evaluate::zot(_cat{}, traits, term[iq], test, trial); + return Evaluate::zot(_cat{}, traits, expr[iq], test, trial); } + /// Calculates `expr[iq] * test * grad(trial)`, maybe partial derivative only. virtual double evalFot1(std::size_t iq, Dune::FieldVector<double,1> const& test, Dune::FieldVector<double,dow> const& grad_trial) const override { - return Evaluate::fot(_cat{}, traits, term[iq], grad_trial, test); + return Evaluate::fot(_cat{}, traits, expr[iq], grad_trial, test); } + /// Calculates `expr[iq] * grad(test) * trial`, maybe partial derivative only. virtual double evalFot2(std::size_t iq, Dune::FieldVector<double,dow> const& grad_test, Dune::FieldVector<double,1> const trial = 1.0) const override { - return Evaluate::fot(_cat{}, traits, term[iq], grad_test, trial); + return Evaluate::fot(_cat{}, traits, expr[iq], grad_test, trial); } + /// Calculates `expr[iq] * grad(test) * grad(trial)`, maybe partial derivatives only. virtual double evalSot(std::size_t iq, Dune::FieldVector<double,dow> const& grad_test, Dune::FieldVector<double,dow> const& grad_trial) const override { - return Evaluate::sot(_cat{}, traits, term[iq], grad_test, grad_trial); + return Evaluate::sot(_cat{}, traits, expr[iq], grad_test, grad_trial); } + /// Returns the polynomial degree of the expression virtual int getDegree(Dune::GeometryType const& t) const override { - return term.getDegree(t); + return expr.getDegree(t); } private: - Term term; + /// A local expression, implementing the Concept \ref Concepts::LocalExpression. + Expression expr; + + /// One of [tag::none, VectorComponent, MatrixComponent] Traits traits; - using value_type = std::decay_t< decltype( std::declval<Term>()[std::size_t(0)] ) >; + using value_type = std::decay_t< decltype( std::declval<Expression>()[std::size_t(0)] ) >; using _cat = ValueCategory_t<value_type>; - - // std::vector<value_type> values; // NOTE: maybe caching is not necessary here, since cached already in the term }; } // end namespace AMDiS diff --git a/dune/amdis/SecondOrderAssembler.hpp b/dune/amdis/SecondOrderAssembler.hpp index 48bd00f9..15c62aa1 100644 --- a/dune/amdis/SecondOrderAssembler.hpp +++ b/dune/amdis/SecondOrderAssembler.hpp @@ -9,11 +9,11 @@ namespace AMDiS { - template <class GridView, class Element> + template <class GridView, class LocalContext> class SecondOrderAssembler - : public Assembler<GridView, Element> + : public Assembler<GridView, LocalContext> { - using Super = Assembler<GridView, Element>; + using Super = Assembler<GridView, LocalContext>; static constexpr int dim = GridView::dimension; static constexpr int dow = GridView::dimensionworld; @@ -23,7 +23,7 @@ namespace AMDiS public: template <class Operator> - SecondOrderAssembler(Operator& op, Element const& element) + SecondOrderAssembler(Operator& op, LocalContext const& element) : Super(op, element, 2) { op.initSot(element, Super::getQuadrature()); diff --git a/dune/amdis/ZeroOrderAssembler.hpp b/dune/amdis/ZeroOrderAssembler.hpp index 6d0e6af8..94327d40 100644 --- a/dune/amdis/ZeroOrderAssembler.hpp +++ b/dune/amdis/ZeroOrderAssembler.hpp @@ -10,11 +10,11 @@ namespace AMDiS { /// Assembler for zero-order operators - template <class GridView, class Element> + template <class GridView, class LocalContext> class ZeroOrderAssembler - : public Assembler<GridView, Element> + : public Assembler<GridView, LocalContext> { - using Super = Assembler<GridView, Element>; + using Super = Assembler<GridView, LocalContext>; static constexpr int dim = GridView::dimension; @@ -23,7 +23,7 @@ namespace AMDiS public: template <class Operator> - ZeroOrderAssembler(Operator& op, Element const& element) + ZeroOrderAssembler(Operator& op, LocalContext const& element) : Super(op, element, 0) { op.initZot(element, Super::getQuadrature()); diff --git a/dune/amdis/common/Concepts.hpp b/dune/amdis/common/Concepts.hpp new file mode 100644 index 00000000..a185e259 --- /dev/null +++ b/dune/amdis/common/Concepts.hpp @@ -0,0 +1,149 @@ +#pragma once + +#include <type_traits> + +#include <dune/functions/common/functionconcepts.hh> + +#include <dune/amdis/common/ConceptsBase.hpp> + +namespace AMDiS +{ + /** + * \defgroup concept Concepts + * \brief Concept definitions + * @{ + **/ + namespace Concepts + { +#ifndef DOXYGEN + namespace Definition + { + // a < b + struct LessThanComparable + { + template <class A, class B> + auto requires_(A&& a, B&& b) -> decltype( a < b ); + }; + + // a + b + struct Addable + { + template <class A, class B> + auto requires_(A&& a, B&& b) -> decltype( + Concepts::valid_expr( + a + b, + b + a + )); + }; + + // a - b + struct Subtractable + { + template <class A, class B> + auto requires_(A&& a, B&& b) -> decltype( a - b ); + }; + + // a * b + struct Multiplicable + { + template <class A, class B> + auto requires_(A&& a, B&& b) -> decltype( + Concepts::valid_expr( + a * b, + b * a + )); + }; + + // a / b + struct Divisible + { + template <class A, class B> + auto requires_(A&& a, B&& b) -> decltype( a / b ); + }; + + // f(args...) + struct Callable + { + template <class F, class... Args> + auto requires_(F&& f, Args&&... args) -> decltype( f(std::forward<Args>(args)...)); + }; + + } // end namespace Definition + + namespace traits + { + template <class A, class B> + struct is_compatible + : std::is_same<std::decay_t<A>, std::decay_t<B>> {}; + + template <class A, class B> + struct is_compatible<Types<A>, Types<B>> + : is_compatible<A,B> {}; + + template <> + struct is_compatible<Types<>, Types<>> : std::true_type {}; + + template <class A0, class... As, class B0, class... Bs> + struct is_compatible<Types<A0,As...>, Types<B0,Bs...>> + : and_t<is_compatible<A0, B0>::value, is_compatible<Types<As...>, Types<Bs...>>::value> {}; + + } // end namespace traits + + +#endif // DOXYGEN + + /// \brief Argument `A` id (implicitly) convertible to arguemnt `B` and vice versa. + template <class A, class B > + constexpr bool Convertible = std::is_convertible<A,B>::value && std::is_convertible<B,A>::value; + + /// Types are the same, up to decay of qualifiers + template <class A, class B> + constexpr bool Compatible = traits::is_compatible<A, B>::value; + + /// \brief Argument A is less-than comparable to type B, i.e. A < B is valid + template <class A, class B = A> + constexpr bool LessThanComparable = models<Definition::LessThanComparable(A, B)>; + + /// \brief Argument A is addable to type B, i.e. A + B is valid + template <class A, class B> + constexpr bool Addable = models<Definition::Addable(A, B)>; + + /// \brief Argument A is subtractable by type B, i.e. A - B is valid + template <class A, class B> + constexpr bool Subtractable = models<Definition::Subtractable(A, B)>; + + /// \brief Argument A is multiplicable by type B, i.e. A * B is valid + template <class A, class B> + constexpr bool Multiplicable = models<Definition::Multiplicable(A, B)>; + + /// \brief Argument A is divisible by type B, i.e. A / B is valid + template <class A, class B> + constexpr bool Divisible = models<Definition::Divisible(A, B)>; + + /// \brief A Collable is a function `F` that can be called with arguments of type `Args...`. + /** + * To be used as follows: `Concepts::Collable<F, Args...>`. Returns true, if + * an instance of `F` can be called by `operator()` with arguments of type + * `Args...`. + **/ + template <class F, class... Args> + constexpr bool Callable = models<Definition::Callable(F, Args...)>; + + /// \brief A Functor is a function `F` with signature `Signature`. + /** + * To be used as follows: `Concepts::Functor<F, R(Args...)>`. Returns true, if + * an instance of `F` can be called by `operator()` with arguments of type + * `Args...` and returns a value of type `R`, i.e. `Signature := R(Args...)`. + **/ + template <class F, class Signature> // F, Signature=Return(Arg) + constexpr bool Functor = Dune::Functions::Concept::isFunction<F, Signature>(); + + /// A predicate is a function that returns a boolean. + template <class F, class... Args> + constexpr bool Predicate = Functor<F, bool(Args...)>; + + /** @} **/ + + } // end namespace Concepts + +} // end namespace AMDiS diff --git a/dune/amdis/common/ConceptsBase.hpp b/dune/amdis/common/ConceptsBase.hpp index 9cfc4e7a..026ce0ba 100644 --- a/dune/amdis/common/ConceptsBase.hpp +++ b/dune/amdis/common/ConceptsBase.hpp @@ -1,173 +1,59 @@ -/** \file concepts_base.hpp */ - #pragma once -// std c++ headers -#include <utility> #include <type_traits> -// AMDiS headers -#include <dune/amdis/common/Mpl.hpp> -#include <dune/amdis/common/Utility.hpp> - -// macro to generate concept-checks -#define HAS_MEMBER_GENERATE(name) \ - template <class F, class... Args> \ - inline auto invoke_member_ ## name(F&& f, Args&&... args) -> \ - decltype(std::forward<F>(f).name (std::forward<Args>(args)...)) { \ - return std::forward<F>(f).name (std::forward<Args>(args)...); \ - } \ - inline no_valid_type invoke_member_ ## name(...); \ - template <class, class> struct has_member_ ## name : std::false_type {}; \ - template <class F, class Return, class... Args> \ - struct has_member_ ## name <F, Return(Args...)> \ - : std::is_convertible< \ - decltype(invoke_member_ ## name(std::declval<F>(), std::declval<Args>()...)), \ - Return > {}; - -#define HAS_MEMBER(name) has_member_ ## name - -#include <dune/functions/common/functionconcepts.hh> +#if defined(DOXYGEN) + #define REQUIRES(...) + #define CONCEPT constexpr +#else + #define REQUIRES(...) std::enable_if_t<__VA_ARGS__ , int>* = nullptr + #define CONCEPT constexpr +#endif namespace AMDiS { - namespace traits + namespace Impl { - // Some type traits to test for type-attributes of a class - // ------------------------------------------------------------------------- - - namespace Impl - { - template <class T, class = typename T::value_type > - std::true_type HasValueTypeImpl(int); - - template <class T> std::false_type HasValueTypeImpl(...); - - - template <class T, class = typename T::size_type > - std::true_type HasSizeTypeImpl(int); - - template <class T> std::false_type HasSizeTypeImpl(...); - - - template <class T, class = typename T::result_type > - std::true_type HasResultTypeImpl(int); - - template <class T> std::false_type HasResultTypeImpl(...); - - } // end namespace Impl - - - template <class T> - using HasValueType = decltype( Impl::HasValueTypeImpl<T>(int{}) ); - - template <class T> - using HasSizeType = decltype( Impl::HasSizeTypeImpl<T>(int{}) ); - - template <class T> - using HasResultType = decltype( Impl::HasResultTypeImpl<T>(int{}) ); - - - // ------------------------------------------------------------------------- + // Workaround for MSVC (problems with alias templates in pack expansion) + template <class... Args> + struct VoidType { using type = void; }; + } + template <class... Args> + using Void_t = typename Impl::VoidType<Args...>::type; - namespace Impl - { - template <class F, class Arg> - inline auto invoke_vector(F&& f, Arg&& arg) - -> decltype( std::forward<F>(f)[std::forward<Arg>(arg)] ); - - inline no_valid_type invoke_vector(...); - - } // end namespace Impl - - template <class, class> - struct IsVector : std::false_type {}; - - template <class F, class Return, class Arg> - struct IsVector<F, Return(Arg)> - : std::is_same< decltype( Impl::invoke_vector(std::declval<std::decay_t<F>>(), std::declval<Arg>()) ), - Return > {}; - - template <class, class> - struct IsVectorWeak : std::false_type {}; - - template <class F, class Return, class Arg> - struct IsVectorWeak<F, Return(Arg)> - : std::is_convertible< decltype( Impl::invoke_vector(std::declval<std::decay_t<F>>(), std::declval<Arg>()) ), - Return > {}; - - - namespace Impl - { - template <class T, class = decltype(&T::push_back())> - std::true_type PushBackImpl(int); - template <class T> std::false_type PushBackImpl(...); - - } // end namespace Impl - - template <class T> - using HasPushBack = decltype(Impl::PushBackImpl<T>(int{})); - - } // end namespace traits - - - /// Namespace that implements basic concept checks namespace Concepts { - /** - * \defgroup Concepts Concepts definitions - * @{ - **/ - - /// A Collable is a function `F` that can be called with arguments of type `Args...`. - /** - * To be used as follows: `Concepts::Collable<F, Args...>`. Returns true, if - * an instance of `F` can be called by `operator()` with arguments of type - * `Args...`. - **/ - template <class F, class... Args> - constexpr bool Callable = Dune::Functions::Concept::isCallable<F, Args...>(); - - - /// A Functor is a function `F` with signature `Signature`. - /** - * To be used as follows: `Concepts::Functor<F, R(Args...)>`. Returns true, if - * an instance of `F` can be called by `operator()` with arguments of type - * `Args...` and returns a value of type `R`, i.e. `Signature := R(Args...)`. - **/ - template <class F, class Signature> // F, Signature=Return(Arg) - constexpr bool Functor = Dune::Functions::Concept::isFunction<F, Signature>(); - + namespace _aux + { + template <class Concept, class = Void_t<>> + struct models + : std::false_type + {}; - /// A predicate is a function that returns a boolean. - template <class F, class... Args> - constexpr bool Predicate = Functor<F, bool(Args...)>; + template <class Concept, class... Ts> + struct models<Concept(Ts...), Void_t< decltype(std::declval<Concept>().requires_(std::declval<Ts>()...)) >> + : std::true_type + {}; + struct valid_expr + { + template <class... Ts> + void operator()(Ts&&...) const; - /// Test whether a type `T` supports vector-access `vector[std::size_t]`. - /** - * To be used as follows: `Concepts::Vector<V>`. Returns true, if `V` has a - * value type, accessible by `Value_t<V>` and an instance of `V` can be called - * by square brackets `operator[]` with an index of type `std::size_t`. The return - * value of this element access should be convertible to `Value_t<V>`. - **/ - template <class T> - constexpr bool Vector = traits::IsVectorWeak<T, Value_t<T>(std::size_t)>::value; +#if defined(__GNUC__) && !defined(__clang__) + template <class... Ts> + void operator()(Ts const&...) const; +#endif + }; + } // end namespace _aux - /// Test whether a type `T` supports matrix-access `matrix(std::size_t, std::size_t)`. - /** - * To be used as follows: `Concepts::Matrix<M>`. Returns true, if `M` has a - * value type, accessible by `Value_t<M>` and an instance of `M` can be called - * by round brackets `operator()` with two indices of type `std::size_t`. The return - * value of this element access should be convertible to `Value_t<V>`. - **/ - template <class T> - constexpr bool Matrix = Functor<T, Value_t<T>(std::size_t, std::size_t)>; + template <class Concept> + constexpr bool models = _aux::models<Concept>::value; - /** @} **/ + constexpr _aux::valid_expr valid_expr = {}; } // end namespace Concepts diff --git a/dune/amdis/common/FieldMatVec.hpp b/dune/amdis/common/FieldMatVec.hpp new file mode 100644 index 00000000..cf955e78 --- /dev/null +++ b/dune/amdis/common/FieldMatVec.hpp @@ -0,0 +1,269 @@ +#pragma once + +#include <dune/common/fmatrix.hh> +#include <dune/common/fvector.hh> + +#include <dune/amdis/Math.hpp> +#include <dune/amdis/common/Concepts.hpp> +#include <dune/amdis/common/Operations.hpp> +#include <dune/amdis/common/ScalarTypes.hpp> + +namespace AMDiS +{ + // some arithmetic operations with Dune::FieldVector + + template <class T, int N, class S, + REQUIRES(Concepts::Arithmetic<S>) > + Dune::FieldVector<T,N> operator*(Dune::FieldVector<T,N> v, S factor) + { + return v *= factor; + } + + template <class S, class T, int N, + REQUIRES(Concepts::Arithmetic<S>) > + Dune::FieldVector<T,N> operator*(S factor, Dune::FieldVector<T,N> v) + { + return v *= factor; + } + + template <class T, int N, class S, + REQUIRES(Concepts::Arithmetic<S>) > + Dune::FieldVector<T,N> operator/(Dune::FieldVector<T,N> v, S factor) + { + return v /= factor; + } + + // some arithmetic operations with Dune::FieldMatrix + + template <class T, int N, int M> + Dune::FieldMatrix<T,N,M> operator+(Dune::FieldMatrix<T,N,M> A, Dune::FieldMatrix<T,N,M> const& B) + { + return A += B; + } + + template <class T, int N, int M> + Dune::FieldMatrix<T,N,M> operator-(Dune::FieldMatrix<T,N,M> A, Dune::FieldMatrix<T,N,M> const& B) + { + return A -= B; + } + + // ---------------------------------------------------------------------------- + + /// Cross-product a 2d-vector = orthogonal vector + template <class T> + Dune::FieldVector<T, 2> cross(Dune::FieldVector<T, 2> const& a) + { + return {{ a[1], -a[0] }}; + } + + /// Cross-product of two vectors (in 3d only) + template <class T> + Dune::FieldVector<T, 3> cross(Dune::FieldVector<T, 3> const& a, Dune::FieldVector<T, 3> const& b) + { + return {{ a[1]*b[2] - a[2]*b[1], + a[2]*b[0] - a[0]*b[2], + a[0]*b[1] - a[1]*b[0] }}; + } + + /// Dot product (vec1^T * vec2) + template <class T, class S, int N> + auto dot(Dune::FieldVector<T,N> const& vec1, Dune::FieldVector<S,N> const& vec2) + { + return vec1.dot(vec2); + } + + // ---------------------------------------------------------------------------- + + namespace Impl + { + template <class T, int N, class Operation> + T accumulate(Dune::FieldVector<T, N> const& x, Operation op) + { + T result = 0; + for (int i = 0; i < N; ++i) + result = op(result, x[i]); + return result; + } + + } // end namespace Impl + + /// Sum of vector entires. + template <class T, int N> + T sum(Dune::FieldVector<T, N> const& x) + { + return Impl::accumulate(x, Operation::plus{}); + } + + + /// Dot-product with the vector itself + template <class T, int N> + auto unary_dot(Dune::FieldVector<T, N> const& x) + { + auto op = [](auto const& a, auto const& b) { return a + sqr(std::abs(b)); }; + return Impl::accumulate(x, op); + } + + /// Maximum over all vector entries + template <class T, int N> + auto max(Dune::FieldVector<T, N> const& x) + { + return Impl::accumulate(x, Operation::maximum{}); + } + + /// Minimum over all vector entries + template <class T, int N> + auto min(Dune::FieldVector<T, N> const& x) + { + return Impl::accumulate(x, Operation::minimum{}); + } + + /// Maximum of the absolute values of vector entries + template <class T, int N> + auto abs_max(Dune::FieldVector<T, N> const& x) + { + return Impl::accumulate(x, Operation::abs_max{}); + } + + /// Minimum of the absolute values of vector entries + template <class T, int N> + auto abs_min(Dune::FieldVector<T, N> const& x) + { + return Impl::accumulate(x, Operation::abs_min{}); + } + + // ---------------------------------------------------------------------------- + + /** \ingroup vector_norms + * \brief The 1-norm of a vector = sum_i |x_i| + **/ + template <class T, int N> + auto one_norm(Dune::FieldVector<T, 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 ) + **/ + template <class T, int N> + auto two_norm(Dune::FieldVector<T, 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) + **/ + template <int p, class T, int N> + auto p_norm(Dune::FieldVector<T, 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 + **/ + template <class T, int N> + auto infty_norm(Dune::FieldVector<T, N> const& x) + { + return abs_max(x); + } + + // ---------------------------------------------------------------------------- + + /// The euklidean distance between two vectors = |lhs-rhs|_2 + template <class T, int N> + T distance(Dune::FieldVector<T, N> const& lhs, Dune::FieldVector<T, N> const& rhs) + { + T result = 0; + for (int i = 0; i < N; ++i) + result += sqr(lhs[i] - rhs[i]); + return std::sqrt(result); + } + + // ---------------------------------------------------------------------------- + + /// Outer product (vec1 * vec2^T) + template <class T, class S, int N, int M, int K> + auto outer(Dune::FieldMatrix<T,N,K> const& vec1, Dune::FieldMatrix<S,M,K> const& vec2) + { + using result_type = Dune::FieldMatrix<decltype( std::declval<T>() * std::declval<S>() ), N, M>; + result_type mat; + for (int i = 0; i < N; ++i) + for (int j = 0; j < M; ++j) + mat[i][j] = vec1[i].dot(vec2[j]); + return mat; + } + + // ---------------------------------------------------------------------------- + + template <class T> + T det(Dune::FieldMatrix<T, 0, 0> const& /*mat*/) + { + return 0; + } + + /// Determinant of a 1x1 matrix + template <class T> + T det(Dune::FieldMatrix<T, 1, 1> const& mat) + { + return mat[0][0]; + } + + /// Determinant of a 2x2 matrix + template <class T> + T det(Dune::FieldMatrix<T, 2, 2> const& mat) + { + return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0]; + } + + /// Determinant of a 3x3 matrix + template <class T> + T det(Dune::FieldMatrix<T, 3, 3> const& mat) + { + return mat[0][0]*mat[1][1]*mat[2][2] + mat[0][1]*mat[1][2]*mat[2][0] + mat[0][2]*mat[1][0]*mat[2][1] + - (mat[0][2]*mat[1][1]*mat[2][0] + mat[0][1]*mat[1][0]*mat[2][2] + mat[0][0]*mat[1][2]*mat[2][1]); + } + + /// Determinant of a NxN matrix + template <class T, int N> + T det(Dune::FieldMatrix<T, N, N> const& mat) + { + return mat.determinant(); + } + + /// Return the inverse of the matrix `mat` + template <class T, int N> + auto inv(Dune::FieldMatrix<T, N, N> mat) + { + mat.invert(); + return mat; + } + + /// Solve the linear system A*x = b + template <class T, int N> + void solve(Dune::FieldMatrix<T, N, N> const& A, Dune::FieldVector<T, N>& x, Dune::FieldVector<T, N> const& b) + { + A.solve(x, b); + } + + + /// Gramian determinant = sqrt( det( DT^T * DF ) ) + template <class T, int N, int M> + T gramian(Dune::FieldMatrix<T,N,M> const& DF) + { + using std::sqrt; + return sqrt( det(outer(DF, DF)) ); + } + + /// Gramian determinant, specialization for 1 column matrices + template <class T, int M> + T gramian(Dune::FieldMatrix<T, 1, M> const& DF) + { + using std::sqrt; + return sqrt(dot(DF[0], DF[0])); + } + +} // end namespace AMDiS diff --git a/dune/amdis/common/Operations.hpp b/dune/amdis/common/Operations.hpp new file mode 100644 index 00000000..77e20d20 --- /dev/null +++ b/dune/amdis/common/Operations.hpp @@ -0,0 +1,215 @@ +#pragma once + +#include <dune/amdis/Math.hpp> +#include <dune/amdis/common/Concepts.hpp> +#include <dune/amdis/common/ScalarTypes.hpp> + +namespace AMDiS +{ + namespace Operation + { + /** \defgroup operations Operations + * \brief Functors, representing unary/binary operations used in expressions. + * @{ + **/ + + struct assign + { + template <class T, class S> + constexpr void operator()(T const& from, S& to) const + { + to = from; + } + }; + + /// Operation that represents A+B + struct plus + { + template <class T, class S> + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return lhs + rhs; + } + }; + + struct plus_assign + { + template <class T, class S> + constexpr void operator()(T const& from, S& to) const + { + to += from; + } + }; + + /// Operation that represents A-B + struct minus + { + template <class T, class S> + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return lhs - rhs; + } + }; + + struct minus_assign + { + template <class T, class S> + constexpr void operator()(T const& from, S& to) const + { + to -= from; + } + }; + + /// Operation that represents A*B + struct times + { + template <class T, class S> + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return lhs * rhs; + } + }; + + struct times_assign + { + template <class T, class S> + constexpr void operator()(T const& from, S& to) const + { + to *= from; + } + }; + + /// Operation that represents A/B + struct divides + { + template <class T, class S> + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return lhs / rhs; + } + }; + + struct divides_assign + { + template <class T, class S> + constexpr void operator()(T const& from, S& to) const + { + to /= from; + } + }; + + /// Operation that represents A-B + struct negate + { + template <class T> + constexpr auto operator()(T const& x) const + { + return -x; + } + }; + + /// Operation that represents max(A,B) + struct maximum + { + template <class T, class S> + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return Math::max(lhs, rhs); + } + }; + + /// Operation that represents min(A,B) + struct minimum + { + template <class T, class S> + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return Math::min(lhs, rhs); + } + }; + + /// Operation that represents max(|A|,|B|) + struct abs_max + { + template <class T, class S> + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return Math::max(std::abs(lhs), std::abs(rhs)); + } + }; + + /// Operation that represents min(|A|,|B|) + struct abs_min + { + template <class T, class S> + constexpr auto operator()(T const& lhs, S const& rhs) const + { + return Math::min(std::abs(lhs), std::abs(rhs)); + } + }; + + /// Operation that represents A*factor + template <class S, class = void> + struct times_scalar; + + template <class S> + struct times_scalar<S, std::enable_if_t<Concepts::Arithmetic<S>> > + { + constexpr explicit times_scalar(S const& scalar) : scalar(scalar) {} + + template <class T, + REQUIRES(Concepts::Multiplicable<T,S>) > + constexpr auto operator()(T const& lhs) const + { + return lhs * scalar; + } + + private: + S scalar; + }; + + /// Operation that represents A/factor + template <class S, class = void> + struct divides_scalar; + + template <class S> + struct divides_scalar<S, std::enable_if_t<Concepts::Arithmetic<S>> > + { + constexpr explicit divides_scalar(S const& scalar) : scalar(scalar) {} + + template <class T, + REQUIRES(Concepts::Divisible<T,S>) > + constexpr auto operator()(T const& lhs) const + { + return lhs / scalar; + } + + private: + S scalar; + }; + + + /// Operation that represents conj(A) or A, if A is complex or not, respectively + template <class T, class = void> + struct hermitian + { + constexpr T operator()(T const& x) const + { + return std::conj(x); + } + }; + + template <class T> + struct hermitian<T, std::enable_if_t<Concepts::Arithmetic<T>> > + { + constexpr T const& operator()(T const& x) const + { + return x; + } + }; + + /** @} **/ + + } // end namespace Operation + +} // end namespace AMDiS diff --git a/dune/amdis/common/Traits.hpp b/dune/amdis/common/Traits.hpp deleted file mode 100644 index f990e96e..00000000 --- a/dune/amdis/common/Traits.hpp +++ /dev/null @@ -1,121 +0,0 @@ -/** \file traits.hpp */ - -#pragma once - -// std c++ headers -#include <type_traits> - -// AMDiS headers -#include <dune/amdis/common/ConceptsBase.hpp> -#include <dune/amdis/common/Mpl.hpp> -#include <dune/amdis/common/Utility.hpp> - -namespace AMDiS -{ - // some traits to test for binary operators on types - namespace traits - { - namespace Impl - { - struct PlusOp - { - template <class T, class U> - auto operator()(T t, U u) const { return t + u; } - }; - - struct MultipliesOp - { - template <class T, class U> - auto operator()(T t, U u) const { return t * u; } - }; - - } // end namespace Impl - - template <class T, class U = T> - using IsAddable = bool_t<Concepts::Callable<Impl::PlusOp, T, U>>; - - template <class T, class U = T> - using IsMultiplicable = bool_t<Concepts::Callable<Impl::MultipliesOp, T, U>>; - - - template <class T> - using IsTriviallyCopyable = std::is_pod<T>; - - // ------------------------------------------------------------------------- - - // larger types - template <class... Ts> - struct LargerType; - - template <class T1, class T2, class... Ts> - struct LargerType<T1, T2, Ts...> - { - using type = typename std::conditional_t< (sizeof(T1) > sizeof(T2)), - LargerType<T1,Ts...>, - LargerType<T2,Ts...> >::type; - }; - - template <class T1, class T2> - struct LargerType<T1, T2> - { - using type = std::conditional_t<(sizeof(T1) > sizeof(T2)), T1, T2>; - }; - - template <class T> struct LargerType<T> - { - using type = T; - }; - - // maximal size type - template <class... Es> - using MaxSizeType = typename LargerType<Size_t<Es>...>::type; - - // ------------------------------------------------------------------------- - - namespace Impl - { - template <class A, class B> - struct is_compatible - : std::is_same<Decay_t<A>, Decay_t<B>> {}; - - template <class A, class B> - struct is_compatible<Types<A>, Types<B>> - : is_compatible<A,B> {}; - - template <> - struct is_compatible<Types<>, Types<>> : std::true_type {}; - - template <class A0, class... As, class B0, class... Bs> - struct is_compatible<Types<A0,As...>, Types<B0,Bs...>> - : and_t<is_compatible<A0, B0>::value, is_compatible<Types<As...>, Types<Bs...>>::value> {}; - - } // end namespace Impl - - template <class A, class B> - using IsCompatible = Impl::is_compatible<A,B>; - - } // end namespace traits - - namespace Concepts - { - /** \addtogroup Concepts - * @{ - **/ - - /// Types support \f$ A + B \f$ - template <class A, class B> - constexpr bool Addable = traits::IsAddable<A, B>::value; - - /// Types support \f$ A \cdot B \f$ - template <class A, class B> - constexpr bool Multiplicable = traits::IsMultiplicable<A, B>::value; - - /// Types are the same, up to decay of qualifiers - template <class A, class B> - constexpr bool Compatible = traits::IsCompatible<A, B>::value; - - /** @} **/ - - } // end namespace Concepts - -} // end namespace AMDiS diff --git a/dune/amdis/terms/CoordsTerm.hpp b/dune/amdis/terms/CoordsTerm.hpp index 87161e1c..14acf991 100644 --- a/dune/amdis/terms/CoordsTerm.hpp +++ b/dune/amdis/terms/CoordsTerm.hpp @@ -10,7 +10,7 @@ #include <dune/istl/bvector.hh> #include <dune/amdis/OperatorTermBase.hpp> -#include <dune/amdis/common/Traits.hpp> +#include <dune/amdis/common/Concepts.hpp> namespace AMDiS { diff --git a/dune/amdis/terms/DOFVectorTerm.hpp b/dune/amdis/terms/DOFVectorTerm.hpp index d72cc1b5..5d5441c8 100644 --- a/dune/amdis/terms/DOFVectorTerm.hpp +++ b/dune/amdis/terms/DOFVectorTerm.hpp @@ -11,7 +11,6 @@ #include <dune/istl/bvector.hh> #include <dune/amdis/OperatorTermBase.hpp> -#include <dune/amdis/common/Traits.hpp> #include <dune/amdis/utility/GetDegree.hpp> namespace AMDiS diff --git a/dune/amdis/terms/FunctorTerm.hpp b/dune/amdis/terms/FunctorTerm.hpp index 5350f5fd..a75dbfc8 100644 --- a/dune/amdis/terms/FunctorTerm.hpp +++ b/dune/amdis/terms/FunctorTerm.hpp @@ -11,7 +11,7 @@ #include <dune/amdis/OperatorTermBase.hpp> #include <dune/amdis/common/IndexSeq.hpp> -#include <dune/amdis/common/Traits.hpp> +#include <dune/amdis/common/Concepts.hpp> namespace AMDiS { @@ -32,7 +32,7 @@ namespace AMDiS std::tuple<Terms...> terms; template <class F, class... Ts, - class = std::enable_if_t<Concepts::Compatible<Functor, F>> > + REQUIRES(Concepts::Compatible<Functor, F>) > FunctorTerm(F&& f, Ts&&... ts) : fct(std::forward<F>(f)) , terms(std::forward<Ts>(ts)...) diff --git a/dune/amdis/terms/TermGenerator.hpp b/dune/amdis/terms/TermGenerator.hpp index d61781ec..044db05a 100644 --- a/dune/amdis/terms/TermGenerator.hpp +++ b/dune/amdis/terms/TermGenerator.hpp @@ -9,7 +9,7 @@ #include <dune/amdis/terms/ConstantTerm.hpp> #include <dune/amdis/terms/CoordsTerm.hpp> -#include <dune/amdis/common/ConceptsBase.hpp> +#include <dune/amdis/common/Concepts.hpp> #include <dune/amdis/common/ValueCategory.hpp> namespace AMDiS diff --git a/dune/amdis/utility/GetEntity.hpp b/dune/amdis/utility/GetEntity.hpp index 2a756eca..9f345a2d 100644 --- a/dune/amdis/utility/GetEntity.hpp +++ b/dune/amdis/utility/GetEntity.hpp @@ -1,18 +1,9 @@ #pragma once +#include <dune/amdis/common/ConceptsBase.hpp> + namespace AMDiS { - namespace Impl - { - // Workaround for MSVC (problems with alias templates in pack expansion) - template <class... Args> - struct VoidType { using type = void; }; - } - - template <class... Args> - using Void_t = typename Impl::VoidType<Args...>::type; - - namespace Impl { template <class E, class = Void_t<>> -- GitLab