Skip to content
GitLab
Commit e54ba68f authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Concepts redesigned and FieldVector/FieldMatrix operations added

parent 0bcc9707
Hide whitespace changes
Inline Side-by-side
dune/amdis/Assembler.hpp
View file @ e54ba68f
......@@ -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:
......
dune/amdis/DirichletBC.hpp
View file @ e54ba68f
......@@ -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))
......
dune/amdis/FirstOrderAssembler.hpp
View file @ e54ba68f
......@@ -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)
......
dune/amdis/Operator.hpp
View file @ e54ba68f
......@@ -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);
......
dune/amdis/OperatorTermBase.hpp
View file @ e54ba68f
......@@ -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
dune/amdis/SecondOrderAssembler.hpp
View file @ e54ba68f
......@@ -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());
......
dune/amdis/ZeroOrderAssembler.hpp
View file @ e54ba68f
......@@ -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());
......
dune/amdis/common/Concepts.hpp 0 → 100644
View file @ e54ba68f
#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
dune/amdis/common/ConceptsBase.hpp
View file @ e54ba68f
/** \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