Externalise VTK writing

src/Makefile.am

@@ -23,7 +23,7 @@ run-one-body-sample-gdb: one-body-sample
 	libtool --mode execute gdb ./one-body-sample
 
 one_body_sample_SOURCES = \
-	one-body-sample.cc assemblers.cc LambertW.cc compute_state.cc
+	one-body-sample.cc assemblers.cc LambertW.cc compute_state.cc vtk.cc
 one_body_sample_CPPFLAGS = \
 	$(AM_CPPFLAGS) -Dsrcdir=\"$(srcdir)\"
 one_body_sample_LDADD = \

src/one-body-sample.cc

@@ -34,7 +34,6 @@
 #include <dune/common/stdstreams.hh>
 #include <dune/common/timer.hh>
 #include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/io/file/vtk/vtkwriter.hh>
 #include <dune/grid/yaspgrid.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
@@ -46,7 +45,6 @@
 #include <dune/fufem/assemblers/transferoperatorassembler.hh>
 #include <dune/fufem/dunepython.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functions/vtkbasisgridfunction.hh>
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 #include <dune/fufem/sharedpointermap.hh>
@@ -66,6 +64,7 @@
 
 #include "assemblers.hh"
 #include "compute_state.hh"
+#include "vtk.hh"
 
 int const dim = 2;
 
@@ -416,6 +415,7 @@ int main(int argc, char *argv[]) {
 
       // Compute von Mises stress and write everything to a file
       if (parset.get<bool>("writeVTK")) {
+        // Compute von Mises stress
         auto const displacement =
             Dune::make_shared<BasisGridFunction<P1Basis, VectorType> const>(
                 p1Basis, u4);
@@ -424,22 +424,9 @@ int main(int argc, char *argv[]) {
         FunctionalAssembler<P0Basis>(p0Basis)
             .assemble(localStressAssembler, vonMisesStress, true);
 
-        auto const displacement_ptr =
-            Dune::make_shared<VTKBasisGridFunction<P1Basis, VectorType> const>(
-                p1Basis, u4, "displacement");
-        auto const state_ptr = Dune::make_shared<
-            VTKBasisGridFunction<P1Basis, SingletonVectorType> const>(
-            p1Basis, *s4_new, "state");
-        auto const vonmises_ptr = Dune::make_shared<
-            VTKBasisGridFunction<P0Basis, SingletonVectorType> const>(
-            p0Basis, vonMisesStress, "stress");
-        Dune::VTKWriter<GridView> writer(leafView);
-        writer.addVertexData(state_ptr);
-        writer.addVertexData(displacement_ptr);
-        writer.addCellData(vonmises_ptr);
-
-        std::string const filename((boost::format("obs%d") % run).str());
-        writer.write(filename.c_str());
+        writeVtk<P1Basis, P0Basis, VectorType, SingletonVectorType, GridView>(
+            p1Basis, u4, *s4_new, p0Basis, vonMisesStress, leafView,
+            (boost::format("obs%d") % run).str());
       }
 
       if (parset.get<bool>("solver.gs.use")) {

src/vtk.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "vtk.hh"
+
+template <class VertexBasis, class CellBasis, class VectorType,
+          class SingletonVectorType, class GridView>
+void writeVtk(VertexBasis vertexBasis, VectorType displacement,
+              SingletonVectorType state, CellBasis cellBasis,
+              SingletonVectorType stress, GridView gridView,
+              std::string filename) {
+  auto const displacement_ptr =
+      Dune::make_shared<VTKBasisGridFunction<VertexBasis, VectorType> const>(
+          vertexBasis, displacement, "displacement");
+  auto const state_ptr = Dune::make_shared<
+      VTKBasisGridFunction<VertexBasis, SingletonVectorType> const>(
+      vertexBasis, state, "state");
+  auto const vonmises_ptr = Dune::make_shared<
+      VTKBasisGridFunction<CellBasis, SingletonVectorType> const>(
+      cellBasis, stress, "stress");
+
+  Dune::VTKWriter<GridView> writer(gridView);
+  writer.addVertexData(state_ptr);
+  writer.addVertexData(displacement_ptr);
+  writer.addCellData(vonmises_ptr);
+
+  writer.write(filename.c_str());
+}
+
+#include "vtk_tmpl.cc"

src/vtk.hh

+#ifndef VTK_HH
+#define VTK_HH
+
+#include <dune/common/shared_ptr.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+#include <dune/fufem/functions/vtkbasisgridfunction.hh>
+
+template <class VertexBasis, class CellBasis, class VectorType,
+          class SingletonVectorType, class GridView>
+void writeVtk(VertexBasis vertexBasis, VectorType displacement,
+              SingletonVectorType state, CellBasis cellBasis,
+              SingletonVectorType stress, GridView gridView,
+              std::string filename);
+
+#endif

src/vtk_tmpl.cc

+#include "vtk.hh"
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/functionspacebases/p0basis.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+
+typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
+
+// {{{ 2D
+typedef Dune::FieldVector<double, 2> SmallVector2;
+typedef Dune::FieldMatrix<double, 2, 2> SmallMatrix2;
+typedef Dune::BlockVector<SmallVector2> VectorType2;
+
+typedef Dune::YaspGrid<2> GridType2;
+typedef GridType2::LeafGridView GridView2;
+typedef P1NodalBasis<GridView2, double> P1Basis2;
+typedef P0Basis<GridView2, double> P0Basis2;
+
+template void writeVtk<P1Basis2, P0Basis2, VectorType2, SingletonVectorType,
+                       GridView2>(P1Basis2 vertexBasis,
+                                  VectorType2 displacement,
+                                  SingletonVectorType state, P0Basis2 cellBasis,
+                                  SingletonVectorType stress,
+                                  GridView2 gridView, std::string filename);
+// }}}
+
+// {{{ 3D
+typedef Dune::FieldVector<double, 3> SmallVector3;
+typedef Dune::FieldMatrix<double, 3, 3> SmallMatrix3;
+typedef Dune::BlockVector<SmallVector3> VectorType3;
+
+typedef Dune::YaspGrid<3> GridType3;
+typedef GridType3::LeafGridView GridView3;
+typedef P1NodalBasis<GridView3, double> P1Basis3;
+typedef P0Basis<GridView3, double> P0Basis3;
+
+template void writeVtk<P1Basis3, P0Basis3, VectorType3, SingletonVectorType,
+                       GridView3>(P1Basis3 vertexBasis,
+                                  VectorType3 displacement,
+                                  SingletonVectorType state, P0Basis3 cellBasis,
+                                  SingletonVectorType stress,
+                                  GridView3 gridView, std::string filename);
+// }}}