multi-body-problem.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef HAVE_IPOPT
#undef HAVE_IPOPT
#endif

#include <atomic>
#include <cmath>
#include <csignal>
#include <exception>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <memory>

#include <dune/common/bitsetvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/function.hh>
#include <dune/common/fvector.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>

#include <dune/grid/common/mcmgmapper.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>

#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/formatstring.hh>
#include <dune/fufem/hdf5/file.hh>

#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>

#include <dune/contact/common/deformedcontinuacomplex.hh>
#include <dune/contact/common/couplingpair.hh>
#include <dune/contact/projections/normalprojection.hh>

#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/globalfriction.hh>

#include "assemblers.hh"
#include "gridselector.hh"
#include "explicitgrid.hh"
#include "explicitvectors.hh"

#include "data-structures/enumparser.hh"
#include "data-structures/enums.hh"
#include "data-structures/contactnetwork.hh"
#include "data-structures/matrices.hh"
#include "data-structures/program_state.hh"

#include "factories/stackedblocksfactory.hh"
#include "factories/threeblocksfactory.hh"

//#include "io/hdf5-levelwriter.hh"
#include "io/hdf5/restart-io.hh"
#include "io/vtk.hh"

#include "multi-body-problem-data/bc.hh"
#include "multi-body-problem-data/mybody.hh"
#include "multi-body-problem-data/grid/mygrids.hh"

#include "spatial-solving/tnnmg/functional.hh"
//#include "spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
#include "spatial-solving/tnnmg/localbisectionsolver.hh"
#include "spatial-solving/contact/nbodyassembler.hh"
#include "spatial-solving/solverfactory.hh"

#include "time-stepping/adaptivetimestepper.hh"
#include "time-stepping/rate.hh"
#include "time-stepping/state.hh"
#include "time-stepping/updaters.hh"

#include "utils/debugutils.hh"
#include "utils/diameter.hh"

// for getcwd
#include <unistd.h>

//#include <tbb/tbb.h> //TODO multi threading preconditioner?
//#include <pthread.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/src/multi-body-problem.cfg", parset);
  Dune::ParameterTreeParser::readINITree(
      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem-%dD.cfg", dims), parset);
  Dune::ParameterTreeParser::readOptions(argc, argv, parset);
  return parset;
}

