brittle-fracture.cc 32.5 KB
Newer Older
1
2
3
4
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>

5
6
#include <fenv.h>

7
8
9
10
11
#include <dune/common/parametertree.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/parametertreeparser.hh>

#include <dune/grid/uggrid.hh>
12
#include <dune/grid/onedgrid.hh>
13
14
15
16
17
18
19
20
#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>

21
#include <dune/functions/functionspacebases/compositebasis.hh>
22
23
24
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
Gräser, Carsten's avatar
Gräser, Carsten committed
25
26
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
27
#include <dune/functions/functionspacebases/transformedindexbasis.hh>
28
#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
Gräser, Carsten's avatar
Gräser, Carsten committed
29
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
30
31
32

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

37
#include <dune/fracture-phasefields/brittlefractureenergynormassembler.hh>
Gräser, Carsten's avatar
Gräser, Carsten committed
38
#include <dune/fracture-phasefields/cracksurfaceassembler.hh>
39
40
#include <dune/fracture-phasefields/elementquantityassembler.hh>

41
#include <dune/solvers/solvers/loopsolver.hh>
42
#include <dune/solvers/solvers/umfpacksolver.hh>
43
44
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
45
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
46
#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
47
#include <dune/solvers/iterationsteps/blockgssteps.hh>
48

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

51
#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
52
#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
53
#include <dune/tnnmg/functionals/sumfunctional.hh>
54
#include <dune/tnnmg/functionals/sumfunctionalconstrainedlinearization.hh>
Sander, Oliver's avatar
Sander, Oliver committed
55
#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
56
57
58
#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
#include <dune/tnnmg/iterationsteps/tnnmgacceleration.hh>
#include <dune/tnnmg/problem-classes/bisection.hh>
59
#include <dune/tnnmg/projections/obstacledefectprojection.hh>
60

61
#include <dune/fracture-phasefields/brittlefracturejetevaluator.hh>
62
63
#include <dune/fracture-phasefields/gradientvaluejetevaluator.hh>
#include <dune/fracture-phasefields/gradienttostrainvaluemap.hh>
64
#include <dune/fracture-phasefields/brittlefracturelocalsolver.hh>
65
#include <dune/fracture-phasefields/brittlefracturepreconditionedlocalsolver.hh>
66
#include <dune/fracture-phasefields/integralfunctional.hh>
67
#include <dune/fracture-phasefields/integralfunctionalconstrainedlinearization.hh>
68
69
#include <dune/fracture-phasefields/damagedelasticenergydensity_decomposed.hh>
#include <dune/fracture-phasefields/damagedelasticenergydensity_nondecomposed.hh>
70
#include <dune/fracture-phasefields/integraldensity.hh>
71

72

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

76
//#define GAUSSSEIDEL
Gräser, Carsten's avatar
Gräser, Carsten committed
77
using namespace Dune;
78
79


