add vtk output to mmdg

dune/mmdg/dg.hh

@@ -38,7 +38,7 @@ public:
   static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);
 
   //constructor
-  DG (const Grid& grid, const GridView& gridView, const Mapper mapper, const Problem& problem) :
+  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),
@@ -72,7 +72,7 @@ public:
   const int dof_;
 
 protected:
-  DG (const Grid& grid, const GridView& gridView, const Mapper mapper,
+  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)
@@ -451,13 +451,13 @@ protected:
   }
 
   //writes the solution to a vtk file
-  void writeVTKoutput () const
+  void writeVTKoutput (const int bulkDOF) 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);
+    Dune::DynamicVector<Scalar> pressure(bulkDOF, 0.0);
+    Dune::DynamicVector<Scalar> exactPressure(bulkDOF, 0.0);
 
     for (const auto& elem : elements(gridView_))
     {
@@ -502,10 +502,17 @@ protected:
     vtkWriter.pwrite("pressureData", "", "data");
   }
 
-
   std::shared_ptr<Matrix> A; //stiffness matrix
   Vector b; //load vector
   Vector d; //solution vector
+
+private:
+  //writes the solution to a vtk file
+  void writeVTKoutput () const
+  {
+    writeVTKoutput(dof_);
+  }
+
 };
 
 #endif

dune/mmdg/mmdg.hh

@@ -9,6 +9,7 @@ class MMDG : public DG<Grid, 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 Scalar = typename Base::Scalar; //NOTE using Base::Scalar
@@ -18,7 +19,9 @@ public:
   using InterfaceBasis = std::array<Coordinate, dim-1>;
   using QRuleInterface = typename Base::QRuleEdge;
   using QRulesInterface = typename Base::QRulesEdge;
-  static constexpr auto iSimplex = Dune::GeometryTypes::simplex(dim-1);
+  using VTKFunction = Dune::VTKFunction<InterfaceGridView>;
+  using P1Function = Dune::NonconformingP1VTKFunction<InterfaceGridView,
+    Dune::DynamicVector<Scalar>>;
 
   //constructor
   MMDG (const Grid& grid, const GridView& gridView, const Mapper& mapper,
@@ -32,6 +35,13 @@ public:
   {
     assembleSLE(mu, xi);
     Base::A->compress();
+
+    Dune::InverseOperatorResult result;
+    Dune::UMFPack<Matrix> solver(*Base::A);
+    solver.apply(Base::d, Base::b, result);
+
+    //write solution to a vtk file
+    writeVTKoutput();
   }
 
   const InterfaceGridView& iGridView_;
@@ -96,8 +106,8 @@ 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());
+      const InterfaceBasis iFrame =
+        interfaceFrame(intersct.centerUnitOuterNormal());
 
       //in the total system of linear equations (SLE) Ad = b,
       //the index iElemIdxSLE refers to the basis function (0, phi_iElem,0)
@@ -558,7 +568,7 @@ private:
   //3d: tau_1 ~ [ normal[1] + normal[2],     tau_2 ~ [ -normal[1],
   //             -normal[0],                            normal[0] + normal[2]
   //             -normal[0]             ],             -normal[1]            ]
-  InterfaceBasis interfaceFrame(const Coordinate& normal)
+  InterfaceBasis interfaceFrame(const Coordinate& normal) const
   {
     InterfaceBasis iBasis;
 
@@ -587,6 +597,71 @@ private:
     return iBasis;
   }
 
