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