brittle-fracture.cc

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>

#include <fenv.h>

#include <dune/common/parametertree.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/parametertreeparser.hh>

#include <dune/grid/uggrid.hh>
#include <dune/grid/onedgrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/io/file/vtk.hh>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/preconditioners.hh>

#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/transformedindexbasis.hh>
#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>

#include <dune/fufem/functions/localbasisderivativefunction.hh>
#include <dune/fufem/assemblers/dunefunctionsoperatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/dunefunctionsl2functionalassembler.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/fufem/mappedmatrix.hh>

#include <dune/fracture-phasefields/brittlefractureenergynormassembler.hh>
#include <dune/fracture-phasefields/cracksurfaceassembler.hh>
#include <dune/fracture-phasefields/elementquantityassembler.hh>

#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/solvers/umfpacksolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>

#include <dune/matrix-vector/transpose.hh>

#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
#include <dune/tnnmg/functionals/sumfunctional.hh>
#include <dune/tnnmg/functionals/sumfunctionalconstrainedlinearization.hh>
#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
#include <dune/tnnmg/iterationsteps/tnnmgacceleration.hh>
#include <dune/tnnmg/problem-classes/bisection.hh>
#include <dune/tnnmg/projections/obstacledefectprojection.hh>

#include <dune/fracture-phasefields/brittlefracturejetevaluator.hh>
#include <dune/fracture-phasefields/gradientvaluejetevaluator.hh>
#include <dune/fracture-phasefields/gradienttostrainvaluemap.hh>
#include <dune/fracture-phasefields/brittlefracturelocalsolver.hh>
#include <dune/fracture-phasefields/brittlefracturepreconditionedlocalsolver.hh>
#include <dune/fracture-phasefields/integralfunctional.hh>
#include <dune/fracture-phasefields/integralfunctionalconstrainedlinearization.hh>
#include <dune/fracture-phasefields/damagedelasticenergydensity_decomposed.hh>
#include <dune/fracture-phasefields/damagedelasticenergydensity_nondecomposed.hh>
#include <dune/fracture-phasefields/integraldensity.hh>


//#include <dune/matrix-vector/densematrixview.hh>
//#include <dune/matrix-vector/densevectorview.hh>

//#define GAUSSSEIDEL
using namespace Dune;


/** \brief Sum of two functionals defined on the set of real numbers
 *
 * \tparam DamageFunctional Expected to model BoxConstrainedQuadraticFunctional
 * \tparam DamagedElasticityFunctional Expected to be smooth but nonquadratic
 */
template <class DamageFunctional, class DamagedElasticityFunctional>
class ScalarSumFunctional
{
public:
  ScalarSumFunctional(const DamageFunctional& damageFunctional,
                      const DamagedElasticityFunctional& damagedElasticityFunctional)
  : damageFunctional_(damageFunctional),
    damagedElasticityFunctional_(damagedElasticityFunctional)
    {}

  Solvers::Interval<double> domain() const
  {
    return Solvers::Interval<double>{damageFunctional_.lowerObstacle(),
                                     damageFunctional_.upperObstacle() };
  }

  Solvers::Interval<double> subDifferential(double z) const
  {
    constexpr double tolerance = 1e-15;

    // quadratic part with obstacle
    Solvers::Interval<double> dJ = z*damageFunctional_.quadraticPart() - damageFunctional_.linearPart();

    // Compute the subdifferential of the obstacle constraint by hand
    // TODO: This is stupid -- the functional should do this for us
    if (z < damageFunctional_.lowerObstacle())
      dJ = -std::numeric_limits<double>::max();

    if (z <= tolerance + damageFunctional_.lowerObstacle())
      dJ[0] = -std::numeric_limits<double>::max();

    if (z >= -tolerance + damageFunctional_.upperObstacle())
      dJ[1] = std::numeric_limits<double>::max();

    if (z > damageFunctional_.upperObstacle())
      dJ = std::numeric_limits<double>::max();

    // Smooth nonquadratic part
    dJ += damagedElasticityFunctional_.subDifferential(z);

    return dJ;
  }

