From 78ab641516ab0a8bac3b15bf284100af3cc07391 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Sat, 11 May 2019 18:59:28 +0200
Subject: [PATCH] Add generalized derivative and local-to-global adapter to
handle the global derivatives uniformly
---
src/amdis/Marker.hpp | 2 +-
src/amdis/common/DerivativeTraits.hpp | 58 ++++
.../gridfunctions/AnalyticGridFunction.hpp | 64 ++--
src/amdis/gridfunctions/CMakeLists.txt | 5 +-
.../gridfunctions/ConstantGridFunction.hpp | 26 +-
.../gridfunctions/CoordsGridFunction.hpp | 5 +-
src/amdis/gridfunctions/DOFVectorView.hpp | 2 +-
src/amdis/gridfunctions/Derivative.hpp | 65 ++++
.../gridfunctions/DerivativeGridFunction.hpp | 93 ++++--
src/amdis/gridfunctions/DiscreteFunction.hpp | 10 +-
.../gridfunctions/DiscreteFunction.inc.hpp | 316 +++++++++++++-----
.../gridfunctions/FunctorGridFunction.hpp | 22 +-
...dFunctionConcepts.hpp => GridFunction.hpp} | 96 +-----
src/amdis/gridfunctions/Order.hpp | 51 +++
.../localoperators/FirstOrderPartialTest.hpp | 30 +-
.../FirstOrderTestPartialTrial.hpp | 35 +-
.../SecondOrderPartialTestPartialTrial.hpp | 54 +--
src/amdis/localoperators/ZeroOrderTest.hpp | 4 +-
src/amdis/utility/LocalBasisCache.hpp | 21 ++
src/amdis/utility/LocalToGlobalAdapter.hpp | 207 ++++++++++++
src/amdis/utility/QuadratureFactory.hpp | 27 +-
test/CMakeLists.txt | 3 +
test/DiscreteFunctionTest.cpp | 4 +-
test/GradientTest.cpp | 91 +++++
test/Tests.hpp | 2 +-
25 files changed, 908 insertions(+), 385 deletions(-)
create mode 100644 src/amdis/common/DerivativeTraits.hpp
create mode 100644 src/amdis/gridfunctions/Derivative.hpp
rename src/amdis/gridfunctions/{GridFunctionConcepts.hpp => GridFunction.hpp} (62%)
create mode 100644 src/amdis/gridfunctions/Order.hpp
create mode 100644 src/amdis/utility/LocalToGlobalAdapter.hpp
create mode 100644 test/GradientTest.cpp
diff --git a/src/amdis/Marker.hpp b/src/amdis/Marker.hpp
index 48f5ccd2..3f99fe7a 100644
--- a/src/amdis/Marker.hpp
+++ b/src/amdis/Marker.hpp
@@ -10,7 +10,7 @@
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/TypeTraits.hpp>
-#include <amdis/gridfunctions/GridFunctionConcepts.hpp>
+#include <amdis/gridfunctions/GridFunction.hpp>
#include <amdis/AdaptInfo.hpp>
#include <amdis/Flag.hpp>
diff --git a/src/amdis/common/DerivativeTraits.hpp b/src/amdis/common/DerivativeTraits.hpp
new file mode 100644
index 00000000..e86e0ff6
--- /dev/null
+++ b/src/amdis/common/DerivativeTraits.hpp
@@ -0,0 +1,58 @@
+#pragma once
+
+#include <dune/functions/common/defaultderivativetraits.hh>
+
+namespace Dune
+{
+ template <class K, int n>
+ class FieldVector;
+
+ template <class K, int n, int m>
+ class FieldMatrix;
+}
+
+namespace AMDiS
+{
+ namespace tag
+ {
+ struct gradient {};
+ struct divergence {};
+ struct partial { std::size_t comp = 0; };
+
+ // register possible types for derivative traits
+ struct derivative_type : gradient, divergence, partial {};
+ }
+
+ template <class Sig, class Type>
+ struct DerivativeTraits;
+
+ template <class R, class D>
+ struct DerivativeTraits<R(D), tag::gradient>
+ : public Dune::Functions::DefaultDerivativeTraits<R(D)>
+ {};
+
+ template <class R, class D>
+ struct DerivativeTraits<R(D), tag::partial>
+ {
+ using Range = R;
+ };
+
+ template <class K, class D>
+ struct DerivativeTraits<K(D), tag::divergence>
+ {
+ // error
+ };
+
+ template <class K, int n, class D>
+ struct DerivativeTraits<Dune::FieldVector<K,n>(D), tag::divergence>
+ {
+ using Range = K;
+ };
+
+ template <class K, int n, int m, class D>
+ struct DerivativeTraits<Dune::FieldMatrix<K,n,m>(D), tag::divergence>
+ {
+ using Range = Dune::FieldVector<K,m>;
+ };
+
+} // end namespace AMDiS
\ No newline at end of file
diff --git a/src/amdis/gridfunctions/AnalyticGridFunction.hpp b/src/amdis/gridfunctions/AnalyticGridFunction.hpp
index 1c48728c..496c6b19 100644
--- a/src/amdis/gridfunctions/AnalyticGridFunction.hpp
+++ b/src/amdis/gridfunctions/AnalyticGridFunction.hpp
@@ -7,7 +7,9 @@
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <amdis/Operations.hpp>
-#include <amdis/gridfunctions/GridFunctionConcepts.hpp>
+#include <amdis/common/DerivativeTraits.hpp>
+#include <amdis/gridfunctions/GridFunction.hpp>
+#include <amdis/gridfunctions/Order.hpp>
namespace AMDiS
{
@@ -97,20 +99,19 @@ namespace AMDiS
* **Requirements:**
* - The functor `F` must fulfill the concept \ref Concepts::HasPartial
**/
- template <class R, class D, class LC, class F>
- auto derivative(AnalyticLocalFunction<R(D),LC,F> const& lf)
+ template <class R, class D, class LC, class F, class Type>
+ auto derivative(AnalyticLocalFunction<R(D),LC,F> const& lf, Type const& type)
{
- static_assert(Concepts::HasPartial<F>,
- "No partial(...,_0) defined for Functor F of AnalyticLocalFunction.");
+ static_assert(Concepts::HasDerivative<F,Type>,
+ "No derivative(F,DerivativeType) defined for Functor F of AnalyticLocalFunction.");
- auto df = partial(lf.fct(), index_<0>);
+ auto df = derivative(lf.fct(), type);
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
- using DerivativeSignature = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range(D);
+ using DerivativeSignature = typename DerivativeTraits<RawSignature,Type>::Range(D);
return AnalyticLocalFunction<DerivativeSignature,LC,decltype(df)>{df};
}
-
/// \class AnalyticGridFunction
/// \brief A Gridfunction that evaluates a function with global coordinates.
/**
@@ -152,7 +153,7 @@ namespace AMDiS
}
/// Return the \ref AnalyticLocalFunction of the AnalyticGridFunction.
- LocalFunction localFunction() const
+ LocalFunction makeLocalFunction() const
{
return {fct_};
}
@@ -171,6 +172,30 @@ namespace AMDiS
EntitySet entitySet_;
};
+
+ 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
+
+ } // end namespace Concepts
+
+
// Creator for the AnalyticGridFunction
template <class Function>
struct GridFunctionCreator<Function, std::enable_if_t<Concepts::CallableDomain<Function>>>
@@ -183,30 +208,21 @@ namespace AMDiS
};
- /// \fn localFunction
- /// \brief Return the LocalFunction of the AnalyticGridFunction.
- /// \relates AnalyticGridfunction
- template <class F, class GV>
- auto localFunction(AnalyticGridFunction<F,GV> const& gf)
- {
- return gf.localFunction();
- }
-
/// \fn derivative
/// \brief Return a GridFunction representing the derivative of a functor.
/**
* \relates AnalyticGridfunction
*
* **Requirements:**
- * - Functor `F` must fulfill the concept \ref Concepts::HasPartial
+ * - Functor `F` must fulfill the concept \ref Concepts::HasDerivative<Type>
**/
- template <class F, class GV>
- auto derivative(AnalyticGridFunction<F,GV> const& gf)
+ template <class F, class GV, class Type>
+ auto derivative(AnalyticGridFunction<F,GV> const& gf, Type const& type)
{
- static_assert(Concepts::HasPartial<F>,
- "No partial(...,_0) defined for Functor of AnalyticLocalFunction.");
+ static_assert(Concepts::HasDerivative<F,Type>,
+ "No derivative(F,DerivativeType) defined for Functor of AnalyticLocalFunction.");
- auto df = partial(gf.fct(), index_<0>);
+ auto df = derivative(gf.fct(), type);
return AnalyticGridFunction<decltype(df), GV>{df, gf.entitySet().gridView()};
}
diff --git a/src/amdis/gridfunctions/CMakeLists.txt b/src/amdis/gridfunctions/CMakeLists.txt
index e1a017f5..e9672a2e 100644
--- a/src/amdis/gridfunctions/CMakeLists.txt
+++ b/src/amdis/gridfunctions/CMakeLists.txt
@@ -4,11 +4,14 @@ install(FILES
AnalyticGridFunction.hpp
ConstantGridFunction.hpp
CoordsGridFunction.hpp
+ Derivative.hpp
DerivativeGridFunction.hpp
+ DerivativeTraits.hpp
DiscreteFunction.hpp
DiscreteFunction.inc.hpp
DOFVectorView.hpp
FunctorGridFunction.hpp
- GridFunctionConcepts.hpp
+ GridFunction.hpp
OperationsGridFunction.hpp
+ Order.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/gridfunctions)
diff --git a/src/amdis/gridfunctions/ConstantGridFunction.hpp b/src/amdis/gridfunctions/ConstantGridFunction.hpp
index 1cfe69d7..96134909 100644
--- a/src/amdis/gridfunctions/ConstantGridFunction.hpp
+++ b/src/amdis/gridfunctions/ConstantGridFunction.hpp
@@ -6,11 +6,11 @@
#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
-#include <dune/functions/common/defaultderivativetraits.hh>
+#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/gridfunctions/AnalyticGridFunction.hpp>
-#include <amdis/gridfunctions/GridFunctionConcepts.hpp>
+#include <amdis/gridfunctions/GridFunction.hpp>
namespace AMDiS
{
@@ -53,10 +53,11 @@ namespace AMDiS
/// \brief Create a \ref ConstantLocalFunction representing the derivative
/// of a constant function, that ist, the value 0.
- auto derivative() const
+ template <class Type>
+ auto makeDerivative(Type const& type) const
{
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
- using DerivativeRange = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range;
+ using DerivativeRange = typename DerivativeTraits<RawSignature,Type>::Range;
DerivativeRange diff(0);
return ConstantLocalFunction<DerivativeRange(D),LC,DerivativeRange>{diff};
}
@@ -71,13 +72,6 @@ namespace AMDiS
T value_;
};
- /// \relates ConstantLocalFunction
- template <class R, class D, class LC, class T>
- auto derivative(ConstantLocalFunction<R(D), LC, T> const& lf)
- {
- return lf.derivative();
- }
-
/// \brief Gridfunction returning a constant value.
/**
@@ -122,7 +116,7 @@ namespace AMDiS
}
/// \brief Create an \ref ConstantLocalFunction with the stores `value_`.
- LocalFunction localFunction() const
+ LocalFunction makeLocalFunction() const
{
return {value_};
}
@@ -133,14 +127,6 @@ namespace AMDiS
};
- /// \relates ConstantGridFunction
- template <class T, class GV>
- auto localFunction(ConstantGridFunction<T,GV> const& gf)
- {
- return gf.localFunction();
- }
-
-
namespace Concepts
{
/** \addtogroup Concepts
diff --git a/src/amdis/gridfunctions/CoordsGridFunction.hpp b/src/amdis/gridfunctions/CoordsGridFunction.hpp
index 28965dff..cfdfd72c 100644
--- a/src/amdis/gridfunctions/CoordsGridFunction.hpp
+++ b/src/amdis/gridfunctions/CoordsGridFunction.hpp
@@ -5,6 +5,7 @@
#include <dune/common/diagonalmatrix.hh>
#include <dune/common/typeutilities.hh>
+#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/gridfunctions/AnalyticGridFunction.hpp>
namespace AMDiS
@@ -48,7 +49,7 @@ namespace AMDiS
return Dune::DiagonalMatrix<T,N>{T(1)};
}
};
- friend Derivative partial(CoordsFunction const& /*f*/, index_t<0>)
+ friend Derivative derivative(CoordsFunction const& /*f*/, tag::gradient)
{
return Derivative{};
}
@@ -105,7 +106,7 @@ namespace AMDiS
int comp_;
};
- friend Derivative partial(Self const& f, index_t<0>)
+ friend Derivative derivative(Self const& f, tag::gradient)
{
return Derivative{f.comp_};
}
diff --git a/src/amdis/gridfunctions/DOFVectorView.hpp b/src/amdis/gridfunctions/DOFVectorView.hpp
index 8a443e2f..aee45435 100644
--- a/src/amdis/gridfunctions/DOFVectorView.hpp
+++ b/src/amdis/gridfunctions/DOFVectorView.hpp
@@ -1,8 +1,8 @@
#pragma once
-#include <amdis/GridFunctions.hpp>
#include <amdis/functions/Interpolate.hpp>
#include <amdis/gridfunctions/DiscreteFunction.hpp>
+#include <amdis/gridfunctions/GridFunction.hpp>
namespace AMDiS
{
diff --git a/src/amdis/gridfunctions/Derivative.hpp b/src/amdis/gridfunctions/Derivative.hpp
new file mode 100644
index 00000000..6516cd6b
--- /dev/null
+++ b/src/amdis/gridfunctions/Derivative.hpp
@@ -0,0 +1,65 @@
+#pragma once
+
+#include <type_traits>
+
+#include <amdis/common/Concepts.hpp>
+#include <amdis/common/DerivativeTraits.hpp>
+#include <amdis/common/Index.hpp>
+
+namespace AMDiS
+{
+ /// The derivative of a localfunction as localfunction itself
+ template <class LocalFunction, class Type,
+ REQUIRES(std::is_convertible<tag::derivative_type, Type>::value),
+ class = void_t<decltype(std::declval<LocalFunction>().makeDerivative(std::declval<Type>()))> >
+ auto derivative(LocalFunction const& lf, Type const& type)
+ {
+ return lf.makeDerivative(type);
+ }
+
+ namespace Concepts
+ {
+ /** \addtogroup Concepts
+ * @{
+ **/
+
+ namespace Definition
+ {
+ struct HasDerivative
+ {
+ template <class F, class T>
+ auto requires_(F&& f, T&& t) -> decltype( derivative(f,t) );
+ };
+
+ struct HasLocalFunctionDerivative
+ {
+ template <class F, class T>
+ auto requires_(F&& f, T&& t) -> decltype( derivative(localFunction(f),t) );
+ };
+
+ struct HasPartial
+ {
+ template <class F, class I>
+ auto requires_(F&& f, I&& i) -> decltype( partial(f, i) );
+ };
+
+ } // end namespace Definition
+
+
+ /// \brief GridFunction GF has free function `derivative(F)`
+ template <class GF, class Type>
+ constexpr bool HasDerivative = models<Definition::HasDerivative(GF,Type)>;
+
+ /// \brief GridFunction GF has free function `derivative(localFunction(F))`
+ template <class GF, class Type>
+ constexpr bool HasLocalFunctionDerivative = models<Definition::HasLocalFunctionDerivative(GF,Type)>;
+
+ /// \brief Functor F has free function `partial(F,_0)`
+ template <class F>
+ constexpr bool HasPartial = models<Definition::HasPartial(F,index_t<0>)>;
+
+ /** @} **/
+
+ } // end namespace Concepts
+
+} // end namespace AMDiS
diff --git a/src/amdis/gridfunctions/DerivativeGridFunction.hpp b/src/amdis/gridfunctions/DerivativeGridFunction.hpp
index 19693c8d..eb96d470 100644
--- a/src/amdis/gridfunctions/DerivativeGridFunction.hpp
+++ b/src/amdis/gridfunctions/DerivativeGridFunction.hpp
@@ -2,40 +2,65 @@
#include <type_traits>
-#include <dune/functions/common/defaultderivativetraits.hh>
#include <dune/grid/utility/hierarchicsearch.hh>
-#include <amdis/gridfunctions/GridFunctionConcepts.hpp>
+#include <amdis/common/Concepts.hpp>
+#include <amdis/common/DerivativeTraits.hpp>
+#include <amdis/gridfunctions/Derivative.hpp>
+#include <amdis/gridfunctions/GridFunction.hpp>
namespace AMDiS
{
+ namespace Impl
+ {
+ template <class LF, class Sig>
+ struct CheckFunctorConcept
+ {
+ static_assert(Concepts::Functor<LF, Sig>, "Derivative of LocalFunction can not be called as a functor.");
+ };
+
+ template <class Traits>
+ struct CheckValidRange
+ {
+ static_assert(!std::is_same<typename Traits::Range, Dune::Functions::InvalidRange>::value, "Invalid Range.");
+ };
+ }
+
/// \class DerivativeGridFunction
/// \brief A Gridfunction that returns the derivative when calling localFunction.
/**
* \ingroup GridFunctions
- * Wrapps the GridFunction so that \ref localFunction returns a LocalFunction
+ * Wraps the GridFunction so that \ref localFunction returns a LocalFunction
* representing the derivative of the LocalFunction of the GridFunction.
*
* \tparam GridFunction The GridFunction that is wrapped.
+ * \tparam The type of derivative, i.e. one of tag::gradient, tag::partial{i}, or tag::divergence
*
* **Requirements:**
* - `GridFunction` models \ref Concepts::GridFunction and the LocalFunction
* has a derivative.
**/
- template <class GridFunction>
+ template <class GridFunction, class Type>
class DerivativeGridFunction
{
using GridFctRange = typename GridFunction::Range;
using GridFctDomain = typename GridFunction::Domain;
using RawSignature = typename Dune::Functions::SignatureTraits<GridFctRange(GridFctDomain)>::RawSignature;
- using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
- using LocalFunction = TYPEOF( derivative(localFunction(std::declval<GridFunction>())) ) ;
+
+ using Traits = DerivativeTraits<RawSignature, Type>;
+ using LocalFunction = TYPEOF( derivative(localFunction(std::declval<GridFunction>()), std::declval<Type>()) ) ;
+
+ using LocalFctRange = typename Traits::Range;
+ using LocalFctDomain = typename GridFunction::EntitySet::LocalCoordinate;
+
+ using _CHECK1_ = Impl::CheckValidRange<Traits>;
+ using _CHECK2_ = Impl::CheckFunctorConcept<LocalFunction, LocalFctRange(LocalFctDomain)>;
enum { hasDerivative = false };
public:
/// The Range of the derivative of the GridFunction
- using Range = typename DerivativeTraits::Range;
+ using Range = typename Traits::Range;
/// The domain of the GridFunction
using Domain = GridFctDomain;
@@ -45,12 +70,10 @@ namespace AMDiS
public:
/// Constructor. Stores a copy of gridFct.
- explicit DerivativeGridFunction(GridFunction const& gridFct)
+ explicit DerivativeGridFunction(GridFunction const& gridFct, Type const& type)
: gridFct_{gridFct}
- {
- static_assert(isValidRange<DerivativeTraits>(),
- "Derivative of GridFunction not defined");
- }
+ , type_{type}
+ {}
/// Evaluate derivative in global coordinates. NOTE: expensive
Range operator()(Domain const& x) const
@@ -65,15 +88,15 @@ namespace AMDiS
auto element = hsearch.findEntity(x);
auto geometry = element.geometry();
- auto localFct = derivative(localFunction(gridFct_));
+ auto localFct = makeLocalFunction();
localFct.bind(element);
return localFct(geometry.local(x));
}
/// Return the derivative-localFunction of the GridFunction.
- friend LocalFunction localFunction(DerivativeGridFunction const& gf)
+ LocalFunction makeLocalFunction() const
{
- return derivative(localFunction(gf.gridFct_));
+ return derivative(localFunction(gridFct_), type_);
}
/// Return the \ref EntitySet of the \ref GridFunction.
@@ -84,6 +107,7 @@ namespace AMDiS
private:
GridFunction gridFct_;
+ Type type_;
};
@@ -100,19 +124,19 @@ namespace AMDiS
* `GridFct::hasDerivative == false`.
* - The localFunction of the `GridFct` models `Concepts::HasDerivative`.
**/
- template <class GridFct,
- class LocalFct = decltype(localFunction(std::declval<GridFct>())),
+ template <class GridFct, class Type,
+ class LocalFct = decltype( localFunction(std::declval<GridFct>()) ),
REQUIRES(not GridFct::hasDerivative)>
- auto derivative(GridFct const& gridFct)
+ auto derivative(GridFct const& gridFct, Type const& type)
{
- static_assert(Concepts::HasDerivative<LocalFct>,
- "derivative(LocalFunction) not defined!");
- return DerivativeGridFunction<GridFct>{gridFct};
+ static_assert(Concepts::HasDerivative<LocalFct,Type>,
+ "derivative(LocalFunction,type) not defined!");
+ return DerivativeGridFunction<GridFct,Type>{gridFct, type};
}
#ifndef DOXYGEN
- template <class Expr>
+ template <class Expr, class Type>
struct DerivativePreGridFunction
{
using Self = DerivativePreGridFunction;
@@ -122,17 +146,18 @@ namespace AMDiS
template <class GridView>
static auto create(Self const& self, GridView const& gridView)
{
- return derivative(makeGridFunction(self.expr_, gridView));
+ return derivative(makeGridFunction(self.expr_, gridView), self.type_);
}
};
Expr expr_;
+ Type type_ = {};
};
namespace Traits
{
- template <class Expr>
- struct IsPreGridFunction<DerivativePreGridFunction<Expr>>
+ template <class Expr,class Type>
+ struct IsPreGridFunction<DerivativePreGridFunction<Expr,Type>>
: std::true_type {};
}
#endif
@@ -143,7 +168,7 @@ namespace AMDiS
/// \relates DerivativeGridFunction
/**
* \ingroup GridFunctions
- * Generates a Gridfunction representing the derivative of a GridFunction.
+ * Generates a Gridfunction representing the gradient of a GridFunction.
* See \ref DerivativeGridFunction.
*
* **Examples:**
@@ -153,7 +178,21 @@ namespace AMDiS
template <class Expr>
auto gradientAtQP(Expr const& expr)
{
- return DerivativePreGridFunction<Expr>{expr};
+ return DerivativePreGridFunction<Expr, tag::gradient>{expr};
+ }
+
+ /// Generates a Gridfunction representing the divergence of a vector-valued GridFunction.
+ template <class Expr>
+ auto divergenceAtQP(Expr const& expr)
+ {
+ return DerivativePreGridFunction<Expr, tag::divergence>{expr};
+ }
+
+ /// Generates a Gridfunction representing the partial derivative of a GridFunction.
+ template <class Expr>
+ auto partialAtQP(Expr const& expr, std::size_t i)
+ {
+ return DerivativePreGridFunction<Expr, tag::partial>{expr, tag::partial{i}};
}
/** @} **/
diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp
index 9323112d..8a9ba447 100644
--- a/src/amdis/gridfunctions/DiscreteFunction.hpp
+++ b/src/amdis/gridfunctions/DiscreteFunction.hpp
@@ -9,7 +9,6 @@
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/typetree/childextraction.hh>
-#include <amdis/GridFunctions.hpp>
#include <amdis/LinearAlgebra.hpp>
#include <amdis/typetree/FiniteElementType.hpp>
#include <amdis/typetree/RangeType.hpp>
@@ -60,7 +59,12 @@ namespace AMDiS
public:
/// A LocalFunction representing the derivative of the DOFVector on a bound element
+ template <class Type>
+ class DerivativeLocalFunctionBase;
+
class GradientLocalFunction;
+ class PartialLocalFunction;
+ class DivergenceLocalFunction;
/// A LocalFunction representing the value the DOFVector on a bound element
class LocalFunction;
@@ -83,9 +87,9 @@ namespace AMDiS
Range operator()(Domain const& x) const;
/// \brief Create a local function for this view on the DOFVector. \relates LocalFunction
- friend LocalFunction localFunction(DiscreteFunction const& self)
+ LocalFunction makeLocalFunction() const
{
- return LocalFunction{self};
+ return LocalFunction{*this};
}
/// \brief Return a \ref Dune::Functions::GridViewEntitySet
diff --git a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
index 28478de7..4f4af74e 100644
--- a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
+++ b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
@@ -1,7 +1,9 @@
#pragma once
+#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/utility/LocalBasisCache.hpp>
+#include <amdis/utility/LocalToGlobalAdapter.hpp>
#include <dune/common/ftraits.hh>
#include <dune/grid/utility/hierarchicsearch.hh>
@@ -55,9 +57,19 @@ public:
Range operator()(Domain const& x) const;
/// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
- friend GradientLocalFunction derivative(LocalFunction const& localFunction)
+ GradientLocalFunction makeDerivative(tag::gradient type) const
{
- return GradientLocalFunction{localFunction.globalFunction_};
+ return GradientLocalFunction{globalFunction_, type};
+ }
+
+ DivergenceLocalFunction makeDerivative(tag::divergence type) const
+ {
+ return DivergenceLocalFunction{globalFunction_, type};
+ }
+
+ PartialLocalFunction makeDerivative(tag::partial type) const
+ {
+ return PartialLocalFunction{globalFunction_, type};
}
/// \brief The \ref polynomialDegree() of the LocalFunctions
@@ -83,16 +95,59 @@ private:
template <class GB, class VT, class TP>
-class DiscreteFunction<GB,VT,TP>::GradientLocalFunction
+typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
+LocalFunction::operator()(Domain const& x) const
+{
+ assert( bound_ );
+ Range y(0);
+
+ auto&& coefficients = *globalFunction_.dofVector_;
+ auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
+ for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
+ {
+ auto localBasisCache = makeNodeCache(node);
+ auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(), x);
+ std::size_t size = node.finiteElement().size();
+
+ // Get range entry associated to this node
+ auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
+
+ for (std::size_t i = 0; i < size; ++i) {
+ auto&& multiIndex = localView_.index(node.localIndex(i));
+
+ // Get coefficient associated to i-th shape function
+ auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
+
+ // Get value of i-th shape function
+ auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
+
+ std::size_t dimC = c.size();
+ std::size_t dimV = v.size();
+ assert(dimC*dimV == std::size_t(re.size()));
+ for(std::size_t j = 0; j < dimC; ++j) {
+ auto&& c_j = c[j];
+ for(std::size_t k = 0; k < dimV; ++k)
+ re[j*dimV + k] += c_j*v[k];
+ }
+ }
+ });
+
+ return y;
+}
+
+
+template <class GB, class VT, class TP>
+ template <class Type>
+class DiscreteFunction<GB,VT,TP>::DerivativeLocalFunctionBase
{
using R = typename DiscreteFunction::Range;
using D = typename DiscreteFunction::Domain;
using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
- using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
+ using Traits = DerivativeTraits<RawSignature, Type>;
public:
using Domain = typename EntitySet::LocalCoordinate;
- using Range = typename DerivativeTraits::Range;
+ using Range = typename Traits::Range;
enum { hasDerivative = false };
@@ -103,15 +158,17 @@ private:
public:
/// Constructor. Stores a copy of the DiscreteFunction.
- GradientLocalFunction(DiscreteFunction const& globalFunction)
+ DerivativeLocalFunctionBase(DiscreteFunction const& globalFunction, Type const& type)
: globalFunction_(globalFunction)
+ , type_(type)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
/// Copy constructor.
- GradientLocalFunction(GradientLocalFunction const& other)
+ DerivativeLocalFunctionBase(DerivativeLocalFunctionBase const& other)
: globalFunction_(other.globalFunction_)
+ , type_(other.type_)
, localView_(globalFunction_.basis()->localView())
, subTree_(&child(localView_.tree(), globalFunction_.treePath()))
{}
@@ -130,9 +187,6 @@ public:
bound_ = false;
}
- /// Evaluate Gradient at bound element in local coordinates
- Range operator()(Domain const& x) const;
-
int order() const
{
assert( bound_ );
@@ -146,8 +200,9 @@ public:
return localView_.element();
}
-private:
+protected:
DiscreteFunction globalFunction_;
+ Type type_;
LocalView localView_;
SubTree const* subTree_;
Dune::Std::optional<Geometry> geometry_;
@@ -156,105 +211,198 @@ private:
template <class GB, class VT, class TP>
-typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
-LocalFunction::operator()(Domain const& x) const
+class DiscreteFunction<GB,VT,TP>::GradientLocalFunction
+ : public DerivativeLocalFunctionBase<tag::gradient>
{
- assert( bound_ );
- Range y(0);
+ using Super = DerivativeLocalFunctionBase<tag::gradient>;
+public:
+ using Range = typename Super::Range;
+ using Domain = typename Super::Domain;
- auto&& coefficients = *globalFunction_.dofVector_;
- auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
- for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
+ using Super::DerivativeLocalFunctionBase;
+
+ /// Evaluate Gradient at bound element in local coordinates
+ Range operator()(Domain const& x) const
{
- auto&& fe = node.finiteElement();
- auto&& localBasis = fe.localBasis();
+ assert( Super::bound_ );
+ Range dy(0);
+
+ auto&& coefficients = *Super::globalFunction_.dofVector_;
+ auto&& nodeToRangeEntry = Super::globalFunction_.nodeToRangeEntry_;
+ for_each_leaf_node(*Super::subTree_, [&,this](auto const& node, auto const& tp)
+ {
+ auto localBasis = makeLocalToGlobalBasisAdapter(node, Super::geometry_.value());
+ auto const& gradients = localBasis.gradientsAt(x);
+
+ // Get range entry associated to this node
+ auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
+
+ for (std::size_t i = 0; i < localBasis.size(); ++i) {
+ auto&& multiIndex = Super::localView_.index(node.localIndex(i));
+
+ // Get coefficient associated to i-th shape function
+ auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
+
+ // Get value of i-th transformed reference gradient
+ auto grad = Dune::Functions::flatVectorView(gradients[i]);
+
+ std::size_t dimC = c.size();
+ std::size_t dimV = grad.size();
+ assert(dimC*dimV == std::size_t(re.size()));
+ for(std::size_t j = 0; j < dimC; ++j) {
+ auto&& c_j = c[j];
+ for(std::size_t k = 0; k < dimV; ++k)
+ re[j*dimV + k] += c_j*grad[k];
+ }
+ }
+ });
- NodeCache<TYPEOF(node)> localBasisCache(localBasis);
- auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(),x);
+ return dy;
+ }
+};
- // Get range entry associated to this node
- auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
- for (std::size_t i = 0; i < localBasis.size(); ++i) {
- auto&& multiIndex = localView_.index(node.localIndex(i));
+template <class GB, class VT, class TP>
+class DiscreteFunction<GB,VT,TP>::DivergenceLocalFunction
+ : public DerivativeLocalFunctionBase<tag::divergence>
+{
+ using Super = DerivativeLocalFunctionBase<tag::divergence>;
- // Get coefficient associated to i-th shape function
- auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
+public:
+ using Range = typename Super::Range;
+ using Domain = typename Super::Domain;
- // Get value of i-th shape function
- auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
+ using Super::DerivativeLocalFunctionBase;
- std::size_t dimC = c.size();
- std::size_t dimV = v.size();
- assert(dimC*dimV == std::size_t(re.size()));
- for(std::size_t j = 0; j < dimC; ++j) {
- auto&& c_j = c[j];
- for(std::size_t k = 0; k < dimV; ++k)
- re[j*dimV + k] += c_j*v[k];
+ /// Evaluate divergence at bound element in local coordinates
+ Range operator()(Domain const& x) const
+ {
+ return evaluate_node(x, bool_t<SubTree::isPower && SubTree::degree() == GridView::dimensionworld>{});
+ }
+
+private:
+ Range evaluate_node(Domain const& x, std::false_type) const
+ {
+ error_exit("Cannot calculate divergence(node) where node is not a power node.");
+ return Range(0);
+ }
+
+ Range evaluate_node(Domain const& x, std::true_type) const
+ {
+ static_assert(Size_v<Range> == 1, "Implemented for scalar coefficients only.");
+
+ assert( Super::bound_ );
+ Range dy(0);
+
+ auto&& coefficients = *Super::globalFunction_.dofVector_;
+ auto&& node = *Super::subTree_;
+
+ auto localBasis = makeLocalToGlobalBasisAdapter(node.child(0), Super::geometry_.value());
+ auto const& gradients = localBasis.gradientsAt(x);
+
+ auto re = Dune::Functions::flatVectorView(dy);
+ assert(int(re.size()) == 1);
+ for (std::size_t i = 0; i < localBasis.size(); ++i) {
+ auto grad = Dune::Functions::flatVectorView(gradients[i]);
+
+ assert(int(grad.size()) == GridView::dimensionworld);
+ for (std::size_t j = 0; j < GridView::dimensionworld; ++j) {
+ auto&& multiIndex = Super::localView_.index(node.child(j).localIndex(i));
+ re[0] += coefficients[multiIndex] * grad[j];
}
}
- });
- return y;
-}
+ return dy;
+ }
+};
template <class GB, class VT, class TP>
-typename DiscreteFunction<GB,VT,TP>::GradientLocalFunction::Range DiscreteFunction<GB,VT,TP>::
-GradientLocalFunction::operator()(Domain const& x) const
+class DiscreteFunction<GB,VT,TP>::PartialLocalFunction
+ : public DerivativeLocalFunctionBase<tag::partial>
{
- assert( bound_ );
- Range dy(0);
+ using Super = DerivativeLocalFunctionBase<tag::partial>;
- auto&& coefficients = *globalFunction_.dofVector_;
- auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
- for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
- {
- // TODO: may DOFVectorView::Range to FieldVector type if necessary
- using LocalDerivativeTraits
- = Dune::Functions::DefaultDerivativeTraits<VT(Domain)>;
- using GradientBlock = typename LocalDerivativeTraits::Range;
-
- auto&& fe = node.finiteElement();
- auto&& localBasis = fe.localBasis();
- std::size_t size = localBasis.size();
-
- NodeCache<TYPEOF(node)> localBasisCache(localBasis);
- auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element().type(),x);
-
- // The transposed inverse Jacobian of the map from the reference element to the element
- auto&& jacobian = geometry_.value().jacobianInverseTransposed(x);
+public:
+ using Range = typename Super::Range;
+ using Domain = typename Super::Domain;
- // Compute the shape function gradients on the real element
- std::vector<GradientBlock> gradients(referenceGradients.size());
- for (std::size_t i = 0; i < gradients.size(); ++i) // J^(-T) * D[phi] -> grad^T
- Dune::MatVec::as_vector(gradients[i]) = jacobian * Dune::MatVec::as_vector(referenceGradients[i]);
+ using Super::DerivativeLocalFunctionBase;
- // Get range entry associated to this node
- auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
+ /// Evaluate partial derivative at bound element in local coordinates
+ Range operator()(Domain const& x) const
+ {
+ assert( Super::bound_ );
+ Range dy(0);
+
+ auto&& coefficients = *Super::globalFunction_.dofVector_;
+ auto&& nodeToRangeEntry = Super::globalFunction_.nodeToRangeEntry_;
+ for_each_leaf_node(*Super::subTree_, [&,this](auto const& node, auto const& tp)
+ {
+ auto localBasis = makeLocalToGlobalBasisAdapter(node, Super::geometry_.value());
+ auto const& partial = localBasis.partialsAt(x, Super::type_.comp);
+
+ // Get range entry associated to this node
+ auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
+
+ for (std::size_t i = 0; i < localBasis.size(); ++i) {
+ auto&& multiIndex = Super::localView_.index(node.localIndex(i));
+
+ // Get coefficient associated to i-th shape function
+ auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
+
+ // Get value of i-th transformed reference partial_derivative
+ auto d_comp = Dune::Functions::flatVectorView(partial[i]);
+
+ std::size_t dimC = c.size();
+ std::size_t dimV = d_comp.size();
+ assert(dimC*dimV == std::size_t(re.size()));
+ for(std::size_t j = 0; j < dimC; ++j) {
+ auto&& c_j = c[j];
+ for(std::size_t k = 0; k < dimV; ++k)
+ re[j*dimV + k] += c_j*d_comp[k];
+ }
+ }
+ });
- for (std::size_t i = 0; i < size; ++i) {
- auto&& multiIndex = localView_.index(node.localIndex(i));
+ return dy;
+ }
- // Get coefficient associated to i-th shape function
- auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
+private:
+ template <class T, class J, class RG>
+ void test_partial(std::vector<T> const& partial, std::size_t comp, J const& jacobian, RG const& referenceGradients) const
+ {
+ for (std::size_t i = 0; i < referenceGradients.size(); ++i) { // J^(-T) * D[phi] -> grad^T
+ auto grad = jacobian * Dune::MatVec::as_vector(referenceGradients[i]);
+ auto grad_ = Dune::Functions::flatVectorView(grad);
- // Get value of i-th transformed reference gradient
- auto grad = Dune::Functions::flatVectorView(gradients[i]);
+ if (std::abs(grad_[comp] - partial[i]) > 1.e-10) {
+ print_error(partial[i], comp, grad, jacobian, referenceGradients[i][0]);
+ error_exit("wrong partial derivatives");
+ }
+ }
+ }
- std::size_t dimC = c.size();
- std::size_t dimV = grad.size();
- assert(dimC*dimV == std::size_t(re.size()));
- for(std::size_t j = 0; j < dimC; ++j) {
- auto&& c_j = c[j];
- for(std::size_t k = 0; k < dimV; ++k)
- re[j*dimV + k] += c_j*grad[k];
+ template <class T, class G, class J, class RG>
+ void print_error(T const& partial_i, std::size_t comp, G const& grad_i, J const& jacobian, RG const& referenceGradient) const
+ {
+ msg(__PRETTY_FUNCTION__);
+ msg("jacobian = {}", jacobian);
+ auto&& jacobian_mat = Dune::MatVec::as_matrix(jacobian);
+ msg("jacobian_ = [");
+ for (std::size_t i = 0; i < jacobian_mat.N(); ++i) {
+ for (std::size_t j = 0; j < jacobian_mat.M(); ++j) {
+ msg_("{} ",jacobian_mat[i][j]);
}
+ msg(";");
}
- });
+ msg("]");
- return dy;
-}
+ msg("ref-grad = {}", referenceGradient);
+ auto grad_ = Dune::Functions::flatVectorView(grad_i);
+ msg("comp = {}, grad_ = ({}, {}), partial = {}", comp, grad_[0], grad_[1], partial_i);
+ }
+};
template <class GB, class VT, class TP>
diff --git a/src/amdis/gridfunctions/FunctorGridFunction.hpp b/src/amdis/gridfunctions/FunctorGridFunction.hpp
index 2054b342..3c6fe3c4 100644
--- a/src/amdis/gridfunctions/FunctorGridFunction.hpp
+++ b/src/amdis/gridfunctions/FunctorGridFunction.hpp
@@ -8,7 +8,9 @@
#include <amdis/Operations.hpp>
#include <amdis/common/ForEach.hpp>
#include <amdis/common/Logical.hpp>
-#include <amdis/gridfunctions/GridFunctionConcepts.hpp>
+#include <amdis/gridfunctions/Derivative.hpp>
+#include <amdis/gridfunctions/GridFunction.hpp>
+#include <amdis/gridfunctions/Order.hpp>
namespace AMDiS
{
@@ -120,9 +122,9 @@ namespace AMDiS
* **Requirements:**
* - The Functor `F` must model `Concepts::HasPartial`
**/
- template <class Sig, class F, class... LFs,
+ template <class Sig, class F, class... LFs, class Type,
REQUIRES(Concepts::HasPartial<F>)>
- auto derivative(FunctorLocalFunction<Sig,F,LFs...> const& lf)
+ auto derivative(FunctorLocalFunction<Sig,F,LFs...> const& lf, Type const& type)
{
auto index_seq = std::index_sequence_for<LFs...>{};
@@ -135,7 +137,7 @@ namespace AMDiS
using std::get;
auto const& lgfs_i = *get<VALUE(_i)>(lf.localFcts());
- return makeFunctorGridFunction(Operation::Multiplies{}, di_f, derivative(lgfs_i));
+ return makeFunctorGridFunction(Operation::Multiplies{}, di_f, derivative(lgfs_i, type));
};
// sum_i [ d_i(f)[lgfs...] * derivative(lgfs_i)
@@ -226,7 +228,7 @@ namespace AMDiS
return std::get<0>(gridFcts_).entitySet();
}
- LocalFunction localFunction() const
+ LocalFunction makeLocalFunction() const
{
return Tools::apply([&](auto const&... gridFcts) { return LocalFunction{fct_, gridFcts...}; }, gridFcts_);
}
@@ -237,16 +239,6 @@ namespace AMDiS
};
- /// \fn localFunction
- /// \brief Creates a LocalFunction from the LocalFunctions of the GridFunctions.
- /// \relates FunctorLocalFunction
- template <class F, class... GFs>
- auto localFunction(FunctorGridFunction<F,GFs...> const& gf)
- {
- return gf.localFunction();
- }
-
-
// Generator function for FunctorGridFunction expressions
template <class Functor, class... GridFcts>
auto makeFunctorGridFunction(Functor const& f, GridFcts const&... gridFcts)
diff --git a/src/amdis/gridfunctions/GridFunctionConcepts.hpp b/src/amdis/gridfunctions/GridFunction.hpp
similarity index 62%
rename from src/amdis/gridfunctions/GridFunctionConcepts.hpp
rename to src/amdis/gridfunctions/GridFunction.hpp
index 14c62890..434ec4c6 100644
--- a/src/amdis/gridfunctions/GridFunctionConcepts.hpp
+++ b/src/amdis/gridfunctions/GridFunction.hpp
@@ -4,16 +4,29 @@
#include <amdis/common/Concepts.hpp>
#include <amdis/common/Logical.hpp>
+#include <amdis/gridfunctions/Derivative.hpp>
+#include <amdis/gridfunctions/Order.hpp>
namespace AMDiS
{
+ // The localFunction of a GridFunction
+ template <class GridFunction,
+ class = void_t<decltype(std::declval<GridFunction>().makeLocalFunction())> >
+ auto localFunction(GridFunction const& gf)
+ {
+ return gf.makeLocalFunction();
+ }
+
+
namespace Traits
{
template <class T>
struct IsPreGridFunction
: std::false_type {};
+
} // end namespace Traits
+
namespace Concepts
{
/** \addtogroup Concepts
@@ -28,36 +41,6 @@ namespace AMDiS
auto requires_(F&& f) -> decltype( localFunction(f) );
};
- struct HasDerivative
- {
- template <class F>
- auto requires_(F&& f) -> decltype( derivative(f) );
- };
-
- struct HasLocalFunctionDerivative
- {
- template <class F>
- auto requires_(F&& f) -> decltype( derivative(localFunction(f)) );
- };
-
- struct HasLocalFunctionOrder
- {
- template <class F>
- auto requires_(F&& f) -> decltype( order(localFunction(f)) );
- };
-
- struct HasPartial
- {
- template <class F>
- auto requires_(F&& f) -> decltype( partial(f, index_<0>) );
- };
-
- struct HasOrder
- {
- template <class F>
- auto requires_(F&& f) -> decltype( order(f) );
- };
-
struct HasGridFunctionTypes
{
template <class GF>
@@ -67,10 +50,6 @@ namespace AMDiS
typename GF::EntitySet >;
};
- template <class F, int dow>
- constexpr bool CallableDow =
- Concepts::Callable<F, Dune::FieldVector<double, dow>>;
-
} // end namespace Definition
@@ -78,38 +57,6 @@ namespace AMDiS
template <class GF>
constexpr bool HasLocalFunction = models<Definition::HasLocalFunction(GF)>;
- /// \brief GridFunction GF has free function `derivative(F)`
- template <class GF>
- constexpr bool HasDerivative = models<Definition::HasDerivative(GF)>;
-
- /// \brief GridFunction GF has free function `derivative(localFunction(F))`
- template <class GF>
- constexpr bool HasLocalFunctionDerivative = models<Definition::HasLocalFunctionDerivative(GF)>;
-
- /// \brief GridFunction GF has free function `order(localFunction(F))`
- template <class GF>
- constexpr bool HasLocalFunctionOrder = models<Definition::HasLocalFunctionOrder(GF)>;
-
- /// \brief Functor F has free function `partial(F,_0)`
- template <class F>
- constexpr bool HasPartial = models<Definition::HasPartial(F)>;
-
- /// \brief LocalFuncion LF has free function `order(F)`
- template <class LF>
- constexpr bool HasOrder = models<Definition::HasOrder(LF)>;
-
-
- /// \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
-
/// \brief GridFunction GF is a Type that has LocalFunction and provides some
/// typedefs for `Domain`, `Range`, and `EntitySet`.
template <class GF>
@@ -127,23 +74,6 @@ namespace AMDiS
} // end namespace Concepts
- template <class Traits,
- bool valid = !std::is_same<typename Traits::Range, Dune::Functions::InvalidRange>::value>
- constexpr bool isValidRange()
- {
- static_assert(valid, "Traits::Range is invalid");
- return valid;
- }
-
-
- // polynomial order of local functions
- template <class LocalFunction,
- class = void_t<decltype(std::declval<LocalFunction>().order())> >
- int order(LocalFunction const& lf)
- {
- return lf.order();
- }
-
#ifndef DOXYGEN
template <class PreGridFct, class = void>
diff --git a/src/amdis/gridfunctions/Order.hpp b/src/amdis/gridfunctions/Order.hpp
new file mode 100644
index 00000000..ff2cef28
--- /dev/null
+++ b/src/amdis/gridfunctions/Order.hpp
@@ -0,0 +1,51 @@
+#pragma once
+
+#include <amdis/common/Concepts.hpp>
+
+namespace AMDiS
+{
+ // polynomial order of local functions
+ template <class LocalFunction,
+ class = void_t<decltype(std::declval<LocalFunction>().order())> >
+ int order(LocalFunction const& lf)
+ {
+ return lf.order();
+ }
+
+
+ namespace Concepts
+ {
+ /** \addtogroup Concepts
+ * @{
+ **/
+
+ namespace Definition
+ {
+ struct HasLocalFunctionOrder
+ {
+ template <class F>
+ auto requires_(F&& f) -> decltype( order(localFunction(f)) );
+ };
+
+ struct HasOrder
+ {
+ template <class F>
+ auto requires_(F&& f) -> decltype( order(f) );
+ };
+
+ } // end namespace Definition
+
+
+ /// \brief GridFunction GF has free function `order(localFunction(F))`
+ template <class GF>
+ constexpr bool HasLocalFunctionOrder = models<Definition::HasLocalFunctionOrder(GF)>;
+
+ /// \brief LocalFuncion LF has free function `order(F)`
+ template <class LF>
+ constexpr bool HasOrder = models<Definition::HasOrder(LF)>;
+
+ /** @} **/
+
+ } // end namespace Concepts
+
+} // end namespace AMDiS
diff --git a/src/amdis/localoperators/FirstOrderPartialTest.hpp b/src/amdis/localoperators/FirstOrderPartialTest.hpp
index 46fcb9f3..eb0192ef 100644
--- a/src/amdis/localoperators/FirstOrderPartialTest.hpp
+++ b/src/amdis/localoperators/FirstOrderPartialTest.hpp
@@ -43,40 +43,22 @@ namespace AMDiS
static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only.");
auto const& quad = this->getQuadratureRule(contextGeo.type(), node);
- auto const& localFE = node.finiteElement();
- std::size_t feSize = localFE.size();
+ auto localBasis = makeLocalToGlobalBasisAdapter(node, contextGeo.geometry());
- using RangeFieldType = typename NodeQuadCache<Node>::Traits::RangeFieldType;
- NodeQuadCache<Node> cache(localFE.localBasis());
-
- auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(contextGeo, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
decltype(auto) local = contextGeo.local(quad[iq].position());
- // The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
- assert(jacobian.M() == CG::dim);
-
// The multiplicative factor in the integral transformation formula
- const auto dx = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
- const auto c = Super::coefficient(local);
-
- // The gradients of the shape functions on the reference element
- auto const& shapeGradients = shapeGradientsCache[iq];
+ auto dx = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
+ dx *= Super::coefficient(local);
// Compute the shape function gradients on the real element
- thread_local std::vector<RangeFieldType> partial;
- partial.resize(shapeGradients.size());
- for (std::size_t i = 0; i < partial.size(); ++i) {
- partial[i] = jacobian[comp_][0] * shapeGradients[i][0][0];
- for (std::size_t j = 1; j < jacobian.M(); ++j)
- partial[i] += jacobian[comp_][j] * shapeGradients[i][0][j];
- }
+ auto const& partial = localBasis.partialsAt(local, comp_);
- for (std::size_t j = 0; j < feSize; ++j) {
+ for (std::size_t j = 0; j < localBasis.size(); ++j) {
const auto local_j = node.localIndex(j);
- elementVector[local_j] += dx * (c * partial[j]);
+ elementVector[local_j] += dx * partial[j];
}
}
}
diff --git a/src/amdis/localoperators/FirstOrderTestPartialTrial.hpp b/src/amdis/localoperators/FirstOrderTestPartialTrial.hpp
index b3d67ff3..08a5c6f5 100644
--- a/src/amdis/localoperators/FirstOrderTestPartialTrial.hpp
+++ b/src/amdis/localoperators/FirstOrderTestPartialTrial.hpp
@@ -3,9 +3,9 @@
#include <type_traits>
#include <amdis/GridFunctionOperator.hpp>
-#include <amdis/utility/LocalBasisCache.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/StaticSize.hpp>
+#include <amdis/utility/LocalToGlobalAdapter.hpp>
namespace AMDiS
{
@@ -51,49 +51,30 @@ namespace AMDiS
"Operator can be applied to Leaf-Nodes only.");
auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
+
auto const& rowLocalFE = rowNode.finiteElement();
- auto const& colLocalFE = colNode.finiteElement();
std::size_t rowSize = rowLocalFE.size();
- std::size_t colSize = colLocalFE.size();
-
- using RangeFieldType = typename NodeQuadCache<CN>::Traits::RangeFieldType;
-
NodeQuadCache<RN> rowCache(rowLocalFE.localBasis());
- NodeQuadCache<CN> colCache(colLocalFE.localBasis());
+ auto colLocalBasis = makeLocalToGlobalBasisAdapter(colNode, contextGeo.geometry());
auto const& shapeValuesCache = rowCache.evaluateFunctionAtQP(contextGeo, quad);
- auto const& shapeGradientsCache = colCache.evaluateJacobianAtQP(contextGeo, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
decltype(auto) local = contextGeo.local(quad[iq].position());
- // The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
- auto&& jacobian_mat = Dune::MatVec::as_matrix(jacobian);
- assert(jacobian.M() == CG::dim);
-
// The multiplicative factor in the integral transformation formula
- const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
- const auto c = Super::coefficient(local);
+ auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
+ factor *= Super::coefficient(local);
// the values of the shape functions on the reference element at the quadrature point
auto const& shapeValues = shapeValuesCache[iq];
- // The gradients of the shape functions on the reference element
- auto const& shapeGradients = shapeGradientsCache[iq];
-
// Compute the shape function gradients on the real element
- thread_local std::vector<RangeFieldType> colPartial;
- colPartial.resize(shapeGradients.size());
- for (std::size_t i = 0; i < colPartial.size(); ++i) {
- colPartial[i] = jacobian_mat[comp_][0] * shapeGradients[i][0][0];
- for (std::size_t j = 1; j < jacobian_mat.M(); ++j)
- colPartial[i] += jacobian_mat[comp_][j] * shapeGradients[i][0][j];
- }
+ auto const& colPartial = colLocalBasis.partialsAt(local, comp_);
- for (std::size_t j = 0; j < colSize; ++j) {
+ for (std::size_t j = 0; j < colLocalBasis.size(); ++j) {
const auto local_j = colNode.localIndex(j);
- const auto value = factor * (c * colPartial[j]);
+ const auto value = factor * colPartial[j];
for (std::size_t i = 0; i < rowSize; ++i) {
const auto local_i = rowNode.localIndex(i);
elementMatrix[local_i][local_j] += value * shapeValues[i];
diff --git a/src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp b/src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp
index 32dd3b3e..1724a6a9 100644
--- a/src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp
+++ b/src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp
@@ -3,9 +3,9 @@
#include <type_traits>
#include <amdis/GridFunctionOperator.hpp>
-#include <amdis/utility/LocalBasisCache.hpp>
#include <amdis/Output.hpp>
#include <amdis/common/StaticSize.hpp>
+#include <amdis/utility/LocalToGlobalAdapter.hpp>
namespace AMDiS
{
@@ -44,60 +44,28 @@ namespace AMDiS
template <class CG, class RN, class CN, class Mat>
void getElementMatrix(CG const& contextGeo, RN const& rowNode, CN const& colNode, Mat& elementMatrix)
{
- static_assert(RN::isLeaf && CN::isLeaf,
- "Operator can be applied to Leaf-Nodes only.");
+ static_assert(RN::isLeaf && CN::isLeaf, "Operator can be applied to Leaf-Nodes only.");
auto const& quad = this->getQuadratureRule(contextGeo.type(), rowNode, colNode);
- auto const& rowLocalFE = rowNode.finiteElement();
- auto const& colLocalFE = colNode.finiteElement();
- std::size_t rowSize = rowLocalFE.size();
- std::size_t colSize = colLocalFE.size();
- using RowFieldType = typename NodeQuadCache<RN>::Traits::RangeFieldType;
- using ColFieldType = typename NodeQuadCache<CN>::Traits::RangeFieldType;
+ auto rowLocalBasis = makeLocalToGlobalBasisAdapter(rowNode, contextGeo.geometry());
+ auto colLocalBasis = makeLocalToGlobalBasisAdapter(colNode, contextGeo.geometry());
- NodeQuadCache<RN> rowCache(rowLocalFE.localBasis());
- NodeQuadCache<CN> colCache(colLocalFE.localBasis());
-
- auto const& rowGradientsCache = rowCache.evaluateJacobianAtQP(contextGeo, quad);
- auto const& colGradientsCache = colCache.evaluateJacobianAtQP(contextGeo, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
decltype(auto) local = contextGeo.local(quad[iq].position());
- // The transposed inverse Jacobian of the map from the reference element to the element
- const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
- auto&& jacobian_mat = Dune::MatVec::as_matrix(jacobian);
-
// The multiplicative factor in the integral transformation formula
- const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
- const auto c = Super::coefficient(local);
-
- // The gradients of the shape functions on the reference element
- auto const& rowShapeGradients = rowGradientsCache[iq];
- auto const& colShapeGradients = colGradientsCache[iq];
-
- // Compute the shape function gradients on the real element
- thread_local std::vector<RowFieldType> rowPartial;
- rowPartial.resize(rowShapeGradients.size());
- for (std::size_t i = 0; i < rowPartial.size(); ++i) {
- rowPartial[i] = jacobian_mat[compTest_][0] * rowShapeGradients[i][0][0];
- for (std::size_t j = 1; j < jacobian_mat.M(); ++j)
- rowPartial[i] += jacobian_mat[compTest_][j] * rowShapeGradients[i][0][j];
- }
+ auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
+ factor *= Super::coefficient(local);
- thread_local std::vector<ColFieldType> colPartial;
- colPartial.resize(colShapeGradients.size());
- for (std::size_t i = 0; i < colPartial.size(); ++i) {
- colPartial[i] = jacobian_mat[compTrial_][0] * colShapeGradients[i][0][0];
- for (std::size_t j = 1; j < jacobian_mat.M(); ++j)
- colPartial[i] += jacobian_mat[compTrial_][j] * colShapeGradients[i][0][j];
- }
+ auto const& rowPartial = rowLocalBasis.partialsAt(local, compTest_);
+ auto const& colPartial = colLocalBasis.partialsAt(local, compTrial_);
- for (std::size_t j = 0; j < colSize; ++j) {
+ for (std::size_t j = 0; j < colLocalBasis.size(); ++j) {
const auto local_j = colNode.localIndex(j);
- const auto value = factor * (c * colPartial[j]);
- for (std::size_t i = 0; i < rowSize; ++i) {
+ const auto value = factor * colPartial[j];
+ for (std::size_t i = 0; i < rowLocalBasis.size(); ++i) {
const auto local_i = rowNode.localIndex(i);
elementMatrix[local_i][local_j] += value * rowPartial[i];
}
diff --git a/src/amdis/localoperators/ZeroOrderTest.hpp b/src/amdis/localoperators/ZeroOrderTest.hpp
index a6bf83bf..a064d44f 100644
--- a/src/amdis/localoperators/ZeroOrderTest.hpp
+++ b/src/amdis/localoperators/ZeroOrderTest.hpp
@@ -51,8 +51,8 @@ namespace AMDiS
decltype(auto) local = contextGeo.local(quad[iq].position());
// The multiplicative factor in the integral transformation formula
- const auto factor = Super::coefficient(local) * contextGeo.integrationElement(quad[iq].position())
- * quad[iq].weight();
+ auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
+ factor *= Super::coefficient(local);
auto const& shapeValues = shapeValuesCache[iq];
for (std::size_t i = 0; i < size; ++i) {
diff --git a/src/amdis/utility/LocalBasisCache.hpp b/src/amdis/utility/LocalBasisCache.hpp
index 1ee17656..3ef62785 100644
--- a/src/amdis/utility/LocalBasisCache.hpp
+++ b/src/amdis/utility/LocalBasisCache.hpp
@@ -68,6 +68,10 @@ namespace AMDiS
: localBasis_(&localBasis)
{}
+ NodeCache(Node const& node)
+ : NodeCache(node.finiteElement().localBasis())
+ {}
+
/// Bind the cache to a local basis for later evaluations.
// NOTE: Should be called before the first evaluation of the local basis function or
// jacobian.
@@ -113,6 +117,12 @@ namespace AMDiS
LocalBasisType const* localBasis_ = nullptr;
};
+ /// Generator function for \ref NodeCache, \relates NodeCache.
+ template <class Node>
+ NodeCache<Node> makeNodeCache(Node const& node)
+ {
+ return NodeCache<Node>{node};
+ }
/// \brief Cache of LocalBasis evaluations and jacobians at quadrature points
/**
@@ -174,6 +184,10 @@ namespace AMDiS
: localBasis_(&localBasis)
{}
+ NodeQuadCache(Node const& node)
+ : NodeQuadCache(node.finiteElement().localBasis())
+ {}
+
/// Bind the cache to a local basis for later evaluations.
// NOTE: Should be called before the first evaluation of the local basis function or
// jacobian.
@@ -223,4 +237,11 @@ namespace AMDiS
LocalBasisType const* localBasis_ = nullptr;
};
+ /// Generator function for \ref NodeQuadCache, \relates NodeQuadCache.
+ template <class Node>
+ NodeQuadCache<Node> makeNodeQuadCache(Node const& node)
+ {
+ return NodeQuadCache<Node>{node};
+ }
+
} // end namespace AMDiS
diff --git a/src/amdis/utility/LocalToGlobalAdapter.hpp b/src/amdis/utility/LocalToGlobalAdapter.hpp
new file mode 100644
index 00000000..b0df3522
--- /dev/null
+++ b/src/amdis/utility/LocalToGlobalAdapter.hpp
@@ -0,0 +1,207 @@
+#pragma once
+
+#include <cstddef>
+#include <vector>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/geometry/type.hh>
+
+#include <amdis/common/DerivativeTraits.hpp>
+#include <amdis/common/FieldMatVec.hpp>
+#include <amdis/utility/LocalBasisCache.hpp>
+
+namespace AMDiS
+{
+ /// Traits class for local-to-global basis adaptors
+ /**
+ * \tparam LocalBasisTraits Traits class of the LocalBasis to be adapted.
+ * \tparam dimGlobal Dimension of the global coordinates,
+ * i.e. Geometry::coorddimension, if the global
+ * coordinates are determined by a Geometry.
+ *
+ * \implements BasisInterface::Traits
+ */
+ template <class LocalBasisTraits, std::size_t dimGlobal>
+ struct LocalToGlobalBasisAdapterTraits
+ {
+ static const std::size_t dimDomainLocal = LocalBasisTraits::dimDomain;
+ static const std::size_t dimDomainGlobal = dimGlobal;
+ static const std::size_t dimRange = LocalBasisTraits::dimRange;
+
+ using DomainField = typename LocalBasisTraits::DomainFieldType;
+ using DomainLocal = typename LocalBasisTraits::DomainType;
+ using DomainGlobal = FieldVector<DomainField, dimDomainGlobal>;
+
+ using RangeField = typename LocalBasisTraits::RangeFieldType;
+ using Range = typename LocalBasisTraits::RangeType;
+
+ using GradientRange = typename DerivativeTraits<Range(DomainGlobal), tag::gradient>::Range;
+ using PartialRange = typename DerivativeTraits<Range(DomainGlobal), tag::partial>::Range;
+ };
+
+ /// \brief Convert a simple (scalar) local basis into a global basis
+ /**
+ * The local basis must be scalar, i.e. LocalBasis::Traits::dimRange must be 1
+ * It's values are not transformed.
+ *
+ * For scalar function \f$f\f$, the gradient is equivalent to the transposed
+ * Jacobian \f$\nabla f|_x = J_f^T(x)\f$. The Jacobian is thus transformed
+ * using
+ * \f[
+ * \nabla f|_{\mu(\hat x)} =
+ * \hat J_\mu^{-T}(\hat x) \cdot \hat\nabla\hat f|_{\hat x}
+ * \f]
+ * Here the hat \f$\hat{\phantom x}\f$ denotes local quantities and
+ * \f$\mu\f$ denotes the local-to-global map of the geometry.
+ *
+ * \tparam BasisNode Type of the typetree node containing a local basis to adopt.
+ * \tparam Geometry Type of the local-to-global transformation.
+ *
+ * NOTE: The adapter implements a caching of local basis evaluations at coordinates.
+ */
+ template <class BasisNode, class Geometry>
+ class LocalToGlobalBasisAdapter
+ {
+ using LocalBasis = typename BasisNode::FiniteElement::Traits::LocalBasisType;
+
+ static_assert(std::is_same<typename LocalBasis::Traits::DomainFieldType, typename Geometry::ctype>::value,
+ "LocalToGlobalBasisAdapter: LocalBasis must use the same ctype as Geometry");
+
+ static_assert(std::size_t(LocalBasis::Traits::dimDomain) == std::size_t(Geometry::mydimension),
+ "LocalToGlobalBasisAdapter: LocalBasis domain dimension must match local dimension of Geometry");
+
+ public:
+ using Cache = NodeCache<BasisNode>;
+ using Traits = LocalToGlobalBasisAdapterTraits<typename LocalBasis::Traits, Geometry::coorddimension>;
+
+ /// \biref Construct a LocalToGlobalBasisAdapter
+ /**
+ * \param node The basis node in the typetree containing the local basis to adopt.
+ * \param geometry The geometry object to use for adaption.
+ *
+ * \note This class stores the references passed here. Any use of this
+ * class after these references have become invalid results in
+ * undefined behaviour. The exception is that the destructor of
+ * this class may still be called.
+ */
+ LocalToGlobalBasisAdapter(BasisNode const& node, Geometry const& geometry)
+ : localBasis_(node.finiteElement().localBasis())
+ , localBasisCache_(localBasis_)
+ , geometry_(geometry)
+ , size_(localBasis_.size())
+ {}
+
+ /// Return the number of local basis functions
+ std::size_t size() const { return size_; }
+
+ /// \brief Return maximum polynomial order of the base function
+ /**
+ * This is to determine the required quadrature order. For an affine
+ * geometry this is the same order as for the local basis. For other
+ * geometries this returns the order of the local basis plus the global
+ * dimension minus 1. The assumption for non-affine geometries is that
+ * they are still multi-linear.
+ */
+ std::size_t order() const
+ {
+ if (geometry_.affine())
+ // affine linear
+ return localBasis_.order();
+ else
+ // assume at most order dim
+ return localBasis_.order() + Traits::dimDomainGlobal - 1;
+ }
+
+ /// Evaluate the local basis functions in the local coordinate `in`
+ void evaluateFunction(typename Traits::DomainLocal const& x,
+ std::vector<typename Traits::Range>& out) const
+ {
+ out = localBasisCache_.evaluateFunction(geometry_.type(), x);
+ }
+
+ /// Evaluate the local basis functions in the local coordinate `in` and
+ /// return the result using a reference to a thread_local (or static) vector.
+ auto const& valuesAt(typename Traits::DomainLocal const& x) const
+ {
+ return localBasisCache_.evaluateFunction(geometry_.type(), x);
+ }
+
+ /// Return the full (global) gradient of the local basis functions in
+ /// the local coordinate `in`
+ void evaluateGradient(typename Traits::DomainLocal const& x,
+ std::vector<typename Traits::GradientRange>& out) const
+ {
+ auto const& localJacobian = localBasisCache_.evaluateJacobian(geometry_.type(), x);
+ auto const& geoJacobian = geometry_.jacobianInverseTransposed(x);
+
+ out.resize(size_);
+ for (std::size_t i = 0; i < size_; ++i)
+ geoJacobian.mv(Dune::MatVec::as_vector(localJacobian[i]),
+ Dune::MatVec::as_vector(out[i]));
+ }
+
+ /// Return the full (global) gradient of the local basis functions in
+ /// the local coordinate `in` and return the result using a reference
+ /// to a thread_local (or static) vector.
+ auto const& gradientsAt(typename Traits::DomainLocal const& x) const
+ {
+ thread_local std::vector<typename Traits::GradientRange> grad;
+ evaluateGradient(x, grad);
+ return grad;
+ }
+
+ /// Return the (global) partial derivative in direction `comp` of the
+ /// local basis functions in the local coordinate `in`
+ void evaluatePartial(typename Traits::DomainLocal const& x,
+ std::size_t comp,
+ std::vector<typename Traits::PartialRange>& out) const
+ {
+ auto const& localJacobian = localBasisCache_.evaluateJacobian(geometry_.type(), x);
+ // NOTE: geoJacobian might be a Dune::DiagonalMatrix with limited interface!
+ auto const& geoJacobian = geometry_.jacobianInverseTransposed(x);
+
+ out.resize(size_);
+ typename Traits::GradientRange grad;
+ auto&& grad_ = Dune::MatVec::as_vector(grad);
+ for (std::size_t i = 0; i < size_; ++i) {
+ geoJacobian.mv(Dune::MatVec::as_vector(localJacobian[i]), grad_);
+ out[i] = grad_[comp];
+ }
+ }
+
+ /// Return the (global) partial derivative in direction `comp` of the
+ /// local basis functions in the local coordinate `in` and return the
+ /// result using a reference to a thread_local (or static) vector.
+ auto const& partialsAt(typename Traits::DomainLocal const& x, std::size_t comp) const
+ {
+ thread_local std::vector<typename Traits::PartialRange> d_comp;
+ evaluatePartial(x, comp, d_comp);
+ return d_comp;
+ }
+
+ private:
+ /// Reference to the local basis bound to the adapter
+ LocalBasis const& localBasis_;
+
+ /// A basis cache for the evaluation at local coordinates
+ Cache localBasisCache_;
+
+ /// Reference to the geometry this adapter is bound to
+ Geometry const& geometry_;
+
+ /// The number of basis functions
+ std::size_t size_;
+ };
+
+
+ /// Generator function for the \ref LocalToGlobalBasisAdapter. \relates LocalToGlobalBasisAdapter.
+ template <class BasisNode, class Geometry>
+ LocalToGlobalBasisAdapter<BasisNode, Geometry>
+ makeLocalToGlobalBasisAdapter(BasisNode const& node, Geometry const& geometry)
+ {
+ return LocalToGlobalBasisAdapter<BasisNode, Geometry>{node, geometry};
+ }
+
+} // namespace AMDiS
diff --git a/src/amdis/utility/QuadratureFactory.hpp b/src/amdis/utility/QuadratureFactory.hpp
index 865c6cd8..cf326b53 100644
--- a/src/amdis/utility/QuadratureFactory.hpp
+++ b/src/amdis/utility/QuadratureFactory.hpp
@@ -5,7 +5,7 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>
-#include <amdis/GridFunctions.hpp>
+#include <amdis/gridfunctions/Order.hpp>
namespace AMDiS
{
@@ -40,7 +40,7 @@ namespace AMDiS
template <class LF>
- using HasLocalFunctionOrder = decltype( AMDiS::order(std::declval<LF>()) );
+ using HasLocalFunctionOrder = decltype( order(std::declval<LF>()) );
/// \brief Factory for quadrature rule, that calculates the coefficient order from
@@ -81,14 +81,6 @@ namespace AMDiS
return QuadFactoryFromLocalFunction<ctype, dim, LocalFunction>{};
}
- template <class ctype, int dim, class LocalFunction>
- auto makeQuadratureFactoryPtr(PreQuadFactoryFromLocalFunction const& /*pre*/)
- {
- using Factory = QuadFactoryFromLocalFunction<ctype, dim, LocalFunction>;
- return std::make_unique<Factory>();
- }
-
-
/// \brief Factory for quadrature rule, that takes to order of the localFunction
/// and a quadrature-type
@@ -137,14 +129,6 @@ namespace AMDiS
return QuadFactoryFromOrder<ctype, dim, LocalFunction>{pre.order, pre.qt};
}
- template <class ctype, int dim, class LocalFunction>
- auto makeQuadratureFactoryPtr(PreQuadFactoryFromOrder const& pre)
- {
- using Factory = QuadFactoryFromOrder<ctype, dim, LocalFunction>;
- return std::make_unique<Factory>(pre.order, pre.qt);
- }
-
-
/// \brief Factory for quadrature rule, that is based on an existing rule
template <class ctype, int dim, class LocalFunction>
@@ -187,11 +171,4 @@ namespace AMDiS
return QuadFactoryFromRule<ctype, dim, LocalFunction>{pre.rule};
}
- template <class ctype, int dim, class LocalFunction, class QuadRule>
- auto makeQuadratureFactoryPtr(PreQuadFactoryFromRule<QuadRule> const& pre)
- {
- using Factory = QuadFactoryFromRule<ctype, dim, LocalFunction>;
- return std::make_unique<Factory>(pre.rule);
- }
-
} // end namespace AMDiS
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index fe265d45..282214b3 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -46,6 +46,9 @@ dune_add_test(SOURCES FiniteElementTypeTest.cpp
dune_add_test(SOURCES FilesystemTest.cpp
LINK_LIBRARIES amdis)
+dune_add_test(SOURCES GradientTest.cpp
+ LINK_LIBRARIES amdis)
+
dune_add_test(SOURCES HybridSizeTest.cpp
LINK_LIBRARIES amdis)
diff --git a/test/DiscreteFunctionTest.cpp b/test/DiscreteFunctionTest.cpp
index e330eaaf..143811dd 100644
--- a/test/DiscreteFunctionTest.cpp
+++ b/test/DiscreteFunctionTest.cpp
@@ -75,9 +75,9 @@ int main(int argc, char** argv)
AMDIS_TEST( comp(U0, U1) );
AMDIS_TEST( comp(U0, U2) );
- auto du0 = derivative(u0);
+ auto du0 = derivative(u0, tag::gradient{});
auto localFct = localFunction(u0);
- auto dlocalFct = derivative(localFct);
+ auto dlocalFct = derivative(localFct, tag::gradient{});
for (auto const& e : elements(prob.gridView()))
{
localFct.bind(e);
diff --git a/test/GradientTest.cpp b/test/GradientTest.cpp
new file mode 100644
index 00000000..17b62282
--- /dev/null
+++ b/test/GradientTest.cpp
@@ -0,0 +1,91 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <iostream>
+
+#include <amdis/AMDiS.hpp>
+#include <amdis/LinearAlgebra.hpp>
+#include <amdis/GridFunctions.hpp>
+#include <amdis/Integrate.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+template <std::size_t p>
+void test()
+{
+ using Grid = Dune::YaspGrid<2>;
+ using Basis1 = LagrangeBasis<Grid,p,p>;
+
+ Grid grid({1.0, 1.0}, {16,16});
+ auto gridView = grid.leafGridView();
+ auto basis = Basis1::create(gridView);
+
+ auto uVector = makeDOFVector(basis);
+
+ auto u = makeDOFVectorView(uVector);
+ auto u_0 = makeDOFVectorView(uVector, 0);
+ auto u_1 = makeDOFVectorView(uVector, 1);
+
+ // eval a functor at global coordinates (at quadrature points)
+ u_0 << evalAtQP([](auto const& x) { return x[0] + x[1]; });
+ u_1 << evalAtQP([](auto const& x) { return -2*x[0] + 3*x[1]; });
+
+ // test gradient
+ AMDIS_TEST_APPROX((std::sqrt(integrate(unary_dot(gradientAtQP(u_0)), gridView))), std::sqrt(2.0));
+ AMDIS_TEST_APPROX((std::sqrt(integrate(unary_dot(gradientAtQP(u_1)), gridView))), std::sqrt(13.0));
+
+ // test divergence
+ AMDIS_TEST_APPROX((integrate(divergenceAtQP(u), gridView)), 4.0);
+
+ // test partial derivative
+ double d0u_0 = integrate(partialAtQP(u_0,0), gridView);
+ double d1u_1 = integrate(partialAtQP(u_1,1), gridView);
+ AMDIS_TEST_APPROX(d0u_0, 1.0);
+ AMDIS_TEST_APPROX(d1u_1, 3.0);
+
+ auto gf_d0u_0 = makeGridFunction(partialAtQP(u_0,0), gridView);
+ auto lf_d0u_0 = localFunction(gf_d0u_0);
+ for (auto const& e : elements(gridView)) {
+ lf_d0u_0.bind(e);
+ double v = lf_d0u_0(e.geometry().center());
+ AMDIS_TEST_APPROX(v, 1.0);
+ lf_d0u_0.unbind();
+ }
+
+ u_0 << evalAtQP([](auto const& x) { return std::sin(x[0]) + std::cos(x[1]); });
+ u_1 << evalAtQP([](auto const& x) { return -std::sin(x[0]) * std::cos(x[1]); });
+
+ AMDIS_TEST_APPROX((integrate(divergenceAtQP(u), gridView)), (integrate(partialAtQP(u_0,0) + partialAtQP(u_1,1), gridView)));
+
+ auto vVector(uVector);
+ auto v = makeDOFVectorView(vVector);
+ auto v_0 = makeDOFVectorView(vVector, 0);
+ auto v_1 = makeDOFVectorView(vVector, 1);
+
+ // test gradient
+ v << evalAtQP([](auto const& x) { return Dune::FieldVector<double,2>{std::cos(x[0]), -std::sin(x[1])}; });
+ AMDIS_TEST((std::sqrt(integrate(unary_dot(v - gradientAtQP(u_0)), gridView))) < 0.05);
+
+ // test divergence
+ v_0 << evalAtQP([](auto const& x) { return std::cos(x[0]) + std::sin(x[0])*std::sin(x[1]); });
+ AMDIS_TEST((integrate(v_0 - divergenceAtQP(u), gridView)) < 0.05);
+
+ // test partial derivative
+ v_0 << evalAtQP([](auto const& x) { return std::cos(x[0]); });
+ v_1 << evalAtQP([](auto const& x) { return std::sin(x[0])*std::sin(x[1]); });
+ AMDIS_TEST((integrate(v_0 - partialAtQP(u_0,0), gridView)) < 0.05);
+ AMDIS_TEST((integrate(v_1 - partialAtQP(u_1,1), gridView)) < 0.05);
+}
+
+int main(int argc, char** argv)
+{
+ Environment env(argc, argv);
+
+ test<1u>();
+ test<2u>();
+ test<3u>();
+
+ return report_errors();
+}
diff --git a/test/Tests.hpp b/test/Tests.hpp
index c44226e5..319e565d 100644
--- a/test/Tests.hpp
+++ b/test/Tests.hpp
@@ -129,7 +129,7 @@ namespace AMDiS
// implementation for ranges
template <class T1, class T2,
- class = decltype(((void)std::declval<T1>().cbegin(), (void)std::declval<T2>().cbegin())) >
+ class = decltype(std::declval<T1>().cbegin(), std::declval<T2>().cbegin()) >
inline void _test_approx(T1 const& expr_range, T2 const& range,
char const* expr_str, char const* value_str, char const* file, size_t line)
{
--
GitLab