From 2cfdb544d37ff0f09de5200982925de9ffc4abdc Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Thu, 4 Feb 2021 15:34:14 +0100
Subject: [PATCH] switch to new io
---
src/foam/foam.cc | 150 +++++++++++------------------------------------
1 file changed, 33 insertions(+), 117 deletions(-)
diff --git a/src/foam/foam.cc b/src/foam/foam.cc
index 46e8bc65..59660df5 100644
--- a/src/foam/foam.cc
+++ b/src/foam/foam.cc
@@ -56,6 +56,7 @@
#include <dune/tectonic/factories/twoblocksfactory.hh>
+#include <dune/tectonic/io/io-handler.hh>
#include <dune/tectonic/io/hdf5-writer.hh>
#include <dune/tectonic/io/hdf5/restart-io.hh>
#include <dune/tectonic/io/vtk.hh>
@@ -88,12 +89,13 @@
std::vector<std::vector<double>> allReductionFactors;
size_t const dims = MY_DIM;
+const std::string sourcePath = "/home/joscha/software/dune/dune-tectonic/src/foam/";
Dune::ParameterTree getParameters(int argc, char *argv[]) {
Dune::ParameterTree parset;
- Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/foam/foam.cfg", parset);
+ Dune::ParameterTreeParser::readINITree(sourcePath + "foam.cfg", parset);
Dune::ParameterTreeParser::readINITree(
- Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/foam/foam-%dD.cfg", dims), parset);
+ Dune::Fufem::formatString(sourcePath + "foam-%dD.cfg", dims), parset);
Dune::ParameterTreeParser::readOptions(argc, argv, parset);
return parset;
}
@@ -112,12 +114,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 + "foam.cfg", outPath, copyOptions);
+ std::filesystem::copy(Dune::Fufem::formatString(sourcePath + "foam-%dD.cfg", dims), outPath, copyOptions);
+ std::filesystem::current_path(outPath);
+
std::ofstream out("foam.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;
@@ -131,8 +143,6 @@ int main(int argc, char *argv[]) {
ContactNetwork& contactNetwork = blocksFactory.contactNetwork();
-
-
const size_t bodyCount = contactNetwork.nBodies();
for (size_t i=0; i<contactNetwork.nLevels(); i++) {
@@ -145,12 +155,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");
}
// ----------------------------
@@ -170,33 +180,13 @@ int main(int argc, char *argv[]) {
using MyProgramState = ProgramState<Vector, ScalarVector>;
MyProgramState programState(nVertices);
- auto const firstRestart = parset.get<size_t>("io.restarts.first");
- auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
- auto const writeRestarts = parset.get<bool>("io.restarts.write");
- auto const writeData = parset.get<bool>("io.data.write");
- bool const handleRestarts = writeRestarts or firstRestart > 0;
-
-
- auto dataFile =
- writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
- auto restartFile = handleRestarts
- ? std::make_unique<HDF5::File>(
- "restarts.h5",
- writeRestarts ? HDF5::Access::READWRITE
- : HDF5::Access::READONLY)
- : nullptr;
+ IOHandler<Assembler, ContactNetwork> ioHandler(parset.sub("io"), contactNetwork);
-
- auto restartIO = handleRestarts
- ? std::make_unique<RestartIO<MyProgramState>>(
- *restartFile, nVertices)
- : nullptr;
-
- if (firstRestart > 0) // automatically adjusts the time and timestep
- restartIO->read(firstRestart, programState);
- else
- programState.setupInitialConditions(parset, contactNetwork);
+ bool restartRead = ioHandler.read(programState);
+ if (!restartRead) {
+ programState.setupInitialConditions(parset, contactNetwork);
+ }
auto& nBodyAssembler = contactNetwork.nBodyAssembler();
for (size_t i=0; i<bodyCount; i++) {
@@ -213,87 +203,9 @@ int main(int argc, char *argv[]) {
auto& globalFriction = contactNetwork.globalFriction();
globalFriction.updateAlpha(programState.alpha);
- using MyVertexBasis = typename Assembler::VertexBasis;
- using MyCellBasis = typename Assembler::CellBasis;
- std::vector<Vector> vertexCoordinates(bodyCount);
- std::vector<const MyVertexBasis* > vertexBases(bodyCount);
- std::vector<const MyCellBasis* > cellBases(bodyCount);
-
- auto& wPatches = blocksFactory.weakPatches();
- std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches(bodyCount);
-
-
- for (size_t i=0; i<bodyCount; i++) {
- const auto& body = contactNetwork.body(i);
- vertexBases[i] = &(body->assembler()->vertexBasis);
- cellBases[i] = &(body->assembler()->cellBasis);
-
- weakPatches[i].resize(1);
- weakPatches[i][0] = wPatches[i].get();
-
- auto& vertexCoords = vertexCoordinates[i];
- vertexCoords.resize(nVertices[i]);
-
- auto hostLeafView = body->grid()->hostGrid().leafGridView();
- Dune::MultipleCodimMultipleGeomTypeMapper<
- LeafGridView, Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout());
- for (auto &&v : vertices(hostLeafView))
- vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
- }
-
- typename ContactNetwork::BoundaryPatches frictionBoundaries;
- contactNetwork.boundaryPatches("friction", frictionBoundaries);
-
- std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches;
- contactNetwork.frictionPatches(frictionPatches);
-
- auto dataWriter =
- writeData ? std::make_unique<
- HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
- *dataFile, vertexCoordinates, vertexBases,
- frictionPatches, weakPatches)
- : nullptr;
-
- const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/home/joscha/software/dune/dune-tectonic/body");
-
IterationRegister iterationCount;
- auto const report = [&](bool initial = false) {
- if (writeData) {
- dataWriter->reportSolution(programState, contactNetwork, globalFriction);
- if (!initial)
- dataWriter->reportIterations(programState, iterationCount);
- dataFile->flush();
- }
-
- if (writeRestarts and !initial and
- programState.timeStep % restartSpacing == 0) {
- restartIO->write(programState);
- restartFile->flush();
- }
-
- if (parset.get<bool>("io.printProgress"))
- std::cout << "timeStep = " << std::setw(6) << programState.timeStep
- << ", time = " << std::setw(12) << programState.relativeTime
- << ", tau = " << std::setw(12) << programState.relativeTau
- << std::endl;
-
- if (parset.get<bool>("io.vtk.write")) {
- std::vector<ScalarVector> stress(bodyCount);
-
- for (size_t i=0; i<bodyCount; i++) {
- const auto& body = contactNetwork.body(i);
- body->assembler()->assembleVonMisesStress(body->data()->getYoungModulus(),
- body->data()->getPoissonRatio(),
- programState.u[i], stress[i]);
-
- }
-
- vtkWriter.write(programState.timeStep, programState.u, programState.v,
- programState.alpha, stress);
- }
- };
- report(true);
+ ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, true);
// -------------------
// Set up TNNMG solver
@@ -398,8 +310,8 @@ int main(int argc, char *argv[]) {
energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
}
- std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance << std::endl;
- std::cout << "must refine: " << (energyNorm > refinementTolerance) << std::endl;
+ //std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance << std::endl;
+ //std::cout << "must refine: " << (energyNorm > refinementTolerance) << std::endl;
return energyNorm > refinementTolerance;
};
@@ -432,10 +344,14 @@ int main(int argc, char *argv[]) {
typename ContactNetwork::ExternalForces externalForces;
contactNetwork.externalForces(externalForces);
+ StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+ stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
+ externalForces, stateEnergyNorms);
+
AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
- adaptiveTimeStepper(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes, current,
+ adaptiveTimeStepper(stepBase, contactNetwork, current,
programState.relativeTime, programState.relativeTau,
- externalForces, stateEnergyNorms, mustRefine);
+ mustRefine);
size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
@@ -460,7 +376,7 @@ int main(int argc, char *argv[]) {
contactNetwork.setDeformation(programState.u);
- report();
+ ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, false);
if (programState.timeStep==timeSteps) {
std::cout << "limit of timeSteps reached!" << std::endl;
--
GitLab