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

Merge branch 'modularize-spectrally-degraded-energy' into 'master'

Modularize the Decomposed::DamagedElasticEnergyDensity class

See merge request !50
parents 9b8a756a 4dadee14
Pipeline #8770 passed with stage
in 4 minutes and 37 seconds
......@@ -19,7 +19,9 @@ install(FILES
integraldensity.hh
integralfunctional.hh
integralfunctionalconstrainedlinearization.hh
miehesplitting.hh
projectednewtonstep.hh
spectraltophysicaltransform.hh
damagedelasticenergydensity_decomposed.hh
damagedelasticenergydensity_nondecomposed.hh
ductilefracturejet.hh
......
......@@ -6,8 +6,6 @@
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/fracture-phasefields/cracksurfacedensity.hh>
namespace Dune {
namespace FracturePhasefields {
namespace NonDecomposed {
......@@ -238,10 +236,6 @@ public:
gradient[size-1] = -2.0 * (1.0 - state[size-1])*f_->psi_0(state);
// linear part of critical energy density for model with elastic phase
if (f_->crackSurfaceDensity_.type() == CrackSurfaceDensity::AT1) // psi_c = g_c / (4*c_w*l)
gradient[size-1] += f_->crackSurfaceDensity_.g_c() / (4.0*f_->crackSurfaceDensity_.c_w()*f_->crackSurfaceDensity_.l());
return gradient;
}
......@@ -286,14 +280,8 @@ public:
// Damage part
if (block==1)
{
gradient[size-1] = -2.0 * (1.0 - jet[size-1])*f_->psi_0(jet);
// linear part of critical energy density for model with elastic phase
if (f_->crackSurfaceDensity_.type() == CrackSurfaceDensity::AT1) // psi_c = g_c / (4*c_w*l)
gradient[size-1] += f_->crackSurfaceDensity_.g_c() / (4.0*f_->crackSurfaceDensity_.c_w()*f_->crackSurfaceDensity_.l());
}
return gradient;
}
......@@ -319,7 +307,6 @@ 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");
......@@ -330,10 +317,6 @@ public:
{
RangeField energy = degradation(state)*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)
energy += crackSurfaceDensity_.g_c() / (4.0*crackSurfaceDensity_.c_w()*crackSurfaceDensity_.l()) * state[size-1];
return energy;
}
......@@ -380,10 +363,6 @@ public:
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())
......@@ -479,10 +458,6 @@ public:
// 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);
}
......@@ -530,7 +505,6 @@ private:
RangeField lambda_;
RangeField mu_;
RangeField k_;
CrackSurfaceDensity crackSurfaceDensity_;
/** \brief The diagonal entries of a symmetric matrix */
static std::array<unsigned int,dim> diagonals()
......
// -*- 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
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FRACTUREPHASEFIELDS_SPECTRALTOPHYSICALTRANSFORM_HH
#define DUNE_FRACTUREPHASEFIELDS_SPECTRALTOPHYSICALTRANSFORM_HH
namespace Dune::FracturePhasefields
{
/** \brief Transform derivates of a spectral function from spectral to physical coordinates
*
* This class contains only static methods.
*/
template <class field_type, int dim>
class SpectralToPhysicalTransform
{
/** \brief Size of a symmetric dim x dim matrix as a vector
*/
static const std::size_t size = dim*(dim+1)/2;
/** \brief Map the (i,j)-index of a (symmetric) matrix to the underlying index in the data vector
*
* The memory layout is v00, v10, v11, v20, v21, v22, v30, ...
*/
static std::size_t linearIndex(std::size_t i, std::size_t j)
{
return (i >= j) ? i*(i+1)/2+j : j*(j+1)/2+i;
}
public:
/** \brief Transform gradient in spectral space to gradient in physical space
*
* The result is written into the gradientPhysical argument, starting at offset.
*
* \tparam Gradient Type of the result vector. Must be a dense vector.
*
* The formula is described in Lemma 3.1 of
*
* Adrian S. Lewis and Hristo S. Sendov: "Twice Differentiable Spectral Functions"
* SIAM J. Matrix Anal. Appl., Vol. 23, No. 2, pp. 368–386, (2001)
*/
template <class Gradient>
static void gradient(const FieldMatrix<field_type,dim,dim>& eigenVectors,
const FieldVector<field_type,dim>& gradientSpectral,
Gradient& gradientPhysical,
std::size_t offset)
{
std::fill(gradientPhysical.begin()+offset, gradientPhysical.begin()+offset+size, 0.0);
for (std::size_t i=0; i<dim; i++)
for (std::size_t j=0; j<=i; j++)
for (std::size_t a=0; a<dim; a++)
gradientPhysical[linearIndex(j,i)+offset] += ((i==j) ? 1 : 2) * eigenVectors[a][i]*gradientSpectral[a]*eigenVectors[a][j];
}
/** \brief Transform Hesse matrix in spectral space to Hesse matrix in physical space
*
* The result is written into the hessianPhysical argument, starting at (offset,offset).
*
* \tparam Hessian Type of the result matrix. Must be a dense matrix.
*
* The formula is derived and described in
*
* Adrian S. Lewis and Hristo S. Sendov: "Twice Differentiable Spectral Functions"
* SIAM J. Matrix Anal. Appl., Vol. 23, No. 2, pp. 368–386, (2001)
*
* In there it is Theorem 3.3.
*/
template <class Hessian>
static void hessian(const FieldVector<field_type,dim>& eigenValues,
const FieldMatrix<field_type,dim,dim>& eigenVectors,
const FieldVector<field_type,dim>& gradientSpectral,
const FieldMatrix<field_type,dim,dim>& hessianSpectral,
Hessian& hessianPhysical,
std::size_t offset)
{
for (std::size_t i=offset; i<offset+size; ++i)
for (std::size_t j=offset; j<offset+size; ++j)
hessianPhysical[i][j] = 0.0;
// psi,{eps eps} in physical space
for (std::size_t i=0;i < dim; ++i)
for (std::size_t j=0;j <= i; ++j)
for (std::size_t k=0;k < dim; ++k)
for (std::size_t l=0;l <= k; ++l)
{
// The result matrix is symmetric: Compute only the lower triangular part
if (linearIndex(i,j) < linearIndex(k,l))
continue;
field_type value(0.0);
for (std::size_t a=0;a < dim; ++a)
for (std::size_t b=0;b < dim; ++b)
{
// In the Lewis/Sendov paper, this is the expression
// $\nabla^2 f(\lambda(A)) [diag \tilde{H}, diag \tilde{H}]$.
value += eigenVectors[a][i] * eigenVectors[a][j] * hessianSpectral[a][b] * eigenVectors[b][k] * eigenVectors[b][l];
// In the Lewis/Sendov paper, this is the expression
// $\langle \mathcal{A}, \tilde{H} \circ \tilde{H} \rangle$.
if (a!=b)
{
const auto eigenBase = (eigenVectors[a][i]*eigenVectors[b][j]*eigenVectors[a][k]*eigenVectors[b][l])
+ (eigenVectors[a][i]*eigenVectors[b][j]*eigenVectors[b][k]*eigenVectors[a][l]);
// for different eigenvalues
if(fabs(eigenValues[b]-eigenValues[a]) > 1e-8)
value += 0.5*(gradientSpectral[b]-gradientSpectral[a])/(eigenValues[b]-eigenValues[a]) * eigenBase;
// for equal eigenvalues
else
value += 0.5*(hessianSpectral[b][b]-hessianSpectral[a][b]) * eigenBase;
}
}
// Off-diagonal entries count twice: They come in identical pairs,
// but we compute only one of them.
auto ijFactor = (i==j) ? 1 : 2;
auto klFactor = (k==l) ? 1 : 2;
value *= ijFactor * klFactor;
// Add new value
hessianPhysical[linearIndex(i,j)][linearIndex(k,l)] += value;
if (linearIndex(i,j) !=linearIndex(k,l)) // Off-diagonal entry?
hessianPhysical[linearIndex(k,l)][linearIndex(i,j)] += value;
}
}
};
} // namespace Dune::FracturePhasefields
#endif
......@@ -517,6 +517,7 @@ int main (int argc, char *argv[]) try
FracturePhasefields::CrackSurfaceAssembler<decltype(blockedBasis)> localCrackAssembler(crackSurfaceDensity);
// Assemble linear term for rhs
// Currently this only contains the linear term of the AT-1 crack surface density functional
auto functionalAssembler = Dune::Fufem::DuneFunctionsFunctionalAssembler{damageBasis};
functionalAssembler.assembleBulk(Functions::istlVectorBackend(crackRHS), localCrackAssembler);
......@@ -639,10 +640,6 @@ int main (int argc, char *argv[]) try
// Create a solver
//////////////////////////////
Vector rhs;
Functions::istlVectorBackend(rhs).resize(blockedBasis);
Functions::istlVectorBackend(rhs) = 0.0;
using BitVector = Solvers::DefaultBitVector_t<Vector>;
// first create a linear iteration step
......@@ -783,7 +780,7 @@ int main (int argc, char *argv[]) try
// Assemble load vector
/////////////////////////////
rhs = 0;
Vector rhs = crackRHS;
xOld = x;
......
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