From 23a60283115327e7eab3248fc84ce0603040deb8 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Tue, 30 Aug 2022 01:27:17 +0200
Subject: [PATCH] Add VTK-writing options

---
 dune/microstructure/prestrainedMaterial.hh | 123 +++++++++++++++++++--
 src/Cell-Problem-New.cc                    |  99 +----------------
 2 files changed, 118 insertions(+), 104 deletions(-)

diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 3a274edc..afd2dca5 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -245,16 +245,6 @@ public:
 
   // }
 
-
-  // -----------------------------------------------------------------
-  // --- write material (grid)functions to VTK
-  void write_materialFunctions()
-  {
-
-
-	return;
-
-  };
     ///////////////////////////////
     // Getter
     ///////////////////////////////
@@ -276,7 +266,120 @@ public:
     //getLameParameters() .. Throw Exception if isotropic_ = false! 
 
 
+  ///////////////////////////////
+  // VTK - Writer
+  ///////////////////////////////
+  // -----------------------------------------------------------------
+  // --- write material (grid)functions to VTK
+  void writeVTKMaterialFunctions(std::array<int, dim> nElements, const int level)
+  {
+        std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
+        using VTKGridType = YaspGrid<dim, 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 = VTKGridType::LeafGridView;
+        const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+        
+        auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+        auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
+
+        std::vector<double> indicatorFunction_CoeffP0;
+        Functions::interpolate(scalarP0FeBasis, indicatorFunction_CoeffP0, indicatorFunction_);
+        auto indicatorFunction_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, indicatorFunction_CoeffP0);
+        
+        std::vector<double> indicatorFunction_CoeffP1;
+        Functions::interpolate(scalarP1FeBasis, indicatorFunction_CoeffP1, indicatorFunction_);
+        auto indicatorFunction_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, indicatorFunction_CoeffP1);
+        
+        VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
+        
+        MaterialVtkWriter.addVertexData(
+            indicatorFunction_P0,
+            VTK::FieldInfo("indicatorFunction_P0", VTK::FieldInfo::Type::scalar, 1));    
+        MaterialVtkWriter.addCellData(
+            indicatorFunction_P1,
+            VTK::FieldInfo("indicatorFunction_P1", VTK::FieldInfo::Type::scalar, 1));    
+        
+        MaterialVtkWriter.write(outputPath + "/MaterialFunctions-level"+ std::to_string(level) );
+        std::cout << "wrote data to file:" + outputPath +"/MaterialFunctions-level" + std::to_string(level) << std::endl;  
+	      return;
+  };
 
+  // -----------------------------------------------------------------
+  // --- write prestrain (grid)functions to VTK
+  // void writeVTKPrestainFunctions(std::array<int, dim> nElements, const int level)
+  // {
+  //       std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
+  //       using VTKGridType = YaspGrid<dim, 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 = VTKGridType::LeafGridView;
+  //       const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+        
+  //       auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+  //       auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
+
+  //       std::vector<double> prestrainFunction_CoeffP0;
+  //       Functions::interpolate(scalarP0FeBasis, prestrainFunction_CoeffP0, prestrain_);
+  //       auto prestrainFunction_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, prestrainFunction_CoeffP0);
+        
+  //       std::vector<double> prestrainFunction_CoeffP1;
+  //       Functions::interpolate(scalarP1FeBasis, prestrainFunction_CoeffP1, prestrain_);
+  //       auto prestrainFunction_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, prestrainFunction_CoeffP1);
+        
+  //       VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
+        
+  //       MaterialVtkWriter.addVertexData(
+  //           prestrainFunction_P0,
+  //           VTK::FieldInfo("prestrainFunction_P0", VTK::FieldInfo::Type::scalar, 1));    
+  //       MaterialVtkWriter.addCellData(
+  //           prestrainFunction_P1,
+  //           VTK::FieldInfo("prestrainFunction_P1", VTK::FieldInfo::Type::scalar, 1));    
+
+  //       MaterialVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
+  //       std::cout << "wrote data to file:" + outputPath + "/PrestrainFunctions-level" + std::to_string(level) << std::endl;  
+	//       return;
+  // };
+
+  // void writeVTKPrestainFunctions(std::array<int, dim> nElements, const int level)
+  // {
+//     using VTKGridType = YaspGrid<dim, 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},{40,40,40});
+//     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 = VTKGridType::LeafGridView;
+//     const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+   
+//     FTKfillerContainer<dim> VTKFiller;
+//     VTKFiller.vtkPrestrainNorm(gridView_VTK, B_Term, "PrestrainBNorm");
+    
+//     // WORKS Too 
+//     VTKFiller.vtkProblemCell(gridView_VTK, B_Term, muLocal,"VTKProblemCell");;
+    
+    
+//     // TEST 
+//     auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+//     auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
+
+//     std::vector<double> B_CoeffP0;
+//     Functions::interpolate(scalarP0FeBasis, B_CoeffP0, B_Term);
+//     auto B_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, B_CoeffP0);
+        
+
+        
+//     VTKWriter<GridView> PrestrainVtkWriter(gridView_VTK);
+         
+//     PrestrainVtkWriter.addCellData(
+//             B_DGBF_P0,
+//             VTK::FieldInfo("B_P0", VTK::FieldInfo::Type::scalar, 1));     
+        
+//     PrestrainVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
+//     std::cout << "wrote data to file:" + outputPath +"/PrestrainFunctions-level" + std::to_string(level) << std::endl; 
+
+//    // };
+  
 
 }; // end class
 
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index 22cd8d03..ad05b62e 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -570,11 +570,6 @@ int main(int argc, char *argv[])
     auto correctorComputer = CorrectorComputer(Basis_CE, material_, muTerm, lambdaTerm, gamma, log, parameterSet);
     correctorComputer.solve();
 