+  void writeVTKoutput () const
+  {
+    //degrees of freedom of bulk and interface grids
+    const int bulkDOF = (dim + 1) * Base::gridView_.size(0);
+    const int interfaceDOF = dim * iGridView_.size(0);
+
+    //=== write bulk solution to file ===
+    Base::writeVTKoutput(bulkDOF);
+
+    //=== write interface solution to file ===
+
+    //storage for pressure data, for each interface element we store the
+    //pressure at the corners of the interface element, the VTKFunction will
+    //create the output using a linear interpolation for each interface element
+    Dune::DynamicVector<Scalar> iPressure(interfaceDOF, 0.0);
+    Dune::DynamicVector<Scalar> exactIPressure(interfaceDOF, 0.0);
+
+    for (const auto& iElem : elements(iGridView_))
+    {
+      const int iElemIdxSLE = bulkDOF + dim * iMapper_.index(iElem);
+      const auto& iGeo = iElem.geometry();
+
+      //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());
+
+      for (int k = 0; k < iGeo.corners(); k++)
+      {
+        if (Base::problem_.hasExactSolution())
+        {
+          exactIPressure[iElemIdxSLE + k] =
+            Base::problem_.exactInterfaceSolution(iGeo.corner(k));
+        }
+
+        //contribution of the basis function
+        // phi_iElem,0 (x) = indicator(iElem);
+        //at the kth corner of iElem
+        iPressure[iElemIdxSLE + k] = Base::d[iElemIdxSLE];
+
+        for (int i = 0; i < dim - 1; i++)
+        {
+          //contribution of the basis function
+          // phi_iElem,i (x) = (x * tau[i]) * indicator(iElem);
+          //at the kth corner of iElem
+          iPressure[iElemIdxSLE + k] +=
+            Base::d[iElemIdxSLE + i + 1] * (iGeo.corner(k) * iFrame[i]);
+        }
+      }
+    }
+
+    Dune::VTKWriter<InterfaceGridView>
+      vtkWriter(iGridView_, Dune::VTK::nonconforming);
+    vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
+      new P1Function(iGridView_, iPressure, "interfacePressure")));
+
+    if (Base::problem_.hasExactSolution())
+    {
+      vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
+        new P1Function(iGridView_, exactIPressure, "exactInterfacePressure")));
+    }
+
+    vtkWriter.pwrite("interfacePressureData", "", "data");
+  }
+
 };
 
 #endif

dune/mmdg/problems/coupleddgproblem.hh

@@ -14,6 +14,12 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar>
     using Base = DGProblem<Vector, Scalar>;
     using Matrix = typename Base::Matrix;
 
+    //the exact solution on the interface at position pos
+    virtual Scalar exactInterfaceSolution (const Vector& pos) const
+    {
+      return Scalar(0.0);
+    }
+
     //interface source term at position pos
     virtual Scalar qInterface (const Vector& pos) const
     {
@@ -26,12 +32,6 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar>
       return Scalar(1.0);
     }
 
-    //(tangential) gradient of the aperture at position pos
-    virtual Vector gradAperture (const Vector& pos) const
-    {
-      return Vector(0.0);
-    }
-
     //tangential permeability tensor of the interface at position pos
     virtual Matrix Kparallel (const Vector& pos) const
     {

src/mmdg.cc

@@ -33,7 +33,7 @@ int main(int argc, char** argv)
     static constexpr int dim = GRIDDIM;
 
     using Grid = Dune::MovingMesh<dim>;
-    using GridFactory = Dune::DGFGridFactory<Grid>;
+    using GridFactory = Dune::GmshGridFactory<Grid, /*implicit=*/false>;
     using GridView = typename Grid::LeafGridView;
     using Mapper = typename Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
     using IGrid = typename Grid::InterfaceGrid;
@@ -53,7 +53,7 @@ int main(int argc, char** argv)
     const double xi = std::stod(pt["xi"]); //coupling parameter
 
     //create a grid from .dgf file
-    GridFactory gridFactory( "grids/cube" + std::to_string(dim) + "d.dgf" );
+    GridFactory gridFactory( "grids/sphere" + std::to_string(dim) + "d.msh" );
     const Grid& grid = *gridFactory.grid();
     const GridView& gridView = grid.leafGridView();
     const IGrid& iGrid = grid.interfaceGrid();