Commit c8db8ccb authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Reference GridFunction

parent 1ba853f5
......@@ -173,44 +173,6 @@ namespace AMDiS
};
namespace Concepts
{
namespace Definition
{
template <class F, int dow>
constexpr bool CallableDow =
Concepts::Callable<F, Dune::FieldVector<double, dow>>;
}
/// \brief Functor F is collable with GlobalCoordinates `F(Dune::FieldVector<double,DOW>)`
#ifndef WORLDDIM
template <class F>
constexpr bool CallableDomain =
Definition::CallableDow<F, 1> || Definition::CallableDow<F, 2> || Definition::CallableDow<F, 3>;
#else
template <class F>
constexpr bool CallableDomain =
Definition::CallableDow<F, WORLDDIM>;
#endif
template <class F>
using CallableDomain_t = bool_t<CallableDomain<F>>;
} // end namespace Concepts
// Creator for the AnalyticGridFunction
template <class Function>
struct GridFunctionCreator<Function, std::enable_if_t<Concepts::CallableDomain<Function>>>
{
template <class GridView>
static auto create(Function const& fct, GridView const& gridView)
{
return AnalyticGridFunction<Function, GridView>{fct, gridView};
}
};
/// \fn derivative
/// \brief Return a GridFunction representing the derivative of a functor.
/**
......@@ -242,6 +204,9 @@ namespace AMDiS
template <class GridView>
static auto create(Self const& self, GridView const& gridView)
{
using Coordinate = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
static_assert(std::is_invocable_v<Function, Coordinate>);
return AnalyticGridFunction<Function, GridView>{self.fct_, gridView};
}
};
......@@ -269,8 +234,7 @@ namespace AMDiS
* - `evalAtQP([](Dune::FieldVector<double, 2> const& x) { return x[0]; })`
* - `evalAtQP(Operation::TwoNorm{})`
**/
template <class Function,
REQUIRES(Concepts::CallableDomain<Function>)>
template <class Function>
auto evalAtQP(Function const& f)
{
return AnalyticPreGridFunction<Function>{f};
......
......@@ -14,46 +14,44 @@
namespace AMDiS
{
namespace Impl
{
template <class T0, class... Ts>
struct DomainType
{
using type = typename T0::Domain;
};
template <class T0, class... Ts>
struct EntitySetType
{
using type = typename T0::EntitySet;
};
} // end namespace Impl
#ifndef DOXYGEN
template <class Signatur, class Functor, class... LocalFunctions>
template <class Signatur, class Element, class Functor, class... LocalFunctions>
class FunctorLocalFunction;
#endif
// implementation
template <class R, class D, class Functor, class... LocalFunctions>
class FunctorLocalFunction<R(D), Functor, LocalFunctions...>
template <class R, class D, class E, class Functor, class... LocalFunctions>
class FunctorLocalFunction<R(D), E, Functor, LocalFunctions...>
{
public:
using Range = R;
using Domain = D;
using Element = E;
enum { hasDerivative = true };
template <class LocalFct>
struct LocalFunctionWrapper
{
template <class GridFct>
LocalFunctionWrapper(GridFct const& gridFct)
private:
// if we get a GridFunction extract its LocalFunction
template <class GridFct, decltype((localFunction(std::declval<GridFct>()),true)) = true>
LocalFunctionWrapper(GridFct const& gridFct, Dune::PriorityTag<2>)
: localFct_{localFunction(gridFct)}
{}
// if we get a LocalFunction store it directly
template <class LF>
LocalFunctionWrapper(LF const& localFct, Dune::PriorityTag<1>)
: localFct_(localFct)
{}
public:
template <class Fct>
LocalFunctionWrapper(Fct const& fct)
: LocalFunctionWrapper{fct, Dune::PriorityTag<5>{}}
{}
LocalFct const& operator*() const { return localFct_; }
LocalFct & operator*() { return localFct_; }
......@@ -69,7 +67,6 @@ namespace AMDiS
{}
/// Calls \ref bind for all localFunctions
template <class Element>
void bind(Element const& element)
{
Ranges::forEach(localFcts_, [&](auto& localFct) {
......@@ -88,7 +85,8 @@ namespace AMDiS
/// Applies the functor \ref fct_ to the evaluated localFunctions
Range operator()(Domain const& x) const
{
return Ranges::apply([&](auto&&... localFct) { return fct_((*localFct)(x)...); }, localFcts_);
return Ranges::apply([&](auto const&... localFct)
{ return fct_((*localFct)(x)...); }, localFcts_);
}
public:
......@@ -108,6 +106,15 @@ namespace AMDiS
std::tuple<LocalFunctionWrapper<LocalFunctions>...> localFcts_;
};
template <class Element, class Functor, class... LocalFunctions>
auto makeFunctorLocalFunction (Functor const& f, LocalFunctions const&... localFcts)
{
using D = typename Element::Geometry::LocalCoordinate;
using R = TYPEOF(f(localFcts(std::declval<D>())...));
return FunctorLocalFunction<R(D), Element, Functor, LocalFunctions...>{f, localFcts...};
}
/// \fn derivative
/// \brief Derivative of the LocalFunction of a FunctorGridFunction, utilizing
......@@ -118,9 +125,9 @@ namespace AMDiS
* **Requirements:**
* - The Functor `F` must model `Concepts::HasPartial`
**/
template <class Sig, class F, class... LFs, class Type,
template <class Sig, class E, class F, class... LFs, class Type,
REQUIRES(Concepts::HasPartial<F>)>
auto derivativeOf(FunctorLocalFunction<Sig,F,LFs...> const& lf, Type const& type)
auto derivativeOf(FunctorLocalFunction<Sig,E,F,LFs...> const& lf, Type const& type)
{
auto index_seq = std::index_sequence_for<LFs...>{};
......@@ -128,21 +135,21 @@ namespace AMDiS
auto term_i = [&](auto const _i)
{
auto di_f = Ranges::apply([&](auto const&... lgfs) {
return makeFunctorGridFunction(partial(lf.fct(), _i), (*lgfs)...);
return makeFunctorLocalFunction<E>(partial(lf.fct(), _i), (*lgfs)...);
}, lf.localFcts());
using std::get;
auto const& lgfs_i = *get<VALUE(_i)>(lf.localFcts());
return makeFunctorGridFunction(Operation::Multiplies{}, di_f, derivativeOf(lgfs_i, type));
return makeFunctorLocalFunction<E>(Operation::Multiplies{}, di_f, derivativeOf(lgfs_i, type));
};
// sum_i [ d_i(f)[lgfs...] * derivativeOf(lgfs_i)
auto gridFct = Ranges::apply([&](auto const... _i)
auto localFct = Ranges::apply([&](auto const... _i)
{
return makeFunctorGridFunction(Operation::Plus{}, term_i(_i)...);
return makeFunctorLocalFunction<E>(Operation::Plus{}, term_i(_i)...);
}, index_seq);
return localFunction(gridFct);
return localFct;
}
......@@ -154,10 +161,10 @@ namespace AMDiS
* - The functor `F` must model `Concepts::HasFunctorOrder`
* - All localFunctions `LFs...` must model `Concepts::Polynomial`
**/
template <class Sig, class F, class... LFs,
template <class Sig, class E, class F, class... LFs,
REQUIRES(Concepts::HasFunctorOrder<F,sizeof...(LFs)>
&& (Concepts::Polynomial<LFs> &&...))>
int order(FunctorLocalFunction<Sig,F,LFs...> const& lf)
int order(FunctorLocalFunction<Sig,E,F,LFs...> const& lf)
{
return Ranges::apply([&lf](auto const&... lgfs) {
return order(lf.fct(), order(*lgfs)...);
......@@ -173,6 +180,8 @@ namespace AMDiS
* locally it is evaluated
* \f$ f(g_0(x), g_1(x), ...) \f$
*
* \tparam Sig Signature of the GridFunction
* \tparam EntitySet The EntitySet this GridFunction is defined on
* \tparam Functor The type of the outer functor `f`
* \tparam GridFunctions... The GridFunction types of `g_i`
*
......@@ -181,18 +190,16 @@ namespace AMDiS
* - `arity(g_i) == arity(g_j) for i != j`
* - `g_i` models concept \ref GridFunction
**/
template <class Functor, class... GridFunctions>
class FunctorGridFunction
template <class Sig, class EntitySet, class Functor, class... GridFunctions>
class FunctorGridFunction;
template <class R, class D, class ES, class Functor, class... GridFunctions>
class FunctorGridFunction<R(D), ES, Functor, GridFunctions...>
{
public:
/// The result type of the functor when applied to the grid-functions
using Range = decltype(std::declval<Functor>()(std::declval<typename GridFunctions::Range>()...));
/// The argument type that can be applied to the grid-functions
using Domain = typename Impl::DomainType<GridFunctions...>::type;
/// The set of entities this grid-function binds to
using EntitySet = typename Impl::EntitySetType<GridFunctions...>::type;
using Range = R;
using Domain = D;
using EntitySet = ES;
enum { hasDerivative = false };
......@@ -200,51 +207,63 @@ namespace AMDiS
template <class GridFct>
using LocalFct = TYPEOF( localFunction(std::declval<GridFct>()) );
using RawRange = remove_cvref_t<Range>;
using LocalDomain = typename EntitySet::LocalCoordinate;
using Element = typename EntitySet::Element;
public:
using LocalFunction = FunctorLocalFunction<RawRange(LocalDomain), Functor, LocalFct<GridFunctions>...>;
using LocalFunction
= FunctorLocalFunction<Range(LocalDomain), Element, Functor, LocalFct<GridFunctions>...>;
/// \brief Constructor. Stores copies of the functor and gridfunctions.
template <class... GridFcts>
explicit FunctorGridFunction(Functor const& fct, GridFcts&&... gridFcts)
: fct_{fct}
FunctorGridFunction(EntitySet const& entitySet, Functor const& fct, GridFcts&&... gridFcts)
: entitySet_{entitySet}
, fct_{fct}
, gridFcts_{FWD(gridFcts)...}
{}
/// Applies the functor to the evaluated gridfunctions
Range operator()(Domain const& x) const
{
return Ranges::apply([&](auto&&... gridFct) { return fct_(gridFct(x)...); }, gridFcts_);
return Ranges::apply([&](auto const&... gridFct)
{ return fct_(gridFct(x)...); }, gridFcts_);
}
/// Return the stored \ref EntitySet of the first GridFunction
EntitySet const& entitySet() const
{
return std::get<0>(gridFcts_).entitySet();
return entitySet_;
}
LocalFunction makeLocalFunction() const
{
return Ranges::apply([&](auto const&... gridFcts) { return LocalFunction{fct_, gridFcts...}; }, gridFcts_);
return Ranges::apply([&](auto const&... gridFcts)
{ return LocalFunction{fct_, gridFcts...}; }, gridFcts_);
}
private:
EntitySet entitySet_;
Functor fct_;
std::tuple<GridFunctions...> gridFcts_;
};
// Generator function for FunctorGridFunction expressions
template <class Functor, class... GridFcts>
auto makeFunctorGridFunction(Functor const& f, GridFcts const&... gridFcts)
template <class Functor, class GridView, class... GridFcts>
auto makeFunctorGridFunction(Functor const& f, GridView const& gridView,
GridFcts const&... gridFcts)
{
static_assert((Concepts::GridFunction<GridFcts> &&...),
static_assert((Concepts::GridFunction<GridFcts> && ...),
"All passed parameters must be GridFunctions.");
static_assert(Concepts::Callable<Functor, typename GridFcts::Range...>,
"Range types of grid functions are not compatible with the functor.");
return FunctorGridFunction<Functor, GridFcts...>{f, gridFcts...};
using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
using Domain = typename EntitySet::GlobalCoordinate;
using Range = TYPEOF(f(gridFcts(std::declval<Domain>())...));
using FGF = FunctorGridFunction<Range(Domain), EntitySet, Functor, GridFcts...>;
return FGF{EntitySet{gridView},f, gridFcts...};
}
......@@ -261,7 +280,7 @@ namespace AMDiS
static auto create(Self const& self, GridView const& gridView)
{
return Ranges::apply([&](auto const&... pgf) {
return makeFunctorGridFunction(self.fct_,
return makeFunctorGridFunction(self.fct_, gridView,
makeGridFunction(pgf, gridView)...);
}, self.preGridFcts_);
}
......
......@@ -96,19 +96,33 @@ namespace AMDiS
struct GridFunctionCreator
: PreGridFct::Creator {};
// forward declaration
template <class Function>
struct AnalyticPreGridFunction;
namespace Impl
{
// Specialization for type that is already a GridFunction
template <class GridFct, class GridView>
GridFct const& makeGridFunctionImpl(GridFct const& gridFct, GridView const& /*gridView*/, std::true_type)
GridFct const& makeGridFunctionImpl(GridFct const& gridFct, GridView const& /*gridView*/, std::true_type, Dune::PriorityTag<2>)
{
return gridFct;
}
// If F is a callable with global coordinates, create an AnalyticGridFunction
template <class F, class GridView,
class Coordinate = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate,
std::enable_if_t<std::is_invocable_v<F, Coordinate>, int> = 0>
auto makeGridFunctionImpl(F const& f, GridView const& gridView, std::false_type, Dune::PriorityTag<1>)
{
AnalyticPreGridFunction<F> preGridFct{f};
return AnalyticPreGridFunction<F>::Creator::create(preGridFct, gridView);
}
// Use the \ref GridFunctionCreator to create a gridFunction from a preGridFunction
template <class PreGridFct, class GridView,
class Creator = GridFunctionCreator<PreGridFct>>
decltype(auto) makeGridFunctionImpl(PreGridFct const& preGridFct, GridView const& gridView, std::false_type)
class Creator = GridFunctionCreator<PreGridFct>>
decltype(auto) makeGridFunctionImpl(PreGridFct const& preGridFct, GridView const& gridView, std::false_type, Dune::PriorityTag<0>)
{
return Creator::create(preGridFct, gridView);
}
......@@ -154,7 +168,7 @@ namespace AMDiS
decltype(auto) makeGridFunction(PreGridFct const& preGridFct, GridView const& gridView)
{
using isGridFct = Concepts::GridFunction_t<PreGridFct>;
return Impl::makeGridFunctionImpl(preGridFct, gridView, isGridFct{});
return Impl::makeGridFunctionImpl(preGridFct, gridView, isGridFct{}, Dune::PriorityTag<5>{});
}
} // end namespace AMDiS
......@@ -160,6 +160,20 @@ namespace AMDiS
return invokeAtQP(Operation::Get_{i}, FWD(value));
}
/// \brief Return a reference-wrapper as GridFunction. \relates FunctorGridFunction
template <class T>
auto constant(T const& value)
{
return invokeAtQP(Operation::Constant{value});
}
/// \brief Return a reference-wrapper as GridFunction. \relates FunctorGridFunction
template <class T>
auto reference(T& value)
{
return invokeAtQP(Operation::Reference{value});
}
// @}
// unary vector operations
......
......@@ -114,6 +114,38 @@ namespace AMDiS
// -------------------------------------------------------------------------
/// Functor representing a constant value
template <class T>
struct Reference
{
constexpr Reference(T& value)
: value_(&value)
{}
template <class... Ts>
constexpr T const& operator()(Ts const&... /*args*/) const
{
return *value_;
}
private:
T* value_;
};
template <class T, class... Int>
constexpr int order(Reference<T> const&, Int... /*orders*/)
{
return 0;
}
template <class T, std::size_t J>
constexpr auto partial(Reference<T> const&, index_t<J>)
{
return Zero{};
}
// -------------------------------------------------------------------------
template <class T0, class... Ts>
inline constexpr decltype(auto) get_element(index_t<0>, T0&& t0, Ts&&... /*ts*/)
......
......@@ -42,6 +42,9 @@ int main(int argc, char** argv)
// Evaluation of the DOFVector (component) at quadrature points
auto expr9 = prob.solution();
auto expr10 = reference(expr4);
auto expr11 = expr6 * expr10;
// ---------------------------------
// derivative of expressions
......@@ -53,6 +56,10 @@ int main(int argc, char** argv)
auto diff8 = gradientOf(expr8);
auto diff9 = gradientOf(expr9);
auto diff10 = gradientOf(expr9 + expr9);
auto diff11 = gradientOf(2 * expr9);
auto diff12 = gradientOf(expr8 * expr9);
// ---------------------------------
u.interpolate(expr1);
......
Markdown is supported
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