static std::atomic<bool> terminationRequested(false);
void handleSignal(int signum) { terminationRequested = true; }

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("../log.txt");
    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;

    // ----------------------
    // set up contact network
    // ----------------------
    StackedBlocksFactory<Grid, Vector> stackedBlocksFactory(parset);
    using ContactNetwork = typename StackedBlocksFactory<Grid, Vector>::ContactNetwork;
    stackedBlocksFactory.build();

    ContactNetwork& contactNetwork = stackedBlocksFactory.contactNetwork();

    /*ThreeBlocksFactory<Grid, Vector> threeBlocksFactory(parset);
    using ContactNetwork = typename ThreeBlocksFactory<Grid, Vector>::ContactNetwork;
    threeBlocksFactory.build();

    ContactNetwork& contactNetwork = threeBlocksFactory.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(), "debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
    }

    // ----------------------------
    // assemble contactNetwork
    // ----------------------------
    contactNetwork.assemble();

    //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());

    // -----------------
    // init input/output
    // -----------------
    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);
    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;


    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);


    //DUNE_THROW(Dune::Exception, "Just need to stop here!");

    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);

    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);

    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]);

      Dune::MultipleCodimMultipleGeomTypeMapper<
          DefLeafGridView, Dune::MCMGVertexLayout> const vertexMapper(body->gridView(), Dune::mcmgVertexLayout());
      for (auto &&v : vertices(body->gridView()))
        vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
    }

    //typename contactNetwork::BoundaryPatches frictionBoundaries;
    //contactNetwork.boundaryPatches("friction", frictionBoundaries);

    /*
    auto dataWriter =
        writeData ? std::make_unique<
                        HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>(
                        *dataFile, vertexCoordinates, vertexBases,
                        frictionBoundaries) //, weakPatches)
                  : nullptr;*/

    const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "/storage/mi/podlesny/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);

    // -------------------
    // Set up TNNMG solver
    // -------------------

    BitVector totalDirichletNodes;
    contactNetwork.totalNodes("dirichlet", totalDirichletNodes);

    /*for (size_t i=0; i<totalDirichletNodes.size(); i++) {
        bool val = false;
        for (size_t d=0; d<dims; d++) {
            val = val || totalDirichletNodes[i][d];
        }

        totalDirichletNodes[i] = val;
        for (size_t d=0; d<dims; d++) {
            totalDirichletNodes[i][d] = val;
        }
    }*/

    print(totalDirichletNodes, "totalDirichletNodes:");

    //using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
    using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
    using NonlinearFactory = SolverFactory<Functional, BitVector>;

    using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
    using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
    using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
                               StateUpdater<ScalarVector, Vector>>;

    BoundaryFunctions velocityDirichletFunctions;
    contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);

    BoundaryNodes dirichletNodes;
    contactNetwork.boundaryNodes("dirichlet", dirichletNodes);

    /*for (size_t i=0; i<dirichletNodes.size(); i++) {
        for (size_t j=0; j<dirichletNodes[i].size(); j++) {
        print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
        }
    }*/

    std::vector<const Dune::BitSetVector<1>*> frictionNodes;
    contactNetwork.frictionNodes(frictionNodes);

    for (size_t i=0; i<frictionNodes.size(); i++) {
        print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
    }

    Updaters current(
        initRateUpdater(
            parset.get<Config::scheme>("timeSteps.scheme"),
            velocityDirichletFunctions,
            dirichletNodes,
            contactNetwork.matrices(),
            programState.u,
            programState.v,
            programState.a),
        initStateUpdater<ScalarVector, Vector>(
            parset.get<Config::stateModel>("boundary.friction.stateModel"),
            programState.alpha,
            nBodyAssembler.getContactCouplings(),
            contactNetwork.couplings())
            );


    auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");

    const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();

    auto const mustRefine = [&](Updaters &coarseUpdater,
                                Updaters &fineUpdater) {

        //return false;
      std::cout << "Step 1" << std::endl;

      std::vector<ScalarVector> coarseAlpha;
      coarseAlpha.resize(bodyCount);
      coarseUpdater.state_->extractAlpha(coarseAlpha);

      print(coarseAlpha, "coarseAlpha:");

      std::vector<ScalarVector> fineAlpha;
      fineAlpha.resize(bodyCount);
      fineUpdater.state_->extractAlpha(fineAlpha);

      print(fineAlpha, "fineAlpha:");

      std::cout << "Step 3" << std::endl;

      ScalarVector::field_type energyNorm = 0;
      for (size_t i=0; i<stateEnergyNorms.size(); i++) {
          std::cout << "for " << i << std::endl;

          std::cout << not stateEnergyNorms[i] << std::endl;

          if (coarseAlpha[i].size()==0 || fineAlpha[i].size()==0)
              continue;

          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;
      return energyNorm > refinementTolerance;
    };


    std::signal(SIGXCPU, handleSignal);
    std::signal(SIGINT, handleSignal);
    std::signal(SIGTERM, handleSignal);

/*
    // set patch preconditioner for linear correction in TNNMG method
    using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
    using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, PatchSolver, Matrix, Vector>;

    const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");

    auto gsStep = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
    PatchSolver patchSolver(gsStep,
                               preconditionerParset.get<size_t>("maximumIterations"),
                               preconditionerParset.get<double>("tolerance"),
                               nullptr,
                               preconditionerParset.get<Solver::VerbosityMode>("verbosity"),
                               false); // absolute error

    Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
    Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<Preconditioner::Mode>("mode"));
    preconditioner.setPatchSolver(patchSolver);
    preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
*/
    // set adaptive time stepper
    typename ContactNetwork::ExternalForces externalForces;
    contactNetwork.externalForces(externalForces);

    AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
        adaptiveTimeStepper(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes, current,
                            programState.relativeTime, programState.relativeTau,
                            externalForces, stateEnergyNorms, mustRefine);

    while (!adaptiveTimeStepper.reachedEnd()) {
      programState.timeStep++;

      //preconditioner.build();
      iterationCount = adaptiveTimeStepper.advance();

      programState.relativeTime = adaptiveTimeStepper.relativeTime_;
      programState.relativeTau = adaptiveTimeStepper.relativeTau_;
      current.rate_->extractDisplacement(programState.u);
      current.rate_->extractVelocity(programState.v);
      current.rate_->extractAcceleration(programState.a);
      current.state_->extractAlpha(programState.alpha);

      print(programState.u, "current u:");
      print(programState.v, "current v:");
      print(programState.a, "current a:");
      print(programState.alpha, "current alpha:");

      contactNetwork.setDeformation(programState.u);

      report();

      if (programState.timeStep==50) {
        std::cout << "limit of timeSteps reached!" << std::endl;
        break; // TODO remove after debugging
      }

      if (terminationRequested) {
        std::cerr << "Terminating prematurely" << std::endl;
        break;
      }
    }


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