96 typedef typename LB::Traits::RangeFieldType Scalar;
97 typedef typename LB::Traits::DomainFieldType Vector;
98 typename F::Traits::RangeType y;
101 fill(out.begin(), out.end(), 0.0);
103 const int qOrder = 3;
104 const QuadratureRule<Scalar,2>& rule1 = QuadratureRules<Scalar,2>::rule(GeometryType(GeometryType::cube,2), qOrder);
106 for (
typename QuadratureRule<Scalar,2>::const_iterator it = rule1.begin();
107 it != rule1.end(); ++it)
109 Dune::FieldVector<Scalar,2> qPos = it->position();
110 typename LB::Traits::DomainType localPos;
113 localPos[1] = qPos[0];
114 localPos[2] = qPos[1];
115 f.evaluate(localPos, y);
116 out[0] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*it->weight()*sign0;
117 out[6] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[0] - 1.0)*it->weight();
118 out[12] += (y[0]*n0[0] + y[1]*n0[1] + y[2]*n0[2])*(2.0*qPos[1] - 1.0)*it->weight();
119 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();
122 localPos[1] = qPos[0];
123 localPos[2] = qPos[1];
124 f.evaluate(localPos, y);
125 out[1] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*it->weight()*sign1;
126 out[7] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[0])*it->weight();
127 out[13] += (y[0]*n1[0] + y[1]*n1[1] + y[2]*n1[2])*(1.0 - 2.0*qPos[1])*it->weight();
128 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();
130 localPos[0] = qPos[0];
132 localPos[2] = qPos[1];
133 f.evaluate(localPos, y);
134 out[2] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*it->weight()*sign2;
135 out[8] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(1.0 - 2.0*qPos[0])*it->weight();
136 out[14] += (y[0]*n2[0] + y[1]*n2[1] + y[2]*n2[2])*(2.0*qPos[1] - 1.0)*it->weight();
137 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();
139 localPos[0] = qPos[0];
141 localPos[2] = qPos[1];
142 f.evaluate(localPos, y);
143 out[3] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*it->weight()*sign3;
144 out[9] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(2.0*qPos[0] - 1.0)*it->weight();
145 out[15] += (y[0]*n3[0] + y[1]*n3[1] + y[2]*n3[2])*(1.0 - 2.0*qPos[1])*it->weight();
146 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();
148 localPos[0] = qPos[0];
149 localPos[1] = qPos[1];
151 f.evaluate(localPos, y);
152 out[4] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*it->weight()*sign4;
153 out[10] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[0])*it->weight();
154 out[16] += (y[0]*n4[0] + y[1]*n4[1] + y[2]*n4[2])*(1.0 - 2.0*qPos[1])*it->weight();
155 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();
157 localPos[0] = qPos[0];
158 localPos[1] = qPos[1];
160 f.evaluate(localPos, y);
161 out[5] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*it->weight()*sign5;
162 out[11] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[0] - 1.0)*it->weight();
163 out[17] += (y[0]*n5[0] + y[1]*n5[1] + y[2]*n5[2])*(2.0*qPos[1] - 1.0)*it->weight();
164 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();
167 const QuadratureRule<Vector,3>& rule2 = QuadratureRules<Vector,3>::rule(GeometryType(GeometryType::cube,3), qOrder);
168 for (
typename QuadratureRule<Vector,3>::const_iterator it = rule2.begin();
169 it != rule2.end(); ++it)
171 FieldVector<double,3> qPos = it->position();
174 out[24] += y[0]*it->weight();
175 out[25] += y[1]*it->weight();
176 out[26] += y[2]*it->weight();
177 out[27] += y[0]*qPos[1]*it->weight();
178 out[28] += y[0]*qPos[2]*it->weight();
179 out[29] += y[1]*qPos[0]*it->weight();
180 out[30] += y[1]*qPos[2]*it->weight();
181 out[31] += y[2]*qPos[0]*it->weight();
182 out[32] += y[2]*qPos[1]*it->weight();
183 out[33] += y[0]*qPos[1]*qPos[2]*it->weight();
184 out[34] += y[1]*qPos[0]*qPos[2]*it->weight();
185 out[35] += y[2]*qPos[0]*qPos[1]*it->weight();