dune-localfunctions  2.8.0
raviartthomas12dlocalinterpolation.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH
4 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH
5 
6 #include <vector>
7 
8 #include <dune/geometry/quadraturerules.hh>
10 
11 namespace Dune
12 {
13 
22  template<class LB>
24  {
25 
26  public:
27 
33  RT12DLocalInterpolation (std::bitset<3> s = 0)
34  {
35  using std::sqrt;
36  for (size_t i=0; i<3; i++)
37  sign_[i] = (s[i]) ? -1.0 : 1.0;
38 
39  n_[0] = { 0.0, -1.0};
40  n_[1] = {-1.0, 0.0};
41  n_[2] = { 1.0/sqrt(2.0), 1.0/sqrt(2.0)};
42 
43  c_ = { 0.5*n_[0][0] - 1.0*n_[0][1],
44  -1.0*n_[1][0] + 0.5*n_[1][1],
45  0.5*n_[2][0] + 0.5*n_[2][1]};
46  }
47 
56  template<typename F, typename C>
57  void interpolate (const F& ff, std::vector<C>& out) const
58  {
59  // f gives v*outer normal at a point on the edge!
60  typedef typename LB::Traits::RangeFieldType Scalar;
61  typedef typename LB::Traits::DomainFieldType Vector;
62 
63  auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
64 
65  out.resize(8);
66  fill(out.begin(), out.end(), 0.0);
67 
68  const int qOrder1 = 4;
69  const auto& rule1 = Dune::QuadratureRules<Scalar,1>::rule(Dune::GeometryTypes::simplex(1), qOrder1);
70 
71  for (auto&& qp : rule1)
72  {
73  Scalar qPos = qp.position();
74  typename LB::Traits::DomainType localPos;
75 
76  localPos = {qPos, 0.0};
77  auto y = f(localPos);
78  out[0] += (y[0]*n_[0][0] + y[1]*n_[0][1])*qp.weight()*sign_[0]/c_[0];
79  out[3] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(2.0*qPos - 1.0)*qp.weight()/c_[0];
80 
81  localPos = {0.0, qPos};
82  y = f(localPos);
83  out[1] += (y[0]*n_[1][0] + y[1]*n_[1][1])*qp.weight()*sign_[1]/c_[1];
84  out[4] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(1.0 - 2.0*qPos)*qp.weight()/c_[1];
85 
86  localPos = {1.0 - qPos, qPos};
87  y = f(localPos);
88  out[2] += (y[0]*n_[2][0] + y[1]*n_[2][1])*qp.weight()*sign_[2]/c_[2];
89  out[5] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(2.0*qPos - 1.0)*qp.weight()/c_[2];
90  }
91 
92  const int qOrder2 = 8;
93  const auto& rule2 = Dune::QuadratureRules<Vector,2>::rule(Dune::GeometryTypes::simplex(2), qOrder2);
94 
95  for (auto&& qp : rule2)
96  {
97  auto qPos = qp.position();
98 
99  auto y = f(qPos);
100  out[6] += y[0]*qp.weight();
101  out[7] += y[1]*qp.weight();
102  }
103  }
104 
105  private:
106  // Edge orientations
107  std::array<typename LB::Traits::RangeFieldType, 3> sign_;
108 
109  // Edge normals
110  std::array<typename LB::Traits::DomainType, 3> n_;
111 
112  std::array<typename LB::Traits::RangeFieldType, 3> c_;
113  };
114 }
115 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH
Definition: bdfmcube.hh:16
First order Raviart-Thomas shape functions on the reference quadrilateral.
Definition: raviartthomas12dlocalinterpolation.hh:24
RT12DLocalInterpolation(std::bitset< 3 > s=0)
Make set number s, where 0 <= s < 8.
Definition: raviartthomas12dlocalinterpolation.hh:33
void interpolate(const F &ff, std::vector< C > &out) const
Interpolate a given function with shape functions.
Definition: raviartthomas12dlocalinterpolation.hh:57