Remove manual implementation of operator* for FieldMatrix

See merge request !56

This is now handled by `operator*` in `FieldMatrix` based on `mtv()`.

Method mv for unblockeddiagonalmatrix

See merge request !54

Remove computation of the nodal weights

See merge request !55

Switch to truncated GS and euclidean norm

I don't see them used anywhere.

Move class ScalarSumFunctional to a separate file

See merge request !53

This is a first step towards making it more general.  There is a copy
of it in ductile-fracture.cc, which I would like to eliminate.

Add damping to local Newton solver

See merge request !52

To avoid numerical instability we in fact allow
that the energy increase, if the relative in crease
is less than threshold 1-e14.

This is now only used in a single place hence having a lambda
for this is overkill.

This ensures energy-monotonicity of the smoother
with local Newton-step

Ensure and track monotonicity

See merge request !51

This is a collection of many interacting patches that help to ensure
monotonicity of the TNNMG iteration and add more diagnostic output to
track the latter:

* Implement `ScalarSumFunctional::operator()` in a way that is
  consistent with the global energy evaluation. Unfortunately
  this means we cannot use the evaluation of the individual
  constrained functionals.
* Implement `derivative(ScalarSumFunctional)`. The functional
  is in general not differentiable, hence we return the minimal
  norm entry of the subdifferential or, equivalently, the projection
  of zero into the subdifferential. Finding a root of the
  subdifferential is equivalent to find a root of this single-valued
  map. Furthermore a subdifferential-aware bisection is equivalent
  to a bisection for the single-valued map.
* To ensure monotonicity of the coarse grid correction add a
  postprocessing step to the line search: After an approximate
  root of the derivative has been found, check if this decreases energy.
  If not, then half the step size until it does. Since we have a descent
  direction this must terminate. This step is necessary because the
  functional is not convex.
* When comparing energies take care to always evaluate it exactly the
  same way, to avoid non-monotonicity due to inconsistent evaluation
  algorithms.
* Track energy decrease of nonlinear smoother and coarse grid
  correction. While the latter is guaranteed by the above changes,
  it may happen that the smoother is not monotone - either because
  of local solves with uncontrolled inexactness (called exact in our terminalogy, arg...)
  or because of round-off-errors.
* Add more diagnostic output to the table:
  * Number of subdifferential evaluations during line search.
    Hint on interpretation: This will be at least 2+#bisectionSteps.
  * Number of functional evaluations during line search.
    Hint on interpretation: This will be 2+#reductionSteps,
    where #reductionSteps is the number of halving steps done
    to ensure monotonicity (see above).
  * Energy after TNNMG step with 16 digits.
  * Energy after nonlinear pre-smoothing with 16 digits.
  * If the nonlinear pre-smoother is not monotone:
    Relative monotonicity violation of the smoother with an additional
    `X` marker for easier grepping.

* Export origin and direction
* Cache and use evaluation matrix times origin
* Implement `operator()`

Modularize the Decomposed::DamagedElasticEnergyDensity class

See merge request !50

Previously, we would use the CrackSurfaceAssembler class to assemble
the linear term of the crack surface density functional that appears
in the AT-1 model, and then throw it away.  The actual computation
would get the term from the DamagedElasticEnergyDensity implementations.

However, conceptually that linear term doesn't really have anything
to do the the damaged elastic energy density.  Therefore, as a
cleanup measure, remove the AT-1 term from the damaged elastic energy
density implementations.

The current Decomposed::DamagedElasticEnergyDensity class does too many
things together:
a) It implements the St.Venant-Kirchhoff density in terms of the strain
   eigenvalues, split into tensile and compressive parts as described
   by Miehe, Welschinger, Hofacker.
b) It degrades the tensile part with a hard-wired degradation function.
c) It transforms the degraded density and its gradient and Hesse matrix
   from spectral to physical coordinates.

This patch modularizes a): It introduces a separate class MieheSplitting
that only computes the positive and negative parts of the St.Venant-
Kirchhoff density in spectral coordinates (and the gradient and
Hesse matrix).

The old version was to pessimistic and led to numerical errors
if the domain was very large (e.g., because all constrained dofs
have been zero).

Patches to speed up TNNMG for materials with a spectral split.

See merge request !48

Instead, simply sum the eigenvalues to get the trace
of the strain matrix.

... and use it in IntegralFunctionalConstrainedLinearization::bind.

This saves another 5% for the model with the spectrally split
energy density.

Previously, computing the gradient and the Hesse matrix of the energy
density where two separate operations.  For the spectrally split
density this was wasteful, because each call involved an expensive
eigenvector decomposition.  This commit introduces a new method
called 'gradientHessian' which returns both first and second
derivatives together.  That way, the eigendecomposition is
done only once, and run-time is saved.

This patch does a lot of refactoring which, to my eyes,
makes the code easier to understand.

Also, it adds documentation for the transformation rules
for spectral functions.  These rules are nontrivial, and we
need to state where they were taken from.

This patch adds the references to a paper by
Lewis and Sendov, which derives and explains the formulas.

This can be considerably more efficient, because the two methods
sometimes share a considerable amount of code.

The prime example for this is the energy density of a degraded
elastic material with a spectral split.  In that situation,
both the gradient and the Hesse matrix implementations need
an eigenvector decomposition of the strain, which is quite
expensive.  With the new code, the eigenvector decomposition
is computed only once and used for both the gradient and
the Hesse matrix.

For now it's in the cc file, because I didn't come up
with a nice name. The new version also exactly documents
what it does.

As the nonlinear smoother is so expensive, doing more than the
canonical single multigrid iteration for the linear correction
problem does not appreciably increase the time per iteration.

In preliminary tests it does decrease the number of iterations
considerably, though.

It is not worth doing this change for the efficiency gain alone
(it is not measurable).  But the new code is simpler, too.

It is never used.

Measurements indicate speedup around 15% for a 3d brittle-fracture
simulation.

Make TNNMG faster

See merge request !47

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

The blocks of the evaluation matrix are sparse.  In particular,
they have a block-diagonal structure.  Previously, the implementation
used a FieldMatrix for these blocks.  Which, in theory, is wasteful
because of the matrix sparsity.

This commit introduces a block-diagonal matrix with compile-time size
for this.

Unfortunately, initial tests do not show any speedup.