one-body-sample.cc

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

#ifndef srcdir
#error srcdir unset
#endif

#ifdef HAVE_IPOPT
#undef HAVE_IPOPT
#endif

#ifndef HAVE_PYTHON
#error Python is required
#endif

#include <exception>
#include <iostream>

#include <boost/format.hpp>

#include <Python.h>

#include <cmath>

#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/common/shared_ptr.hh>
#include <dune/common/stdstreams.hh>
#include <dune/common/timer.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>

#include <dune/fufem/assemblers/functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/transferoperatorassembler.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functionspacebases/p0basis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/sharedpointermap.hh>
#include <dune/solvers/common/numproc.hh> // Solver::FULL
#include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>

#include <dune/tectonic/myblockproblem.hh>
#include <dune/tectonic/myconvexproblem.hh>
#include <dune/tectonic/mydirectionalconvexfunction.hh>

#include "assemblers.hh"
#include "compute_state.hh"
#include "print.hh"
#include "vtk.hh"

int const dim = 2;

template <class GridView>
void setup_boundary(GridView const &gridView,
                    Dune::BitSetVector<dim> &ignoreNodes,
                    Dune::BitSetVector<1> &neumannNodes,
                    Dune::BitSetVector<1> &frictionalNodes) {
  typedef typename GridView::template Codim<dim>::Iterator VertexLeafIterator;

  typedef Dune::MultipleCodimMultipleGeomTypeMapper<
      GridView, Dune::MCMGVertexLayout> VertexMapper;
  VertexMapper const myVertexMapper(gridView);

  size_t dirichlet_nodes = 0;
  size_t neumann_nodes = 0;
  size_t frictional_nodes = 0;
  for (VertexLeafIterator it = gridView.template begin<dim>();
       it != gridView.template end<dim>(); ++it) {
    assert(it->geometry().corners() == 1);
    Dune::FieldVector<double, dim> const coordinates = it->geometry().corner(0);
    if (coordinates[1] == 1) {
      ++dirichlet_nodes;
      size_t const id = myVertexMapper.map(*it);
      ignoreNodes[id] = true;
    } else if (coordinates[1] == 0) {
      ++frictional_nodes;
      size_t const id = myVertexMapper.map(*it);
      frictionalNodes[id] = true;
      ignoreNodes[id][1] = true; // Zero displacement in direction y
    } else if (coordinates[0] == 0 || coordinates[0] == 1) {
      ++neumann_nodes;
      size_t const id = myVertexMapper.map(*it);
      neumannNodes[id] = true;
    }
  }
  std::cout << "Number of Neumann nodes: " << neumann_nodes << std::endl;
  std::cout << "Number of Dirichlet nodes: " << dirichlet_nodes << std::endl;
  std::cout << "Number of frictional nodes: " << dirichlet_nodes << std::endl;
}

