3#ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALINTERPOLATION_HH
4#define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALINTERPOLATION_HH
12#include <dune/common/fvector.hh>
13#include <dune/common/math.hh>
14#include <dune/common/rangeutilities.hh>
15#include <dune/common/typetraits.hh>
17#include <dune/geometry/quadraturerules.hh>
34 template<
class D,
class R,
unsigned int dim,
unsigned int order>
37 using DomainType = FieldVector<D, dim>;
38 using FaceDomainType = FieldVector<D, dim-1>;
39 using RangeType = FieldVector<R, dim>;
40 using DomainFieldType = D;
41 using RangeFieldType = R;
43 static constexpr std::size_t numFaces = 2*dim;
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);
54 inline static auto legendre(
unsigned int i,
const DomainFieldType& x )
66 return 20*x*x*x-30*x*x+12*x-1;
68 return ((2*i-1)*(2*x-1)*legendre(x,i-1) - (i-1)*legendre(x,i-2)) / i;
81 template<std::size_t d, std::size_t kMax = std::numeric_limits<std::size_t>::max()>
82 constexpr inline static auto unrank (std::size_t i)
83 -> std::array<std::size_t, d>
85 assert( i < binomial(d+kMax, kMax));
87 std::array<std::size_t, d> mi{};
91 std::size_t k = 0, b = 1;
92 for(;k <= kMax && b <= i; ++k, b = binomial(d+k-1, k))
96 for(; p < d && i > 0; ++p)
98 std::size_t m = 0, c = 1;
99 for(; m <= k && c <= i; ++m, c = binomial(d-p+m-2, m))
116 inline static auto embed (
unsigned int face,
const FaceDomainType& x )
122 for (
auto i : range(dim))
135 std::fill(sign_.begin(), sign_.end(), 1.0);
145 for (
auto i : range(numFaces))
147 sign_[i] = s[i] ? -1 : 1;
148 normal_[i][i/2] = i%2 ? 1 : -1;
161 template<
class F,
class C>
164 out.resize(numFaces*faceDofs + interiorDofs);
165 std::fill(out.begin(),out.end(), 0.0);
167 for(
auto i : range(numFaces))
184 template<
class F,
class C>
185 void trace (
unsigned int face,
const F& f, C& out)
const
187 assert(out.size() >= numFaces*faceDofs + interiorDofs);
188 assert(face < numFaces);
190 const auto o = face*faceDofs;
191 const auto& n = normal_[face];
192 const auto& s = sign_[face];
193 const auto& c = n[face/2];
195 const auto& rule = QuadratureRules<DomainFieldType, dim-1>::rule(GeometryTypes::cube(dim-1), order+(order-1)+5);
196 for (
const auto& qp : rule)
198 const auto& x = qp.position();
199 auto y = f(embed(face, x));
201 for (
auto i : range(faceDofs))
203 auto mi = unrank<dim-1,order-1>(i);
205 RangeFieldType phi = 1.;
206 for (
auto j : range(dim-1))
207 phi *= legendre(mi[j], x[j]);
209 out[o+i] += (y*n) * phi * qp.weight() * (i%2 ? c : s);
223 template<
class F,
class C>
226 assert(out.size() >= numFaces*faceDofs + interiorDofs);
228 const auto o = numFaces*faceDofs;
230 const auto& rule = QuadratureRules<DomainFieldType, dim>::rule(GeometryTypes::cube(dim), order+std::max((
int)order-2,(
int)0)+5);
231 for(
const auto& qp : rule)
233 const auto& x = qp.position();
236 for (
auto i : range(interiorDofs/dim))
238 auto mi = unrank<dim,order-2>(i);
240 RangeFieldType phi = 1.;
241 for (
auto j : range(dim))
242 phi *= legendre(mi[j], x[j]);
244 for (
auto j : range(dim))
245 out[o+dim*i+j] += y[j]* phi * qp.weight();
251 std::array<RangeFieldType, numFaces> sign_;
252 std::array<DomainType, numFaces> normal_;
256 template<
class D,
class R,
unsigned int dim,
unsigned int order>
257 constexpr std::size_t BDFMCubeLocalInterpolation<D, R, dim, order>::numFaces;
259 template<
class D,
class R,
unsigned int dim,
unsigned int order>
260 constexpr unsigned int BDFMCubeLocalInterpolation<D, R, dim, order>::interiorDofs;
262 template<
class D,
class R,
unsigned int dim,
unsigned int order>
263 constexpr unsigned int BDFMCubeLocalInterpolation<D, R, dim, order>::faceDofs;
267 template<
class D,
class R,
unsigned int dim>
268 class BDFMCubeLocalInterpolation<D, R, dim, 0>
270 static_assert(AlwaysFalse<D>::value,
271 "`BDFMCubeLocalCoefficients` not defined for order 0.");
Definition bdfmcube.hh:18
Interpolation for Brezzi-Douglas-Fortin-Marini shape functions on cubes.
Definition brezzidouglasfortinmarini/cube/localinterpolation.hh:36
BDFMCubeLocalInterpolation(std::bitset< numFaces > s)
Make set number s, where 0 <= s < 2^(2*dim)
Definition brezzidouglasfortinmarini/cube/localinterpolation.hh:143
BDFMCubeLocalInterpolation()
Standard constructor.
Definition brezzidouglasfortinmarini/cube/localinterpolation.hh:133
void interpolate(const F &f, C &out) const
Interpolate a given function with shape functions.
Definition brezzidouglasfortinmarini/cube/localinterpolation.hh:162
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 brezzidouglasfortinmarini/cube/localinterpolation.hh:185
void interior(const F &f, C &out) const
Interpolate a given function with shape functions in the interior of the reference element.
Definition brezzidouglasfortinmarini/cube/localinterpolation.hh:224