[bugfix] use nonforming vtk writer in dg.hh

dune/mmdg/dg.hh

@@ -26,13 +26,13 @@ public:
 
   const void operator() (const Scalar K, const Scalar mu) const
   {
-    //NOTE: what is an appropriate sparse matrix type?
+    //NOTE: what is an appropriate sparse matrix type? -> BCRS
     const int dof = (1 + dim)*gridView_.size(0); //degrees of freedom
     using Matrix = Dune::DynamicMatrix<Scalar>;
     using Vector = Dune::DynamicVector<Scalar>;
     Matrix A(dof, dof, 0.0); //stiffness matrix
     Vector b(dof, 0.0); //load vector
-    Vector d;
+    Vector d(dof, 0.0);
 
     //NOTE:
 //    const int order = 3; //order of the quadrature rule
@@ -311,6 +311,8 @@ public:
       }
     }
 
+    std::cout << A << std::endl;
+
     //NOTE: what would be an appropiate solver here?
     A.solve(d, b);
 
@@ -348,11 +350,28 @@ public:
       }
     }
 
+    Dune::DynamicVector<Scalar>
+      pressureTEST(gridView_.size(0), 0.0);
+
+    for (const auto& elem : elements(gridView_))
+    {
+      const int elemIdx = mapper_.index(elem);
+      const int elemIdxSLE = (dim + 1)*elemIdx;
+      const auto& center = elem.geometry().center();
+      pressureTEST[elemIdx] = d[elemIdxSLE];
+
+      for (int i = 0; i < dim; i++)
+      {
+        pressureTEST[elemIdx] += d[elemIdxSLE + i + 1] * center[i];
+      }
+    }
+
+
     //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)
@@ -366,7 +385,7 @@ public:
       solution *= -16 / pow(M_PI, 4);
       return solution;
     };
-*/
+
 
     Dune::DynamicVector<double> exactPressure(gridView_.size(0));
     for (const auto& elem : elements(gridView_))
@@ -379,10 +398,10 @@ public:
     //create output directory if necessary
     mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
 
-    auto vtkWriter =
-      std::make_shared<Dune::VTKWriter<GridView>>(gridView_);
-    Dune::VTKSequenceWriter<GridView>
-      vtkSqWriter(vtkWriter, "pressure", "data", "data");
+    Dune::VTKWriter<GridView> vtkWriter(gridView_, Dune::VTK::nonconforming);
+  //    std::make_shared<Dune::VTKWriter<GridView>>(gridView_);
+//    Dune::VTKSequenceWriter<GridView>
+//      vtkSqWriter(vtkWriter, "pressure", "data", "data");
 
     using VTKFunction = Dune::VTKFunction<GridView>;
     using P1Function = Dune::NonconformingP1VTKFunction<GridView,
@@ -391,10 +410,11 @@ public:
       new P1Function(gridView_, pressure, "pressure");
 
     //NOTE: why does this only work with a sequence writer?
-    vtkSqWriter.addVertexData(
+    vtkWriter.addVertexData(
       std::shared_ptr<const VTKFunction>(vtkOperator) );
-    vtkSqWriter.addCellData(exactPressure, "exactPressure");
-    vtkSqWriter.write(0.0);
+    vtkWriter.addCellData(exactPressure, "exactPressure");
+    vtkWriter.addCellData(pressureTEST, "pressureTEST");
+    vtkWriter.write("linearpressure");
   }
 
   const GridView& gridView_;