src/uniform-grid-writer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cmath>
+#include <exception>
+#include <iostream>
+
+#include <dune/common/exceptions.hh>
+// #include <dune/common/parametertree.hh>
+// #include <dune/common/parametertreeparser.hh>
+
+#include "assemblers.hh"
+#include "diameter.hh"
+#include "gridselector.hh"
+#include "one-body-problem-data/mygrid.hh"
+#include "vtk.hh"
+
+size_t const dims = MY_DIM;
+size_t const refinements = 5; // FIXME?
+
+int main(int argc, char *argv[]) {
+  try {
+    // Dune::ParameterTree parset;
+    // Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+
+    using GridView = Grid::LeafGridView;
+    using MyAssembler = MyAssembler<GridView, dims>;
+
+    GridConstructor<Grid> gridConstructor;
+    auto grid = gridConstructor.getGrid();
+
+    // refine uniformly!
+    for (size_t refinement = 0; refinement < refinements; ++refinement)
+      grid->globalRefine(1);
+
+    double minDiameter = std::numeric_limits<double>::infinity();
+    double maxDiameter = 0.0;
+    for (auto &&e : elements(grid->leafGridView())) {
+      auto const geometry = e.geometry();
+      auto const diam = diameter(geometry);
+      minDiameter = std::min(minDiameter, diam);
+      maxDiameter = std::max(maxDiameter, diam);
+    }
+    std::cout << "min diameter: " << minDiameter << std::endl;
+    std::cout << "max diameter: " << maxDiameter << std::endl;
+
+    auto const leafView = grid->leafGridView();
+    auto const leafVertexCount = leafView.size(dims);
+
+    std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
+
+    MyAssembler const myAssembler(leafView);
+    MyVTKWriter<typename MyAssembler::VertexBasis,
+                typename MyAssembler::CellBasis> const
+        vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+    vtkWriter.writeGrid();
+  } catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  }
+}

src/vtk.cc

@@ -32,10 +32,10 @@ void MyVTKWriter<VertexBasis, CellBasis>::write(
           vertexBasis, v, "velocity");
   writer.addVertexData(velocityPointer);
 
-  auto const logStatePointer =
+  auto const AlphaPointer =
       std::make_shared<VTKBasisGridFunction<VertexBasis, ScalarVector> const>(
-          vertexBasis, alpha, "logState");
-  writer.addVertexData(logStatePointer);
+          vertexBasis, alpha, "Alpha");
+  writer.addVertexData(AlphaPointer);
 
   auto const stressPointer =
       std::make_shared<VTKBasisGridFunction<CellBasis, ScalarVector> const>(
@@ -43,7 +43,16 @@ void MyVTKWriter<VertexBasis, CellBasis>::write(
   writer.addCellData(stressPointer);
 
   std::string const filename = prefix + std::to_string(record);
-  writer.write(filename.c_str());
+  writer.write(filename.c_str(), Dune::VTK::appendedraw);
+}
+
+template <class VertexBasis, class CellBasis>
+void MyVTKWriter<VertexBasis, CellBasis>::writeGrid() const {
+  Dune::VTKWriter<typename VertexBasis::GridView> writer(
+      vertexBasis.getGridView());
+
+  std::string const filename = prefix + "_grid";
+  writer.write(filename.c_str(), Dune::VTK::appendedraw);
 }
 
 #include "vtk_tmpl.cc"

src/vtk.hh

@@ -15,5 +15,7 @@ template <class VertexBasis, class CellBasis> class MyVTKWriter {
   template <class Vector, class ScalarVector>
   void write(size_t record, Vector const &u, Vector const &v,
              ScalarVector const &alpha, ScalarVector const &stress) const;
+
+  void writeGrid() const;
 };
 #endif

src/vtk_tmpl.cc

-#ifndef DIM
-#error DIM unset
+#ifndef MY_DIM
+#error MY_DIM unset
 #endif
 
 #include "explicitgrid.hh"