Commit 32e77c49 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'issue/local_basis_cache' into 'develop'

Local basis cache depending on node type

See merge request !53
parents d68cbeaf 7465e2c2
Pipeline #1677 passed with stage
in 48 minutes and 26 seconds
#pragma once #pragma once
#include <tuple>
#include <vector> #include <vector>
#include <dune/common/hash.hh> #include <dune/common/hash.hh>
...@@ -7,99 +8,188 @@ ...@@ -7,99 +8,188 @@
#include <dune/geometry/type.hh> #include <dune/geometry/type.hh>
#include <amdis/utility/ConcurrentCache.hpp> #include <amdis/utility/ConcurrentCache.hpp>
#include <amdis/utility/FiniteElementType.hpp>
namespace AMDiS namespace AMDiS
{ {
template <class LocalBasisType> /// \brief Cache of LocalBasis evaluations and jacobians at local coordinates
class LocalBasisCache /**
* Caching is done using the ConcurrentCache data structure with a key that
* depends on the element type and local coordinates inside the element.
* Two methods are provided for evaluation of local basis functions and local
* basis jacobians at the local coordinate.
*
* \tparam Node A leaf basis node in a basis tree
**/
template <class Node>
class NodeCache
{ {
static_assert(Node::isLeaf, "");
using LocalBasisType = typename Node::FiniteElement::Traits::LocalBasisType;
using DomainType = typename LocalBasisType::Traits::DomainType;
public: public:
using Traits = typename LocalBasisType::Traits; using Traits = typename LocalBasisType::Traits;
using DomainType = typename Traits::DomainType; using ShapeValues = std::vector<typename Traits::RangeType>;
using RangeType = typename Traits::RangeType; using ShapeGradients = std::vector<typename Traits::JacobianType>;
using RangeFieldType = typename Traits::RangeFieldType;
using JacobianType = typename Traits::JacobianType;
using QuadratureRule = Dune::QuadratureRule<RangeFieldType,Traits::dimDomain>;
using QuadratureRules = Dune::QuadratureRules<RangeFieldType,Traits::dimDomain>;
private: private:
// Pair of GeometryType and local coordinates
/// Pair of GeometryType and quadrature order and size struct Key
template <class Tag>
struct QuadKey
{ {
unsigned int id; // topologyId unsigned int id; // topologyId
int p; // quadrature order DomainType local; // local coordinate
std::size_t size; // number of quadrature points
struct hasher struct hasher
{ {
std::size_t operator()(QuadKey const& t) const std::size_t operator()(Key const& t) const
{ {
std::size_t seed = 0; std::size_t seed = 0;
Dune::hash_combine(seed, t.id); Dune::hash_combine(seed, t.id);
Dune::hash_combine(seed, t.p); Dune::hash_range(seed, t.local.begin(), t.local.end());
Dune::hash_combine(seed, t.size);
return seed; return seed;
} }
}; };
friend bool operator==(QuadKey const& lhs, QuadKey const& rhs) friend bool operator==(Key const& lhs, Key const& rhs)
{ {
return std::tie(lhs.id,lhs.p,lhs.size) == std::tie(rhs.id,rhs.p,rhs.size); return std::tie(lhs.id,lhs.local) == std::tie(rhs.id,rhs.local);
} }
}; };
/// Pair of GeometryType and local coordinates public:
template <class Tag> /// Create a cache not associated to a local basis.
struct CoordsKey /// Requires a call to \ref bind() before the first evaluation of the cache.
NodeCache() = default;
/// Associate the cache to a localBasis.
NodeCache(LocalBasisType const& localBasis)
: localBasis_(&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.
// NOTE: Binding is only necessary, if the default constructor is used and
// the cache is not associated with a localBasis in the constructor.
void bind(LocalBasisType const& localBasis)
{
localBasis_ = &localBasis;
}
/// Evaluate local basis functions at local coordinates
ShapeValues const& evaluateFunction(Dune::GeometryType const& type, DomainType const& local) const
{
Key key{type.id(), local};
return shapeValues_.get(key, [&](Key const&)
{
ShapeValues data;
localBasis_->evaluateFunction(local, data);
return data;
});
}
/// Evaluate local basis jacobians at local coordinates
ShapeGradients const& evaluateJacobian(Dune::GeometryType const& type, DomainType const& local) const
{
Key key{type.id(), local};
return shapeGradients_.get(key, [&](Key const&)
{
ShapeGradients data;
localBasis_->evaluateJacobian(local, data);
return data;
});
}
private:
template <class Value>
using Cache = ConcurrentCache<Key, Value, ThreadLocalPolicy,
std::unordered_map<Key, Value, typename Key::hasher>>;
Cache<ShapeValues> shapeValues_;
Cache<ShapeGradients> shapeGradients_;
LocalBasisType const* localBasis_ = nullptr;
};
/// \brief Cache of LocalBasis evaluations and jacobians at quadrature points
/**
* Caching is done using the ConcurrentCache data structure with a key that
* depends on the element type and order and size of a quadrature rule.
* Two methods are provided for evaluation of local basis functions and local
* basis jacobians at all quadrature points. A vector of values is returned.
*
* \tparam Node A leaf basis node in a basis tree
**/
template <class Node>
class NodeQuadCache
{
static_assert(Node::isLeaf, "");
using LocalBasisType = typename Node::FiniteElement::Traits::LocalBasisType;
public:
using Traits = typename LocalBasisType::Traits;
using ShapeValues = std::vector<typename Traits::RangeType>;
using ShapeGradients = std::vector<typename Traits::JacobianType>;
using ShapeValuesContainer = std::vector<ShapeValues>;
using ShapeGradientsContainer = std::vector<ShapeGradients>;
private:
// Pair of GeometryType and quadrature order and size
struct Key
{ {
unsigned int id; // topologyId unsigned int id; // topologyId
DomainType local; // local coordinate int p; // quadrature order
std::size_t size; // number of quadrature points
struct hasher struct hasher
{ {
std::size_t operator()(CoordsKey const& t) const std::size_t operator()(Key const& t) const
{ {
std::size_t seed = 0; std::size_t seed = 0;
Dune::hash_combine(seed, t.id); Dune::hash_combine(seed, t.id);
Dune::hash_range(seed, t.local.begin(), t.local.end()); Dune::hash_combine(seed, t.p);
Dune::hash_combine(seed, t.size);
return seed; return seed;
} }
}; };
friend bool operator==(CoordsKey const& lhs, CoordsKey const& rhs) friend bool operator==(Key const& lhs, Key const& rhs)
{ {
return std::tie(lhs.id,lhs.local) == std::tie(rhs.id,rhs.local); return std::tie(lhs.id,lhs.p,lhs.size) == std::tie(rhs.id,rhs.p,rhs.size);
} }
}; };
using LocalValuesKey = CoordsKey<struct ValuesTag>;
using LocalGradientsKey = CoordsKey<struct GradientsTag>;
using ValuesKey = QuadKey<struct ValuesTag>;
using GradientsKey = QuadKey<struct GradientsTag>;
public: public:
using ShapeValues = std::vector<RangeType>; /// Create a cache not associated to a local basis.
using ShapeGradients = std::vector<JacobianType>; /// Requires a call to \ref bind() before the first evaluation of the cache.
NodeQuadCache() = default;
using ShapeValuesContainer = std::vector<ShapeValues>;
using ShapeGradientsContainer = std::vector<ShapeGradients>;
/// Associate the cache to a localBasis.
LocalBasisCache(LocalBasisType const& localBasis) NodeQuadCache(LocalBasisType const& localBasis)
: localBasis_(&localBasis) : localBasis_(&localBasis)
{} {}
// evaluate local basis-function at all quadrature points /// Bind the cache to a local basis for later evaluations.
template <class LocalContext, int d> // NOTE: Should be called before the first evaluation of the local basis function or
ShapeValuesContainer const& evaluateFunctionAtQp(LocalContext const& context, Dune::QuadratureRule<RangeFieldType,d> const& quad) const // jacobian.
// NOTE: Binding is only necessary, if the default constructor is used and
// the cache is not associated with a localBasis in the constructor.
void bind(LocalBasisType const& localBasis)
{
localBasis_ = &localBasis;
}
/// Evaluate local basis functions at all quadrature points of a quadrature rule.
template <class LocalContext, class QuadratureRule>
ShapeValuesContainer const& evaluateFunctionAtQP(LocalContext const& context, QuadratureRule const& quad) const
{ {
ValuesKey key{context.type().id(),quad.order(),quad.size()}; Key key{context.type().id(), quad.order(), quad.size()};
return ShapeValuesContainerCache::get(key, [&](ValuesKey const&) return shapeValuesContainer_.get(key, [&](Key const&)
{ {
ShapeValuesContainer data(quad.size()); ShapeValuesContainer data(quad.size());
for (std::size_t iq = 0; iq < quad.size(); ++iq) for (std::size_t iq = 0; iq < quad.size(); ++iq)
...@@ -108,12 +198,12 @@ namespace AMDiS ...@@ -108,12 +198,12 @@ namespace AMDiS
}); });
} }
// evaluate local basis-gradients at all quadrature points /// Evaluate local basis jacobians at all quadrature points of a quadrature rule.
template <class LocalContext, int d> template <class LocalContext, class QuadratureRule>
ShapeGradientsContainer const& evaluateJacobianAtQp(LocalContext const& context, Dune::QuadratureRule<RangeFieldType,d> const& quad) const ShapeGradientsContainer const& evaluateJacobianAtQP(LocalContext const& context, QuadratureRule const& quad) const
{ {
GradientsKey key{context.type().id(),quad.order(),quad.size()}; Key key{context.type().id(), quad.order(), quad.size()};
return ShapeGradientsContainerCache::get(key, [&](GradientsKey const&) return shapeGradientsContainer_.get(key, [&](Key const&)
{ {
ShapeGradientsContainer data(quad.size()); ShapeGradientsContainer data(quad.size());
for (std::size_t iq = 0; iq < quad.size(); ++iq) for (std::size_t iq = 0; iq < quad.size(); ++iq)
...@@ -122,45 +212,15 @@ namespace AMDiS ...@@ -122,45 +212,15 @@ namespace AMDiS
}); });
} }
// evaluate local basis-function at point
template <class Element>
ShapeValues const& evaluateFunction(Element const& element, DomainType const& local) const
{
LocalValuesKey key{element.type().id(),local};
return ShapeValuesCache::get(key, [&](LocalValuesKey const&)
{
ShapeValues data;
localBasis_->evaluateFunction(local, data);
return data;
});
}
// evaluate local basis-gradients at point
template <class Element>
ShapeGradients const& evaluateJacobian(Element const& element, DomainType const& local) const
{
LocalGradientsKey key{element.type().id(),local};
return ShapeGradientsCache::get(key, [&](LocalGradientsKey const&)
{
ShapeGradients data;
localBasis_->evaluateJacobian(local, data);
return data;
});
}
private: private:
using ShapeValuesCache = ConcurrentCache<LocalValuesKey, ShapeValues, ThreadLocalPolicy, template <class Value>
std::unordered_map<LocalValuesKey, ShapeValues, typename LocalValuesKey::hasher>>; using Cache = ConcurrentCache<Key, Value, ThreadLocalPolicy,
using ShapeGradientsCache = ConcurrentCache<LocalGradientsKey, ShapeGradients, ThreadLocalPolicy, std::unordered_map<Key, Value, typename Key::hasher>>;
std::unordered_map<LocalGradientsKey, ShapeGradients, typename LocalGradientsKey::hasher>>;
using ShapeValuesContainerCache = ConcurrentCache<ValuesKey, ShapeValuesContainer, ThreadLocalPolicy, Cache<ShapeValuesContainer> shapeValuesContainer_;
std::unordered_map<ValuesKey, ShapeValuesContainer, typename ValuesKey::hasher>>; Cache<ShapeGradientsContainer> shapeGradientsContainer_;
using ShapeGradientsContainerCache = ConcurrentCache<GradientsKey, ShapeGradientsContainer, ThreadLocalPolicy,
std::unordered_map<GradientsKey, ShapeGradientsContainer, typename GradientsKey::hasher>>;
LocalBasisType const* localBasis_; LocalBasisType const* localBasis_ = nullptr;
}; };
} // end namespace AMDiS } // end namespace AMDiS
#pragma once
#include <vector>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>
#include <amdis/LocalBasisCache.hpp>
#include <amdis/common/TupleUtility.hpp>
#include <amdis/utility/ConcurrentCache.hpp>
namespace AMDiS
{
namespace Impl
{
/// Pair of GeometryType and quadrature order and size
template <class LocalBasisType, class Tag>
struct EvaluatorCacheKey
{
using type = std::tuple<HashableGeometryType,typename LocalBasisType::Traits::DomainType>;
type value;
friend std::size_t hash_value(EvaluatorCacheKey const& t)
{
std::size_t seed = hash_value(std::get<0>(t.value));
Dune::hash_range(seed, std::get<1>(t.value).begin(), std::get<1>(t.value).end());
return seed;
}
friend bool operator==(EvaluatorCacheKey const& lhs, EvaluatorCacheKey const& rhs)
{
return lhs.value == rhs.value;
}
};
}
} // end namespace AMDiS
DUNE_DEFINE_HASH(DUNE_HASH_TEMPLATE_ARGS(class LocalBasisType, class Tag),DUNE_HASH_TYPE(AMDiS::Impl::EvaluatorCacheKey<LocalBasisType,Tag>))
namespace AMDiS
{
template <class LocalBasisType>
class LocalBasisEvaluatorCache
{
public:
using Traits = typename LocalBasisType::Traits;
using DomainType = typename Traits::DomainType;
using RangeType = typename Traits::RangeType;
using RangeFieldType = typename Traits::RangeFieldType;
using JacobianType = typename Traits::JacobianType;
using QuadratureRule = Dune::QuadratureRule<RangeFieldType,Traits::dimDomain>;
using QuadratureRules = Dune::QuadratureRules<RangeFieldType,Traits::dimDomain>;
private:
using ValuesKey = Impl::EvaluatorCacheKey<LocalBasisType, struct ValuesTag>;
using GradientsKey = Impl::EvaluatorCacheKey<LocalBasisType, struct GradientsTag>;
using Policy = tag::thread_local_policy;
public:
using ShapeValues = std::vector<RangeType>;
using ShapeGradients = std::vector<JacobianType>;
LocalBasisEvaluatorCache(LocalBasisType const& localBasis)
: localBasis_(&localBasis)
{}
template <class LocalContext>
ShapeValues const& evaluateFunction(LocalContext const& context, DomainType const& local)
{
ValuesKey key{typename ValuesKey::type{context.type(),local}};
return ShapeValuesCache::get(key, [&](ShapeValues* data, ValuesKey const&)
{
localBasis_->evaluateFunction(local, *data);
});
}
template <class LocalContext>
ShapeGradients const& evaluateJacobian(LocalContext const& context, DomainType const& local)
{
GradientsKey key{typename GradientsKey::type{context.type(),local}};
return ShapeGradientsCache::get(key, [&](ShapeGradients* data, GradientsKey const&)
{
localBasis_->evaluateJacobian(local, *data);
});
}
private:
using ShapeValuesCache = ConcurrentCache<ValuesKey, ShapeValues, Policy>;
using ShapeGradientsCache = ConcurrentCache<GradientsKey, ShapeGradients, Policy>;
LocalBasisType const* localBasis_;
};
} // end namespace AMDiS
...@@ -58,8 +58,7 @@ namespace AMDiS ...@@ -58,8 +58,7 @@ namespace AMDiS
static_assert(std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>{}, static_assert(std::is_same<FiniteElementType_t<RowNode>, FiniteElementType_t<ColNode>>{},
"Galerkin operator requires equal finite elements for test and trial space." ); "Galerkin operator requires equal finite elements for test and trial space." );
using LocalBasisType = typename FiniteElementType_t<RowNode>::Traits::LocalBasisType; using RangeFieldType = typename NodeQuadCache<RowNode>::Traits::RangeFieldType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType;
auto localFctA = localFunction(gridFctA_); localFctA.bind(context.element()); auto localFctA = localFunction(gridFctA_); localFctA.bind(context.element());
auto localFctB = localFunction(gridFctB_); localFctB.bind(context.element()); auto localFctB = localFunction(gridFctB_); localFctB.bind(context.element());
...@@ -75,9 +74,9 @@ namespace AMDiS ...@@ -75,9 +74,9 @@ namespace AMDiS
using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::mydimension>; using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::mydimension>;
auto const& quad = QuadratureRules::rule(context.type(), quadDeg); auto const& quad = QuadratureRules::rule(context.type(), quadDeg);
LocalBasisCache<LocalBasisType> cache(localFE.localBasis()); NodeQuadCache<RowNode> cache(localFE.localBasis());
auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad); auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
auto const& shapeValuesCache = cache.evaluateFunctionAtQp(context, quad); auto const& shapeValuesCache = cache.evaluateFunctionAtQP(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) { for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
...@@ -146,8 +145,6 @@ namespace AMDiS ...@@ -146,8 +145,6 @@ namespace AMDiS
static_assert(Node::isLeaf, static_assert(Node::isLeaf,
"Operator can be applied to Leaf-Nodes only." ); "Operator can be applied to Leaf-Nodes only." );
using LocalBasisType = typename FiniteElementType_t<Node>::Traits::LocalBasisType;
auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element()); auto localFctF = localFunction(gridFctF_); localFctF.bind(context.element());
auto const& localFE = node.finiteElement(); auto const& localFE = node.finiteElement();
...@@ -158,8 +155,8 @@ namespace AMDiS ...@@ -158,8 +155,8 @@ namespace AMDiS
using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::dimension>; using QuadratureRules = Dune::QuadratureRules<typename Context::Geometry::ctype, Context::LocalContext::dimension>;
auto const& quad = QuadratureRules::rule(context.type(), quad_order); auto const& quad = QuadratureRules::rule(context.type(), quad_order);
LocalBasisCache<LocalBasisType> cache(localFE.localBasis()); NodeQuadCache<Node> cache(localFE.localBasis());
auto const& shapeValuesCache = cache.evaluateFunctionAtQp(context, quad); auto const& shapeValuesCache = cache.evaluateFunctionAtQP(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) { for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
......
...@@ -50,12 +50,10 @@ namespace AMDiS ...@@ -50,12 +50,10 @@ namespace AMDiS
auto const& localFE = node.child(0).finiteElement(); auto const& localFE = node.child(0).finiteElement();
std::size_t feSize = localFE.size(); std::size_t feSize = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType; using RangeFieldType = typename NodeQuadCache<typename Node::ChildType>::Traits::RangeFieldType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType; NodeQuadCache<typename Node::ChildType> cache(localFE.localBasis());
LocalBasisCache<LocalBasisType> cache(localFE.localBasis()); auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) { for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(quad[iq].position()); decltype(auto) local = context.local(quad[iq].position());
......
...@@ -48,12 +48,10 @@ namespace AMDiS ...@@ -48,12 +48,10 @@ namespace AMDiS
auto const& localFE = node.finiteElement(); auto const& localFE = node.finiteElement();
std::size_t feSize = localFE.size(); std::size_t feSize = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType; using RangeFieldType = typename NodeQuadCache<Node>::Traits::RangeFieldType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType; NodeQuadCache<Node> cache(localFE.localBasis());
LocalBasisCache<LocalBasisType> cache(localFE.localBasis()); auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) { for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(quad[iq].position()); decltype(auto) local = context.local(quad[iq].position());
......
...@@ -49,12 +49,10 @@ namespace AMDiS ...@@ -49,12 +49,10 @@ namespace AMDiS
auto const& localFE = node.finiteElement(); auto const& localFE = node.finiteElement();
std::size_t feSize = localFE.size(); std::size_t feSize = localFE.size();
using LocalBasisType = typename std::decay_t<decltype(localFE)>::Traits::LocalBasisType; using RangeFieldType = typename NodeQuadCache<Node>::Traits::RangeFieldType;
using RangeFieldType = typename LocalBasisType::Traits::RangeFieldType; NodeQuadCache<Node> cache(localFE.localBasis());
LocalBasisCache<LocalBasisType> cache(localFE.localBasis()); auto const& shapeGradientsCache = cache.evaluateJacobianAtQP(context, quad);
auto const& shapeGradientsCache = cache.evaluateJacobianAtQp(context, quad);
for (std::size_t iq = 0; iq < quad.size(); ++iq) { for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element // Position of the current quadrature point in the reference element
decltype(auto) local = context.local(quad[iq].position()); decltype(auto) local = context.local(quad[iq].position());
......
...@@ -54,16 +54,13 @@ namespace AMDiS ...@@ -54,16 +54,13 @@ namespace AMDiS
std::size_t rowSize = rowLocalFE.size();