1#ifndef DUNE_COPASI_LOCAL_FINITE_ELEMENT_CACHE_HH
2#define DUNE_COPASI_LOCAL_FINITE_ELEMENT_CACHE_HH
4#include <dune-copasi-config.hh>
8#include <dune/geometry/quadraturerules.hh>
9#include <dune/geometry/referenceelements.hh>
11#include <dune/common/exceptions.hh>
12#include <dune/common/hash.hh>
26template<
class LocalBasisTraits>
29 using Domain =
typename LocalBasisTraits::DomainType;
30 using DomainField =
typename LocalBasisTraits::DomainFieldType;
31 static constexpr int dimDomain = LocalBasisTraits::dimDomain;
33 using Range =
typename LocalBasisTraits::RangeType;
34 using RangeField =
typename LocalBasisTraits::RangeFieldType;
35 using Jacobian =
typename LocalBasisTraits::JacobianType;
41 Key(std::type_index fem_id =
typeid(
void),
42 void const* quad_id =
nullptr,
63 std::size_t hash = 233;
93 const Dune::QuadratureRule<DomainField, dim>& quad_rule,
94 auto&& quad_proj)
const noexcept
96 const auto& fem_id =
typeid(finite_element.localBasis());
97 const auto quad_id = &quad_rule;
99 quad_proj(referenceElement<DomainField, dim>(quad_rule.type()).position(0, 0));
100 return { fem_id, quad_id, proj_id };
106 void bind(
const auto& finite_element,
107 const Dune::QuadratureRule<DomainField, dimDomain>& quad_rule,
110 bind(finite_element, quad_rule, std::identity{}, force);
130 template<
int dim,
typename ProjQuad>
131 requires std::is_invocable_r_v<Domain, ProjQuad, FieldVector<DomainField, dim>>
132 void bind(
const auto& finite_element,
133 const Dune::QuadratureRule<DomainField, dim>& quad_rule,
134 ProjQuad&& quad_proj,
140 auto new_key =
key(finite_element, quad_rule, std::forward<ProjQuad>(quad_proj));
141 if (std::exchange(_key, new_key) == new_key) {
146 _fe_size = finite_element.size();
147 _quad_size = quad_rule.size();
149 auto rit = (_key ==
Key{}) ? _range.end() : _range.find(_key);
150 if (rit != _range.end()) {
151 _bound_range = rit->second.data();
153 _bound_range = cache_evaluate(finite_element.localBasis(), quad_rule, quad_proj, _key);
156 auto jit = (_key ==
Key{}) ? _jacobian.end() : _jacobian.find(_key);
157 if (jit != _jacobian.end()) {
158 _bound_jacobian = jit->second.data();
160 _bound_jacobian = cache_jacobian(finite_element.localBasis(), quad_rule, quad_proj, _key);
169 assert(std::exchange(_key,
Key{}) !=
Key{});
170 assert(std::exchange(_bound_range,
nullptr));
171 assert(std::exchange(_bound_jacobian,
nullptr));
183 assert(_bound_range);
184 assert(position_id < _quad_size);
185 return { _bound_range + _fe_size * position_id, _fe_size };
197 assert(_bound_jacobian);
198 assert(position_id < _quad_size);
199 return { _bound_jacobian + _fe_size * position_id, _fe_size };
203 Range* cache_evaluate(
const auto& basis,
const auto& quad_rule,
const auto& projection, Key
key)
205 auto [it, inserted] = _range.emplace(
key, std::vector<Range>{});
206 assert(inserted or
key == Key{});
207 auto& range = it->second;
209 for (
const auto& quad : quad_rule) {
210 basis.evaluateFunction(projection(quad.position()), _quad_range);
211 range.insert(end(range), begin(_quad_range), end(_quad_range));
213 assert(range.size() == _fe_size * _quad_size);
217 Jacobian* cache_jacobian(
const auto& basis,
218 const auto& quad_rule,
219 const auto& projection,
222 auto [it, inserted] = _jacobian.emplace(
key, std::vector<Jacobian>{});
223 assert(inserted or
key == Key{});
224 auto& jacobian = it->second;
226 for (
const auto& quad : quad_rule) {
227 basis.evaluateJacobian(projection(quad.position()), _quad_jacobian);
228 jacobian.insert(end(jacobian), begin(_quad_jacobian), end(_quad_jacobian));
230 assert(jacobian.size() == _fe_size * _quad_size);
231 return jacobian.data();
235 std::size_t _fe_size = 0, _quad_size = 0;
236 Range
const* _bound_range =
nullptr;
237 Jacobian
const* _bound_jacobian =
nullptr;
238 std::unordered_map<Key, std::vector<Range>,
typename Key::Hash> _range;
239 std::unordered_map<Key, std::vector<Jacobian>,
typename Key::Hash> _jacobian;
240 std::vector<Range> _quad_range;
241 std::vector<Jacobian> _quad_jacobian;
This class describes a local basis cache.
Definition: local_basis_cache.hh:28
void bind(const auto &finite_element, const Dune::QuadratureRule< DomainField, dimDomain > &quad_rule, bool force=false)
Binds a finite element to this cache.
Definition: local_basis_cache.hh:106
std::span< const Jacobian > evaluateJacobian(std::size_t position_id) const
Evaluate jacobian with the bound finite element and quadrature rule.
Definition: local_basis_cache.hh:195
LocalBasisCache()
Constructs a new instance with empty values.
Definition: local_basis_cache.hh:75
void unbind() noexcept
Unbind finite element from this cache.
Definition: local_basis_cache.hh:167
void bind(const auto &finite_element, const Dune::QuadratureRule< DomainField, dim > &quad_rule, ProjQuad &&quad_proj, bool force=false)
Binds a finite element to this cache.
Definition: local_basis_cache.hh:132
Key key(const auto &finite_element, const Dune::QuadratureRule< DomainField, dim > &quad_rule, auto &&quad_proj) const noexcept
Key to match previews applications of this cache.
Definition: local_basis_cache.hh:92
std::span< const Range > evaluateFunction(std::size_t position_id) const
Evaluate function with the bound finite element and quadrature rule.
Definition: local_basis_cache.hh:181
Definition: axis_names.hh:7
Definition: local_basis_cache.hh:57
std::size_t operator()(const Key &key) const
Definition: local_basis_cache.hh:61
std::size_t result_type
Definition: local_basis_cache.hh:59
Key for the local basis cache.
Definition: local_basis_cache.hh:40
Domain _proj_id
Definition: local_basis_cache.hh:52
void const * _quad_id
Definition: local_basis_cache.hh:51
std::type_index _fem_id
Definition: local_basis_cache.hh:50
bool operator==(const Key &) const =default
Key(std::type_index fem_id=typeid(void), void const *quad_id=nullptr, Domain proj_id={})
Definition: local_basis_cache.hh:41