  const DamageFunctional& damageFunctional_;
  const DamagedElasticityFunctional& damagedElasticityFunctional_;

};

/* Class to project a correction onto the admissible set.
 * All we have is bound constraints.  Therefore, in an ideal world, we would
 * simply use the ObstacleDefectProjection class instead of this one.  However,
 * our functional is a sum of two terms, and only the first one has the
 * lowerObstacle/upperObstacle members.  Therefore we need this little shim.
 */
struct DamageDefectProjection
{
  template <class Functional, class Vector>
  constexpr void operator()(const Functional& f, const Vector& x, Vector& v)
  {
    TNNMG::ObstacleDefectProjection projection;
    projection(std::get<0>(f.functions()), x, v);
  }
};


int main (int argc, char *argv[]) try
{
  feenableexcept(FE_INVALID);
  // Init MPI, if present
  MPIHelper::instance(argc,argv);

  if (argc < 3)
    DUNE_THROW(Exception, "Usage: ./brittle-fracture <python path> <python module without extension>");

  // Start Python interpreter
  Python::start();
  auto pyMain = Python::main();
  pyMain.runStream()
      << std::endl << "import math"
      << std::endl << "import sys"
      << std::endl << "sys.path.append('" << argv[1] << "')"
      << std::endl;
  auto pyModule = pyMain.import(argv[2]);

  // parse data file
  ParameterTree parameterSet;
  pyModule.get("parameterSet").toC(parameterSet);

  // read possible further parameters from the command line
  ParameterTreeParser::readOptions(argc, argv, parameterSet);

  // Print all parameters, to make them appear in the log file
  std::cout << "Executable: brittle-fracture, with parameters:" << std::endl;
  parameterSet.report();

  int numLevels         = parameterSet.get<int>("numLevels");
  int numTimeSteps      = parameterSet.get<int>("numTimeSteps");

  // material parameters
  const ParameterTree materialParameters = parameterSet.sub("materialParameters");
  double k              = materialParameters.get<double>("k");      // Residual elastic stiffness

  // We would like to be able to do
  // auto crackSurfaceDensity = parameterSet.get<FracturePhasefields::CrackSurfaceDensity>("materialParameters");
  // but this does not work so far.
  auto crackSurfaceDensity = FracturePhasefields::parseCrackSurfaceDensity(materialParameters);

  // solver settings
  double tolerance      = parameterSet.get<double>("tolerance");
  int maxIterations     = parameterSet.get<int>("maxIterations");
  int nu1               = parameterSet.get<int>("nu1");       // number of pre-smoothing steps
  int nu2               = parameterSet.get<int>("nu2");       // number of post-smoothing steps
#ifdef GAUSSSEIDEL
  double baseTolerance  = parameterSet.get<double>("baseTolerance");
#endif
  std::string linearCorrection = parameterSet.get<std::string>("linearCorrection");

  std::string resultPath = parameterSet.get("resultPath", "");

  int materialModelType  = parameterSet.get<int>("materialModelType");    //0 = small strain non decomposed
                                                                          //1 = small strain decomposed

  int outerIterations = parameterSet.get("outerIterations",1);
  int displacementIterations = parameterSet.get("displacementIterations",1);

  if (outerIterations == -1 && displacementIterations == 1)    // Displacement and damage are smoothed together
    std::cout << "using smoother 1" << std::endl;              // via operator-splitting.  The displacement substep
                                                               // applies only one Newton iteration

  if (outerIterations == 1 && displacementIterations == 1)     // Displacement and damage are smoothed separately.
    std::cout << "using smoother 2" << std::endl;              // Displacement only gets one Newton Step.

  if (outerIterations == -1 && displacementIterations == -1)   // Displacement and damage are smoothed together
    std::cout << "using smoother 3" << std::endl;              // via operator-splitting. The displacement substep
                                                               // is solved exactly by a Newton method.

  std::string localSolverType = parameterSet.get<std::string>("localSolver.type", "");


  ///////////////////////////////
  //   Generate the grid
  ///////////////////////////////

  // The grid dimension
#ifndef BRITTLE_FRACTURE_DIMENSION
  const int dim = 2;
#else
  const int dim = BRITTLE_FRACTURE_DIMENSION;
#endif

  // The grid type (dim=1: OneDGrid; dim>1:UGGrid)
  using GridType = std::conditional<dim==1,OneDGrid,UGGrid<dim> >::type;

  std::shared_ptr<GridType> grid;
  FieldVector<double,dim> lower(0), upper(1);
  std::array<unsigned int,dim> elements;

  if (parameterSet.get<bool>("structuredGrid"))
  {
    lower = parameterSet.get<FieldVector<double,dim> >("lower");
    upper = parameterSet.get<FieldVector<double,dim> >("upper");

    elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
    grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
  }
  else
  {
    std::string path     = parameterSet.get<std::string>("path");
    std::string gridFile = parameterSet.get<std::string>("gridFile");
    grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
  }

  typedef GridType::LeafGridView GridView;
  using Coordinate = Dune::FieldVector<double, dim>;

  GridView gridView = grid->leafGridView();

  // refine uniformly
  grid->globalRefine(numLevels-1);

  // The number of damage values
  const size_t nDamageComponents = 1;


  // number of components
  static const std::size_t N=dim+1;

  using namespace Dune::Indices;
  using namespace Dune::Functions::BasisBuilder;

  // Create an index transformation that transforms the raw blocked indices
  // to implement leaf blocking for composite nodes
  auto leafBlockingTransformation = Experimental::indexTransformation(
      [](auto& multiIndex, const auto& basis) {
        if (multiIndex[0] == 0)
          multiIndex = {multiIndex[1], multiIndex[2]};
        else if (multiIndex[0] == 1)
          multiIndex = {multiIndex[1], dim};
      },
      [](const auto& prefix, const auto& basis) -> std::size_t {
        if (prefix.size()==0)
          return basis.size(Dune::ReservedVector<std::size_t,1>({0}));
        if (prefix.size()==1)
          return dim+1;
        return 0;
      },
      _2, _2);  // Min/max index size

  auto blockedBasis = makeBasis(
    grid->leafGridView(),
    Experimental::transformIndices(
      composite(
        power<dim>(lagrange<1>()),
        lagrange<1>()
      ),
    leafBlockingTransformation));

  auto deformationBasis = Functions::subspaceBasis(blockedBasis, _0);
  auto damageBasis = Functions::subspaceBasis(blockedBasis, _1);

  // Vector marking all degrees of freedom used as Dirichlet nodes
  Dune::BlockVector<Dune::FieldVector<bool,dim+1>> dirichletNodes(gridView.size(dim));

  // Compute dirichlet nodes by evaluating a function given in python
  auto dirichletIndicatorFunction = Python::make_function<Dune::TupleVector<Dune::FieldVector<double,dim>,double>>(pyModule.get("dirichlet_indicator"));
  Functions::interpolate(blockedBasis, dirichletNodes, dirichletIndicatorFunction);

  //////////////////////////////////////////////////////////////////////////////////
  //   Assemble matrix and rhs for crack surface density
  //////////////////////////////////////////////////////////////////////////////////

  using Vector = BlockVector<FieldVector<double,N> >;
  using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,N,N> >;

  auto matrix = Matrix();
  auto crackRHS = Vector();

  {
    FracturePhasefields::CrackSurfaceAssembler<decltype(blockedBasis)> localCrackAssembler(crackSurfaceDensity);

    // Assemble linear term for rhs
    auto functionalAssembler = Dune::Fufem::DuneFunctionsFunctionalAssembler{damageBasis};
    functionalAssembler.assembleBulk(Functions::istlVectorBackend(crackRHS), localCrackAssembler);

    // Assemble bilinear term for matrix
    // To be sure that the pattern can hold the Hessian for a linearization,
    // we assemble the pattern and entries independently. While the former
    // is assembled for the whole basis, the latter only considers the
    // damage components.
    auto matrixBackend = Dune::Fufem::istlMatrixBackend(matrix);
    {
      // Assemble full basis
      auto operatorAssembler = Dune::Fufem::DuneFunctionsOperatorAssembler{blockedBasis, blockedBasis};
      auto patternBuilder = matrixBackend.patternBuilder();
      operatorAssembler.assembleBulkPattern(patternBuilder);
      patternBuilder.setupMatrix();
      matrixBackend.assign(0);
    }
    {
      // Assemble entries only for the damage subspace
      auto operatorAssembler = Dune::Fufem::DuneFunctionsOperatorAssembler{damageBasis, damageBasis};
      operatorAssembler.assembleBulkEntries(Fufem::istlMatrixBackend(matrix), localCrackAssembler);
    }
  }

  //////////////////////////////////////////////////////////////////////////////////
  //   Compute an H^1/L^2 matrix to measure the convergence with
  //////////////////////////////////////////////////////////////////////////////////

  auto energyNormMatrix = Matrix();
  FracturePhasefields::BrittleFractureEnergyNormAssembler<decltype(blockedBasis)>
      localEnergyNormAssembler(blockedBasis,
                               crackSurfaceDensity,
                               parameterSet.sub("materialParameters"),
                               k);   // Scaling for the elasticity part
  auto energyNormAssembler = Fufem::duneFunctionsOperatorAssembler(blockedBasis, blockedBasis);

  energyNormAssembler.assembleBulk(Fufem::istlMatrixBackend(energyNormMatrix),
                                   localEnergyNormAssembler);

  //////////////////////////////////////////////////////////////////////////////////
  //   Compute a vector that contains a the integrals over each shape function
  //////////////////////////////////////////////////////////////////////////////////

  auto scalarBasis = makeBasis(grid->leafGridView(), lagrange<1>());
  using ScalarLocalFE = decltype(scalarBasis)::LocalView::Tree::FiniteElement;

  auto constantOne = Functions::makeAnalyticGridViewFunction([](auto x){return 1.0;}, gridView);

  auto nodalWeightsAssembler = Fufem::makeDuneFunctionsL2FunctionalAssembler<GridView,ScalarLocalFE>(constantOne,
                                                                                                     QuadratureRuleKey(dim,0));

  BlockVector<double> nodalWeights;

  auto damageFunctionalAssembler = Fufem::DuneFunctionsFunctionalAssembler(scalarBasis);
  damageFunctionalAssembler.assembleBulk(Functions::istlVectorBackend(nodalWeights),
                                         nodalWeightsAssembler);

  //////////////////////////////////////////////////////////////
  //  Create nonlinear coupling functional
  //////////////////////////////////////////////////////////////

  // number of quadrature points per element
  constexpr std::size_t maxQuadPointsSelect[3]={2,4,8};
  static const std::size_t maxQuadPoints=maxQuadPointsSelect[dim-1];

  // quadrature order
  static const std::size_t quadOrder=2;

  // Infrastructure for precomputing things at each quadrature point
  FracturePhasefields::GradientValueJetEvaluator<decltype(scalarBasis)> gradientValueJetEvaluator;

  using SmallJet = decltype(gradientValueJetEvaluator)::Jet;
  auto gradientValueAssembler = elementQuantityAssembler<SmallJet, maxQuadPoints>(scalarBasis,
                                                                                  gradientValueJetEvaluator,
                                                                                  quadOrder);

  // Evaluate all quadrature weights on the grid
  using GlobalQuadratureWeightVector = decltype(gradientValueAssembler)::WeightVector;
  auto globalQuadratureWeightVector = gradientValueAssembler.assembleWeightVector();

  // Evaluate shape function values and gradients at each quadrature node
  using GradientValueMatrix = BCRSMatrix<FieldMatrix<double,SmallJet::size(),1> >;
  auto gradientValueMatrix = gradientValueAssembler.assembleEvaluationMatrix<GradientValueMatrix>();

  // Shape function strains and values at each quadrature node
  // These are not evaluated and stored, but computed on the fly from the shape function gradients and values.
  FracturePhasefields::GradientToStrainValueMap<GradientValueMatrix,dim, false> gradientToStrainValueMap;
  using StrainValueMatrix = MappedMatrix<decltype(gradientToStrainValueMap)>;
  StrainValueMatrix strainValueMatrix(gradientToStrainValueMap, gradientValueMatrix);

  // The transpose of the matrix with strains and values
  // Again only the gradients and values are actually stored, and the strains are computed as they are accessed.
  using TransposedGradientValueMatrix = Dune::MatrixVector::Transposed<GradientValueMatrix>;
  auto transposedGradientValueMatrix = TransposedGradientValueMatrix{};
  Dune::MatrixVector::transpose(gradientValueMatrix, transposedGradientValueMatrix);

  FracturePhasefields::GradientToStrainValueMap<TransposedGradientValueMatrix,dim,true> transposedMap;
  MappedMatrix<decltype(transposedMap)> transposedStrainValueMatrix(transposedMap, transposedGradientValueMatrix);

  //////////////////////////////////////////////////////////////
  //  Create initial iterate
  //////////////////////////////////////////////////////////////

  using Vector = BlockVector<FieldVector<double,dim+nDamageComponents> >;

  // Initial iterate
  Vector x;
  Functions::istlVectorBackend(x).resize(blockedBasis);
  Functions::istlVectorBackend(x) = 0.0;

  
  // Compute initial value by evaluating a function given in python
  auto initial_damage = Python::make_function<double>(pyModule.get("initial_damage"));
  Functions::interpolate(damageBasis, x, initial_damage);

//  PythonFunction<Coordinate, Coordinate> boundary_deformation(pyModule.get("boundary_deformation"));
//  Functions::interpolate(deformationBasis, x, boundary_deformation);

  //////////////////////////////
  //   Create a solver
  //////////////////////////////

  Vector rhs;
  Functions::istlVectorBackend(rhs).resize(blockedBasis);
  Functions::istlVectorBackend(rhs) = 0.0;

  using BitVector = Solvers::DefaultBitVector_t<Vector>;

  // first create a linear iteration step
  // Gauss-Seidel is the base solver

#ifdef GAUSSSEIDEL
  TruncatedBlockGSStep<Matrix,Vector> linearBaseSolverStep;

  EnergyNorm<Matrix,Vector> baseEnergyNorm(linearBaseSolverStep);
  Solvers::LoopSolver<Vector> linearBaseSolver(linearBaseSolverStep,
                                          1000,
                                          baseTolerance,
                                          baseEnergyNorm,
                                          Solver::QUIET);
#else
  // The matrix of the linear correction problem is symmetric,
  // but possibly indefinite.
  Solvers::UMFPackSolver<Matrix,Vector> linearBaseSolver;
#endif

  std::shared_ptr<LinearIterationStep<Matrix,Vector> > linearIterationStep;

  if (linearCorrection == "geometricMultigrid")
  {
    auto localLinearSolver = Dune::Solvers::BlockGS::LocalSolvers::direct();
    using LinearSmoother = Dune::Solvers::BlockGSStep<std::decay_t<decltype(localLinearSolver)>,Matrix,Vector, BitVector>;

//    using LinearSmoother = TruncatedBlockGSStep<Matrix, Vector>;

    // Make pre and postsmoothers
    auto linearSmoother  = std::make_shared<LinearSmoother>(localLinearSolver);

    auto multigridStep = std::make_shared<MultigridStep<Matrix, Vector> >();

    multigridStep->setMGType(1, nu1, nu2);
    multigridStep->setBaseSolver(linearBaseSolver);
    multigridStep->setSmoother(linearSmoother);

    multigridStep->mgTransfer_.resize(numLevels-1);
    for (size_t i=0; i<multigridStep->mgTransfer_.size(); i++)
    {
      CompressedMultigridTransfer<Vector>* newTransfer = new CompressedMultigridTransfer<Vector>;
      newTransfer->setup(*grid, i, i+1);
      multigridStep->mgTransfer_[i] = Dune::stackobject_to_shared_ptr(*newTransfer);
    }

    linearIterationStep = multigridStep;
  } else
    DUNE_THROW(Exception, "linearCorrection has to be 'geometricMultigrid'");

  // Hack: for the time being, we need the nodal weights in a vector of the same type as all the other vectors.
  // Otherwise, the dissipation functional implementation is not able to correctly restrict the entries.
  Vector coefficients(x.size());
  for (size_t i=0; i<nodalWeights.size(); i++)
    coefficients[i] = nodalWeights[i];

  // Define ignore
  BitVector ignore;
  ignore.resize(dirichletNodes.size());
  for(std::size_t k=0; k<dirichletNodes.size(); ++k)
    for(std::size_t j=0; j<dirichletNodes[k].size(); ++j)
      ignore[k][j] = dirichletNodes[k][j];


  /////////////////////////////
  // Definde elasticity model
  /////////////////////////////

  using Jet = FieldVector<double, dim*(dim+1)/2 + 1>;
  using LocalEnergy = FracturePhasefields::IntegralDensity<Jet>;
   //  small strain non decomposed (default)
       using LocalEnergyNonDecomposed = FracturePhasefields::NonDecomposed::DamagedElasticEnergyDensity<dim>;
  auto localEnergy = LocalEnergy(LocalEnergyNonDecomposed(materialParameters));

   switch (materialModelType)
     {
     case 0:   //  small strain non decomposed (default)
       std::cout << "using small strain non decomposed model" << std::endl;
       break;

     case 1:   //  small strain decomposed
       std::cout << "using small strain decomposed model" << std::endl;
       using LocalEnergyDecomposed = FracturePhasefields::Decomposed::DamagedElasticEnergyDensity<dim>;
       localEnergy = LocalEnergy(LocalEnergyDecomposed(materialParameters));
       break;

    default:
      DUNE_THROW(NotImplemented, "Unknown material model type");
   }

  ////////////////////////////////////////////////////////
  //  Time loop
  ////////////////////////////////////////////////////////

  Vector xOld;

  // Time-dependent Dirichlet boundary value function
   auto boundaryDeformation = Python::make_function<Coordinate>(pyModule.get("boundary_deformation"));

  // Open log files for various measurements
  std::ofstream iterationsFile(resultPath + "iterations.dat", std::ofstream::out);
  std::ofstream reacForceFile(resultPath + "load-displacement.dat", std::ofstream::out);
  std::ofstream wallTimeFile(resultPath + "wall-time.dat", std::ofstream::out);

  for (int i=0; i<numTimeSteps; i++)
  {
    // Construct the local solver / the smoother
    using NewtonLocalSolver = FracturePhasefields::BrittleFractureLocalSolver<dim>;
    using PreconditionedLocalSolver = std::vector<FracturePhasefields::BrittleFracturePreconditionedLocalSolver<dim>>;
    auto localSolverVariant = std::variant<std::monostate, NewtonLocalSolver, PreconditionedLocalSolver>{};

    if (localSolverType == "exact")
      localSolverVariant = FracturePhasefields::BrittleFractureLocalSolver<dim>(outerIterations,displacementIterations);
    if (localSolverType == "preconditioned")
    {
      // There is a separate local solver object for each Lagrange point.
      // Each local solver object gets the diagonal matrix block of the linear elasticity matrix.
      PreconditionedLocalSolver preconditionedLocalSolvers(energyNormMatrix.N());

      for (std::size_t j=0; j<energyNormMatrix.N(); j++)
      {
        for (int ii=0; ii<dim; ii++)
          for (int jj=0; jj<dim; jj++)
            // (1+k) times the elasticity matrix is an upper model.
            // The energyNormMatrix contains the elasticity matrix _scaled with k_.
            // Therefore, divide by k here.
            preconditionedLocalSolvers[j].displacementHessian()[ii][jj] = ((1+k)/k)*energyNormMatrix[j][j][ii][jj];
      }

      localSolverVariant = preconditionedLocalSolvers;
    }

    std::visit(overload([&](std::monostate& impl) {}, [&](const std::monostate& impl) {}, [&](auto&& localSolver) {

    std::cout << "  ------  time step " << i << "  ------" << std::endl;

    /////////////////////////////
    //  Assemble load vector
    /////////////////////////////

    rhs = 0;

    xOld = x;
    
    // Dirichlet boundary value function at the current time
    auto currentBoundaryDeformation = [&boundaryDeformation, &i](const FieldVector<double,dim>& x)
    {
        return boundaryDeformation(x, i);
    };

    Functions::interpolate(deformationBasis,
                           x,
                           currentBoundaryDeformation,
                           dirichletNodes);

    Vector lowerObstacle(x.size());
    Vector upperObstacle(x.size());

    for (std::size_t j=0; j<x.size(); j++)
    {
      lowerObstacle[j] = -std::numeric_limits<double>::infinity();
      upperObstacle[j] =  std::numeric_limits<double>::infinity();
      lowerObstacle[j][dim] = xOld[j][dim];
      upperObstacle[j][dim] = 1.0;
    }

    /////////////////////////////
    //  Compute solution
    /////////////////////////////

    using QuadraticPartFunctional = TNNMG::BoxConstrainedQuadraticFunctional<Matrix,Vector,Vector,Vector,double>;
    auto quadraticPartFunctional = QuadraticPartFunctional(matrix,
                                                           rhs,
                                                           lowerObstacle,
                                                           upperObstacle);

    using DissipationFunctional = FracturePhasefields::IntegralFunctional<Vector, StrainValueMatrix, decltype(transposedStrainValueMatrix), Jet, GlobalQuadratureWeightVector, LocalEnergy, double>;
    auto dissipationFunctional = DissipationFunctional(strainValueMatrix, transposedStrainValueMatrix, globalQuadratureWeightVector, localEnergy);

    using Functional = TNNMG::SumFunctional<QuadraticPartFunctional,DissipationFunctional>;
    auto J = Functional(quadraticPartFunctional,dissipationFunctional);

    auto bisectionSolver = [](auto& xi, const auto& restriction, const auto& ignoreI) {

        ScalarSumFunctional<decltype(std::get<0>(restriction.functions())),
                          decltype(std::get<1>(restriction.functions()))> sumFunctional(std::get<0>(restriction.functions()), std::get<1>(restriction.functions()));
      Bisection bisection;

      xi = 1;
      int bisectionSteps;
      xi = bisection.minimize(sumFunctional, xi, xi, bisectionSteps);
    };

    using QuadraticPartLinearization = TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<QuadraticPartFunctional, BitVector>;

    QuadraticPartLinearization quadraticPartLinearization(quadraticPartFunctional);
    quadraticPartLinearization.enableGradientAwareTruncation(parameterSet.get<bool>("gradientAwareTruncation"));

    using DissipationLinearization = FracturePhasefields::IntegralFunctionalConstrainedLinearization<DissipationFunctional,
                                                                                                     Matrix,
                                                                                                     BitVector>;
    DissipationLinearization dissipationLinearization(dissipationFunctional);


    using Linearization = TNNMG::SumFunctionalConstrainedLinearization<Functional,
                                                                       BitVector,
                                                                       QuadraticPartLinearization,
                                                                       DissipationLinearization>;

    auto sumLinearization = std::make_shared<Linearization>(quadraticPartLinearization,
                                                            dissipationLinearization);

    using DefectProjection = DamageDefectProjection;
    using LineSearchSolver = decltype(bisectionSolver);

    using NonlinearSmoother = TNNMG::NonlinearGSStep<Functional, decltype(localSolver), BitVector>;

    auto nonlinearSmoother = std::make_shared<NonlinearSmoother>(J, x, localSolver);

    using Step = TNNMG::TNNMGStep<Functional,
                                  BitVector,
                                  Linearization,
                                  DefectProjection,
                                  LineSearchSolver>;

    auto step = Step(J, x, nonlinearSmoother, linearIterationStep, 1, DefectProjection(), bisectionSolver);
    step.setLinearization(sumLinearization);

    auto norm = EnergyNorm<Matrix, Vector>(energyNormMatrix);
    Solvers::LoopSolver<Vector> solver(step, maxIterations, tolerance, norm, Solver::FULL);

    step.setPreSmoothingSteps(nu1);

    step.setIgnore(ignore);

    solver.addCriterion(
            [&](){
                return Dune::formatString("   % 12d", step.linearization().truncated().count());
            },
            "   truncated   ");

    solver.addCriterion(
            [&](){
                return Dune::formatString("   % 12.5e", step.lastDampingFactor());
            },
            "   step size     ");

    solver.addCriterion(
            [&](){
                return Dune::formatString("   % 12.5e", J(x));
            },
            "   energy      ");

    solver.preprocess();
    Timer timer;

    solver.solve();

    std::cout << "Solving the spatial problem took " << timer.elapsed() << " seconds." << std::endl;

    // Write number of iterations and wall time to log files
    iterationsFile << solver.iterationCount()+1 << std::endl;
    wallTimeFile << timer.elapsed() << std::endl;


    /////////////////////////////
    //  Output result
    /////////////////////////////

    VTKWriter<GridView> vtkWriter(grid->leafGridView());

    // Write the displacement field
    auto deformationFunction = Functions::makeDiscreteGlobalBasisFunction<Coordinate>(deformationBasis, x);
    vtkWriter.addVertexData(deformationFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));

    // Write the damage field
    auto damageFunction = Functions::makeDiscreteGlobalBasisFunction<double>(damageBasis, x);
    vtkWriter.addVertexData(damageFunction, VTK::FieldInfo("damage", VTK::FieldInfo::Type::scalar, 1));

    // Actually write the file
    vtkWriter.write(resultPath + "brittle-fracture-result-" + std::to_string(i));

    ////////////////
    // Nodal force
    ////////////////
    
    // Vector marking all nodes from which the reaction forces are needed
    Dune::BlockVector<Dune::FieldVector<bool,dim+1>> reacNodes(gridView.size(dim));
    auto reacIndicatorFunction = Python::make_function<TupleVector<FieldVector<double,dim>, double>>(pyModule.get("reac_indicator"));
    Functions::interpolate(blockedBasis, reacNodes, reacIndicatorFunction);

    BitVector reacIgnore;
    reacIgnore.resize(reacNodes.size());
    for(std::size_t k=0; k<reacNodes.size(); ++k)
      for(std::size_t j=0; j<reacNodes[k].size(); ++j)
        reacIgnore[k][j] = 0;

    // Get Stiffness Matrix
    Linearization linearization(J,reacIgnore);
    linearization.bind(x);   //  x ist die Stelle, an der linearisiert wird (u & d)
    const auto& stiffness = linearization.hessian();
    //const auto& stress = linearization.negativeGradient();

    // Get only mechanical DOF for mechanical forces, damage = 0
    auto xmech=x;
    for(std::size_t k=0; k<xmech.size(); ++k)
        xmech[k][dim] = 0;

    // Compute dual reaction Forces F=K*x
    Vector force;
    Functions::istlVectorBackend(force).resize(blockedBasis);
    Functions::istlVectorBackend(force) = 0.0;

    stiffness.mv(xmech,force);

    // Sum up Force
    double sumForce = 0;
    double disp = 0;
    for(std::size_t k=0; k<reacNodes.size(); ++k)
      for(std::size_t j=0; j<dim; ++j) //sum up only mechanical forces
        if (reacNodes[k][j])           //get sumForce for nodes that are set in the input file
          {
            sumForce = sumForce + force[k][j];
            disp = x[k][j];
          }

    // write force into file
    reacForceFile << disp << "  " << sumForce << std::endl;
  }), localSolverVariant);
  }


  iterationsFile.close();
  reacForceFile.close();
  wallTimeFile.close();

} catch (Exception& e)
{
  std::cout << e.what() << std::endl;
}