#ifndef SRC_IO_HANDLER_HH #define SRC_IO_HANDLER_HH #include <filesystem> #include <memory> #include <dune/common/parametertree.hh> #include <dune/fufem/boundarypatch.hh> #include <dune/fufem/geometry/convexpolyhedron.hh> //#include <dune/tectonic/explicitgrid.hh> #include <dune/fufem/hdf5/file.hh> #include "../utils/geocoordinate.hh" #include "../time-stepping/adaptivetimestepper.hh" #include "hdf5/restart-io.hh" #include "hdf5/iteration-writer.hh" #include "hdf5/time-writer.hh" #include "hdf5-bodywriter.hh" #include "hdf5-writer.hh" #include "vtk.hh" template <class Assembler, class ContactNetwork, class ProgramState> class IOHandler { private: using Vector = typename Assembler::Vector; using LocalVector = typename Vector::block_type; using MyVertexBasis = typename Assembler::VertexBasis; using MyCellBasis = typename Assembler::CellBasis; using HDF5Writer = HDF5Writer<ProgramState, MyVertexBasis, typename Assembler::GV>; const size_t bodyCount_; std::vector<size_t> nVertices_; std::vector<Vector> vertexCoordinates_; std::vector<const MyVertexBasis* > vertexBases_; std::vector<const MyCellBasis* > cellBases_; std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches_; std::vector<std::vector<const ConvexPolyhedron<LocalVector>*>> weakPatches_; std::unique_ptr<HDF5::File> dataFile_; std::unique_ptr<HDF5Writer> dataWriter_; bool writeVTK_; bool writeData_; bool writeRestarts_; size_t restartIdx_; size_t restartSpacing_; bool printProgress_; public: IOHandler(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) : bodyCount_(contactNetwork.nBodies()), nVertices_(bodyCount_), vertexCoordinates_(bodyCount_), vertexBases_(bodyCount_), cellBases_(bodyCount_), weakPatches_(bodyCount_), writeVTK_(parset.get<bool>("vtk.write")), writeData_(parset.get<bool>("data.write")), writeRestarts_(parset.get<bool>("restarts.write")), restartIdx_(parset.get<size_t>("restarts.first")), restartSpacing_(parset.get<size_t>("restarts.spacing")), printProgress_(parset.get<bool>("printProgress")) { for (size_t i=0; i<bodyCount_; i++) { nVertices_[i] = contactNetwork.body(i)->nVertices(); } for (size_t i=0; i<bodyCount_; i++) { const auto& body = contactNetwork.body(i); vertexBases_[i] = &(body->assembler()->vertexBasis); cellBases_[i] = &(body->assembler()->cellBasis); auto& vertexCoords = vertexCoordinates_[i]; vertexCoords.resize(nVertices_[i]); auto hostLeafView = body->grid()->hostGrid().leafGridView(); Dune::MultipleCodimMultipleGeomTypeMapper< decltype(hostLeafView), Dune::MCMGVertexLayout> const vertexMapper(hostLeafView, Dune::mcmgVertexLayout()); 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"); } contactNetwork.frictionPatches(frictionPatches_); dataFile_ = std::make_unique<HDF5::File>("output.h5"); dataWriter_ = std::make_unique<HDF5Writer>(*dataFile_, vertexCoordinates_, vertexBases_, frictionPatches_); } template <class GlobalFriction> void write(const ProgramState& programState, const ContactNetwork& contactNetwork, const GlobalFriction& friction, const IterationRegister& iterationCount, bool initial = false) { if (writeData_) { writeData(programState, contactNetwork, friction, iterationCount, initial); } if (writeVTK_) { writeVTK(programState, contactNetwork); } if (writeRestarts_) { writeRestarts(programState); } if (printProgress_) { std::cout << "timeStep = " << std::setw(6) << programState.timeStep << ", time = " << std::setw(12) << programState.relativeTime << ", tau = " << std::setw(12) << programState.relativeTau << std::endl; } } bool read(ProgramState& programState) { using ADS = RestartIO<ProgramState>; if (restartIdx_ > 0) {// automatically adjusts the time and timestep auto restartFile = std::make_unique<HDF5::File>("restarts.h5", HDF5::Access::READONLY); auto restartIO = std::make_unique<ADS>(*restartFile, nVertices_); restartIO->read(restartIdx_, programState); return true; } else { return false; } } private: void writeVTK(const ProgramState& programState, const ContactNetwork& contactNetwork) { using ScalarVector = typename Assembler::ScalarVector; 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]); } const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases_, vertexBases_, "iterates/"); vtkWriter.write(programState.timeStep, programState.u, programState.v, programState.alpha, stress); } void writeRestarts(const ProgramState& programState) { if (programState.timeStep % restartSpacing_ == 0) { auto restartFile = std::make_unique<HDF5::File>("restarts.h5", HDF5::Access::READWRITE); auto restartIO = std::make_unique<RestartIO<ProgramState>>(*restartFile, nVertices_); restartIO->write(programState); restartFile->flush(); } } template <class GlobalFriction> void writeData(const ProgramState& programState, const ContactNetwork& contactNetwork, const GlobalFriction& friction, const IterationRegister& iterationCount, bool initial) { dataWriter_->reportSolution(programState, contactNetwork, friction); if (!initial) { dataWriter_->reportIterations(programState, iterationCount); } else { dataWriter_->reportWeightedNormalStress(programState); } dataFile_->flush(); } }; #endif