hdf5test.cc

#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;
  }
}