Dune Core Modules (2.8.0)

localinterpolation.hh
1#ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALINTERPOLATION_HH
2#define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALINTERPOLATION_HH
3
4#include <algorithm>
5#include <array>
6#include <bitset>
7#include <vector>
8#include <limits>
9
11#include <dune/common/math.hh>
14
16
17#include <dune/localfunctions/common/localinterpolation.hh>
18
19
20namespace Dune
21{
22
34 template<class D, class R, unsigned int dim, unsigned int order>
36 {
38 using FaceDomainType = FieldVector<D, dim-1>;
40 using DomainFieldType = D;
41 using RangeFieldType = R;
42
43 static constexpr std::size_t numFaces = 2*dim;
44
45 static constexpr unsigned int interiorDofs = dim*binomial(dim+order-2, order-2);
46 static constexpr unsigned int faceDofs = binomial(dim+order-2, order-1);
47
54 inline static auto legendre( unsigned int i, const DomainFieldType& x )
55 -> RangeFieldType
56 {
57 switch( i )
58 {
59 case 0:
60 return 1;
61 case 1:
62 return 2*x-1;
63 case 2:
64 return 6*x*x-6*x+1;
65 case 3:
66 return 20*x*x*x-30*x*x+12*x-1;
67 default:
68 return ((2*i-1)*(2*x-1)*legendre(x,i-1) - (i-1)*legendre(x,i-2)) / i;
69 }
70 }
71
82 constexpr inline static auto unrank (std::size_t i)
83 -> std::array<std::size_t, d>
84 {
85 assert( i < binomial(d+kMax, kMax));
86
87 std::array<std::size_t, d> mi{};
88 if (i == 0)
89 return mi;
90
91 std::size_t k = 0, b = 1;
92 for(;k <= kMax && b <= i; ++k, b = binomial(d+k-1, k))
93 i -= b;
94
95 std::size_t p = 0;
96 for(; p < d && i > 0; ++p)
97 {
98 std::size_t m = 0, c = 1;
99 for(; m <= k && c <= i; ++m, c = binomial(d-p+m-2, m))
100 i -= c;
101
102 mi[p] = k-m;
103 k = m;
104 }
105 mi[p] = k;
106
107 return mi;
108 }
109
116 inline static auto embed (unsigned int face, const FaceDomainType& x )
117 -> DomainType
118 {
119 DomainType y;
120
121 std::size_t j = 0;
122 for (auto i : range(dim))
123 if (i == face/2)
124 y[i] = face%2;
125 else
126 y[i] = x[j++];
127
128 return y;
129 }
130
131 public:
134 {
135 std::fill(sign_.begin(), sign_.end(), 1.0);
136 }
137
143 BDFMCubeLocalInterpolation (std::bitset<numFaces> s)
144 {
145 for (auto i : range(numFaces))
146 {
147 sign_[i] = s[i] ? -1 : 1;
148 normal_[i][i/2] = i%2 ? 1 : -1;
149 }
150 }
151
161 template<class F, class C>
162 void interpolate (const F& ff, C& out) const
163 {
164 out.resize(numFaces*faceDofs + interiorDofs);
165 std::fill(out.begin(),out.end(), 0.0);
166
167 auto&& f = Impl::makeFunctionWithCallOperator<DomainType>(ff);
168
169 for(auto i : range(numFaces))
170 trace(i, f, out);
171
172 interior(f, out);
173 }
174
186 template<class F, class C>
187 void trace (unsigned int face, const F& f, C& out) const
188 {
189 assert(out.size() >= numFaces*faceDofs + interiorDofs);
190 assert(face < numFaces);
191
192 const auto o = face*faceDofs;
193 const auto& n = normal_[face];
194 const auto& s = sign_[face];
195 const auto& c = n[face/2];
196
197 const auto& rule = QuadratureRules<DomainFieldType, dim-1>::rule(GeometryTypes::cube(dim-1), order+(order-1)+5);
198 for (const auto& qp : rule)
199 {
200 const auto& x = qp.position();
201 auto y = f(embed(face, x));
202
203 for (auto i : range(faceDofs))
204 {
205 auto mi = unrank<dim-1,order-1>(i);
206
207 RangeFieldType phi = 1.;
208 for (auto j : range(dim-1))
209 phi *= legendre(mi[j], x[j]);
210
211 out[o+i] += (y*n) * phi * qp.weight() * (i%2 ? c : s);
212 }
213 }
214 }
215
225 template<class F, class C>
226 void interior (const F& f, C& out) const
227 {
228 assert(out.size() >= numFaces*faceDofs + interiorDofs);
229
230 const auto o = numFaces*faceDofs;
231
232 const auto& rule = QuadratureRules<DomainFieldType, dim>::rule(GeometryTypes::cube(dim), order+std::max((int)order-2,(int)0)+5);
233 for(const auto& qp : rule)
234 {
235 const auto& x = qp.position();
236 auto y = f(x);
237
238 for (auto i : range(interiorDofs/dim))
239 {
240 auto mi = unrank<dim,order-2>(i);
241
242 RangeFieldType phi = 1.;
243 for (auto j : range(dim))
244 phi *= legendre(mi[j], x[j]);
245
246 for (auto j : range(dim))
247 out[o+dim*i+j] += y[j]* phi * qp.weight();
248 }
249 }
250 }
251
252 private:
253 std::array<RangeFieldType, numFaces> sign_;
254 std::array<DomainType, numFaces> normal_;
255 };
256
257
258 template<class D, class R, unsigned int dim, unsigned int order>
259 constexpr std::size_t BDFMCubeLocalInterpolation<D, R, dim, order>::numFaces;
260
261 template<class D, class R, unsigned int dim, unsigned int order>
262 constexpr unsigned int BDFMCubeLocalInterpolation<D, R, dim, order>::interiorDofs;
263
264 template<class D, class R, unsigned int dim, unsigned int order>
265 constexpr unsigned int BDFMCubeLocalInterpolation<D, R, dim, order>::faceDofs;
266
267
268#ifndef DOXYGEN
269 template<class D, class R, unsigned int dim>
270 class BDFMCubeLocalInterpolation<D, R, dim, 0>
271 {
272 static_assert(AlwaysFalse<D>::value,
273 "`BDFMCubeLocalCoefficients` not defined for order 0.");
274 };
275#endif //#ifndef DOXYGEN
276
277} // namespace Dune
278
279#endif // #ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALINTERPOLATION_HH
Interpolation for Brezzi-Douglas-Fortin-Marini shape functions on cubes.
Definition: localinterpolation.hh:36
void interpolate(const F &ff, C &out) const
Interpolate a given function with shape functions.
Definition: localinterpolation.hh:162
BDFMCubeLocalInterpolation(std::bitset< numFaces > s)
Make set number s, where 0 <= s < 2^(2*dim)
Definition: localinterpolation.hh:143
BDFMCubeLocalInterpolation()
Standard constructor.
Definition: localinterpolation.hh:133
void trace(unsigned int face, const F &f, C &out) const
Interpolate a given function with shape functions on a face of the reference element.
Definition: localinterpolation.hh:187
void interior(const F &f, C &out) const
Interpolate a given function with shape functions in the interior of the reference element.
Definition: localinterpolation.hh:226
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:280
Traits for type conversions and type information.
Implements a vector constructed from a given type representing a field and a compile-time given size.
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:470
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:11
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:128
Utilities for reduction like operations on ranges.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)