#ifdef HAVE_CONFIG_H #include "config.h" #endif #ifdef HAVE_IPOPT #undef HAVE_IPOPT #endif #define MY_DIM 2 #include <atomic> #include <cmath> #include <csignal> #include <exception> #include <fstream> #include <iostream> #include <iomanip> #include <memory> #include <dune/common/exceptions.hh> #include <dune/common/parallel/mpihelper.hh> #include <dune/common/parametertree.hh> #include <dune/common/parametertreeparser.hh> #include <dune/fufem/formatstring.hh> #include <dune/fufem/hdf5/file.hh> #include <dune/tectonic/assemblers.hh> #include <dune/tectonic/gridselector.hh> #include <dune/tectonic/explicitgrid.hh> #include <dune/tectonic/explicitvectors.hh> #include <dune/tectonic/data-structures/enumparser.hh> #include <dune/tectonic/data-structures/enums.hh> #include <dune/tectonic/data-structures/network/contactnetwork.hh> #include <dune/tectonic/data-structures/program_state.hh> #include <dune/tectonic/factories/strikeslipfactory.hh> #include <dune/tectonic/io/vtk.hh> #include <dune/tectonic/io/hdf5-writer.hh> #include <dune/tectonic/utils/geocoordinate.hh> #include <dune/tectonic/utils/debugutils.hh> // for getcwd #include <unistd.h> size_t const dims = MY_DIM; Dune::ParameterTree getParameters(int argc, char *argv[]) { Dune::ParameterTree parset; Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/hdf5test/hdf5test.cfg", parset); Dune::ParameterTreeParser::readINITree( Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/dune/tectonic/tests/hdf5test/hdf5test-%dD.cfg", dims), parset); Dune::ParameterTreeParser::readOptions(argc, argv, parset); return parset; } int main(int argc, char *argv[]) { try { Dune::MPIHelper::instance(argc, argv); char buffer[256]; char *val = getcwd(buffer, sizeof(buffer)); if (val) { std::cout << buffer << std::endl; std::cout << argv[0] << std::endl; } std::ofstream out("hdf5test.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); // ---------------------- // set up contact network // ---------------------- using BlocksFactory = StrikeSlipFactory<Grid, Vector>; BlocksFactory blocksFactory(parset); using ContactNetwork = typename BlocksFactory::ContactNetwork; blocksFactory.build(); auto& contactNetwork = blocksFactory.contactNetwork(); const size_t bodyCount = contactNetwork.nBodies(); /* for (size_t i=0; i<contactNetwork.nLevels(); i++) { // printDofLocation(contactNetwork.body(i)->gridView()); //Vector def(contactNetwork.deformedGrids()[i]->size(dims)); //def = 1; //deformedGridComplex.setDeformation(def, i); 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)); } }*/ for (size_t i=0; i<bodyCount; i++) { writeToVTK(contactNetwork.body(i)->gridView(), "", "initial_body_" + std::to_string(i) + "_leaf"); } // ---------------------------- // assemble contactNetwork // ---------------------------- contactNetwork.assemble(); //printMortarBasis<Vector>(contactNetwork.nBodyAssembler()); std::vector<size_t> nVertices(bodyCount); for (size_t i=0; i<bodyCount; i++) { nVertices[i] = contactNetwork.body(i)->nVertices(); } using MyProgramState = ProgramState<Vector, ScalarVector>; MyProgramState programState(nVertices); programState.setupInitialConditions(parset, contactNetwork); auto& nBodyAssembler = contactNetwork.nBodyAssembler(); for (size_t i=0; i<bodyCount; i++) { contactNetwork.body(i)->setDeformation(programState.u[i]); } nBodyAssembler.assembleTransferOperator(); nBodyAssembler.assembleObstacle(); // ------------------------ // assemble global friction // ------------------------ contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress); auto& globalFriction = contactNetwork.globalFriction(); globalFriction.updateAlpha(programState.alpha); // ----------------- // init input/output // ----------------- using Assembler = MyAssembler<DefLeafGridView, dims>; 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()); } std::vector<const typename ContactNetwork::BoundaryPatch*> frictionPatches; contactNetwork.frictionPatches(frictionPatches); auto const writeData = parset.get<bool>("io.data.write"); auto dataFile = writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr; auto dataWriter = writeData ? std::make_unique< HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>( *dataFile, vertexCoordinates, vertexBases, frictionPatches, weakPatches) : nullptr; IterationRegister iterationCount; auto const report = [&](bool initial = false) { if (writeData) { dataWriter->reportSolution(programState, contactNetwork, globalFriction); if (!initial) dataWriter->reportIterations(programState, iterationCount); dataFile->flush(); } }; report(true); std::cout.rdbuf(coutbuf); //reset to standard output again } catch (Dune::Exception &e) { Dune::derr << "Dune reported error: " << e << std::endl; } catch (std::exception &e) { std::cerr << "Standard exception: " << e.what() << std::endl; } }