Commit 36f1462e authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Allow to compute local gradient and Hesse matrix together

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.
parent 30bc76c2
......@@ -25,8 +25,7 @@ class BrittleFractureLocalSolver
auto negativeGradient = derivative(quadraticFunctional)(x);
negativeGradient *= -1;
auto hessian = quadraticFunctional.quadraticPart();
negativeGradient -= nonlinearFunctional.gradient(x, 0);
nonlinearFunctional.addHessian(x, hessian, 0, 0);
nonlinearFunctional.addNegativeGradientHessian(x, negativeGradient, hessian, 0, 0);
// extract upper left block
FieldMatrix<double, dim, dim> H00(0);
......
......@@ -204,6 +204,60 @@ public:
}
}
/** \brief Add only a sub-block of the Hessian
*
* Block decomposition of the Hessian depends on the integral density
*
* \param row Block row of the gradient vector and Hesse matrix block to be added
* \param col Block column of the Hesse matrix block to be added
*/
template<class Hessian>
void addNegativeGradientHessian(const Vector& v, Vector& negativeGradient, Hessian& hessian,
size_t row, size_t col) const
{
const auto& BT = f_->originalFunctional().transposedEvaluationMatrix();
const auto& weights = quadratureWeights();
const auto& densityGradient = derivative(density());
const auto& densityHessian = derivative(derivative(density()));
for (auto&& [BT_ij, j] : sparseRange(BT[index()]))
{
auto jet_j = f_->jet()[j];
BT_ij.umtv(v, jet_j);
// This is the gradient part
auto g_j = densityGradient(jet_j, row);
BT_ij.usmv(-1.0 * weights[j], g_j, negativeGradient);
// This is the Hessian part
auto Hj = densityHessian(jet_j, row, col);
Hj *= weights[j];
// compute BT_ij * Hj * (BT_ij)^T
if constexpr (Std::is_detected_v<Impl::isUnblockedDiagonalMatrix_t, decltype(BT_ij)>)
BT_ij.leftRightMultiply(Hj, hessian, row, col);
else
{
auto rowRange = densityHessian.getBlockRange(row);
auto colRange = densityHessian.getBlockRange(col);
for (size_t k : rowRange)
{
for (size_t l : colRange)
{
auto&& BT_ij_k = BT_ij[k];
auto&& BT_ij_l = BT_ij[l];
for (auto&& [Hj_r, r] : sparseRange(Hj))
for (auto&& [Hj_rs, s] : sparseRange(Hj_r))
hessian[k][l] += BT_ij_k[r] * Hj_rs * BT_ij_l[s];
}
}
}
}
}
protected:
const GlobalFunctional* f_;
......
......@@ -76,6 +76,17 @@ public:
addHessian(x, hessian, 1, 1);
}
void addNegativeGradientHessian(const FieldVector<double,size>& x,
FieldVector<double,size>& negativeGradient,
FieldMatrix<double,size,size>& hessian,
std::size_t row, std::size_t col) const
{
auto grad = gradient(x, row);
negativeGradient -= grad;
addHessian(x, hessian, row, col);
}
FieldVector<double,size> center_;
double p_;
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment