write output to new directory structure

dune/tectonic/io/io-handler.hh

 #ifndef SRC_IO_HANDLER_HH
 #define SRC_IO_HANDLER_HH
 
+#include <filesystem>
 
 #include <dune/common/parametertree.hh>
 
@@ -82,6 +83,11 @@ class IOHandler {
           for (auto &&v : vertices(hostLeafView))
             vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
         }
+
+        if (writeVTK_) {
+            if (!std::filesystem::is_directory("iterates/"))
+                std::filesystem::create_directory("iterates");
+        }
     }
 
     template <class ProgramState, class GlobalFriction>
@@ -135,7 +141,7 @@ class IOHandler {
 
         }
 
-        const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases_, vertexBases_, "../debug_print/vtk/");
+        const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases_, vertexBases_, "iterates/");
         vtkWriter.write(programState.timeStep, programState.u, programState.v,
                         programState.alpha, stress);
     }

dune/tectonic/time-stepping/adaptivetimestepper.cc

@@ -78,7 +78,7 @@ int AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::coarsen(
 
       setDeformation(current_);
       auto C_Step = Step(stepBase_, current_, currentNBodyAssembler, relativeTime_, 2 * relativeTau_, iterationRegister_);
-      C_Step.run(Step::Mode::newThread); //newThread
+      C_Step.run(Step::Mode::sameThread); //newThread
 
       //updateReductionFactors(reductionFactors);
       //std::cout << "AdaptiveTimeStepper C computed!" << std::endl << std::endl;
@@ -92,7 +92,7 @@ int AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::coarsen(
       //auto R2_linearSolver = makeLinearSolver();
       auto&& nBodyAssembler = step(currentNBodyAssembler);
       auto R2_Step = Step(stepBase_, R1_.updaters, nBodyAssembler, relativeTime_ + relativeTau_, relativeTau_, iterationRegister_);
-      R2_Step.run(Step::Mode::newThread); //newThread
+      R2_Step.run(Step::Mode::sameThread); //newThread
 
       //updateReductionFactors(reductionFactors);
       //std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
@@ -195,7 +195,7 @@ int AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::determin
     setDeformation(current_);
     //auto C_linearSolver = makeLinearSolver();
     auto C_Step = Step(stepBase_, current_, currentNBodyAssembler, relativeTime_, 2 * relativeTau_, iterationRegister_);
-    C_Step.run(Step::Mode::newThread); // newThread
+    C_Step.run(Step::Mode::sameThread); // newThread
     //updateReductionFactors(reductionFactors);
     //std::cout << "AdaptiveTimeStepper C computed!" << std::endl << std::endl;
 
@@ -203,7 +203,7 @@ int AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::determin
     //auto R2_linearSolver = makeLinearSolver();
     auto&& nBodyAssembler = step(currentNBodyAssembler);
     auto R2_Step = Step(stepBase_, R1_.updaters, nBodyAssembler, relativeTime_ + relativeTau_, relativeTau_, iterationRegister_);
-    R2_Step.run(Step::Mode::newThread); //newThread
+    R2_Step.run(Step::Mode::sameThread); //newThread
 
     //updateReductionFactors(reductionFactors);
     //std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
@@ -221,7 +221,7 @@ int AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::determin
         setDeformation(current_);
         //auto F1_linearSolver = makeLinearSolver();
         auto F1_Step = Step(stepBase_, current_, currentNBodyAssembler, relativeTime_, relativeTau_ / 2.0, iterationRegister_);
-        F1_Step.run(Step::Mode::newThread); //newThread
+        F1_Step.run(Step::Mode::sameThread); //newThread
         //updateReductionFactors(reductionFactors);
         //std::cout << "AdaptiveTimeStepper F1 computed!" << std::endl << std::endl;
 

src/strikeslip/strikeslip-2D.cfg

 # -*- mode:conf -*-
 [body0]
-smallestDiameter = 0.02 # 2e-3 [m]
+smallestDiameter = 0.006 # 2e-3 [m]
 
 [body1]
-smallestDiameter = 0.02  # 2e-3 [m]
+smallestDiameter = 0.006  # 2e-3 [m]
 
 [timeSteps]
 refinementTolerance = 1e-3# 1e-5

src/strikeslip/strikeslip.cc

@@ -14,6 +14,7 @@
 #include <iostream>
 #include <iomanip>
 #include <memory>
+#include <filesystem>
 
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/exceptions.hh>
@@ -88,13 +89,14 @@
 #include <dune/tectonic/utils/reductionfactors.hh>
 std::vector<std::vector<double>> allReductionFactors;
 
-size_t const dims = MY_DIM;
+const size_t dims = MY_DIM;
+const std::string sourcePath = "/home/joscha/software/dune/dune-tectonic/src/strikeslip/";
 
 Dune::ParameterTree getParameters(int argc, char *argv[]) {
   Dune::ParameterTree parset;
-  Dune::ParameterTreeParser::readINITree("/home/joscha/software/dune/dune-tectonic/src/strikeslip/strikeslip.cfg", parset);
+  Dune::ParameterTreeParser::readINITree(sourcePath + "strikeslip.cfg", parset);
   Dune::ParameterTreeParser::readINITree(
-      Dune::Fufem::formatString("/home/joscha/software/dune/dune-tectonic/src/strikeslip/strikeslip-%dD.cfg", dims), parset);
+      Dune::Fufem::formatString(sourcePath + "strikeslip-%dD.cfg", dims), parset);
   Dune::ParameterTreeParser::readOptions(argc, argv, parset);
   return parset;
 }
@@ -113,12 +115,22 @@ int main(int argc, char *argv[]) {
         std::cout << argv[0] << std::endl;
     }
 
+    auto const parset = getParameters(argc, argv);
+
+    auto outPath = std::filesystem::current_path();
+    outPath +=  "/output/" + parset.get<std::string>("outPath");
+    if (!std::filesystem::is_directory(outPath))
+        std::filesystem::create_directories(outPath);
+
+    const auto copyOptions = std::filesystem::copy_options::overwrite_existing;
+    std::filesystem::copy(sourcePath + "strikeslip.cfg", outPath, copyOptions);
+    std::filesystem::copy(Dune::Fufem::formatString(sourcePath + "strikeslip-%dD.cfg", dims), outPath, copyOptions);
+    std::filesystem::current_path(outPath);
+
     std::ofstream out("strikeslip.log");
     std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
     std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
 
-    auto const parset = getParameters(argc, argv);
-
     using Assembler = MyAssembler<DefLeafGridView, dims>;
     using field_type = Matrix::field_type;
 
@@ -144,12 +156,12 @@ int main(int argc, char *argv[]) {
         const auto& level = *contactNetwork.level(i);
 
         for (size_t j=0; j<level.nBodies(); j++) {
-            writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+            //writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
         }
     }
 
     for (size_t i=0; i<bodyCount; i++) {
-        writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
+        //writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
     }
 
     // ----------------------------

src/strikeslip/strikeslip.cfg

+outPath = test  # output written to ./output/strikeslip/outPath
+
 # -*- mode:conf -*-
 gravity         = 9.81     # [m/s^2]
 
@@ -73,7 +75,7 @@ relativeTau = 2e-2 # 1e-6
 
 [timeSteps]
 scheme = newmark
-timeSteps = 5
+timeSteps = 15
 
 [u0.solver]
 maximumIterations = 100