37 for (
size_t i=0; i<3; i++)
38 sign_[i] = (s[i]) ? -1.0 : 1.0;
42 n_[2] = { 1.0/sqrt(2.0), 1.0/sqrt(2.0)};
44 c_ = { 0.5*n_[0][0] - 1.0*n_[0][1],
45 -1.0*n_[1][0] + 0.5*n_[1][1],
46 0.5*n_[2][0] + 0.5*n_[2][1]};
61 typedef typename LB::Traits::RangeFieldType Scalar;
62 typedef typename LB::Traits::DomainFieldType Vector;
65 fill(out.begin(), out.end(), 0.0);
67 const int qOrder1 = 4;
68 const auto& rule1 = Dune::QuadratureRules<Scalar,1>::rule(Dune::GeometryTypes::simplex(1), qOrder1);
70 for (
auto&& qp : rule1)
72 Scalar qPos = qp.position();
73 typename LB::Traits::DomainType localPos;
75 localPos = {qPos, 0.0};
77 out[0] += (y[0]*n_[0][0] + y[1]*n_[0][1])*qp.weight()*sign_[0]/c_[0];
78 out[3] += (y[0]*n_[0][0] + y[1]*n_[0][1])*(2.0*qPos - 1.0)*qp.weight()/c_[0];
80 localPos = {0.0, qPos};
82 out[1] += (y[0]*n_[1][0] + y[1]*n_[1][1])*qp.weight()*sign_[1]/c_[1];
83 out[4] += (y[0]*n_[1][0] + y[1]*n_[1][1])*(1.0 - 2.0*qPos)*qp.weight()/c_[1];
85 localPos = {1.0 - qPos, qPos};
87 out[2] += (y[0]*n_[2][0] + y[1]*n_[2][1])*qp.weight()*sign_[2]/c_[2];
88 out[5] += (y[0]*n_[2][0] + y[1]*n_[2][1])*(2.0*qPos - 1.0)*qp.weight()/c_[2];
91 const int qOrder2 = 8;
92 const auto& rule2 = Dune::QuadratureRules<Vector,2>::rule(Dune::GeometryTypes::simplex(2), qOrder2);
94 for (
auto&& qp : rule2)
96 auto qPos = qp.position();
99 out[6] += y[0]*qp.weight();
100 out[7] += y[1]*qp.weight();