Moving to CHOLMOD was a mistake:  After all, the brittle-fracture
increment functional is non-convex, and hence second derivatives
can be indefinite.

Rather, the only the gradients of scalar shape functions, and compute
the strains on the fly.  This makes the evaluation matrices quite
a bit smaller (roughly by a factor dim).

This seems to reduce the number of TNNMG iterations in a few cases.

In the first version of the manuscript, we used to control convergence
with an energy norm where the elasticity part was scaled with the
residual stiffness 'k'.  As far as I remember it this wasn't really
intentional, but the matrix for that norm just "happened to be
available".

When we started to implement the operator-splitting algorithm we wanted
to consolidate the convergence criterion, and we replaced the scaled
energy norm by the true energy norm.  As it turned out, this
overestimates the error in some situations where there is a complete
crack.  This is particularly noticeable for the models with spectrally
split energy, for reasons that we don't quite understand.

This patch puts the scaling back into the norm.  From a theoretical
point of view the scaled energy norm is not worse than the unscaled
one, but the iteration numbers look much better this way.

That's too much magic for my taste:  If no material model type
is provided in the parameter file and the command line options,
then throw an error.

This allows to get rid of the PDElab assembler.

That factor got lost when I rewrote the preconditioned smoother.
The +k is very small, but it is required for the convergence
theory.  Without it, the quadratic model used by the preconditioner
is not dominating.

They are symmetric after all.

This does not lead to a measurable speedup for the benchmark problems
in the Gräser/Kienle/Sander paper, but it avoids having to explain
why we are using a nonsymmetric solver for symmetric problems.

Previously the code contained the computation of the residual stiffness
twice: Once in the damaged-elastic-density-classes and once in
BrittleFractureLocalOperator.  That was confusing and error-prone:
One had to take good care to make sure that the residual stiffness
was computed only once.

This patch simplifies the situation by removing the residual stiffness
from the BrittleFractureLocalOperator class.  It is now always
computed by the damaged elastic energy density.  This also makes more
sense conceptually, and it should influence the run-time either.

The preconditioned local solver needs the local block of the
linear elasticity matrix, for use as its quadratic model.
Previously, this matrix was extracted from the global
Hesse matrix, with some tricks to account for the scaling
of that global matrix with the residual stiffness k.

With this patch, there is now a dedicated local solver object
for each vertex, and each objects stores the correct local
matrix.  This simplifies the code a little.

This is the data to be shown in the paper.

Previously the reaction force while would simply go to the
current directory.

I usually archive the complete shell output of program runs
used for a paper.  With this patch this output includes all
parameters of the run.

IMO this makes automated parameter studies easier.

The BrittleFractureLocalOperator class still used the old
formula with 'psi_c', which lead to subtle bugs in various
situations.  In other parts of this module it has already
been replaced by the CrackSurfaceDensity class.
This commit completes the conversion of the module.

Warning: As some of the psi_c settings were computed
incorrectly in the parameter files, this patch may slightly
alter your problem specifications!

Having the residual stiffness 'k' in two variables k and k_local
was a constant source of confusion.  This patch eliminates
the k_local variable, and makes sure the correct values
are used.

This leads to a chance compared to previous behavior.  Previously,
the same norm was used for convergence testing, with the difference
that the displacement part was scaled with k (i.e., something
very small).  The new norm gives equal importance to displacement
and damage parts.

This is one step towards getting rid of the dune-pdelab dependency.

Because it is not independent from g_c, and that is just
very confusing.

This class is problem-independent, and we therefore move it
to dune-tnnmg.  With this patch, the new version from dune-tnnmg
is used right away.

For compatibility reasons the old target still works.

In calls to Functions::interpolate it is not necessary anymore.

In cases where it is necessary, it should be (and is) replaced
by Functions::istlVectorBackend.

It was always supposed to be in dune-functions anyway, and now it is
there.

This requires that the python functions range
has a suitable blocking structure. Furthermore you
need https://git.imp.fu-berlin.de/agnumpde/dune-fufem/-/merge_requests/82
in order to convert from python to Dune::TupleVector