From 6841ef4ed67361d0e7732bd310a51a5e00cf3697 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Tue, 10 Oct 2023 15:28:36 +0200
Subject: [PATCH] update VTK-writers and add option for subsamplingRefinement

---
 dune/microstructure/CorrectorComputer.hh      |  6 ++++-
 dune/microstructure/prestrainedMaterial.hh    | 27 ++++++-------------
 experiment/material_orthotropic.py            |  3 +++
 experiment/parametrized_laminate.py           |  3 +++
 .../isotrop_orthotrop_rotation.py             |  3 +++
 .../wood-bilayer-variant/wood_inclusion.py    |  3 +++
 .../wood_upper_laminated.py                   |  3 +++
 .../wood-bilayer/wood_european_beech.py       |  3 +++
 src/Cell-Problem.cc                           |  2 +-
 9 files changed, 32 insertions(+), 21 deletions(-)

diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 3e9b4e47..20c824c7 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 193cb5c8..2cc7fcfe 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 0eb37bb3..be1611ca 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 0fa4bbef..0d82204c 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 2afba12d..6a5b26f6 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 15fec129..074270fd 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 a46aef60..987afb71 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 3f321222..c8c326ab 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 a86425e3..f2897022 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)... 
-- 
GitLab