#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/stackedblocksfactory.hh>
#include <dune/tectonic/factories/threeblocksfactory.hh>
#include <dune/tectonic/io/io-handler.hh>
#include <dune/tectonic/io/hdf5-writer.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/stepbase.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;
#include <dune/tectonic/utils/reductionfactors.hh>
std::vector<std::vector<double>> allReductionFactors;
Dune::ParameterTree getParameters(int argc, char *argv[]) {
Dune::ParameterTree parset;
Dune::ParameterTreeParser::readINITree("/home/joscha/software/dune/dune-tectonic/src/multi-body-problem/multi-body-problem.cfg", parset);
Dune::ParameterTreeParser::readINITree(
Dune::Fufem::formatString("/home/joscha/software/dune/dune-tectonic/src/multi-body-problem/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("multi-body-problem.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
// ----------------------
using BlocksFactory = StackedBlocksFactory<Grid, Vector>;
//using BlocksFactory = ThreeBlocksFactory<Grid, Vector>;
BlocksFactory blocksFactory(parset);
using ContactNetwork = typename BlocksFactory::ContactNetwork;
blocksFactory.build();
ContactNetwork& 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++) {
printDofLocation(contactNetwork.body(i)->gridView());
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);
IOHandler<Assembler, ContactNetwork> ioHandler(parset.sub("io"), contactNetwork);
bool restartRead = ioHandler.read(programState);
if (!restartRead) {
programState.setupInitialConditions(parset, contactNetwork);
}
//print(programState.u, "u:");
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);
IterationRegister iterationCount;
ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, true);
//DUNE_THROW(Dune::Exception, "Just need to stop here!");
// -------------------
// 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);
StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
externalForces, stateEnergyNorms);
AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
adaptiveTimeStepper(stepBase, contactNetwork, current,
programState.relativeTime, programState.relativeTau,
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);
ioHandler.write(programState, contactNetwork, globalFriction, iterationCount, false);
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;
}
}
// output reduction factors
std::vector<double> avgReductionFactors(allReductionFactors.size());
//std::cout << "allReductionFactors.size: " << allReductionFactors.size() << std::endl;
for (size_t i=0; i<avgReductionFactors.size(); i++) {
avgReductionFactors[i] = 1;
size_t c = 0;
//std::cout << "allReductionFactors[i].size: " << allReductionFactors[i].size() << std::endl;
//print(allReductionFactors[i], "all reduction factors: ");
for (size_t j=0; j<allReductionFactors[i].size(); j++) {
avgReductionFactors[i] *= allReductionFactors[i][j];
c++;
}
avgReductionFactors[i] = std::pow(avgReductionFactors[i], 1.0/((double)c));
}
print(avgReductionFactors, "average reduction factors: ");
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;
}
}