From 523c9fb7c3350472783a5a0ae160f9c229772fc1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
 <maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Fri, 20 Mar 2020 12:35:06 +0100
Subject: [PATCH] use quadratureFactor instead of weight*intElem in dg

---
 dune/mmdg/dg.hh | 42 ++++++++++++++++++++----------------------
 1 file changed, 20 insertions(+), 22 deletions(-)

diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index f352e23..0b8059c 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -111,7 +111,7 @@ public:
   const Mapper& mapper_; //element mapper for gridView_
   const Problem& problem_; //the DG problem to be solved
 
-  const int dof_;
+  const int dof_; //degrees of freedom
 
 protected:
   DG (const Grid& grid, const GridView& gridView, const Mapper& mapper,
@@ -156,25 +156,24 @@ protected:
       //basis function phi_elem,i
       const int elemIdxSLE = (dim + 1)*elemIdx;
 
-      //fill load vector b using a quadrature rule (bulk term)
+      //fill load vector b using a quadrature rule
       for (const auto& ip : rule_q)
       {
         //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 quadratureFactor =
+          ip.weight() * geo.integrationElement(qpLocal);
         const Scalar qEvaluation(problem_.q(qpGlobal));
 
         //quadrature for int_elem q*phi_elem,0 dV
-        b[elemIdxSLE] += weight * integrationElem * qEvaluation;
+        b[elemIdxSLE] += quadratureFactor * qEvaluation;
 
         //quadrature for int_elem q*phi_elem,i dV
         for (int i = 0; i < dim; i++)
         {
-          b[elemIdxSLE + i + 1] +=
-            weight * integrationElem * qpGlobal[i] * qEvaluation;
+          b[elemIdxSLE + i + 1] += quadratureFactor * qpGlobal[i] * qEvaluation;
         }
       }
 
@@ -188,15 +187,15 @@ protected:
         const auto& qpLocal = ip.position();
         const auto& qpGlobal = geo.global(qpLocal);
 
-        const Scalar weight(ip.weight());
-        const Scalar integrationElem(geo.integrationElement(qpLocal));
+        const Scalar quadratureFactor =
+          ip.weight() * geo.integrationElement(qpLocal);
 
         for (int i = 0; i < dim; i++)
         {
           for (int j = 0; j <= i; j++)
           {
             A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
-              weight * integrationElem * problem_.K(qpGlobal)[j][i];
+              quadratureFactor * problem_.K(qpGlobal)[j][i];
           }
         }
       }
@@ -240,19 +239,19 @@ protected:
           const auto& qpLocal = ip.position();
           const auto& qpGlobal = intersctGeo.global(qpLocal);
 
-          const Scalar weight(ip.weight());
-          const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
+          const Scalar quadratureFactor =
+            ip.weight() * intersctGeo.integrationElement(qpLocal);
 
           for (int i = 0; i < dim; i++)
           {
             //use midpoint rule for exact evaluation of
-            // int_intersct x[i] ds
+            // int_intersct x[i] ds //TODO
             linearIntegrals[i] = intersctVol * intersctCenter[i];
 
             for (int j = 0; j <= i; j++)
             {
               quadraticIntregrals[i][j] +=
-                weight * integrationElem * qpGlobal[i] * qpGlobal[j];
+                quadratureFactor * qpGlobal[i] * qpGlobal[j];
               }
             }
           }
@@ -273,15 +272,14 @@ protected:
           const auto& qpLocal = ip.position();
           const auto& qpGlobal = intersctGeo.global(qpLocal);
 
-          const Scalar weight(ip.weight());
-          const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
+          const Scalar quadratureFactor =
+            ip.weight() * intersctGeo.integrationElement(qpLocal);
           const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
             problem_.K(qpGlobal);
 
           for (int i = 0; i < dim; i++)
           {
-            KIntegrals[i] +=
-              weight * integrationElem * (KEvaluation[i] * normal);
+            KIntegrals[i] += quadratureFactor * (KEvaluation[i] * normal);
           }
         }
 
@@ -292,8 +290,8 @@ protected:
           const auto& qpLocal = ip.position();
           const auto& qpGlobal = intersctGeo.global(qpLocal);
 
-          const Scalar weight(ip.weight());
-          const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
+          const Scalar quadratureFactor =
+            ip.weight() * intersctGeo.integrationElement(qpLocal);
           const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
             problem_.K(qpGlobal);
 
@@ -301,7 +299,7 @@ protected:
           {
             for (int j = 0; j < dim; j++)
             {
-              KIntegralsX[i][j] += weight * integrationElem *
+              KIntegralsX[i][j] += quadratureFactor *
                 qpGlobal[j] * (KEvaluation[i] * normal);
             }
           }
@@ -315,7 +313,7 @@ protected:
         {
           //index of the neighboring element
           const int neighborIdx = mapper_.index(intersct.outside());
-          const int neighborIdxSLE = (dim + 1)*neighborIdx;
+          const int neighborIdxSLE = (dim + 1) * neighborIdx;
 
           for (int i = 0; i < dim; i++)
           { //we use the relations
-- 
GitLab