diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 3e9b4e47b8afa565d5fe714ebab7d096dee9d3ed..20c824c7be01912ea2eab74cf9f0afdc746bf775 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -5,6 +5,7 @@
 
 #include <dune/common/parametertree.hh>
 // #include <dune/common/float_cmp.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
 #include <dune/istl/matrixindexset.hh>
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/gridfunctions/gridviewfunction.hh> 
@@ -1208,7 +1209,10 @@ public:
       std::string vtkOutputName = parameterSet_.get("outputPath", "../../outputs") + "/" + baseName; 
       std::cout << "vtkOutputName:" << vtkOutputName << std::endl;
 
-      Dune::VTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView());
+      int subsamplingRefinement = parameterSet_.get<int>("subsamplingRefinement", 2);
+      Dune::SubsamplingVTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView(), Dune::refinementLevels(subsamplingRefinement));
+      // Dune::VTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView());
+
       vtkWriter.addVertexData(
           Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_1_),
           Dune::VTK::FieldInfo("Corrector phi_1 level"+ std::to_string(level) , Dune::VTK::FieldInfo::Type::vector, dim));
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 193cb5c867aae89de87adbe4acad133b50e5116b..2cc7fcfea294054b0e38bdb29710636b36c564d6 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -2,6 +2,7 @@
 #define DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
 
 #include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
 #include <dune/grid/yaspgrid.hh>
 #include <dune/functions/gridfunctions/gridviewfunction.hh>
 #include <dune/functions/gridfunctions/gridviewentityset.hh>
@@ -597,35 +598,23 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
   ///////////////////////////////
   // -----------------------------------------------------------------
   // --- write material (grid)functions to VTK
-  void writeVTKMaterialFunctions(std::array<int, dim> nElements, const int level)
+  // void writeVTKMaterialFunctions(std::array<int, dim> nElements, const int level)
+  void writeVTKMaterialFunctions(const int level)
   {
         std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
-        using VTKGridType = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<double, dim> >;
-        // VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
-        VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},nElements);
-        
-        using GridViewVTK = typename VTKGridType::LeafGridView;
-        const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
-        
-        auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
-        auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
 
+        int subsamplingRefinement = parameterSet_.get<int>("subsamplingRefinement", 2);
+
+        auto scalarP0FeBasis = makeBasis(gridView_,lagrange<0>());
         std::vector<double> indicatorFunction_CoeffP0;
         Dune::Functions::interpolate(scalarP0FeBasis, indicatorFunction_CoeffP0, indicatorFunction_);
         auto indicatorFunction_P0 = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, indicatorFunction_CoeffP0);
         
-        std::vector<double> indicatorFunction_CoeffP1;
-        Dune::Functions::interpolate(scalarP1FeBasis, indicatorFunction_CoeffP1, indicatorFunction_);
-        auto indicatorFunction_P1 = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, indicatorFunction_CoeffP1);
-        
-        Dune::VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
+        Dune::SubsamplingVTKWriter<GridView> MaterialVtkWriter(gridView_, Dune::refinementLevels(subsamplingRefinement));
         
         MaterialVtkWriter.addCellData(
             indicatorFunction_P0,
-            Dune::VTK::FieldInfo("indicatorFunction_P0", Dune::VTK::FieldInfo::Type::scalar, 1));
-        MaterialVtkWriter.addVertexData(
-            indicatorFunction_P1,
-            Dune::VTK::FieldInfo("indicatorFunction_P1", Dune::VTK::FieldInfo::Type::scalar, 1));
+            Dune::VTK::FieldInfo("materialIndicatorFunction_P0", Dune::VTK::FieldInfo::Type::scalar, 1));
 
         std::string baseName = parameterSet_.get("baseName", "CellProblem-result");
         MaterialVtkWriter.write(outputPath + "/" + baseName + "_MaterialFunction-level"+ std::to_string(level) );
diff --git a/experiment/material_orthotropic.py b/experiment/material_orthotropic.py
index 0eb37bb3bb5e98d0ae046705711e172434e01839..be1611ca1492f21bc192282ad1f747be3ffed52c 100644
--- a/experiment/material_orthotropic.py
+++ b/experiment/material_orthotropic.py
@@ -111,6 +111,9 @@ parameterSet.print_debug = 0  #(default=false)
 # --- Write Correctos to VTK-File:  
 parameterSet.write_VTK = 1
 
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 2
+
 # --- (Optional output) L2Error, integral mean: 
 #parameterSet.write_L2Error = 1
 #parameterSet.write_IntegralMean = 1      
