Commit e06e72c6 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Local solver: Compute density gradient and Hessian in one go

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.
parent 389a7f58
......@@ -55,6 +55,8 @@ public:
*/
class DamagedElasticEnergyDensityGradient
{
friend class DamagedElasticEnergyDensity<dim>;
public:
using Range = Domain;
......@@ -209,6 +211,8 @@ public:
class DamagedElasticEnergyDensityHessian : public DamagedElasticEnergyDensityGradient
{
friend class DamagedElasticEnergyDensity<dim>;
public:
using Range = FieldMatrix<RangeField,size,size>;
using Gradient = DamagedElasticEnergyDensityGradient;
......@@ -455,6 +459,81 @@ public:
return energy;
}
auto gradientHessian(const FieldVector<RangeField, size>& state, size_t row, size_t col) const
{
if (row!=col) {
std::cerr << "row != col not implemented!" << std::endl;
std::terminate();
}
// Gradient part
std::tuple<typename DamagedElasticEnergyDensityGradient::Range, typename DamagedElasticEnergyDensityHessian::Range> result;
auto& gradient = std::get<0>(result);
auto& hessian = std::get<1>(result);
// Displacement block
if (row==0)
{
const auto [eigenValues, eigenVectors] = eigenVectorDecomposition(state);
// Compute gradient
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);
// Compute Hesse matrix
// First, 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);
}
// Damage block
if (row==1)
{
// Gradient part
// get tensile and compressive part of elastic energy density
FieldVector<double, dim> eigen = eigenValues(state);
const auto [psiP, psiM] = psiSpec(state,eigen);
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());
// Hesse matrix part
hessian[size-1][size-1] = 2.0*psiP;
}
return result;
}
private:
/** \brief Compute eigenvalues of the strain
......
......@@ -338,6 +338,80 @@ public:
return energy;
}
auto gradientHessian(const FieldVector<RangeField, size>& state, size_t row, size_t col) const
{
if (row!=col) {
std::cerr << "row != col not implemented!" << std::endl;
std::terminate();
}
// Gradient part
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);
// Displacement part
if (row==0)
{
// Gradient part
for (auto&& diagonal : diagonals())
gradient[diagonal] = degradationValue * lambda_ * trEps;
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);
}
for (std::size_t i=0; i<size-1; i++)
gradient[i] += degradationValue * 2.0 * mu_ * state[i] * II[i];
// Hesse matrix part
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;
}
// Damage part
if (row==1)
{
// Gradient part
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());
// Hesse matrix part
hessian[size-1][size-1] = 2.0 * psi_0(state);
}
return result;
}
private:
RangeField traceEps(const FieldVector<RangeField, size>& stateVector) const
{
......
......@@ -204,11 +204,17 @@ public:
return (*pimpl_)(jet);
}
auto gradientHessian(const Jet& jet, size_t row, size_t col) const
{
return pimpl_->gradientHessian(jet, row, col);
}
private:
struct Concept
{
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, size_t row, size_t col) const = 0;
virtual IntegralDensityGradient<Jet> derivative_() const = 0;
virtual Concept* clone() const = 0;
};
......@@ -218,6 +224,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, 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 ); }
T object;
......
......@@ -226,12 +226,9 @@ public:
auto jet_j = f_->jet()[j];
BT_ij.umtv(v, jet_j);
// This is the gradient part
auto g_j = densityGradient(jet_j, row);
// Get gradient and Hesse matrix of the density
auto [g_j, Hj] = density().gradientHessian(jet_j, row, col);
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
......
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