LocalBasisCache.hpp 5.45 KB
Newer Older
1
2
3
4
#pragma once

#include <vector>

5
#include <dune/common/hash.hh>
6
7
8
9
10
11
12
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>

#include <amdis/utility/ConcurrentCache.hpp>

namespace AMDiS
{
13
14
  template <class LocalBasisType>
  class LocalBasisCache
15
  {
16
17
  public:
    using Traits = typename LocalBasisType::Traits;
18

19
20
21
22
    using DomainType = typename Traits::DomainType;
    using RangeType = typename Traits::RangeType;
    using RangeFieldType = typename Traits::RangeFieldType;
    using JacobianType = typename Traits::JacobianType;
23

24
25
    using QuadratureRule = Dune::QuadratureRule<RangeFieldType,Traits::dimDomain>;
    using QuadratureRules = Dune::QuadratureRules<RangeFieldType,Traits::dimDomain>;
26

27
  private:
28
29

    /// Pair of GeometryType and quadrature order and size
30
31
    template <class Tag>
    struct QuadKey
32
    {
33
34
35
      unsigned int id;    // topologyId
      int p;              // quadrature order
      std::size_t size;   // number of quadrature points
36

37
      struct hasher
38
      {
39
40
41
42
43
44
45
46
47
48
49
50
51
        std::size_t operator()(QuadKey const& t) const
        {
          std::size_t seed = 0;
          Dune::hash_combine(seed, t.id);
          Dune::hash_combine(seed, t.p);
          Dune::hash_combine(seed, t.size);
          return seed;
        }
      };

      friend bool operator==(QuadKey const& lhs, QuadKey const& rhs)
      {
        return std::tie(lhs.id,lhs.p,lhs.size) == std::tie(rhs.id,rhs.p,rhs.size);
52
      }
53
54
55
56
57
58
59
60
    };

    /// Pair of GeometryType and local coordinates
    template <class Tag>
    struct CoordsKey
    {
      unsigned int id;    // topologyId
      DomainType local;   // local coordinate
61

62
63
64
65
66
67
68
69
70
71
72
73
      struct hasher
      {
        std::size_t operator()(CoordsKey const& t) const
        {
          std::size_t seed = 0;
          Dune::hash_combine(seed, t.id);
          Dune::hash_range(seed, t.local.begin(), t.local.end());
          return seed;
        }
      };

      friend bool operator==(CoordsKey const& lhs, CoordsKey const& rhs)
74
      {
75
        return std::tie(lhs.id,lhs.local) == std::tie(rhs.id,rhs.local);
76
77
78
      }
    };

79
80
    using LocalValuesKey = CoordsKey<struct ValuesTag>;
    using LocalGradientsKey = CoordsKey<struct GradientsTag>;
81

82
83
    using ValuesKey = QuadKey<struct ValuesTag>;
    using GradientsKey = QuadKey<struct GradientsTag>;
84

85
86
87
  public:
    using ShapeValues = std::vector<RangeType>;
    using ShapeGradients = std::vector<JacobianType>;
88

89
90
    using ShapeValuesContainer = std::vector<ShapeValues>;
    using ShapeGradientsContainer = std::vector<ShapeGradients>;
91
92
93
94
95
96


    LocalBasisCache(LocalBasisType const& localBasis)
      : localBasis_(&localBasis)
    {}

97
98
99
    // evaluate local basis-function at all quadrature points
    template <class LocalContext>
    ShapeValuesContainer const& evaluateFunctionAtQp(LocalContext const& context, QuadratureRule const& quad) const
100
    {
101
102
      ValuesKey key{context.type().id(),quad.order(),quad.size()};
      return ShapeValuesContainerCache::get(key, [&](ValuesKey const&)
103
      {
104
        ShapeValuesContainer data(quad.size());
105
        for (std::size_t iq = 0; iq < quad.size(); ++iq)
106
107
          localBasis_->evaluateFunction(context.local(quad[iq].position()), data[iq]);
        return data;
108
109
110
      });
    }

111
112
113
    // evaluate local basis-gradients at all quadrature points
    template <class LocalContext>
    ShapeGradientsContainer const& evaluateJacobianAtQp(LocalContext const& context, QuadratureRule const& quad) const
114
    {
115
116
      GradientsKey key{context.type().id(),quad.order(),quad.size()};
      return ShapeGradientsContainerCache::get(key, [&](GradientsKey const&)
117
      {
118
        ShapeGradientsContainer data(quad.size());
119
        for (std::size_t iq = 0; iq < quad.size(); ++iq)
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
          localBasis_->evaluateJacobian(context.local(quad[iq].position()), data[iq]);
        return data;
      });
    }


    // 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;
149
150
151
152
      });
    }

  private:
153
154
155
156
157
158
159
160
161
    using ShapeValuesCache = ConcurrentCache<LocalValuesKey, ShapeValues, ThreadLocalPolicy,
      std::unordered_map<LocalValuesKey, ShapeValues, typename LocalValuesKey::hasher>>;
    using ShapeGradientsCache = ConcurrentCache<LocalGradientsKey, ShapeGradients, ThreadLocalPolicy,
      std::unordered_map<LocalGradientsKey, ShapeGradients, typename LocalGradientsKey::hasher>>;

    using ShapeValuesContainerCache = ConcurrentCache<ValuesKey, ShapeValuesContainer, ThreadLocalPolicy,
      std::unordered_map<ValuesKey, ShapeValuesContainer, typename ValuesKey::hasher>>;
    using ShapeGradientsContainerCache = ConcurrentCache<GradientsKey, ShapeGradientsContainer, ThreadLocalPolicy,
      std::unordered_map<GradientsKey, ShapeGradientsContainer, typename GradientsKey::hasher>>;
162
163
164
165
166

    LocalBasisType const* localBasis_;
  };

} // end namespace AMDiS