From 5880f105e9b69d7a3d4b2012c54534ece2c5b0fe 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 18:28:27 +0100
Subject: [PATCH] remove template parameter Grid

---
 dune/mmdg/dg.hh        | 18 +++++++-----------
 dune/mmdg/mmdg.hh      | 16 ++++++++--------
 src/convergenceTest.py |  6 +++---
 src/dg.cc              | 13 ++++++-------
 src/mmdg.cc            |  4 ++--
 5 files changed, 26 insertions(+), 31 deletions(-)

diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 0b8059c..f0f3495 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -17,7 +17,7 @@
 
 #include <dune/mmdg/nonconformingp1vtkfunction.hh>
 
-template<class Grid, class GridView, class Mapper, class Problem>
+template<class GridView, class Mapper, class Problem>
 class DG
 {
 public:
@@ -37,10 +37,8 @@ public:
   using QRulesEdge = Dune::QuadratureRules<Scalar, dim-1>;
 
   //constructor
-  DG (const Grid& grid, const GridView& gridView, const Mapper& mapper, const Problem& problem) :
-  //NOTE: why does this not work?
-    //grid_(grid), gridView_(grid.leafGridView()), mapper_(*(new Mapper(gridView_, Dune::mcmgElementLayout()))), problem_(problem), dof_((1 + dim) * gridView_.size(0))
-    grid_(grid), gridView_(gridView), mapper_(mapper), problem_(problem),
+  DG (const GridView& gridView, const Mapper& mapper, const Problem& problem) :
+    gridView_(gridView), mapper_(mapper), problem_(problem),
     dof_((1 + dim) * gridView.size(0))
     {
       //initialize stiffness matrix A, load vector b and solution vector d
@@ -106,7 +104,6 @@ public:
     return std::sqrt(error);
   }
 
-  const Grid& grid_;
   const GridView& gridView_;
   const Mapper& mapper_; //element mapper for gridView_
   const Problem& problem_; //the DG problem to be solved
@@ -114,10 +111,9 @@ public:
   const int dof_; //degrees of freedom
 
 protected:
-  DG (const Grid& grid, const GridView& gridView, const Mapper& mapper,
-    const Problem& problem, const int dof) :
-    grid_(grid), gridView_(gridView), mapper_(mapper), problem_(problem),
-    dof_(dof)
+  DG (const GridView& gridView, const Mapper& mapper, const Problem& problem,
+    const int dof) :
+    gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof)
     {
       //initialize stiffness matrix A, load vector b and solution vector d
       A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
@@ -205,7 +201,7 @@ protected:
       {
         if constexpr (dim != 1)
         {
-          if (grid_.isInterface(intersct))
+          if (gridView_.grid().isInterface(intersct))
           {
             continue;
           }
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 62d2d0b..ccaef4d 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -3,15 +3,15 @@
 
 #include <dune/mmdg/dg.hh>
 
-template<class Grid, class GridView, class Mapper, class InterfaceGridView,
+template<class GridView, class Mapper, class InterfaceGridView,
   class InterfaceMapper, class Problem>
-class MMDG : public DG<Grid, GridView, Mapper, Problem>
+class MMDG : public DG<GridView, Mapper, Problem>
 {
 public:
   static constexpr int dim = GridView::dimension;
   static constexpr auto iSimplex = Dune::GeometryTypes::simplex(dim-1);
 
-  using Base = DG<Grid, GridView, Mapper, Problem>;
+  using Base = DG<GridView, Mapper, Problem>;
   using Scalar = typename Base::Scalar; //NOTE using Base::Scalar
   using Matrix = typename Base::Matrix;
   using Vector = typename Base::Vector;
@@ -24,10 +24,10 @@ public:
     Dune::DynamicVector<Scalar>>;
 
   //constructor
-  MMDG (const Grid& grid, const GridView& gridView, const Mapper& mapper,
+  MMDG (const GridView& gridView, const Mapper& mapper,
     const InterfaceGridView& iGridView, const InterfaceMapper& iMapper,
     const Problem& problem) :
-    Base(grid, gridView, mapper, problem,
+    Base(gridView, mapper, problem,
       (1 + dim) * gridView.size(0) + dim * iGridView.size(0)),
     iGridView_(iGridView), iMapper_(iMapper) {}
 
@@ -59,7 +59,7 @@ public:
       const QRuleInterface& qrule =
         QRulesInterface::rule(iSimplex, quadratureOrder);
       const InterfaceBasis iFrame = interfaceFrame(
-        Base::grid_.asIntersection(iElem).centerUnitOuterNormal());
+        Base::gridView_.grid().asIntersection(iElem).centerUnitOuterNormal());
 
       //determine L2-error on the interface element iElem using a
       //quadrature rule
@@ -147,7 +147,7 @@ private:
     {
       const int iElemIdx = iMapper_.index(iElem);
       const auto& iGeo = iElem.geometry();
-      const auto& intersct = Base::grid_.asIntersection(iElem);
+      const auto& intersct = Base::gridView_.grid().asIntersection(iElem);
 
       //get basis of the interface coordinate system for evaluation of
       //basis function phi_iElem,i with i=1,...,dim-1
@@ -667,7 +667,7 @@ private:
       //get basis of the interface coordinate system for evaluation of
       //basis function phi_iElem,i with i=1,...,dim-1
       const InterfaceBasis iFrame = interfaceFrame(
-        Base::grid_.asIntersection(iElem).centerUnitOuterNormal());
+        Base::gridView_.grid().asIntersection(iElem).centerUnitOuterNormal());
 
       for (int k = 0; k < iGeo.corners(); k++)
       {
diff --git a/src/convergenceTest.py b/src/convergenceTest.py
index 2048858..92ddfa4 100644
--- a/src/convergenceTest.py
+++ b/src/convergenceTest.py
@@ -31,7 +31,7 @@ def writeDGF(N, dim):
     file.close()
 
 
-# reads error, run time and number of elements, i.e. the output from the
+# reads error, run time and maximum edge length, i.e. the output from the
 # excuted discontinuous Galerkin scheme
 def readData():
     file = open("convergenceData", "r")
@@ -54,12 +54,12 @@ if not exists("plots"):
     mkdir("plots")
 
 # number of runs from which the run times are taken as mean value
-numberOfMeans = 3
+numberOfMeans = 1
 
 # numbers of grid elements per direction
 numElemsPerDir = [1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130]
 
-dimensions = [1,2]
+dimensions = [1, 2]
 
 for dim in dimensions:
     # storage for errors, run times and largestEdgeLength
diff --git a/src/dg.cc b/src/dg.cc
index 60d8760..744b48c 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -44,7 +44,7 @@ int main(int argc, char** argv)
     using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
     using Coordinate = Dune::FieldVector<double, dim>;
     using Problem = DGProblem<Coordinate>;
-    using DG = DG<Grid, GridView, Mapper, Problem>;
+    using DG = DG<GridView, Mapper, Problem>;
 
     //get parameters from file "parameterDG.ini"
     Dune::ParameterTree pt;
@@ -83,8 +83,7 @@ int main(int argc, char** argv)
         "parameterDG.ini (use 'poisson' or TODO)");
     }
 
-    //DG<GridView, Mapper, Problem> dg(gridView, mapper, *problem);
-    DG dg(grid, gridView, mapper, *problem);
+    DG dg(gridView, mapper, *problem);
 
     dg(mu);
 
@@ -92,9 +91,6 @@ int main(int argc, char** argv)
     const double error = dg.computeL2error();
     const double timeTotal = ClockHelper::elapsedRuntime(timeInitial);
 
-    //write error, runtime and number of grid elements to file
-    std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc);
-
     //determine maximum edge length
     double largestEdgeLength = 0.;
     for (const auto& edge : edges(gridView))
@@ -102,9 +98,12 @@ int main(int argc, char** argv)
       largestEdgeLength = std::max(largestEdgeLength, edge.geometry().volume());
     }
 
+    //write error, runtime and maximum edge length
+    std::ofstream errorFile("convergenceData", std::ios::out | std::ios::trunc);
+
     if (errorFile.is_open())
     {
-      errorFile << error << "\n" << timeTotal << "\n" << largestEdgeLength;//gridView.size(0);
+      errorFile << error << "\n" << timeTotal << "\n" << largestEdgeLength;
       errorFile.close();
     }
     else
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 47c2e83..2ed9097 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -44,7 +44,7 @@ int main(int argc, char** argv)
       typename Dune::MultipleCodimMultipleGeomTypeMapper<IGridView>;
     using Coordinate = Dune::FieldVector<double, dim>;
     using Problem = CoupledDGProblem<Coordinate>;
-    using MMDG = MMDG<Grid, GridView, Mapper, IGridView, IMapper, Problem>;
+    using MMDG = MMDG<GridView, Mapper, IGridView, IMapper, Problem>;
 
     //get parameters from file "parameterDG.ini"
     Dune::ParameterTree pt;
@@ -92,7 +92,7 @@ int main(int argc, char** argv)
     const Mapper mapper(gridView, Dune::mcmgElementLayout());
     const IMapper iMapper(iGridView, Dune::mcmgElementLayout());
 
-    MMDG mmdg(grid, gridView, mapper, iGridView, iMapper, *problem);
+    MMDG mmdg(gridView, mapper, iGridView, iMapper, *problem);
     mmdg(mu, xi);
 
     //print total run time
-- 
GitLab