int main(int argc, char *argv[]) {
  try {
    typedef SharedPointerMap<std::string, Dune::VirtualFunction<double, double>>
    FunctionMap;
    FunctionMap functions;
    {
      Python::start();

      Python::run("import sys");
      Python::run("sys.path.append('" srcdir "')");

      Python::import("one-body-sample_neumann")
          .get("Functions")
          .toC<FunctionMap::Base>(functions);
    }

    Dune::ParameterTree parset;
    Dune::ParameterTreeParser::readINITree(srcdir "/one-body-sample.parset",
                                           parset); // FIXME
    Dune::ParameterTreeParser::readOptions(argc, argv, parset);

    Dune::Timer timer;

    typedef Dune::FieldVector<double, dim> SmallVector;
    typedef Dune::FieldMatrix<double, dim, dim> SmallMatrix;
    typedef Dune::BCRSMatrix<SmallMatrix> OperatorType;
    typedef Dune::BlockVector<SmallVector> VectorType;
    typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;

    auto const E = parset.get<double>("body.E");
    auto const nu = parset.get<double>("body.nu");
    auto const solver_tolerance = parset.get<double>("solver.tolerance");

    auto const refinements = parset.get<size_t>("grid.refinements");
    size_t const levels = refinements + 1;

    auto const verbose = parset.get<bool>("verbose");
    Solver::VerbosityMode const verbosity =
        verbose ? Solver::FULL : Solver::QUIET;

    // {{{ Set up grid
    typedef Dune::YaspGrid<dim> GridType;
    Dune::FieldVector<double, dim> const end_points(
        1); // nth dimension (zero-indexed) goes from 0 to end_points[n]
    auto grid = Dune::make_shared<GridType>(
        end_points,
        Dune::FieldVector<int, dim>(1), // number of elements in each direction
        Dune::FieldVector<bool, dim>(false), // non-periodic in each direction
        0);                                  // zero overlap (whatever that is)
    grid->globalRefine(refinements);
    size_t const finestSize = grid->size(grid->maxLevel(), dim);

    typedef GridType::LeafGridView GridView;
    GridView const leafView = grid->leafView();
    // }}}

    // Set up bases
    typedef P0Basis<GridView, double> P0Basis;
    typedef P1NodalBasis<GridView, double> P1Basis;
    P0Basis const p0Basis(leafView);
    P1Basis const p1Basis(leafView);

    // Assemble elastic force on the body
    StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
                               P1Basis::LocalFiniteElement> const
    localStiffness(E, nu);
    OperatorType stiffnessMatrix;
    {
      timer.reset();
      OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
          .assemble(localStiffness, stiffnessMatrix);
      std::cout << "Assembled stiffness matrix in " << timer.elapsed() << "s"
                << std::endl;
    }
    EnergyNorm<OperatorType, VectorType> energyNorm(stiffnessMatrix);

    // Set up the boundary
    Dune::BitSetVector<dim> ignoreNodes(grid->size(grid->maxLevel(), dim),
                                        false);
    Dune::BitSetVector<1> neumannNodes(grid->size(grid->maxLevel(), dim),
                                       false);
    Dune::BitSetVector<1> frictionalNodes(grid->size(grid->maxLevel(), dim),
                                          false);
    setup_boundary(leafView, ignoreNodes, neumannNodes, frictionalNodes);

    typedef MyConvexProblem<OperatorType, VectorType> MyConvexProblemType;
    typedef MyBlockProblem<MyConvexProblemType> MyBlockProblemType;
    GenericNonlinearGS<MyBlockProblemType> nonlinearGSStep;
    nonlinearGSStep.ignoreNodes_ = &ignoreNodes;

    VectorType u4(finestSize);
    u4 = 0.0; // Has to be zero!
    VectorType u5 = u4;

    SingletonVectorType s4_old(finestSize);
    s4_old = parset.get<double>("boundary.friction.state.initial");
    SingletonVectorType s5_old = s4_old;

    VectorType u4_diff(finestSize);
    u4_diff = 0.0; // Has to be zero!
    VectorType u5_diff = u4_diff;

    auto s4_new = Dune::make_shared<SingletonVectorType>(finestSize);
    *s4_new = s4_old;
    auto s5_new = Dune::make_shared<SingletonVectorType>(finestSize);
    *s5_new = s5_old;

    SingletonVectorType vonMisesStress;

    VectorType b4;
    VectorType b5;

    auto const nodalIntegrals =
        assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
            leafView, p1Basis, frictionalNodes);
    // {{{ Set up TNNMG solver
    // linear iteration step components
    TruncatedBlockGSStep<OperatorType, VectorType> linearBaseSolverStep;
    EnergyNorm<OperatorType, VectorType> baseEnergyNorm(linearBaseSolverStep);
    LoopSolver<VectorType> linearBaseSolver(
        &linearBaseSolverStep,
        parset.get<int>("solver.tnnmg.linear.maxiterations"), solver_tolerance,
        &baseEnergyNorm, Solver::QUIET);
    TruncatedBlockGSStep<OperatorType, VectorType> linearPresmoother;
    TruncatedBlockGSStep<OperatorType, VectorType> linearPostsmoother;

    // linear iteration step
    MultigridStep<OperatorType, VectorType> linearIterationStep;
    linearIterationStep.setNumberOfLevels(levels);
    linearIterationStep.setMGType(parset.get<int>("solver.tnnmg.linear.mu"),
                                  parset.get<int>("solver.tnnmg.linear.nu1"),
                                  parset.get<int>("solver.tnnmg.linear.nu2"));
    linearIterationStep.basesolver_ = &linearBaseSolver;
    linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);

    // transfer operators
    std::vector<CompressedMultigridTransfer<VectorType> *> transferOperators(
        refinements);
    for (auto &x : transferOperators)
      x = new CompressedMultigridTransfer<VectorType>;
    TransferOperatorAssembler<GridType>(*grid)
        .assembleOperatorPointerHierarchy(transferOperators);
    linearIterationStep.setTransferOperators(transferOperators);

    // tnnmg iteration step
    typedef GenericNonlinearGS<MyBlockProblemType> NonlinearSmootherType;
    NonlinearSmootherType nonlinearSmoother;
    TruncatedNonsmoothNewtonMultigrid<MyBlockProblemType, NonlinearSmootherType>
    multigridStep(linearIterationStep, nonlinearSmoother);
    multigridStep.setSmoothingSteps(parset.get<int>("solver.tnnmg.main.nu1"),
                                    parset.get<int>("solver.tnnmg.main.mu"),
                                    parset.get<int>("solver.tnnmg.main.nu2"));
    multigridStep.ignoreNodes_ = &ignoreNodes;
    // }}}

    std::fstream octave_writer("data", std::fstream::out);
    timer.reset();
    auto const timesteps = parset.get<size_t>("timesteps");
    octave_writer << "# name: A" << std::endl << "# type: matrix" << std::endl
                  << "# rows: " << timesteps << std::endl << "# columns: 3"
                  << std::endl;
    double const h = 1.0 / timesteps;
    for (size_t run = 1; run <= timesteps; ++run) {
      if (parset.get<bool>("printProgress")) {
        std::cout << ".";
        std::cout.flush();
      }

      if (parset.get<bool>("solver.tnnmg.use")) {
        assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
            leafView, p1Basis, neumannNodes, b4,
            functions.get("sampleFunction"), h * run);
        stiffnessMatrix.mmv(u4, b4);
        for (int state_fpi = 0;
             state_fpi < parset.get<int>("solver.tnnmg.fixed_point_iterations");
             ++state_fpi) {
          auto myGlobalNonlinearity =
              assemble_nonlinearity<VectorType, OperatorType>(
                  finestSize, parset, nodalIntegrals, s4_new, h);
          MyConvexProblemType const myConvexProblem(stiffnessMatrix,
                                                    *myGlobalNonlinearity, b4);

          MyBlockProblemType myBlockProblem(parset, myConvexProblem);
          multigridStep.setProblem(u4_diff, myBlockProblem);

          LoopSolver<VectorType> overallSolver(
              &multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
              solver_tolerance, &energyNorm, verbosity);
          overallSolver.solve();

          if (!parset.get<bool>("boundary.friction.state.evolve"))
            continue;

          for (size_t i = 0; i < frictionalNodes.size(); ++i) {
            if (frictionalNodes[i][0]) {
              double const L = parset.get<double>("boundary.friction.ruina.L");
              double const unorm = u4_diff[i].two_norm();

              (*s4_new)[i] = compute_state_update(h, unorm, L, s4_old[i]);
            }
          }
        }

        if (parset.get<bool>("printEvolution"))
          print_evolution<SingletonVectorType, VectorType>(
              frictionalNodes, *s4_new, u4, functions.get("sampleFunction"),
              run, h * run, octave_writer);
      }

      u4 += u4_diff;
      s4_old = *s4_new;

      if (parset.get<bool>("benchmarks.fpi.enable")) {
        assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
            leafView, p1Basis, neumannNodes, b5,
            functions.get("sampleFunction"), h * run);
        stiffnessMatrix.mmv(u5, b5);
        for (int state_fpi = 0;
             state_fpi < parset.get<int>("benchmarks.fpi.iterations");
             ++state_fpi) {
          auto myGlobalNonlinearity =
              assemble_nonlinearity<VectorType, OperatorType>(
                  finestSize, parset, nodalIntegrals, s5_new, h);
          MyConvexProblemType const myConvexProblem(stiffnessMatrix,
                                                    *myGlobalNonlinearity, b5);

          MyBlockProblemType myBlockProblem(parset, myConvexProblem);
          multigridStep.setProblem(u5_diff, myBlockProblem);

          LoopSolver<VectorType> overallSolver(
              &multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
              solver_tolerance, &energyNorm, verbosity);
          overallSolver.solve();

          if (!parset.get<bool>("boundary.friction.state.evolve"))
            continue;

          for (size_t i = 0; i < frictionalNodes.size(); ++i) {
            if (frictionalNodes[i][0]) {
              double const L = parset.get<double>("boundary.friction.ruina.L");
              double const unorm = u5_diff[i].two_norm();

              (*s5_new)[i] = compute_state_update(h, unorm, L, s5_old[i]);
            }
          }
        }
      }

      u5 += u5_diff;
      s5_old = *s5_new;

      if (parset.get<bool>("benchmarks.fpi.enable"))
        std::cout << energyNorm.diff(u5, u4) / energyNorm(u5) << std::endl;

      // Compute von Mises stress and write everything to a file
      if (parset.get<bool>("writeVTK")) {
        // Compute von Mises stress
        VonMisesStressAssembler<GridType> localStressAssembler(
            E, nu,
            Dune::make_shared<BasisGridFunction<P1Basis, VectorType> const>(
                p1Basis, u4));
        FunctionalAssembler<P0Basis>(p0Basis)
            .assemble(localStressAssembler, vonMisesStress, true);

        writeVtk<P1Basis, P0Basis, VectorType, SingletonVectorType, GridView>(
            p1Basis, u4, *s4_new, p0Basis, vonMisesStress, leafView,
            (boost::format("obs%d") % run).str());
      }
    }
    std::cout << std::endl;
    std::cout << "Making " << timesteps << " time steps took "
              << timer.elapsed() << "s" << std::endl;

    // Clean up after the TNNMG solver
    for (auto &x : transferOperators)
      delete x;

    octave_writer.close();

    if (parset.get<bool>("printFrictionalBoundary")) {
      // Print displacement on frictional boundary
      boost::format const formatter("u4[%02d] = %+3e, "
                                    "%|40t|u5[%02d] = %+3e");
      for (size_t i = 0; i < frictionalNodes.size(); ++i)
        if (frictionalNodes[i][0])
          std::cout << boost::format(formatter) % i % u4[i] % i % u5[i]
                    << std::endl;
    }
    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;
  }
}