5#ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH
6#define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH
10#include <dune/geometry/quadraturerules.hh>
38 for (
size_t i=0; i<3; i++)
39 sign_[i] = (s[i]) ? -1.0 : 1.0;
43 n_[2] = { 1.0/sqrt(2.0), 1.0/sqrt(2.0)};
45 c_ = { 0.5*n_[0][0] - 1.0*n_[0][1],
46 -1.0*n_[1][0] + 0.5*n_[1][1],
47 0.5*n_[2][0] + 0.5*n_[2][1]};
58 template<
typename F,
typename C>
62 typedef typename LB::Traits::RangeFieldType Scalar;
63 typedef typename LB::Traits::DomainFieldType Vector;
65 auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
68 fill(out.begin(), out.end(), 0.0);
70 const int qOrder1 = 4;
71 const auto& rule1 = Dune::QuadratureRules<Scalar,1>::rule(Dune::GeometryTypes::simplex(1), qOrder1);
73 for (
auto&& qp : rule1)
75 Scalar qPos = qp.position();
76 typename LB::Traits::DomainType localPos;
78 localPos = {qPos, 0.0};
80 out[0] += (y[0]*n_[0][0] + y[1]*n_[0][1])*qp.weight()*sign_[0]/c_[0];
81 out[3] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(2.0*qPos - 1.0)*qp.weight()/c_[0];
83 localPos = {0.0, qPos};
85 out[1] += (y[0]*n_[1][0] + y[1]*n_[1][1])*qp.weight()*sign_[1]/c_[1];
86 out[4] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(1.0 - 2.0*qPos)*qp.weight()/c_[1];
88 localPos = {1.0 - qPos, qPos};
90 out[2] += (y[0]*n_[2][0] + y[1]*n_[2][1])*qp.weight()*sign_[2]/c_[2];
91 out[5] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(2.0*qPos - 1.0)*qp.weight()/c_[2];
94 const int qOrder2 = 8;
95 const auto& rule2 = Dune::QuadratureRules<Vector,2>::rule(Dune::GeometryTypes::simplex(2), qOrder2);
97 for (
auto&& qp : rule2)
99 auto qPos = qp.position();
102 out[6] += y[0]*qp.weight();
103 out[7] += y[1]*qp.weight();
109 std::array<typename LB::Traits::RangeFieldType, 3> sign_;
112 std::array<typename LB::Traits::DomainType, 3> n_;
114 std::array<typename LB::Traits::RangeFieldType, 3> c_;
Definition: bdfmcube.hh:18
First order Raviart-Thomas shape functions on the reference quadrilateral.
Definition: raviartthomas12dlocalinterpolation.hh:26
RT12DLocalInterpolation(std::bitset< 3 > s=0)
Make set number s, where 0 <= s < 8.
Definition: raviartthomas12dlocalinterpolation.hh:35
void interpolate(const F &ff, std::vector< C > &out) const
Interpolate a given function with shape functions.
Definition: raviartthomas12dlocalinterpolation.hh:59