dune-localfunctions  2.2.1
raviartthomas1q3dlocalinterpolation.hh
Go to the documentation of this file.
1 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH
2 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH
3 
4 #include <vector>
5 
6 #include <dune/geometry/quadraturerules.hh>
7 
8 namespace Dune
9 {
18  template<class LB>
20  {
21 
22 public:
25  {
26  sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
27  }
28 
34  RT1Q3DLocalInterpolation (unsigned int s)
35  {
36  sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
37  if (s & 1)
38  {
39  sign0 = -1.0;
40  }
41  if (s & 2)
42  {
43  sign1 = -1.0;
44  }
45  if (s & 4)
46  {
47  sign2 = -1.0;
48  }
49  if (s & 8)
50  {
51  sign3 = -1.0;
52  }
53  if (s & 16)
54  {
55  sign4 = -1.0;
56  }
57  if (s & 32)
58  {
59  sign5 = -1.0;
60  }
61 
62  n0[0] = -1.0;
63  n0[1] = 0.0;
64  n0[2] = 0.0;
65  n1[0] = 1.0;
66  n1[1] = 0.0;
67  n1[2] = 0.0;
68  n2[0] = 0.0;
69  n2[1] = -1.0;
70  n2[2] = 0.0;
71  n3[0] = 0.0;
72  n3[1] = 1.0;
73  n3[2] = 0.0;
74  n4[0] = 0.0;
75  n4[1] = 0.0;
76  n4[2] = -1.0;
77  n5[0] = 0.0;
78  n5[1] = 0.0;
79  n5[2] = 1.0;
80  }
81 
90  template<class F, class C>
91  void interpolate (const F& f, std::vector<C>& out) const
92  {
93  // f gives v*outer normal at a point on the edge!
94  typedef typename LB::Traits::RangeFieldType Scalar;
95  typedef typename LB::Traits::DomainFieldType Vector;
96  typename F::Traits::RangeType y;
97 
98  out.resize(36);
99  fill(out.begin(), out.end(), 0.0);
100 
101  const int qOrder = 3;
102  const QuadratureRule<Scalar,2>& rule1 = QuadratureRules<Scalar,2>::rule(GeometryType(GeometryType::cube,2), qOrder);
103 
104  for (typename QuadratureRule<Scalar,2>::const_iterator it = rule1.begin();
105  it != rule1.end(); ++it)
106  {
107  Dune::FieldVector<Scalar,2> qPos = it->position();
108  typename LB::Traits::DomainType localPos;
109 
110  localPos[0] = 0.0;
111  localPos[1] = qPos[0];
112  localPos[2] = qPos[1];
113  f.evaluate(localPos, y);
114  out[0] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*it->weight()*sign0;
115  out[6] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*it->weight();
116  out[12] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[1] - 1.0)*it->weight();
117  out[18] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight();
118 
119  localPos[0] = 1.0;
120  localPos[1] = qPos[0];
121  localPos[2] = qPos[1];
122  f.evaluate(localPos, y);
123  out[1] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*it->weight()*sign1;
124  out[7] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*it->weight();
125  out[13] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[1])*it->weight();
126  out[19] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight();
127 
128  localPos[0] = qPos[0];
129  localPos[1] = 0.0;
130  localPos[2] = qPos[1];
131  f.evaluate(localPos, y);
132  out[2] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*it->weight()*sign2;
133  out[8] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*it->weight();
134  out[14] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(2.0*qPos[1] - 1.0)*it->weight();
135  out[20] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight();
136 
137  localPos[0] = qPos[0];
138  localPos[1] = 1.0;
139  localPos[2] = qPos[1];
140  f.evaluate(localPos, y);
141  out[3] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*it->weight()*sign3;
142  out[9] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*it->weight();
143  out[15] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(1.0 - 2.0*qPos[1])*it->weight();
144  out[21] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight();
145 
146  localPos[0] = qPos[0];
147  localPos[1] = qPos[1];
148  localPos[2] = 0.0;
149  f.evaluate(localPos, y);
150  out[4] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*it->weight()*sign4;
151  out[10] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*it->weight();
152  out[16] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[1])*it->weight();
153  out[22] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*(2.0*qPos[1] - 1.0)*it->weight();
154 
155  localPos[0] = qPos[0];
156  localPos[1] = qPos[1];
157  localPos[2] = 1.0;
158  f.evaluate(localPos, y);
159  out[5] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*it->weight()*sign5;
160  out[11] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*it->weight();
161  out[17] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[1] - 1.0)*it->weight();
162  out[23] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*(2.0*qPos[1] - 1.0)*it->weight();
163  }
164 
165  const QuadratureRule<Vector,3>& rule2 = QuadratureRules<Vector,3>::rule(GeometryType(GeometryType::cube,3), qOrder);
166  for (typename QuadratureRule<Vector,3>::const_iterator it = rule2.begin();
167  it != rule2.end(); ++it)
168  {
169  FieldVector<double,3> qPos = it->position();
170 
171  f.evaluate(qPos, y);
172  out[24] += y[0]*it->weight();
173  out[25] += y[1]*it->weight();
174  out[26] += y[2]*it->weight();
175  out[27] += y[0]*qPos[1]*it->weight();
176  out[28] += y[0]*qPos[2]*it->weight();
177  out[29] += y[1]*qPos[0]*it->weight();
178  out[30] += y[1]*qPos[2]*it->weight();
179  out[31] += y[2]*qPos[0]*it->weight();
180  out[32] += y[2]*qPos[1]*it->weight();
181  out[33] += y[0]*qPos[1]*qPos[2]*it->weight();
182  out[34] += y[1]*qPos[0]*qPos[2]*it->weight();
183  out[35] += y[2]*qPos[0]*qPos[1]*it->weight();
184  }
185  }
186 
187 private:
188  typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3, sign4, sign5;
189  typename LB::Traits::DomainType n0, n1, n2, n3, n4, n5;
190  };
191 }
192 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH