Commit 7e59fd36 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Add non-blocked Density::gradientHessian method

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

This saves another 5% for the model with the spectrally split
energy density.
parent d5d09ccd
......@@ -459,6 +459,81 @@ public:
return energy;
}
/** \brief Compute gradient and Hesse matrix together
*
* \return std::tuple of FieldVector for the gradient and FieldMatrix for the Hesse matrix
*/
auto gradientHessian(const FieldVector<RangeField, size>& state) const
{
std::tuple<typename DamagedElasticEnergyDensityGradient::Range, typename DamagedElasticEnergyDensityHessian::Range> result;
auto& gradient = std::get<0>(result);
auto& hessian = std::get<1>(result);
const auto [eigenValues, eigenVectors] = eigenVectorDecomposition(state);
// Compute gradient
// Displacement block
FieldVector<RangeField,dim> sigP, sigM;
const auto gradientSpectral = derivative(*this).damagedSigSpec(state,eigenValues,sigP,sigM);
// transform sigma in physical space
gradient = DamagedElasticEnergyDensityGradient::spectralToPhysical(eigenVectors, gradientSpectral);
// Damage block
const auto [psiP, psiM] = psiSpec(state,eigenValues);
gradient[size-1] = -2.0 * (1.0 - state[size-1])*psiP;
// linear part of critical energy density for model with elastic phase
if (crackSurfaceDensity_.type() == CrackSurfaceDensity::AT1)
gradient[size-1] += crackSurfaceDensity_.g_c() / (4.0*crackSurfaceDensity_.c_w()*crackSurfaceDensity_.l());
// Compute Hesse matrix
// Displacement-displacement block
// The Hesse matrix in spectral space
FieldMatrix<double, dim,dim> hessianSpecP;
FieldMatrix<double, dim,dim> hessianSpecM;
const double trEps = traceEps(state);
if (trEps>0) {
hessianSpecP = lambda_;
hessianSpecM = 0;
} else {
hessianSpecP = 0.0;
hessianSpecM = lambda_;
}
for (std::size_t i=0;i < dim; ++i)
{
if (eigenValues[i]>0 )
hessianSpecP[i][i] += 2.0*mu_;
else
hessianSpecM[i][i] += 2.0*mu_;
}
// degrade only tensile part
FieldMatrix<double, dim,dim> hessianSpec = hessianSpecP * degradation(state)
+ hessianSpecM * (1.0 + k_);
hessian = DamagedElasticEnergyDensityHessian::spectralToPhysical(eigenValues, eigenVectors, gradientSpectral, hessianSpec);
// Displacement-damage and damage-displacement blocks
auto undamagedGradientP = DamagedElasticEnergyDensityGradient::spectralToPhysical(eigenVectors,sigP);
for (std::size_t i=0;i < size-1; ++i)
hessian[i][size-1] = hessian[size-1][i] = -2.0 * (1 - state[size-1])*undamagedGradientP[i];
// Damage-damage block
hessian[size-1][size-1] = 2.0*psiP;
return result;
}
/** \brief Compute a block of gradient and Hesse matrix together
*
* \return std::tuple of FieldVector for the gradient and FieldMatrix for the Hesse matrix
*/
auto gradientHessian(const FieldVector<RangeField, size>& state, size_t row, size_t col) const
{
if (row!=col) {
......
......@@ -337,6 +337,85 @@ public:
return energy;
}
/** \brief Compute gradient and Hesse matrix together
*
* \return std::tuple of FieldVector for the gradient and FieldMatrix for the Hesse matrix
*/
auto gradientHessian(const FieldVector<RangeField, size>& state) const
{
std::tuple<typename DamagedElasticEnergyDensityGradient::Range, typename DamagedElasticEnergyDensityHessian::Range> result;
auto& gradient = std::get<0>(result);
auto& hessian = std::get<1>(result);
gradient = 0.0;
hessian = 0.0;
const double degradationValue = degradation(state);
const double trEps = traceEps(state);
FieldVector<double,size-1> II;
switch (dim)
{
case 1:
II = {1.0};
break;
case 2:
II = {1.0, 2.0, 1.0};
break;
case 3:
II = {1.0, 2.0, 1.0, 2.0, 2.0, 1.0};
break;
default:
DUNE_THROW(RangeError, "Not implemented for dim = " << dim);
}
// Compute gradient
for (auto&& diagonal : diagonals())
gradient[diagonal] = degradationValue * lambda_ * trEps;
for (std::size_t i=0; i<size-1; i++)
gradient[i] += degradationValue * 2.0 * mu_ * state[i] * II[i];
gradient[size-1] = -2.0 * (1.0 - state[size-1]) * psi_0(state);
// linear part of critical energy density for model with elastic phase
if (crackSurfaceDensity_.type() == CrackSurfaceDensity::AT1) // psi_c = g_c / (4*c_w*l)
gradient[size-1] += crackSurfaceDensity_.g_c() / (4.0*crackSurfaceDensity_.c_w()*crackSurfaceDensity_.l());
// Compute Hesse matrix
// Displacement-displacement coupling
for (auto&& d1 : diagonals())
for (auto&& d2: diagonals())
hessian[d1][d2] += lambda_;
for (std::size_t i=0; i<size-1; ++i)
hessian[i][i] += 2.0 * mu_ * II[i];
hessian *= degradationValue;
// Displacement-damage coupling
FieldVector<RangeField, size> stresses(0.0);
for (auto&& diagonal : diagonals())
stresses[diagonal] = lambda_*trEps;
for (std::size_t i=0; i<size-1; i++)
stresses[i] += 2.0 * mu_*state[i];
for (std::size_t i=0; i<size-1; ++i)
hessian[i][size-1] = hessian[size-1][i] = -2.0*(1 - state[size-1])*stresses[i] * II[i];
// Damage-damage coupling
hessian[size-1][size-1] = 2.0 * psi_0(state);
return result;
}
/** \brief Compute a block of gradient and Hesse matrix together
*
* \return std::tuple of FieldVector for the gradient and FieldMatrix for the Hesse matrix
*/
auto gradientHessian(const FieldVector<RangeField, size>& state, size_t row, size_t col) const
{
if (row!=col) {
......
......@@ -204,6 +204,12 @@ public:
return (*pimpl_)(jet);
}
auto gradientHessian(const Jet& jet) const
{
return pimpl_->gradientHessian(jet);
}
auto gradientHessian(const Jet& jet, size_t row, size_t col) const
{
return pimpl_->gradientHessian(jet, row, col);
......@@ -214,6 +220,7 @@ private:
{
virtual ~Concept() {}
virtual RangeField operator() (const Jet& jet) const = 0;
virtual std::tuple<typename IntegralDensityGradient<Jet>::Range, typename IntegralDensityHessian<Jet>::Range> gradientHessian(const Jet& jet) const = 0;
virtual std::tuple<typename IntegralDensityGradient<Jet>::Range, typename IntegralDensityHessian<Jet>::Range> gradientHessian(const Jet& jet, size_t row, size_t col) const = 0;
virtual IntegralDensityGradient<Jet> derivative_() const = 0;
virtual Concept* clone() const = 0;
......@@ -224,6 +231,7 @@ private:
{
Model( T const& value ) : object( value ) {}
RangeField operator() (const Jet& jet) const override {return object(jet); }
std::tuple<typename IntegralDensityGradient<Jet>::Range, typename IntegralDensityHessian<Jet>::Range> gradientHessian(const Jet& jet) const override {return object.gradientHessian(jet);}
std::tuple<typename IntegralDensityGradient<Jet>::Range, typename IntegralDensityHessian<Jet>::Range> gradientHessian(const Jet& jet, size_t row, size_t col) const override {return object.gradientHessian(jet,row,col);}
IntegralDensityGradient<Jet> derivative_() const override {return IntegralDensityGradient<Jet>(derivative(object)); }
Concept* clone() const override{ return new Model( object ); }
......
......@@ -72,23 +72,25 @@ namespace Dune::FracturePhasefields {
*/
void bind(const Vector& x)
{
// Compute first derivate of the functional
negativeGradient_ = f_.gradient(x);
negativeGradient_ *= -1;
typename LocalEnergy::Domain BiXv;
// Compute second derivative of the functional
// Compute first and second derivate of the functional
Solvers::resizeInitializeZero(negativeGradient_, x);
hessian_ = 0;
typename LocalEnergy::Domain BiXv;
for (std::size_t i=0; i<f_.evaluationMatrix().N(); ++i)
{
BiXv = 0.0;
BiXv = 0.0;
addRowXVector(f_.evaluationMatrix()[i], x, BiXv);
auto densityHessian = derivative(derivative(f_.density()));
auto DDgamma = densityHessian(BiXv);
auto [Dgamma, DDgamma] = f_.density().gradientHessian(BiXv);
Dgamma *= -1.0 * f_.quadratureWeights()[i];
DDgamma *= f_.quadratureWeights()[i];
addRowTransposedXVector(f_.evaluationMatrix()[i],
Dgamma,
negativeGradient_);
addRowTransposedXMatrixXRow(f_.evaluationMatrix()[i],
DDgamma,
hessian_);
......@@ -156,6 +158,14 @@ namespace Dune::FracturePhasefields {
entry.umv(x[i], y);
}
static void addRowTransposedXVector(const typename EvaluationMatrix::row_type& evaluationMatrixRow,
const typename LocalEnergy::Domain& x, // Actually: the range of the gradient
Vector& y)
{
for (auto [entry, i] : sparseRange(evaluationMatrixRow))
entry.umtv(x, y[i]);
}
static void addRowTransposedXMatrixXRow(const typename EvaluationMatrix::row_type& row,
const FunctionalHessian& X,
Matrix& Y)
......
......@@ -99,6 +99,17 @@ public:
return jet.two_norm2();
}
auto gradientHessian(const Jet& jet) const
{
std::tuple<FieldVector<double,dim>, ScaledIdentityMatrix<double,dim> > result;
std::get<0>(result) = derivative(*this)(jet);
std::get<1>(result) = derivative(derivative(*this))(jet);
return result;
}
friend auto derivative(const DirichletEnergyDensity& f) {
return DirichletEnergyDensityGradient(f);
}
......
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