dune-localfunctions  2.2.1
raviartthomas12dlocalinterpolation.hh
Go to the documentation of this file.
1 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH
2 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH
3 
4 #include <vector>
5 
6 #include <dune/geometry/quadraturerules.hh>
7 
8 namespace Dune
9 {
10 
19  template<class LB>
21  {
22 
23 public:
26  {
27  sign0 = sign1 = sign2 = 1.0;
28  }
29 
35  RT12DLocalInterpolation (unsigned int s)
36  {
37  sign0 = sign1 = sign2 = 1.0;
38  if (s & 1)
39  {
40  sign0 = -1.0;
41  }
42  if (s & 2)
43  {
44  sign1 = -1.0;
45  }
46  if (s & 4)
47  {
48  sign2 = -1.0;
49  }
50  n0[0] = 0.0;
51  n0[1] = -1.0;
52  n1[0] = -1.0;
53  n1[1] = 0.0;
54  n2[0] = 1.0/sqrt(2.0);
55  n2[1] = 1.0/sqrt(2.0);
56  c0 = 0.5*n0[0] - 1.0*n0[1];
57  c1 = -1.0*n1[0] + 0.5*n1[1];
58  c2 = 0.5*n2[0] + 0.5*n2[1];
59  }
60 
69  template<typename F, typename C>
70  void interpolate (const F& f, std::vector<C>& out) const
71  {
72  // f gives v*outer normal at a point on the edge!
73  typedef typename LB::Traits::RangeFieldType Scalar;
74  typedef typename LB::Traits::DomainFieldType Vector;
75  typename F::Traits::RangeType y;
76 
77  out.resize(8);
78  fill(out.begin(), out.end(), 0.0);
79 
80  const int qOrder1 = 4;
81  const Dune::QuadratureRule<Scalar,1>& rule1 = Dune::QuadratureRules<Scalar,1>::rule(Dune::GeometryType(Dune::GeometryType::simplex,1), qOrder1);
82 
83  for (typename Dune::QuadratureRule<Scalar,1>::const_iterator it = rule1.begin();
84  it != rule1.end(); ++it)
85  {
86  Scalar qPos = it->position();
87  typename LB::Traits::DomainType localPos;
88 
89  localPos[0] = qPos;
90  localPos[1] = 0.0;
91  f.evaluate(localPos, y);
92  out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0/c0;
93  out[3] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight()/c0;
94 
95  localPos[0] = 0.0;
96  localPos[1] = qPos;
97  f.evaluate(localPos, y);
98  out[1] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1/c1;
99  out[4] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight()/c1;
100 
101  localPos[0] = 1.0 - qPos;
102  localPos[1] = qPos;
103  f.evaluate(localPos, y);
104  out[2] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2/c2;
105  out[5] += (y[0]*n2[0] + y[1]*n2[1])*(2.0*qPos - 1.0)*it->weight()/c2;
106  }
107 
108  const int qOrder2 = 8;
109  const Dune::QuadratureRule<Vector,2>& rule2 = Dune::QuadratureRules<Vector,2>::rule(Dune::GeometryType(Dune::GeometryType::simplex,2), qOrder2);
110 
111  for (typename Dune::QuadratureRule<Vector,2>::const_iterator it = rule2.begin();
112  it != rule2.end(); ++it)
113  {
114  Dune::FieldVector<double,2> qPos = it->position();
115 
116  f.evaluate(qPos, y);
117  out[6] += y[0]*it->weight();
118  out[7] += y[1]*it->weight();
119  }
120  }
121 
122 private:
123  typename LB::Traits::RangeFieldType sign0,sign1,sign2;
124  typename LB::Traits::DomainType n0,n1,n2;
125  typename LB::Traits::RangeFieldType c0,c1,c2;
126  };
127 }
128 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS12DLOCALINTERPOLATION_HH