1#ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALINTERPOLATION_HH
2#define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALINTERPOLATION_HH
17#include <dune/localfunctions/common/localinterpolation.hh>
34 template<
class D,
class R,
unsigned int dim,
unsigned int order>
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;
82 constexpr inline static auto unrank (std::size_t i)
83 -> std::array<std::size_t, d>
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 auto&& f = Impl::makeFunctionWithCallOperator<DomainType>(ff);
169 for(
auto i : range(numFaces))
186 template<
class F,
class C>
187 void trace (
unsigned int face,
const F& f, C& out)
const
189 assert(out.size() >= numFaces*faceDofs + interiorDofs);
190 assert(face < numFaces);
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];
198 for (
const auto& qp : rule)
200 const auto& x = qp.position();
201 auto y = f(embed(face, x));
203 for (
auto i : range(faceDofs))
205 auto mi = unrank<dim-1,order-1>(i);
207 RangeFieldType phi = 1.;
208 for (
auto j : range(dim-1))
209 phi *= legendre(mi[j], x[j]);
211 out[o+i] += (y*n) * phi * qp.weight() * (i%2 ? c : s);
225 template<
class F,
class C>
228 assert(out.size() >= numFaces*faceDofs + interiorDofs);
230 const auto o = numFaces*faceDofs;
233 for(
const auto& qp : rule)
235 const auto& x = qp.position();
238 for (
auto i : range(interiorDofs/dim))
240 auto mi = unrank<dim,order-2>(i);
242 RangeFieldType phi = 1.;
243 for (
auto j : range(dim))
244 phi *= legendre(mi[j], x[j]);
246 for (
auto j : range(dim))
247 out[o+dim*i+j] += y[j]* phi * qp.weight();
253 std::array<RangeFieldType, numFaces> sign_;
254 std::array<DomainType, numFaces> normal_;
258 template<
class D,
class R,
unsigned int dim,
unsigned int order>
259 constexpr std::size_t BDFMCubeLocalInterpolation<D, R, dim, order>::numFaces;
261 template<
class D,
class R,
unsigned int dim,
unsigned int order>
262 constexpr unsigned int BDFMCubeLocalInterpolation<D, R, dim, order>::interiorDofs;
264 template<
class D,
class R,
unsigned int dim,
unsigned int order>
265 constexpr unsigned int BDFMCubeLocalInterpolation<D, R, dim, order>::faceDofs;
269 template<
class D,
class R,
unsigned int dim>
270 class BDFMCubeLocalInterpolation<D, R, dim, 0>
272 static_assert(AlwaysFalse<D>::value,
273 "`BDFMCubeLocalCoefficients` not defined for order 0.");
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.