use quadrature rule for quadratic integrals in dg.hh

dune/mmdg/dg.hh

@@ -78,6 +78,7 @@ private:
       QRulesEdge::rule(edge, problem_.quadratureOrder_K());
     const QRuleEdge& rule_Kplus1_Edge =
       QRulesEdge::rule(edge, problem_.quadratureOrder_K() + 1);
+    const QRuleEdge& rule_2ndOrder = QRulesEdge::rule(edge, 2);
 
     //we use the basis
     // phi_elem,0 (x) = indicator(elem);
@@ -147,11 +148,6 @@ private:
         const double intersctVol = intersctGeo.volume();
         const auto& intersctCenter = intersctGeo.center();
 
-        //TODO: quadrature rule cannot be used for dim = 1!
-  //        const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
-  //          Dune::QuadratureRules<double,dim-1>::rule(
-  //          intersct.type(), 2);
-
         //create storage for multiply used integral values
         //linearIntegrals[i] = int_intersct x[i] ds
         Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
@@ -168,37 +164,38 @@ private:
         //KIntegralsX[i][j] = int_intersct x[j]*(K[:,i]*normal) ds
         Dune::FieldMatrix<Scalar, dim, dim> KIntegralsX(0.0);
 
-        //fill vector linearIntegrals and matrix quadraticIntregrals //NOTE
-        for (int i = 0; i < dim; i++)
+        //fill vector linearIntegrals and matrix quadraticIntregrals
+        for (const auto& ip : rule_2ndOrder)
         {
-          //use midpoint rule for exact evaluation of
-          // int_intersct x[i] ds
-          linearIntegrals[i] = intersctVol * intersctCenter[i];
+          //quadrature point in local and global coordinates
+          const auto& qpLocal = ip.position();
+          const auto& qpGlobal = intersctGeo.global(qpLocal);
 
-          for (int j = 0; j <= i; j++)
+          const Scalar weight(ip.weight());
+          const Scalar integrationElem(intersctGeo.integrationElement(qpLocal));
+
+          for (int i = 0; i < dim; i++)
           {
-            const auto& leftCorner = intersctGeo.corner(0);
-            const auto& rightCorner = intersctGeo.corner(1);
-
-            quadraticIntregrals[i][j] = intersctVol / 3 *
-              ( leftCorner[i] * leftCorner[j] + rightCorner[i] * rightCorner[j]
-              + 0.5 * (leftCorner[i] * rightCorner[j] +
-                leftCorner[j] * rightCorner[i]) );
-  /*
-            //use second order quadrature rule for exact evaluation of
-            // int_intersct x_i*x_j ds
-            //NOTE: use Simpson's rule instead manually?
-            for (const auto& ip : secondOrderRule)
+            //use midpoint rule for exact evaluation of
+            // int_intersct x[i] ds
+            linearIntegrals[i] = intersctVol * intersctCenter[i];
+
+            for (int j = 0; j <= i; j++)
             {
-              const auto& qp = ip.position();
-              quadraticIntregrals[i][j] += //NOTE: volume necessary?
-                ip.weight() * qp[i] * qp[j] * intersctVol;
+              quadraticIntregrals[i][j] +=
+                weight * integrationElem * qpGlobal[i] * qpGlobal[j];
+              }
+            }
+          }
+
+          //quadraticIntregrals is symmetric
+          for (int i = 0; i < dim; i++)
+          {
+            for (int j = 0; j < i; j++)
+            {
+              quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
             }
-  */
-            //NOTE: probably unnecessary
-            quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
           }
-        }
 
         //fill vector KIntegrals using a quadrature rule
         for (const auto& ip : rule_K_Edge)