-
-
-//////////////////////////////////////////////////
-   
-
     //--- check Correctors (options):
     if(parameterSet.get<bool>("write_L2Error", false))
          correctorComputer.computeNorms();
@@ -595,6 +590,11 @@ int main(int argc, char *argv[])
     auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material_);
     effectiveQuantitiesComputer.computeEffectiveQuantities();
 
+    //--- write material indicator function to VTK
+    if (write_materialFunctions)
+    {
+        material_.writeVTKMaterialFunctions(nElements,level);
+    }
 
 
 
@@ -652,95 +652,6 @@ int main(int argc, char *argv[])
 
 
 
-
-   if (write_materialFunctions)
-   {
-        using VTKGridType = YaspGrid<dim, 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},{40,40,40});
-        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 = VTKGridType::LeafGridView;
-        const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
-        
-        auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
-        auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
-
-        std::vector<double> mu_CoeffP0;
-        Functions::interpolate(scalarP0FeBasis, mu_CoeffP0, muTerm);
-        auto mu_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, mu_CoeffP0);
-        
-        std::vector<double> mu_CoeffP1;
-        Functions::interpolate(scalarP1FeBasis, mu_CoeffP1, muTerm);
-        auto mu_DGBF_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, mu_CoeffP1);
-        
-        
-        std::vector<double> lambda_CoeffP0;
-        Functions::interpolate(scalarP0FeBasis, lambda_CoeffP0, lambdaTerm);
-        auto lambda_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, lambda_CoeffP0);
-        
-        std::vector<double> lambda_CoeffP1;
-        Functions::interpolate(scalarP1FeBasis, lambda_CoeffP1, lambdaTerm);
-        auto lambda_DGBF_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, lambda_CoeffP1);
-        
-        VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
-        
-        MaterialVtkWriter.addVertexData(
-            mu_DGBF_P1,
-            VTK::FieldInfo("mu_P1", VTK::FieldInfo::Type::scalar, 1));    
-        MaterialVtkWriter.addCellData(
-            mu_DGBF_P0,
-            VTK::FieldInfo("mu_P0", VTK::FieldInfo::Type::scalar, 1));    
-        MaterialVtkWriter.addVertexData(
-            lambda_DGBF_P1,
-            VTK::FieldInfo("lambda_P1", VTK::FieldInfo::Type::scalar, 1));    
-        MaterialVtkWriter.addCellData(
-            lambda_DGBF_P0,
-            VTK::FieldInfo("lambda_P0", VTK::FieldInfo::Type::scalar, 1));    
-        
-        MaterialVtkWriter.write(outputPath + "/MaterialFunctions-level"+ std::to_string(level) );
-        std::cout << "wrote data to file:" + outputPath +"/MaterialFunctions-level" + std::to_string(level) << std::endl;  
-        
-   }
-  
-  
-//    if (write_prestrainFunctions)
-//    {
-//     using VTKGridType = YaspGrid<dim, 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},{40,40,40});
-//     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 = VTKGridType::LeafGridView;
-//     const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
-   
-//     FTKfillerContainer<dim> VTKFiller;
-//     VTKFiller.vtkPrestrainNorm(gridView_VTK, B_Term, "PrestrainBNorm");
-    
-//     // WORKS Too 
-//     VTKFiller.vtkProblemCell(gridView_VTK, B_Term, muLocal,"VTKProblemCell");;
-    
-    
-//     // TEST 
-//     auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
-//     auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
-
-//     std::vector<double> B_CoeffP0;
-//     Functions::interpolate(scalarP0FeBasis, B_CoeffP0, B_Term);
-//     auto B_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, B_CoeffP0);
-        
-
-        
-//     VTKWriter<GridView> PrestrainVtkWriter(gridView_VTK);
-         
-//     PrestrainVtkWriter.addCellData(
-//             B_DGBF_P0,
-//             VTK::FieldInfo("B_P0", VTK::FieldInfo::Type::scalar, 1));     
-        
-//     PrestrainVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
-//     std::cout << "wrote data to file:" + outputPath +"/PrestrainFunctions-level" + std::to_string(level) << std::endl; 
-
-//    }
-  
   
 
 
-- 
GitLab