From d74c68d23dca531d1e51ddcb10e760df80a0bea4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
 <maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Sat, 18 Jan 2020 12:31:24 +0100
Subject: [PATCH] transfer output to separate method, add check
 hasExactSolution() in dg.hh

---
 dune/mmdg/dg.hh | 101 +++++++++++++++++++++++++++---------------------
 1 file changed, 57 insertions(+), 44 deletions(-)

diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 6722715..4ee978c 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -20,6 +20,9 @@ public:
   static constexpr int dim = GridView::dimension;
   using Matrix = Dune::DynamicMatrix<Scalar>; //NOTE: what is an appropriate sparse matrix type? -> BCRS
   using Vector = Dune::DynamicVector<Scalar>;
+  using VTKFunction = Dune::VTKFunction<GridView>;
+  using P1Function = Dune::NonconformingP1VTKFunction<GridView,
+    Dune::DynamicVector<double>>;
 
   //constructor
   DG (const GridView& gridView, const Mapper& mapper,
@@ -40,50 +43,8 @@ public:
     //NOTE: what would be an appropiate solver here?
     A.solve(d, b);
 
-    //storage for pressure data, for each element we store the
-    //pressure at the corners of the element, the VTKFunction will
-    //create the output using a linear interpolation for each element
-    Dune::DynamicVector<Scalar> pressure(dof, 0.0);
-    Dune::DynamicVector<double> exactPressure(dof, 0.0);
-
-    for (const auto& elem : elements(gridView_))
-    {
-      const int elemIdxSLE = (dim + 1)*mapper_.index(elem);
-      const auto& geo = elem.geometry();
-
-      for (int k = 0; k < geo.corners(); k++)
-      {
-        exactPressure[elemIdxSLE + k] = problem_.exactSolution(geo.corner(k));
-
-        //contribution of the basis function
-        // phi_elem,0 (x) = indicator(elem);
-        //at the kth corner of elem
-        pressure[elemIdxSLE + k] = d[elemIdxSLE];
-
-        for (int i = 0; i < dim; i++)
-        {
-          //contribution of the basis function
-          // phi_elem,i (x) = x[i]*indicator(elem);
-          //at the kth corner of elem
-          pressure[elemIdxSLE + k] += d[elemIdxSLE + i + 1] * geo.corner(k)[i];
-        }
-      }
-    }
-
-    //create output directory if necessary
-    mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
-
-    Dune::VTKWriter<GridView> vtkWriter(gridView_, Dune::VTK::nonconforming);
-
-    using VTKFunction = Dune::VTKFunction<GridView>;
-    using P1Function = Dune::NonconformingP1VTKFunction<GridView,
-      Dune::DynamicVector<double>>;
-
-    vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
-      new P1Function(gridView_, pressure, "pressure")));
-    vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
-      new P1Function(gridView_, exactPressure, "exactPressure")));
-    vtkWriter.pwrite("pressureData", "", "data");
+    //write solution to a vtk file
+    writeVTKoutput();
   }
 
   const GridView& gridView_;
@@ -381,6 +342,58 @@ private:
             std::cout << i << ", " << j << std::endl;
   }
 
+  //writes the solution to a vtk file
+  void writeVTKoutput () const
+  {
+    //storage for pressure data, for each element we store the
+    //pressure at the corners of the element, the VTKFunction will
+    //create the output using a linear interpolation for each element
+    Dune::DynamicVector<Scalar> pressure(dof, 0.0);
+    Dune::DynamicVector<Scalar> exactPressure(dof, 0.0);
+
+    for (const auto& elem : elements(gridView_))
+    {
+      const int elemIdxSLE = (dim + 1)*mapper_.index(elem);
+      const auto& geo = elem.geometry();
+
+      for (int k = 0; k < geo.corners(); k++)
+      {
+        if (problem_.hasExactSolution())
+        {
+          exactPressure[elemIdxSLE + k] = problem_.exactSolution(geo.corner(k));
+        }
+
+        //contribution of the basis function
+        // phi_elem,0 (x) = indicator(elem);
+        //at the kth corner of elem
+        pressure[elemIdxSLE + k] = d[elemIdxSLE];
+
+        for (int i = 0; i < dim; i++)
+        {
+          //contribution of the basis function
+          // phi_elem,i (x) = x[i]*indicator(elem);
+          //at the kth corner of elem
+          pressure[elemIdxSLE + k] += d[elemIdxSLE + i + 1] * geo.corner(k)[i];
+        }
+      }
+    }
+
+    //create output directory if necessary
+    mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+
+    Dune::VTKWriter<GridView> vtkWriter(gridView_, Dune::VTK::nonconforming);
+    vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
+      new P1Function(gridView_, pressure, "pressure")));
+
+    if (problem_.hasExactSolution())
+    {
+      vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
+        new P1Function(gridView_, exactPressure, "exactPressure")));
+    }
+
+    vtkWriter.pwrite("pressureData", "", "data");
+  }
+
 
   Matrix A; //stiffness matrix
   Vector b; //load vector
-- 
GitLab