1 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH
2 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH
6 #include <dune/geometry/quadraturerules.hh>
26 sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
36 sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
90 template<
class F,
class C>
94 typedef typename LB::Traits::RangeFieldType Scalar;
95 typedef typename LB::Traits::DomainFieldType Vector;
96 typename F::Traits::RangeType y;
99 fill(out.begin(), out.end(), 0.0);
101 const int qOrder = 3;
102 const QuadratureRule<Scalar,2>& rule1 = QuadratureRules<Scalar,2>::rule(GeometryType(GeometryType::cube,2), qOrder);
104 for (
typename QuadratureRule<Scalar,2>::const_iterator it = rule1.begin();
105 it != rule1.end(); ++it)
107 Dune::FieldVector<Scalar,2> qPos = it->position();
108 typename LB::Traits::DomainType localPos;
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();
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();
128 localPos[0] = qPos[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();
137 localPos[0] = qPos[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();
146 localPos[0] = qPos[0];
147 localPos[1] = qPos[1];
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();
155 localPos[0] = qPos[0];
156 localPos[1] = qPos[1];
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();
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)
169 FieldVector<double,3> qPos = it->position();
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();
188 typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3, sign4, sign5;
189 typename LB::Traits::DomainType n0, n1, n2, n3, n4, n5;
192 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1Q3DLOCALINTERPOLATION_HH