80
81
82
83
84
85
/** \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>
86
87
88
class ScalarSumFunctional
{
public:
89
90
91
92
  ScalarSumFunctional(const DamageFunctional& damageFunctional,
                      const DamagedElasticityFunctional& damagedElasticityFunctional)
  : damageFunctional_(damageFunctional),
    damagedElasticityFunctional_(damagedElasticityFunctional)
93
94
    {}

95
  Solvers::Interval<double> domain() const
96
  {
97
98
    return Solvers::Interval<double>{damageFunctional_.lowerObstacle(),
                                     damageFunctional_.upperObstacle() };
99
100
  }

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

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

108
109
110
111
112
113
114
115
116
117
118
119
120
121
    // 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();

122
123
    // Smooth nonquadratic part
    dJ += damagedElasticityFunctional_.subDifferential(z);
124
125

    return dJ;
126
127
  }

128
129
  const DamageFunctional& damageFunctional_;
  const DamagedElasticityFunctional& damagedElasticityFunctional_;
130
131
132

};

133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
/* 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);
  }
};

149

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

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

159
160
  // Start Python interpreter
  Python::start();
161
162
163
164
165
166
167
  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]);
168

169
170
  // parse data file
  ParameterTree parameterSet;
171
  pyModule.get("parameterSet").toC(parameterSet);
172

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

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

180
181
  int numLevels         = parameterSet.get<int>("numLevels");
  int numTimeSteps      = parameterSet.get<int>("numTimeSteps");
182
183

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

187
188
189
190
191
  // 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);

192
193
194
195
196
197
198
199
200
201
202
203
  // 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", "");

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

207
208
209
210
211
212
213
214
215
216
217
218
219
220
  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.

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


224
225
226
227
228
  ///////////////////////////////
  //   Generate the grid
  ///////////////////////////////

  // The grid dimension
229
#ifndef BRITTLE_FRACTURE_DIMENSION
230
  const int dim = 2;
231
232
233
#else
  const int dim = BRITTLE_FRACTURE_DIMENSION;
#endif
234

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

238
  std::shared_ptr<GridType> grid;
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
  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");
254
    grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
255
256
257
  }

  typedef GridType::LeafGridView GridView;
Gräser, Carsten's avatar
Gräser, Carsten committed
258
259
  using Coordinate = Dune::FieldVector<double, dim>;

260
261
262
263
264
  GridView gridView = grid->leafGridView();

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

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

268
269
270
271
272
273
274

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

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

275
276
277
278
279
280
281
282
283
284
285
  // 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)
286
          return basis.size(Dune::ReservedVector<std::size_t,1>({0}));
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
        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);
304

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

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

Gräser, Carsten's avatar
Gräser, Carsten committed
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
  //////////////////////////////////////////////////////////////////////////////////
  //   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);
    }
  }

350
351
352
353
354
355
  //////////////////////////////////////////////////////////////////////////////////
  //   Compute an H^1/L^2 matrix to measure the convergence with
  //////////////////////////////////////////////////////////////////////////////////

  auto energyNormMatrix = Matrix();
  FracturePhasefields::BrittleFractureEnergyNormAssembler<decltype(blockedBasis)>
356
357
358
359
      localEnergyNormAssembler(blockedBasis,
                               crackSurfaceDensity,
                               parameterSet.sub("materialParameters"),
                               k);   // Scaling for the elasticity part
360
361
362
363
364
  auto energyNormAssembler = Fufem::duneFunctionsOperatorAssembler(blockedBasis, blockedBasis);

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

365
366
367
368
  //////////////////////////////////////////////////////////////////////////////////
  //   Compute a vector that contains a the integrals over each shape function
  //////////////////////////////////////////////////////////////////////////////////

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

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

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

377
  BlockVector<double> nodalWeights;
378

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

383
384
385
386
387
  //////////////////////////////////////////////////////////////
  //  Create nonlinear coupling functional
  //////////////////////////////////////////////////////////////

  // number of quadrature points per element
Sander, Oliver's avatar
Sander, Oliver committed
388
  constexpr std::size_t maxQuadPointsSelect[3]={2,4,8};
Daniel Kienle's avatar
Daniel Kienle committed
389
390
  static const std::size_t maxQuadPoints=maxQuadPointsSelect[dim-1];

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

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

397
398
399
400
  using SmallJet = decltype(gradientValueJetEvaluator)::Jet;
  auto gradientValueAssembler = elementQuantityAssembler<SmallJet, maxQuadPoints>(scalarBasis,
                                                                                  gradientValueJetEvaluator,
                                                                                  quadOrder);
401

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

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

410
411
412
413
414
  // 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);
415

416
417
418
419
420
421
422
423
  // 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);
424

425
426
427
428
  //////////////////////////////////////////////////////////////
  //  Create initial iterate
  //////////////////////////////////////////////////////////////

429
  using Vector = BlockVector<FieldVector<double,dim+nDamageComponents> >;
430
431
432

  // Initial iterate
  Vector x;
Gräser, Carsten's avatar
Gräser, Carsten committed
433
434
  Functions::istlVectorBackend(x).resize(blockedBasis);
  Functions::istlVectorBackend(x) = 0.0;
Gräser, Carsten's avatar
Gräser, Carsten committed
435
436

  
437
  // Compute initial value by evaluating a function given in python
438
  auto initial_damage = Python::make_function<double>(pyModule.get("initial_damage"));
439
  Functions::interpolate(damageBasis, x, initial_damage);
Gräser, Carsten's avatar
Gräser, Carsten committed
440
441

//  PythonFunction<Coordinate, Coordinate> boundary_deformation(pyModule.get("boundary_deformation"));
442
//  Functions::interpolate(deformationBasis, x, boundary_deformation);
Gräser, Carsten's avatar
Gräser, Carsten committed
443

444
445
446
447
  //////////////////////////////
  //   Create a solver
  //////////////////////////////

Gräser, Carsten's avatar
Gräser, Carsten committed
448
449
450
  Vector rhs;
  Functions::istlVectorBackend(rhs).resize(blockedBasis);
  Functions::istlVectorBackend(rhs) = 0.0;
451

452
453
  using BitVector = Solvers::DefaultBitVector_t<Vector>;

454
455
456
457
  // first create a linear iteration step
  // Gauss-Seidel is the base solver

#ifdef GAUSSSEIDEL
458
  TruncatedBlockGSStep<Matrix,Vector> linearBaseSolverStep;
459

460
  EnergyNorm<Matrix,Vector> baseEnergyNorm(linearBaseSolverStep);
461
  Solvers::LoopSolver<Vector> linearBaseSolver(linearBaseSolverStep,
462
463
                                          1000,
                                          baseTolerance,
464
                                          baseEnergyNorm,
465
466
                                          Solver::QUIET);
#else
467
468
469
  // The matrix of the linear correction problem is symmetric,
  // but possibly indefinite.
  Solvers::UMFPackSolver<Matrix,Vector> linearBaseSolver;
470
471
#endif

472
  std::shared_ptr<LinearIterationStep<Matrix,Vector> > linearIterationStep;
473
474
475

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

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

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

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

    multigridStep->setMGType(1, nu1, nu2);
Daniel Kienle's avatar
Daniel Kienle committed
487
    multigridStep->setBaseSolver(linearBaseSolver);
488
489
490
491
492
    multigridStep->setSmoother(linearSmoother);

    multigridStep->mgTransfer_.resize(numLevels-1);
    for (size_t i=0; i<multigridStep->mgTransfer_.size(); i++)
    {
493
      CompressedMultigridTransfer<Vector>* newTransfer = new CompressedMultigridTransfer<Vector>;
494
      newTransfer->setup(*grid, i, i+1);
Daniel Kienle's avatar
Daniel Kienle committed
495
      multigridStep->mgTransfer_[i] = Dune::stackobject_to_shared_ptr(*newTransfer);
496
497
498
499
500
501
502
503
    }

    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.
Gräser, Carsten's avatar
Gräser, Carsten committed
504
  Vector coefficients(x.size());
505
  for (size_t i=0; i<nodalWeights.size(); i++)
506
    coefficients[i] = nodalWeights[i];
507

508
509
510
511
512
513
514
  // 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];

515
516
517
518
519

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

520
  using Jet = FieldVector<double, dim*(dim+1)/2 + 1>;
521
  using LocalEnergy = FracturePhasefields::IntegralDensity<Jet>;
522
523
   //  small strain non decomposed (default)
       using LocalEnergyNonDecomposed = FracturePhasefields::NonDecomposed::DamagedElasticEnergyDensity<dim>;
524
  auto localEnergy = LocalEnergy(LocalEnergyNonDecomposed(materialParameters));
Sander, Oliver's avatar
Sander, Oliver committed
525

526
527
528
529
530
531
532
533
534
   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>;
535
       localEnergy = LocalEnergy(LocalEnergyDecomposed(materialParameters));
536
537
       break;

538
539
    default:
      DUNE_THROW(NotImplemented, "Unknown material model type");
540
541
   }

542
543
544
545
  ////////////////////////////////////////////////////////
  //  Time loop
  ////////////////////////////////////////////////////////

Gräser, Carsten's avatar
Gräser, Carsten committed
546
  Vector xOld;
547
548
549
550

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

551
552
553
554
555
  // 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);

556
557
  for (int i=0; i<numTimeSteps; i++)
  {
558
559
    // Construct the local solver / the smoother
    using NewtonLocalSolver = FracturePhasefields::BrittleFractureLocalSolver<dim>;
560
    using PreconditionedLocalSolver = std::vector<FracturePhasefields::BrittleFracturePreconditionedLocalSolver<dim>>;
561
562
563
564
    auto localSolverVariant = std::variant<std::monostate, NewtonLocalSolver, PreconditionedLocalSolver>{};

    if (localSolverType == "exact")
      localSolverVariant = FracturePhasefields::BrittleFractureLocalSolver<dim>(outerIterations,displacementIterations);
565
566
567
568
569
570
571
572
573
574
    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++)
575
576
577
578
            // (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];
579
580
581
582
      }

      localSolverVariant = preconditionedLocalSolvers;
    }
583
584
585

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

586
587
588
589
590
591
    std::cout << "  ------  time step " << i << "  ------" << std::endl;

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

Gräser, Carsten's avatar
Gräser, Carsten committed
592
    rhs = 0;
593

Gräser, Carsten's avatar
Gräser, Carsten committed
594
    xOld = x;
595
    
596
597
    // Dirichlet boundary value function at the current time
    auto currentBoundaryDeformation = [&boundaryDeformation, &i](const FieldVector<double,dim>& x)
Sander, Oliver's avatar
Sander, Oliver committed
598
    {
599
        return boundaryDeformation(x, i);
600
    };
Daniel Kienle's avatar
Daniel Kienle committed
601

602
    Functions::interpolate(deformationBasis,
603
                           x,
604
                           currentBoundaryDeformation,
605
                           dirichletNodes);
Daniel Kienle's avatar
Daniel Kienle committed
606

607
608
    Vector lowerObstacle(x.size());
    Vector upperObstacle(x.size());
Gräser, Carsten's avatar
Gräser, Carsten committed
609

610
    for (std::size_t j=0; j<x.size(); j++)
Gräser, Carsten's avatar
Gräser, Carsten committed
611
    {
612
613
614
615
      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;
Gräser, Carsten's avatar
Gräser, Carsten committed
616
    }
617
618
619
620
621

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

622
    using QuadraticPartFunctional = TNNMG::BoxConstrainedQuadraticFunctional<Matrix,Vector,Vector,Vector,double>;
623
624
    auto quadraticPartFunctional = QuadraticPartFunctional(matrix,
                                                           rhs,
625
626
                                                           lowerObstacle,
                                                           upperObstacle);
627

628
629
    using DissipationFunctional = FracturePhasefields::IntegralFunctional<Vector, StrainValueMatrix, decltype(transposedStrainValueMatrix), Jet, GlobalQuadratureWeightVector, LocalEnergy, double>;
    auto dissipationFunctional = DissipationFunctional(strainValueMatrix, transposedStrainValueMatrix, globalQuadratureWeightVector, localEnergy);
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644

    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);
    };

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

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

650
651
652
    using DissipationLinearization = FracturePhasefields::IntegralFunctionalConstrainedLinearization<DissipationFunctional,
                                                                                                     Matrix,
                                                                                                     BitVector>;
653
654
    DissipationLinearization dissipationLinearization(dissipationFunctional);

655

656
657
658
659
    using Linearization = TNNMG::SumFunctionalConstrainedLinearization<Functional,
                                                                       BitVector,
                                                                       QuadraticPartLinearization,
                                                                       DissipationLinearization>;
660
661
662
663

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

664
    using DefectProjection = DamageDefectProjection;
665
666
    using LineSearchSolver = decltype(bisectionSolver);

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

Gräser, Carsten's avatar
Gräser, Carsten committed
669
    auto nonlinearSmoother = std::make_shared<NonlinearSmoother>(J, x, localSolver);
670

Sander, Oliver's avatar
Sander, Oliver committed
671
672
673
674
675
676
    using Step = TNNMG::TNNMGStep<Functional,
                                  BitVector,
                                  Linearization,
                                  DefectProjection,
                                  LineSearchSolver>;

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

680
    auto norm = EnergyNorm<Matrix, Vector>(energyNormMatrix);
681
    Solvers::LoopSolver<Vector> solver(step, maxIterations, tolerance, norm, Solver::FULL);
682
683
684
685
686
687
688

    step.setPreSmoothingSteps(nu1);

    step.setIgnore(ignore);

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

    solver.addCriterion(
            [&](){
Sander, Oliver's avatar
Sander, Oliver committed
695
                return Dune::formatString("   % 12.5e", step.lastDampingFactor());
696
697
698
699
700
            },
            "   step size     ");

    solver.addCriterion(
            [&](){
Gräser, Carsten's avatar
Gräser, Carsten committed
701
                return Dune::formatString("   % 12.5e", J(x));
702
703
704
705
            },
            "   energy      ");

    solver.preprocess();
706
    Timer timer;
707
708
709

    solver.solve();

710
711
712
713
714
    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;
715

716

717
718
719
720
721
722
    /////////////////////////////
    //  Output result
    /////////////////////////////

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

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

727
    // Write the damage field
728
    auto damageFunction = Functions::makeDiscreteGlobalBasisFunction<double>(damageBasis, x);
729
    vtkWriter.addVertexData(damageFunction, VTK::FieldInfo("damage", VTK::FieldInfo::Type::scalar, 1));
Gräser, Carsten's avatar
Gräser, Carsten committed
730

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

Daniel Kienle's avatar
Daniel Kienle committed
734
735
736
    ////////////////
    // Nodal force
    ////////////////
737
    
738
    // Vector marking all nodes from which the reaction forces are needed
739
    Dune::BlockVector<Dune::FieldVector<bool,dim+1>> reacNodes(gridView.size(dim));
740
    auto reacIndicatorFunction = Python::make_function<TupleVector<FieldVector<double,dim>, double>>(pyModule.get("reac_indicator"));
741
    Functions::interpolate(blockedBasis, reacNodes, reacIndicatorFunction);
Daniel Kienle's avatar
Daniel Kienle committed
742
743
744
745
746
747
748
749
750

    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);
751
    linearization.bind(x);   //  x ist die Stelle, an der linearisiert wird (u & d)
Daniel Kienle's avatar
Daniel Kienle committed
752
753
754
    const auto& stiffness = linearization.hessian();
    //const auto& stress = linearization.negativeGradient();

755
756
757
758
759
    // 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;

Daniel Kienle's avatar
Daniel Kienle committed
760
761
    // Compute dual reaction Forces F=K*x
    Vector force;
Gräser, Carsten's avatar
Gräser, Carsten committed
762
763
    Functions::istlVectorBackend(force).resize(blockedBasis);
    Functions::istlVectorBackend(force) = 0.0;
764
765

    stiffness.mv(xmech,force);
Daniel Kienle's avatar
Daniel Kienle committed
766
767

    // Sum up Force
768
    double sumForce = 0;
769
    double disp = 0;
Daniel Kienle's avatar
Daniel Kienle committed
770
    for(std::size_t k=0; k<reacNodes.size(); ++k)
771
772
      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
773
774
775
776
          {
            sumForce = sumForce + force[k][j];
            disp = x[k][j];
          }
777

Daniel Kienle's avatar
Daniel Kienle committed
778
    // write force into file
779
    reacForceFile << disp << "  " << sumForce << std::endl;
780
  }), localSolverVariant);
781
782
  }

783
784
785
786
787

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

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