dune-localfunctions  2.2.1
raviartthomas2q2dlocalinterpolation.hh
Go to the documentation of this file.
1 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS2Q2DLOCALINTERPOLATION_HH
2 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS2Q2DLOCALINTERPOLATION_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 = sign3 = 1.0;
28  }
29 
35  RT2Q2DLocalInterpolation (unsigned int s)
36  {
37  sign0 = sign1 = sign2 = sign3 = 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  if (s & 8)
51  {
52  sign3 *= -1.0;
53  }
54 
55  n0[0] = -1.0;
56  n0[1] = 0.0;
57  n1[0] = 1.0;
58  n1[1] = 0.0;
59  n2[0] = 0.0;
60  n2[1] = -1.0;
61  n3[0] = 0.0;
62  n3[1] = 1.0;
63  }
64 
73  template<typename F, typename C>
74  void interpolate (const F& f, std::vector<C>& out) const
75  {
76  // f gives v*outer normal at a point on the edge!
77  typedef typename LB::Traits::RangeFieldType Scalar;
78  typedef typename LB::Traits::DomainFieldType Vector;
79  typename F::Traits::RangeType y;
80 
81  out.resize(24);
82  fill(out.begin(), out.end(), 0.0);
83 
84  const int qOrder = 6;
85  const QuadratureRule<Scalar,1>& rule = QuadratureRules<Scalar,1>::rule(GeometryType(GeometryType::cube,1), qOrder);
86 
87  for (typename QuadratureRule<Scalar,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
88  {
89  Scalar qPos = it->position();
90  typename LB::Traits::DomainType localPos;
91 
92  localPos[0] = 0.0;
93  localPos[1] = qPos;
94  f.evaluate(localPos, y);
95  out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0;
96  out[1] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight();
97  out[2] += (y[0]*n0[0] + y[1]*n0[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign0;
98 
99  localPos[0] = 1.0;
100  localPos[1] = qPos;
101  f.evaluate(localPos, y);
102  out[3] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1;
103  out[4] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight();
104  out[5] += (y[0]*n1[0] + y[1]*n1[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign1;
105 
106  localPos[0] = qPos;
107  localPos[1] = 0.0;
108  f.evaluate(localPos, y);
109  out[6] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2;
110  out[7] += (y[0]*n2[0] + y[1]*n2[1])*(1.0 - 2.0*qPos)*it->weight();
111  out[8] += (y[0]*n2[0] + y[1]*n2[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign2;
112 
113  localPos[0] = qPos;
114  localPos[1] = 1.0;
115  f.evaluate(localPos, y);
116  out[9] += (y[0]*n3[0] + y[1]*n3[1])*it->weight()*sign3;
117  out[10] += (y[0]*n3[0] + y[1]*n3[1])*(2.0*qPos - 1.0)*it->weight();
118  out[11] += (y[0]*n3[0] + y[1]*n3[1])*(6.0*qPos*qPos - 6.0*qPos + 1.0)*it->weight()*sign3;
119  }
120 
121  const QuadratureRule<Vector,2>& rule2 = QuadratureRules<Vector,2>::rule(GeometryType(GeometryType::cube,2), qOrder);
122 
123  for (typename QuadratureRule<Vector,2>::const_iterator it = rule2.begin();
124  it != rule2.end(); ++it)
125  {
126  FieldVector<double,2> qPos = it->position();
127 
128  f.evaluate(qPos, y);
129  out[12] += y[0]*it->weight();
130  out[13] += y[1]*it->weight();
131  out[14] += y[0]*qPos[0]*it->weight();
132  out[15] += y[1]*qPos[0]*it->weight();
133  out[16] += y[0]*qPos[1]*it->weight();
134  out[17] += y[1]*qPos[1]*it->weight();
135  out[18] += y[0]*qPos[0]*qPos[1]*it->weight();
136  out[19] += y[1]*qPos[0]*qPos[1]*it->weight();
137  out[20] += y[0]*qPos[1]*qPos[1]*it->weight();
138  out[21] += y[1]*qPos[0]*qPos[0]*it->weight();
139  out[22] += y[0]*qPos[0]*qPos[1]*qPos[1]*it->weight();
140  out[23] += y[1]*qPos[0]*qPos[0]*qPos[1]*it->weight();
141  }
142  }
143 
144 private:
145  typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3;
146  typename LB::Traits::DomainType n0, n1, n2, n3;
147  };
148 }
149 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS2Q2DLOCALINTERPOLATION_HH