diff --git a/experiment/parametrized_laminate.py b/experiment/parametrized_laminate.py
index 0fa4bbefefd917476a58fb1df8f0e8ba0cc1077e..0d82204c61ee7b7a1440b27064588bff42c4a9f6 100644
--- a/experiment/parametrized_laminate.py
+++ b/experiment/parametrized_laminate.py
@@ -114,6 +114,9 @@ parameterSet.print_debug = 0  #(default=false)
 # --- Write Correctos to VTK-File:  
 parameterSet.write_VTK = 1
 
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 2
+
 # --- (Optional output) L2Error, integral mean: 
 #parameterSet.write_L2Error = 1
 #parameterSet.write_IntegralMean = 1      
diff --git a/experiment/rotation-test/isotrop_orthotrop_rotation.py b/experiment/rotation-test/isotrop_orthotrop_rotation.py
index 2afba12d28ebabe250fd53f1a0914f7db879f5d3..6a5b26f68f92ebc29e864cd992904a4d1c2f4bf5 100644
--- a/experiment/rotation-test/isotrop_orthotrop_rotation.py
+++ b/experiment/rotation-test/isotrop_orthotrop_rotation.py
@@ -138,6 +138,9 @@ parameterSet.print_debug = 0  #(default=false)
 # --- Write Correctos to VTK-File:  
 parameterSet.write_VTK = 1
 
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 2
+
 # --- (Optional output) L2Error, integral mean: 
 #parameterSet.write_L2Error = 1
 #parameterSet.write_IntegralMean = 1      
diff --git a/experiment/wood-bilayer-variant/wood_inclusion.py b/experiment/wood-bilayer-variant/wood_inclusion.py
index 15fec129064376fc81557f45c98286a3f0b3a0a3..074270fd04eff563902c7e698f95186dc2a78ae1 100644
--- a/experiment/wood-bilayer-variant/wood_inclusion.py
+++ b/experiment/wood-bilayer-variant/wood_inclusion.py
@@ -212,6 +212,9 @@ parameterSet.print_debug = 0  #(default=false)
 # --- Write Correctos to VTK-File:  
 parameterSet.write_VTK = 1
 
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 2
+
 # --- (Optional output) L2Error, integral mean: 
 #parameterSet.write_L2Error = 1
 #parameterSet.write_IntegralMean = 1      
diff --git a/experiment/wood-bilayer-variant/wood_upper_laminated.py b/experiment/wood-bilayer-variant/wood_upper_laminated.py
index a46aef6084e89e9e091dd47a4b4a0d195048f8b0..987afb71a73058110bc1200930a3516f4e124d97 100644
--- a/experiment/wood-bilayer-variant/wood_upper_laminated.py
+++ b/experiment/wood-bilayer-variant/wood_upper_laminated.py
@@ -220,6 +220,9 @@ parameterSet.print_debug = 0  #(default=false)
 # --- Write Correctos to VTK-File:  
 parameterSet.write_VTK = 1
 
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 2
+
 # --- (Optional output) L2Error, integral mean: 
 #parameterSet.write_L2Error = 1
 #parameterSet.write_IntegralMean = 1      
diff --git a/experiment/wood-bilayer/wood_european_beech.py b/experiment/wood-bilayer/wood_european_beech.py
index 3f321222bfff470fd6ef46e789b23d28f841fc3f..c8c326abc408137f9cc750e17d2cab2ba5faa154 100644
--- a/experiment/wood-bilayer/wood_european_beech.py
+++ b/experiment/wood-bilayer/wood_european_beech.py
@@ -208,6 +208,9 @@ parameterSet.print_debug = 0  #(default=false)
 # --- Write Correctos to VTK-File:  
 parameterSet.write_VTK = 1
 
+# The grid can be refined several times for a higher resolution in the VTK-file.
+parameterSet.subsamplingRefinement = 2
+
 # --- (Optional output) L2Error, integral mean: 
 #parameterSet.write_L2Error = 1
 #parameterSet.write_IntegralMean = 1      
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index a86425e33e501a316ba4f7e76a70e941f7946fcf..f28970223029293d4c8a9d16115bbf11cbea3778 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -281,7 +281,7 @@ int main(int argc, char *argv[])
     //--- write material indicator function to VTK
     if (parameterSet.get<bool>("write_materialFunctions", false))
     {
-        material.writeVTKMaterialFunctions(nElements,level);
+        material.writeVTKMaterialFunctions(level);
     }
 
     //--- TEST:: Compute Qeff without using the orthogonality (75)...