diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh index 28f71d811225743be7edee4d36453d7fb2fbd3b3..014257a2d65ddcffbe86c1f178539fdc4d200348 100644 --- a/dune/mmdg/dg.hh +++ b/dune/mmdg/dg.hh @@ -26,7 +26,7 @@ public: using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>; using VTKFunction = Dune::VTKFunction<GridView>; using P1Function = Dune::NonconformingP1VTKFunction<GridView, - Dune::DynamicVector<double>>; + Dune::DynamicVector<Scalar>>; //constructor DG (const GridView& gridView, const Mapper& mapper, @@ -63,6 +63,11 @@ private: //assemble stiffness matrix A and load vector b void assembleSLE (const Scalar K, const Scalar mu) { + //quadrature rule + const Dune::QuadratureRule<Scalar, dim>& rule = + Dune::QuadratureRules<Scalar, dim>::rule( + Dune::GeometryTypes::simplex(dim), problem_.quadratureOrder()); + //we use the basis // phi_elem,0 (x) = indicator(elem); // phi_elem,i (x) = x[i]*indicator(elem); @@ -72,43 +77,34 @@ private: const int elemIdx = mapper_.index(elem); const auto& geo = elem.geometry(); const double elemVol = geo.volume(); - const auto& center = geo.center(); //in the system of linear equations (SLE) Ad = b, //the index elemIdxSLE refers to the basis function phi_elem,0 //and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the //basis function phi_elem,i const int elemIdxSLE = (dim + 1)*elemIdx; -/* - //TODO: can be done outside of the loop? - const Dune::QuadratureRule<double,dim>& rule = - Dune::QuadratureRules<double,dim>::rule(geo.type(), problem_.quadratureOrder()); - // Dune::FieldVector<double, dim+1> update(0.0); - //NOTE: how are quadrature rules in Dune applied correctly? + //fill load vector b using a quadrature rule for (const auto& ip : rule) { - const auto& qp = ip.position(); - //NOTE: is the volume of the element taken into account - //automatically? - const double weight = ip.weight(); + //quadrature point in local and global coordinates + const auto& qpLocal = ip.position(); + const auto& qpGlobal = geo.global(qpLocal); + + const Scalar weight(ip.weight()); + const Scalar integrationElem(geo.integrationElement(qpLocal)); + const Scalar qEvalution(problem_.q(qpGlobal)); //quadrature for int_elem q*phi_elem,0 dV - b[elemIdxSLE] += weight * problem_.q(geo.global(qp)) * geo.integrationElement(qp); + b[elemIdxSLE] += weight * qEvalution * integrationElem; //quadrature for int_elem q*phi_elem,i dV for (int i = 0; i < dim; i++) { - b[elemIdxSLE + i + 1] += weight * qp[i] * problem_.q(geo.global(qp)) * geo.integrationElement(qp); + b[elemIdxSLE + i + 1] += + weight * qpGlobal[i] * qEvalution * integrationElem; } } -*/ - //NOTE: makeshift solution for source term q = -1 - b[elemIdxSLE] += -elemVol; - for (int i = 0; i < dim; i++) - { - b[elemIdxSLE + i + 1] += -elemVol * center[i]; - } for (int i = 0; i < dim; i++) {