Commit 4494e98d authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Move the split function in spectral coordinates into a separate file

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).
parent 9b8a756a
......@@ -19,6 +19,7 @@ install(FILES
integraldensity.hh
integralfunctional.hh
integralfunctionalconstrainedlinearization.hh
miehesplitting.hh
projectednewtonstep.hh
damagedelasticenergydensity_decomposed.hh
damagedelasticenergydensity_nondecomposed.hh
......
......@@ -3,12 +3,11 @@
#ifndef DUNE_FRACTUREPHASEFIELDS_DAMAGEDELASTICENERGYDENSITY_DECOMPOSED_HH
#define DUNE_FRACTUREPHASEFIELDS_DAMAGEDELASTICENERGYDENSITY_DECOMPOSED_HH
#include <numeric>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/fracture-phasefields/cracksurfacedensity.hh>
#include <dune/fracture-phasefields/miehesplitting.hh>
namespace Dune::FracturePhasefields {
namespace Decomposed {
......@@ -78,7 +77,7 @@ public:
Range damagedGradient = spectralToPhysical(eigenVectors,sigSpec);
// get tensile and compressive part of elastic energy density
const auto [psiP, psiM] = f_->psiSpec(state,eigenValues);
const auto [psiP, psiM] = f_->mieheSplitting_(eigenValues);
// psi_{d}
damagedGradient[size-1] = -2.0 * (1.0 - state[size-1])*psiP;
......@@ -117,7 +116,7 @@ public:
{
// get tensile and compressive part of elastic energy density
FieldVector<double, dim> eigen = eigenValues(jet);
const auto [psiP, psiM] = f_->psiSpec(jet,eigen);
const auto [psiP, psiM] = f_->mieheSplitting_(eigen);
gradient[size-1] = -2.0 * (1.0 - jet[size-1])*psiP;
......@@ -146,29 +145,6 @@ public:
typedef FieldVector<double, dim> DomainDim;
typedef FieldMatrix<double, dim,dim> MatrixDimDim;
std::tuple<DomainDim,DomainDim> sigSpec(const Domain& state, const DomainDim &eigenValues) const
{
// mechanical part of gradient psi,{eps} in spectral space
const double trEps = std::accumulate(eigenValues.begin(), eigenValues.end(), 0.0);
DomainDim sigP(0);
DomainDim sigM(0);
for (std::size_t i=0; i<dim; i++)
{
if (trEps>0)
sigP[i] += f_->lambda_*trEps;
else
sigM[i] += f_->lambda_*trEps;
if (eigenValues[i]>0)
sigP[i] += 2.0*f_->mu_*eigenValues[i];
else
sigM[i] += 2.0*f_->mu_*eigenValues[i];
}
return std::make_tuple(sigP,sigM);
}
DomainDim damagedSigSpec(const Domain& state, const DomainDim &eigenValues) const
{
// eigenvalues/eigenvectors of eps and tensile/compressive part of sigma
......@@ -179,7 +155,7 @@ public:
DomainDim damagedSigSpec(const Domain& state, const DomainDim &eigenValues, DomainDim &sigP, DomainDim &sigM) const
{
// eigenvalues/eigenvectors of eps and tensile/compressive part of sigma
auto [sigPlus,sigMinus] = sigSpec(state,eigenValues);
auto [sigPlus,sigMinus] = f_->mieheSplitting_.gradient(eigenValues);
sigP = sigPlus;
sigM = sigMinus;
......@@ -234,26 +210,8 @@ public:
// damaged stress in spectral space
auto damagedSigSpec = Gradient::damagedSigSpec(state,eigen, sigP, sigM);
// psi,{eps eps} in spectral space
FieldMatrix<double, dim,dim> hessianSpecP;
FieldMatrix<double, dim,dim> hessianSpecM;
const double trEps = std::accumulate(eigen.begin(), eigen.end(), 0.0);
if (trEps>0) {
hessianSpecP = f__->lambda_;
hessianSpecM = 0;
} else {
hessianSpecP = 0.0;
hessianSpecM = f__->lambda_;
}
for (std::size_t i=0;i < dim; ++i)
{
if (eigen[i]>0 )
hessianSpecP[i][i] += 2.0*f__->mu_;
else
hessianSpecM[i][i] += 2.0*f__->mu_;
}
// Positive and negative parts of the Hesse matrix in spectral space
const auto [hessianSpecP, hessianSpecM] = f__->mieheSplitting_.hessian(eigen);
// degrade only tensile part
FieldMatrix<double, dim,dim> hessianSpec = hessianSpecP * f__->degradation(state)
......@@ -269,7 +227,7 @@ public:
// psi,{d d}
// get tensile part of elastic energy density
const auto [psiP, psiM] = f__->psiSpec(state,eigen);
const auto [psiP, psiM] = f__->mieheSplitting_(eigen);
result[size-1][size-1] = 2.0*psiP;
......@@ -289,26 +247,8 @@ public:
// damaged stress in spectral space
auto damagedSigSpec = Gradient::damagedSigSpec(jet,eigen);
// psi,{eps eps} in spectral space
FieldMatrix<double, dim,dim> hessianSpecP;
FieldMatrix<double, dim,dim> hessianSpecM;
const double trEps = std::accumulate(eigen.begin(), eigen.end(), 0.0);
if (trEps>0) {
hessianSpecP = f__->lambda_;
hessianSpecM = 0;
} else {
hessianSpecP = 0.0;
hessianSpecM = f__->lambda_;
}
for (std::size_t i=0;i < dim; ++i)
{
if (eigen[i]>0 )
hessianSpecP[i][i] += 2.0*f__->mu_;
else
hessianSpecM[i][i] += 2.0*f__->mu_;
}
// Positive and negative parts of the Hesse matrix in spectral space
const auto [hessianSpecP, hessianSpecM] = f__->mieheSplitting_.hessian(eigen);
// degrade only tensile part
FieldMatrix<double, dim,dim> hessianSpec = hessianSpecP * f__->degradation(jet)
......@@ -322,7 +262,7 @@ public:
{
const auto [eigen, eigenV] = DamagedElasticEnergyDensity::eigenVectorDecomposition(jet);
auto [sigP,sigM] = Gradient::sigSpec(jet,eigen);
auto [sigP,sigM] = f__->mieheSplitting_.gradient(eigen);
auto undamagedGradientP = Gradient::spectralToPhysical(eigenV,sigP);
for (std::size_t i=0;i < size-1; ++i)
......@@ -333,7 +273,7 @@ public:
{
const auto [eigen, eigenV] = DamagedElasticEnergyDensity::eigenVectorDecomposition(jet);
auto [sigP,sigM] = Gradient::sigSpec(jet,eigen);
auto [sigP,sigM] = f__->mieheSplitting_.gradient(eigen);
auto undamagedGradientP = Gradient::spectralToPhysical(eigenV,sigP);
for (std::size_t i=0;i < size-1; ++i)
......@@ -345,7 +285,7 @@ public:
{
// get tensile part of elastic energy density
FieldVector<double, dim> eigen = eigenValues(jet);
const auto [psiP, psiM] = f__->psiSpec(jet, eigen);
const auto [psiP, psiM] = f__->mieheSplitting_(eigen);
result[size-1][size-1] = 2.0*psiP;
}
......@@ -438,25 +378,21 @@ public:
};
DamagedElasticEnergyDensity(const ParameterTree& parameters)
: crackSurfaceDensity_(parameters.get<CrackSurfaceDensity::Type>("CrackSurfaceDensityType"), parameters.get<double>("g_c"), parameters.get<double>("l"))
{
lambda_ = parameters.get<double>("lambda");
mu_ = parameters.get<double>("mu");
k_ = parameters.get<double>("k");
}
: crackSurfaceDensity_(parameters.get<CrackSurfaceDensity::Type>("CrackSurfaceDensityType"), parameters.get<double>("g_c"), parameters.get<double>("l")),
mieheSplitting_(parameters),
k_(parameters.get<double>("k"))
{}
RangeField operator () (const FieldVector<RangeField, size>& state) const
{
// get tensile and compressive part of elastic energy density
FieldVector<double, dim> eigenVal = eigenValues(state);
const auto [psiP, psiM] = psiSpec(state, eigenVal);
const auto [psiP, psiM] = mieheSplitting_(eigenValues(state));
// degrade only tensile part
// Degrade only the tensile part
RangeField energy = degradation(state)*psiP + (1.0 + k_)*psiM;
// linear part of critical energy density for model with elastic phase
if (crackSurfaceDensity_.type() == CrackSurfaceDensity::AT1) // psi_c = g_c / (4*c_w*l)
energy += crackSurfaceDensity_.g_c() / (4.0*crackSurfaceDensity_.c_w()*crackSurfaceDensity_.l()) * state[size-1];
energy += crackSurfaceDensity_.g_c() / (4.0*crackSurfaceDensity_.c_w()*crackSurfaceDensity_.l()) * state.back();
return energy;
}
......@@ -483,7 +419,7 @@ public:
gradient = DamagedElasticEnergyDensityGradient::spectralToPhysical(eigenVectors, gradientSpectral);
// Damage block
const auto [psiP, psiM] = psiSpec(state,eigenValues);
const auto [psiP, psiM] = mieheSplitting_(eigenValues);
gradient[size-1] = -2.0 * (1.0 - state[size-1])*psiP;
......@@ -494,26 +430,8 @@ public:
// 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 = std::accumulate(eigenValues.begin(), eigenValues.end(), 0.0);
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_;
}
// Positive and negative parts of the Hesse matrix in spectral space
const auto [hessianSpecP, hessianSpecM] = mieheSplitting_.hessian(eigenValues);
// degrade only tensile part
FieldMatrix<double, dim,dim> hessianSpec = hessianSpecP * degradation(state)
......@@ -562,26 +480,8 @@ public:
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 = std::accumulate(eigenValues.begin(), eigenValues.end(), 0.0);
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_;
}
// Positive and negative parts of the Hesse matrix in spectral space
const auto [hessianSpecP, hessianSpecM] = mieheSplitting_.hessian(eigenValues);
// degrade only tensile part
FieldMatrix<double, dim,dim> hessianSpec = hessianSpecP * degradation(state)
......@@ -596,7 +496,7 @@ public:
// Gradient part
// get tensile and compressive part of elastic energy density
FieldVector<double, dim> eigen = eigenValues(state);
const auto [psiP, psiM] = psiSpec(state,eigen);
const auto [psiP, psiM] = mieheSplitting_(eigen);
gradient[size-1] = -2.0 * (1.0 - state[size-1])*psiP;
......@@ -660,38 +560,14 @@ private:
return ((1.0 - state[size-1])*(1.0-state[size-1]) + k_);
}
std::array<RangeField,2> psiSpec(const FieldVector<RangeField, size>& stateVector,
const FieldVector<RangeField, dim>& eigenValues) const
{
const double trEps = std::accumulate(eigenValues.begin(), eigenValues.end(), 0.0);
RangeField psiP = 0.0;
RangeField psiM = 0.0;
if (trEps>0)
psiP = 0.5*lambda_*trEps*trEps;
else
psiM = 0.5*lambda_*trEps*trEps;
for (std::size_t i = 0; i<dim; i++)
{
if (eigenValues[i]>0)
psiP += mu_*eigenValues[i]*eigenValues[i];
else
psiM += mu_*eigenValues[i]*eigenValues[i];
}
return {psiP, psiM};
}
friend auto derivative(const DamagedElasticEnergyDensity& f) {
return DamagedElasticEnergyDensityGradient(f);
}
private:
CrackSurfaceDensity crackSurfaceDensity_;
RangeField lambda_;
RangeField mu_;
RangeField k_;
const CrackSurfaceDensity crackSurfaceDensity_;
const MieheSplitting<RangeField,dim> mieheSplitting_;
const RangeField k_;
};
/** \brief Compute the positive (i.e., tensile) part of the elastic energy,
......@@ -710,7 +586,7 @@ public:
double operator()(const FieldVector<double,jetSize>& jet)
{
const auto eigen = DamagedElasticEnergyDensity<dim>::eigenValues(jet);
const auto [psiPlus, psiMinus] = damagedElasticEnergyDensity_.psiSpec(jet,eigen);
const auto [psiPlus, psiMinus] = damagedElasticEnergyDensity_.mieheSplitting_(eigen);
return psiPlus;
};
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FRACTUREPHASEFIELDS_MIEHESPLITTING_HH
#define DUNE_FRACTUREPHASEFIELDS_MIEHESPLITTING_HH
#include <array>
#include <numeric>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/parametertree.hh>
namespace Dune::FracturePhasefields
{
/** \brief The St.Venant-Kirchhoff energy density with Miehe's split into tensile and compressive parts
*
* The formula is described in
*
* C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture:
* Variational principles and multi-field FE implementations. Int. J. Numer. Meth. Engng, 83:1273–1311, 2010.
* DOI: https://doi.org/10.1002/nme.2861
*
* This class implements the energy in terms of the strain eigenvalues.
*
* \tparam field_type Number type used for strain components and energy density values
* \tparam dim The space dimension
*/
template <class field_type, int dim>
class MieheSplitting
{
public:
/** \brief Constructor from a set of material parameters
*/
MieheSplitting(const ParameterTree& parameters)
: lambda_(parameters.get<double>("lambda")),
mu_(parameters.get<double>("mu"))
{}
/** \brief Tensile and compressive parts of the elastic energy density
*/
std::array<field_type,2> operator() (const FieldVector<field_type, dim>& eigenValues) const
{
// Tensile and compressive parts of the elastic energy density
std::array<field_type,2> result = {0, 0};
field_type& psiP = result[0];
field_type& psiM = result[1];
// The trace of the strain matrix
const double trEps = std::accumulate(eigenValues.begin(), eigenValues.end(), 0.0);
if (trEps>0)
psiP = 0.5*lambda_*trEps*trEps;
else
psiM = 0.5*lambda_*trEps*trEps;
for (std::size_t i = 0; i<dim; i++)
{
if (eigenValues[i]>0)
psiP += mu_*eigenValues[i]*eigenValues[i];
else
psiM += mu_*eigenValues[i]*eigenValues[i];
}
return result;
}
/** \brief Gradients of the tensile and compressive parts of the energy density
*/
std::array<FieldVector<field_type, dim>,2> gradient(const FieldVector<field_type, dim>& eigenValues) const
{
std::array<FieldVector<field_type, dim>,2> result = {FieldVector<field_type, dim>(0.0), FieldVector<field_type, dim>(0.0)};
auto& gradientPlus = result[0];
auto& gradientMinus = result[1];
const double trEps = std::accumulate(eigenValues.begin(), eigenValues.end(), 0.0);
for (std::size_t i=0; i<dim; i++)
{
if (trEps>0)
gradientPlus[i] += lambda_ * trEps;
else
gradientMinus[i] += lambda_ * trEps;
if (eigenValues[i]>0)
gradientPlus[i] += 2.0 * mu_ * eigenValues[i];
else
gradientMinus[i] += 2.0 * mu_ * eigenValues[i];
}
return result;
}
/** \brief Hesse matrix of the tensile and compressive parts of the energy density
*/
std::array<FieldMatrix<field_type, dim, dim>,2> hessian(const FieldVector<field_type, dim>& eigenValues) const
{
std::array<FieldMatrix<field_type, dim>,2> result = {FieldMatrix<field_type,dim,dim>(0.0), FieldMatrix<field_type,dim,dim>(0.0)};
auto& hessianPlus = result[0];
auto& hessianMinus = result[1];
const double trEps = std::accumulate(eigenValues.begin(), eigenValues.end(), 0.0);
if (trEps>0) {
hessianPlus = lambda_;
hessianMinus = 0;
} else {
hessianPlus = 0.0;
hessianMinus = lambda_;
}
for (std::size_t i=0;i < dim; ++i)
{
if (eigenValues[i]>0 )
hessianPlus[i][i] += 2.0 * mu_;
else
hessianMinus[i][i] += 2.0 * mu_;
}
return result;
}
private:
// The Lamé parameters of the material
const field_type lambda_;
const field_type mu_;
};
} // namespace Dune::FracturePhasefields
#endif
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