foam.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/geometry/convexpolyhedron.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/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/matrices.hh>
#include <dune/tectonic/data-structures/program_state.hh>
#include <dune/tectonic/data-structures/friction/globalfriction.hh>

#include <dune/tectonic/factories/twoblocksfactory.hh>

#include <dune/tectonic/io/hdf5-levelwriter.hh>
#include <dune/tectonic/io/hdf5/restart-io.hh>
#include <dune/tectonic/io/vtk.hh>

#include <dune/tectonic/problem-data/bc.hh>
#include <dune/tectonic/problem-data/mybody.hh>
#include <dune/tectonic/problem-data/grid/mygrids.hh>

#include <dune/tectonic/spatial-solving/tnnmg/functional.hh>
//#include <dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh>
#include <dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh>
#include <dune/tectonic/spatial-solving/solverfactory.hh>

#include <dune/tectonic/time-stepping/adaptivetimestepper.hh>
#include <dune/tectonic/time-stepping/rate.hh>
#include <dune/tectonic/time-stepping/state.hh>
#include <dune/tectonic/time-stepping/updaters.hh>

#include <dune/tectonic/utils/debugutils.hh>
#include <dune/tectonic/utils/diameter.hh>
#include <dune/tectonic/utils/geocoordinate.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/foam/foam.cfg", parset);
  Dune::ParameterTreeParser::readINITree(
      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/foam/foam-%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("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;

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

    ContactNetwork& contactNetwork = twoBlocksFactory.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);

    auto& wPatches = twoBlocksFactory.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]);

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

    size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");

    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);
      globalFriction.updateAlpha(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==timeSteps) {
        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;
  }
}