3 #ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1_CUBE3D_LOCALINTERPOLATION_HH
4 #define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1_CUBE3D_LOCALINTERPOLATION_HH
8 #include <dune/geometry/quadraturerules.hh>
29 sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
39 sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
93 template<
class F,
class C>
97 typedef typename LB::Traits::RangeFieldType Scalar;
98 typedef typename LB::Traits::DomainFieldType Vector;
100 auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
103 fill(out.begin(), out.end(), 0.0);
105 const int qOrder = 3;
106 const QuadratureRule<Scalar,2>& rule1 = QuadratureRules<Scalar,2>::rule(GeometryTypes::cube(2), qOrder);
108 for (
typename QuadratureRule<Scalar,2>::const_iterator it = rule1.begin();
109 it != rule1.end(); ++it)
111 Dune::FieldVector<Scalar,2> qPos = it->position();
112 typename LB::Traits::DomainType localPos;
115 localPos[1] = qPos[0];
116 localPos[2] = qPos[1];
117 auto y = f(localPos);
118 out[0] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*it->weight()*sign0;
119 out[6] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*it->weight();
120 out[12] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[1] - 1.0)*it->weight();
121 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();
124 localPos[1] = qPos[0];
125 localPos[2] = qPos[1];
127 out[1] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*it->weight()*sign1;
128 out[7] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*it->weight();
129 out[13] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[1])*it->weight();
130 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();
132 localPos[0] = qPos[0];
134 localPos[2] = qPos[1];
136 out[2] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*it->weight()*sign2;
137 out[8] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*it->weight();
138 out[14] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(2.0*qPos[1] - 1.0)*it->weight();
139 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();
141 localPos[0] = qPos[0];
143 localPos[2] = qPos[1];
145 out[3] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*it->weight()*sign3;
146 out[9] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*it->weight();
147 out[15] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(1.0 - 2.0*qPos[1])*it->weight();
148 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();
150 localPos[0] = qPos[0];
151 localPos[1] = qPos[1];
154 out[4] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*it->weight()*sign4;
155 out[10] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*it->weight();
156 out[16] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[1])*it->weight();
157 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();
159 localPos[0] = qPos[0];
160 localPos[1] = qPos[1];
163 out[5] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*it->weight()*sign5;
164 out[11] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*it->weight();
165 out[17] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[1] - 1.0)*it->weight();
166 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();
169 const QuadratureRule<Vector,3>& rule2 = QuadratureRules<Vector,3>::rule(GeometryTypes::cube(3), qOrder);
170 for (
typename QuadratureRule<Vector,3>::const_iterator it = rule2.begin();
171 it != rule2.end(); ++it)
173 FieldVector<double,3> qPos = it->position();
176 out[24] += y[0]*it->weight();
177 out[25] += y[1]*it->weight();
178 out[26] += y[2]*it->weight();
179 out[27] += y[0]*qPos[1]*it->weight();
180 out[28] += y[0]*qPos[2]*it->weight();
181 out[29] += y[1]*qPos[0]*it->weight();
182 out[30] += y[1]*qPos[2]*it->weight();
183 out[31] += y[2]*qPos[0]*it->weight();
184 out[32] += y[2]*qPos[1]*it->weight();
185 out[33] += y[0]*qPos[1]*qPos[2]*it->weight();
186 out[34] += y[1]*qPos[0]*qPos[2]*it->weight();
187 out[35] += y[2]*qPos[0]*qPos[1]*it->weight();
192 typename LB::Traits::RangeFieldType sign0, sign1, sign2, sign3, sign4, sign5;
193 typename LB::Traits::DomainType n0, n1, n2, n3, n4, n5;
196 #endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1_CUBE3D_LOCALINTERPOLATION_HH