use quadrature rule for load vector 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++)
       {