[WIP] add debug output

dune/mmdg/dg.hh

@@ -36,6 +36,17 @@ public:
 
     //NOTE:
 //    const int order = 3; //order of the quadrature rule
+    for(const auto& elem : elements(gridView_))
+    {
+      const auto& geo = elem.geometry();
+      std::cout << mapper_.index(elem) << "\t";
+      for (int k = 0; k < geo.corners(); k++)
+      {
+        std::cout << geo.corner(k) << "\t";
+      }
+      std::cout << "\n\n";
+    }
+
 
     //we use the basis
     // phi_elem,0 (x) = indicator(elem);
@@ -48,6 +59,13 @@ public:
       const double elemVol = geo.volume();
       const auto& center = geo.center();
 
+      std::cout << "===============================\nElement " << elemIdx;
+      for (int k = 0; k < geo.corners(); k++)
+      {
+        std::cout << "\t" << geo.corner(k);
+      }
+      std::cout << "\n\n";
+
       //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
@@ -99,6 +117,13 @@ public:
         const double intersctVol = intersctGeo.volume();
         const auto& intersctCenter = intersctGeo.center();
 
+        std::cout << "+++++++++++++++++++++++++++\nIntersection ";
+        for (int k = 0; k < intersctGeo.corners(); k++)
+        {
+          std::cout << "\t" << intersctGeo.corner(k);
+        }
+        std::cout << "\n\n";
+
         //TODO: quadrature rule cannot be used for dim = 1!
 //        const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
 //          Dune::QuadratureRules<double,dim-1>::rule(
@@ -121,11 +146,22 @@ public:
           {
             const auto& leftCorner = intersctGeo.corner(0);
             const auto& rightCorner = intersctGeo.corner(1);
+            std::cout << "left corner: " << leftCorner << "\tright corner: " << rightCorner << "\n";
+
             quadraticIntregrals[i][j] = intersctVol / 3 *
               ( leftCorner[i] * leftCorner[j] +
                 rightCorner[i] * rightCorner[j]
               + 0.5 * (leftCorner[i] * rightCorner[j] +
                 leftCorner[j] * rightCorner[i]) );
+
+            std::cout << i << ", " << j << "\t" << leftCorner[i] * leftCorner[j] << "\t" <<
+              rightCorner[i] * rightCorner[j] << "\t" <<
+              0.5 * (leftCorner[i] * rightCorner[j] +
+              leftCorner[j] * rightCorner[i]) << "\t" << intersctVol / 3 *
+                ( leftCorner[i] * leftCorner[j] +
+                  rightCorner[i] * rightCorner[j]
+                + 0.5 * (leftCorner[i] * rightCorner[j] +
+                  leftCorner[j] * rightCorner[i]) ) << "\t" << quadraticIntregrals[i][j] << "\n\n";
   /*
             //use second order quadrature rule for exact evaluation of
             // int_intersct x_i*x_j ds
@@ -141,6 +177,8 @@ public:
             quadraticIntregrals[i][j] = quadraticIntregrals[j][i];
           }
         }
+        std::cout << "linearIntegrals:\n" << linearIntegrals << "\n\n";
+        std::cout << "quadraticIntregrals:\n" << quadraticIntregrals << "\n\n";
 
         //exact evaluation of
         // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
@@ -150,6 +188,14 @@ public:
 
         if (intersct.neighbor()) //intersct has neighboring element
         {
+          std::cout << "------------------------------\nNeighbor\t";
+          for (int k = 0; k < intersct.outside().geometry().corners(); k++)
+          {
+            std::cout << "\t" << intersct.outside().geometry().corner(k);
+          }
+          std::cout << "\n\n";
+
+
           //index of the neighboring element
           const int neighborIdx = mapper_.index(intersct.outside());
           const int neighborIdxSLE = (dim + 1)*neighborIdx;
@@ -298,6 +344,7 @@ public:
         }
         else //boundary facet
         {
+          std::cout << "------------------------------\nBoundary\n\n";
           for (int i = 0; i < dim; i++)
           { //we use the relations
             // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
@@ -364,10 +411,10 @@ public:
             std::cout << i << ", " << j << std::endl;
 
     //NOTE: dim = 1
-    auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);};
+//    auto exactSolution = [](const auto& pos) {return 0.5*pos[0]*(pos[0] - 1);};
 
     //NOTE: dim = 2
-/*    auto exactSolution = [](const auto& pos)
+    auto exactSolution = [](const auto& pos)
       {
         double solution = 0.0;
         for (int m = 1; m <= 21; m = m + 2)
@@ -381,7 +428,7 @@ public:
         solution *= -16 / pow(M_PI, 4);
         return solution;
       };
-*/
+