#include <Python.h> #ifdef HAVE_CONFIG_H #include "config.h" #endif #ifndef HAVE_PYTHON #error Python is required #endif #ifdef HAVE_IPOPT #undef HAVE_IPOPT #endif #include <cmath> #include <exception> #include <fstream> #include <iostream> #include <iomanip> #include <boost/filesystem/operations.hpp> #include <boost/filesystem/path.hpp> #include <boost/format.hpp> #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/parametertree.hh> #include <dune/common/parametertreeparser.hh> #include <dune/grid/common/mcmgmapper.hh> #include <dune/istl/bcrsmatrix.hh> #include <dune/istl/bvector.hh> #pragma clang diagnostic push #pragma clang diagnostic ignored "-Wdeprecated-declarations" #include <dune/fufem/boundarypatch.hh> #pragma clang diagnostic pop #include <dune/fufem/dunepython.hh> #include <dune/fufem/functions/basisgridfunction.hh> #include <dune/fufem/sharedpointermap.hh> #include <dune/solvers/norms/energynorm.hh> #include <dune/solvers/solvers/loopsolver.hh> #include <dune/solvers/solvers/solver.hh> #include <dune/tnnmg/nonlinearities/zerononlinearity.hh> #include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh> #include <dune/tnnmg/problem-classes/convexproblem.hh> #include <dune/tectonic/myblockproblem.hh> #include <dune/tectonic/globalfriction.hh> #include "adaptivetimestepper.hh" // FIXME #include "assemblers.hh" #include "distances.hh" #include "enumparser.hh" #include "enums.hh" #include "friction_writer.hh" #include "gridselector.hh" #include "matrices.hh" #include "sand-wedge-data/mybody.hh" #include "sand-wedge-data/mygeometry.hh" #include "sand-wedge-data/myglobalfrictiondata.hh" #include "sand-wedge-data/mygrid.hh" #include "sand-wedge-data/special_writer.hh" #include "sand-wedge-data/weakpatch.hh" #include "solverfactory.hh" #include "state.hh" #include "timestepping.hh" #include "vtk.hh" size_t const dims = MY_DIM; void initPython(std::string dataDirectory) { Python::start(); Python::run("import sys"); Python::run(str(boost::format("sys.path.append('%s')") % dataDirectory)); } int main(int argc, char *argv[]) { try { auto const dataDirectory = boost::filesystem::system_complete(boost::filesystem::path(argv[0])) .parent_path() / boost::filesystem::path("sand-wedge-data"); Dune::ParameterTree parset; Dune::ParameterTreeParser::readINITree( (dataDirectory / boost::filesystem::path("parset.cfg")).string(), parset); Dune::ParameterTreeParser::readOptions(argc, argv, parset); MyGeometry::render(); MyGeometry::write(); using GridView = Grid::LeafGridView; using MyAssembler = MyAssembler<GridView, dims>; using Matrix = MyAssembler::Matrix; using LocalMatrix = Matrix::block_type; using Vector = MyAssembler::Vector; using LocalVector = Vector::block_type; using ScalarMatrix = MyAssembler::ScalarMatrix; using ScalarVector = MyAssembler::ScalarVector; using LocalScalarVector = ScalarVector::block_type; auto weakPatch = getWeakPatch<LocalVector>(); // {{{ Set up grid GridConstructor<Grid> gridConstructor; auto grid = gridConstructor.getGrid(); refine(*grid, weakPatch, parset.get<double>("boundary.friction.smallestDiameter")); double minDiameter = std::numeric_limits<double>::infinity(); double maxDiameter = 0.0; for (auto it = grid->template leafbegin<0>(); it != grid->template leafend<0>(); ++it) { auto const geometry = it->geometry(); auto const diam = diameter(geometry); minDiameter = std::min(minDiameter, diam); maxDiameter = std::max(maxDiameter, diam); } std::cout << "min diameter: " << minDiameter << std::endl; std::cout << "max diameter: " << maxDiameter << std::endl; GridView const leafView = grid->leafGridView(); size_t const leafVertexCount = leafView.size(dims); std::cout << "Number of DOFs: " << leafVertexCount << std::endl; auto myFaces = gridConstructor.constructFaces(leafView); // Neumann boundary BoundaryPatch<GridView> const neumannBoundary(leafView); // Frictional Boundary BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower; Dune::BitSetVector<1> frictionalNodes(leafVertexCount); frictionalBoundary.getVertices(frictionalNodes); // Surface BoundaryPatch<GridView> const &surface = myFaces.upper; Dune::BitSetVector<1> surfaceNodes(leafVertexCount); surface.getVertices(surfaceNodes); // Dirichlet Boundary Dune::BitSetVector<dims> noNodes(leafVertexCount); Dune::BitSetVector<dims> dirichletNodes(leafVertexCount); for (size_t i = 0; i < leafVertexCount; ++i) { if (myFaces.right.containsVertex(i)) dirichletNodes[i][0] = true; if (myFaces.lower.containsVertex(i)) dirichletNodes[i][1] = true; #if MY_DIM == 3 if (myFaces.front.containsVertex(i) || myFaces.back.containsVertex(i)) dirichletNodes[i][2] = true; #endif } // Set up functions for time-dependent boundary conditions using Function = Dune::VirtualFunction<double, double>; using FunctionMap = SharedPointerMap<std::string, Function>; FunctionMap functions; { initPython(dataDirectory.string()); Python::import("boundaryconditions") .get("Functions") .toC<typename FunctionMap::Base>(functions); } auto const &velocityDirichletFunction = functions.get("velocityDirichletCondition"), &neumannFunction = functions.get("neumannCondition"); std::fstream timeStepWriter("timeSteps", std::fstream::out); timeStepWriter << std::fixed << std::setprecision(10); auto const reportTimeStep = [&](double _relativeTime, double _relativeTau) { timeStepWriter << _relativeTime << " " << _relativeTau << std::endl; }; MyAssembler myAssembler(leafView); MyBody<dims> const body(parset); Matrices<Matrix> matrices; myAssembler.assembleElasticity(body.getYoungModulus(), body.getPoissonRatio(), matrices.elasticity); myAssembler.assembleViscosity(body.getShearViscosityField(), body.getBulkViscosityField(), matrices.damping); myAssembler.assembleMass(body.getDensityField(), matrices.mass); ScalarMatrix frictionalBoundaryMass; myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary, frictionalBoundaryMass); EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm( frictionalBoundaryMass); // Assemble forces Vector gravityFunctional; myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional); // Problem formulation: right-hand side std::function<void(double, Vector &)> computeExternalForces = [&](double _relativeTime, Vector &_ell) { myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction, _relativeTime); _ell += gravityFunctional; }; // helper auto const solveLinearProblem = [&]( Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix, Vector const &_rhs, Vector &_x, Dune::ParameterTree const &_localParset) { using LinearFactory = SolverFactory< dims, BlockNonlinearTNNMGProblem<ConvexProblem< ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>, Grid>; ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity; LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME *grid, _dirichletNodes); typename LinearFactory::ConvexProblem convexProblem( 1.0, _matrix, zeroNonlinearity, _rhs, _x); typename LinearFactory::BlockProblem problem(parset, convexProblem); auto multigridStep = factory.getStep(); multigridStep->setProblem(_x, problem); EnergyNorm<Matrix, Vector> const norm(_matrix); LoopSolver<Vector> solver( multigridStep.get(), _localParset.get<size_t>("maximumIterations"), _localParset.get<double>("tolerance"), &norm, _localParset.get<Solver::VerbosityMode>("verbosity"), false); // absolute error solver.preprocess(); solver.solve(); }; // {{{ Initial conditions Vector ell0(leafVertexCount); computeExternalForces(0.0, ell0); // Start from a situation of minimal stress, which is // automatically attained in the case [v = 0 = a]. Assuming // dPhi(v = 0) = 0, we thus only have to solve Au = ell0 Vector u_initial(leafVertexCount); solveLinearProblem(dirichletNodes, matrices.elasticity, ell0, u_initial, parset.sub("u0.solver")); ScalarVector normalStress(leafVertexCount); myAssembler.assembleNormalStress(frictionalBoundary, normalStress, body.getYoungModulus(), body.getPoissonRatio(), u_initial); ScalarVector alpha_initial(leafVertexCount); alpha_initial = parset.get<double>("boundary.friction.initialAlpha"); // Start from zero velocity Vector v_initial(leafVertexCount); v_initial = 0.0; // Compute the acceleration from Ma = ell0 - Au (without Dirichlet // constraints), again assuming dPhi(v = 0) = 0 Vector a_initial(leafVertexCount); Vector accelerationRHS = ell0; Arithmetic::subtractProduct(accelerationRHS, matrices.elasticity, u_initial); solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a_initial, parset.sub("a0.solver")); // }}} MyGlobalFrictionData<LocalVector> frictionInfo( parset.sub("boundary.friction"), weakPatch); auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity( parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), frictionalBoundary, frictionInfo, normalStress); myGlobalFriction->updateAlpha(alpha_initial); Vector vertexCoordinates(leafVertexCount); { Dune::MultipleCodimMultipleGeomTypeMapper< GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView); for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it) { auto const geometry = it->geometry(); assert(geometry.corners() == 1); vertexCoordinates[vertexMapper.index(*it)] = geometry.corner(0); } } FrictionWriter<ScalarVector, Vector> frictionWriter( vertexCoordinates, frictionalNodes, "friction", MyGeometry::horizontalProjection); BoundaryWriter<ScalarVector, Vector> verticalSurfaceWriter( vertexCoordinates, surfaceNodes, "verticalSurface", MyGeometry::verticalProjection); BoundaryWriter<ScalarVector, Vector> horizontalSurfaceWriter( vertexCoordinates, surfaceNodes, "horizontalSurface", MyGeometry::horizontalProjection); SpecialWriter<GridView, dims> specialVelocityWriter("specialVelocities", leafView); SpecialWriter<GridView, dims> specialDisplacementWriter( "specialDisplacements", leafView); auto const report = [&](Vector const &_u, Vector const &_v, ScalarVector const &_alpha) { horizontalSurfaceWriter.writeKinetics(_u, _v); verticalSurfaceWriter.writeKinetics(_u, _v); ScalarVector c; myGlobalFriction->coefficientOfFriction(_v, c); frictionWriter.writeKinetics(_u, _v); frictionWriter.writeOther(c, _alpha); BasisGridFunction<typename MyAssembler::VertexBasis, Vector> velocity( myAssembler.vertexBasis, _v); BasisGridFunction<typename MyAssembler::VertexBasis, Vector> displacement( myAssembler.vertexBasis, _u); specialVelocityWriter.write(velocity); specialDisplacementWriter.write(displacement); }; report(u_initial, v_initial, alpha_initial); MyVTKWriter<typename MyAssembler::VertexBasis, typename MyAssembler::CellBasis> const vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs"); if (parset.get<bool>("io.writeVTK")) { ScalarVector stress; myAssembler.assembleVonMisesStress( body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress); vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress); } // Set up TNNMG solver using NonlinearFactory = SolverFactory< dims, MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>, Grid>; NonlinearFactory factory(parset.sub("solver.tnnmg"), *grid, dirichletNodes); using UpdaterPair = std::pair< std::shared_ptr<StateUpdater<ScalarVector, Vector>>, std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dims>>>; UpdaterPair current( initStateUpdater<ScalarVector, Vector>( parset.get<Config::stateModel>("boundary.friction.stateModel"), alpha_initial, frictionalNodes, parset.get<double>("boundary.friction.L"), parset.get<double>("boundary.friction.V0")), initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"), velocityDirichletFunction, dirichletNodes, matrices, u_initial, v_initial, a_initial)); auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance"); auto const mustRefine = [&](UpdaterPair &coarseUpdater, UpdaterPair &fineUpdater) { ScalarVector coarseAlpha; coarseUpdater.first->extractAlpha(coarseAlpha); ScalarVector fineAlpha; fineUpdater.first->extractAlpha(fineAlpha); return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance; }; size_t timeStep = 1; AdaptiveTimeStepper<NonlinearFactory, UpdaterPair, EnergyNorm<ScalarMatrix, ScalarVector>> adaptiveTimeStepper(factory, parset, myGlobalFriction, current, computeExternalForces, stateEnergyNorm, mustRefine); while (!adaptiveTimeStepper.reachedEnd()) { adaptiveTimeStepper.advance(); reportTimeStep(adaptiveTimeStepper.getRelativeTime(), adaptiveTimeStepper.getRelativeTau()); Vector u, v; ScalarVector alpha; current.second->extractDisplacement(u); current.second->extractVelocity(v); current.first->extractAlpha(alpha); report(u, v, alpha); if (parset.get<bool>("io.writeVTK")) { ScalarVector stress; myAssembler.assembleVonMisesStress(body.getYoungModulus(), body.getPoissonRatio(), u, stress); vtkWriter.write(timeStep, u, v, alpha, stress); } timeStep++; } timeStepWriter.close(); Python::stop(); } catch (Dune::Exception &e) { Dune::derr << "Dune reported error: " << e << std::endl; } catch (std::exception &e) { std::cerr << "Standard exception: " << e.what() << std::endl; } }