Commit 8e497a1a authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Move transform from spectral to physical coordinates into separate file

parent 4494e98d
......@@ -21,6 +21,7 @@ install(FILES
integralfunctionalconstrainedlinearization.hh
miehesplitting.hh
projectednewtonstep.hh
spectraltophysicaltransform.hh
damagedelasticenergydensity_decomposed.hh
damagedelasticenergydensity_nondecomposed.hh
ductilefracturejet.hh
......
......@@ -8,6 +8,7 @@
#include <dune/fracture-phasefields/cracksurfacedensity.hh>
#include <dune/fracture-phasefields/miehesplitting.hh>
#include <dune/fracture-phasefields/spectraltophysicaltransform.hh>
namespace Dune::FracturePhasefields {
namespace Decomposed {
......@@ -39,11 +40,6 @@ class DamagedElasticEnergyDensity
};
}
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:
typedef double RangeField;
......@@ -74,7 +70,8 @@ public:
auto sigSpec = damagedSigSpec(state,eigenValues,sigP,sigM);
// transform sigma in physical space
Range damagedGradient = spectralToPhysical(eigenVectors,sigSpec);
Range damagedGradient;
SpectralToPhysicalTransform<RangeField,dim>::gradient(eigenVectors,sigSpec, damagedGradient, 0);
// get tensile and compressive part of elastic energy density
const auto [psiP, psiM] = f_->mieheSplitting_(eigenValues);
......@@ -108,7 +105,7 @@ public:
auto sigSpec = damagedSigSpec(jet,eigenValues,sigP,sigM);
// transform sigma in physical space
return spectralToPhysical(eigenVectors,sigSpec);
SpectralToPhysicalTransform<RangeField,dim>::gradient(eigenVectors,sigSpec,gradient,0);
}
// psi_{d}
......@@ -165,20 +162,6 @@ public:
return (sigP*degradationValue) + (sigM*(1.0 + f_->k_));
}
/** \brief Transform gradient in spectral space to gradient in physical space
*/
static Range spectralToPhysical(const MatrixDimDim &eigenVectors, const DomainDim &gradientSpec)
{
// mechanical part of gradient psi,{eps} in physical space
Range gradient(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++)
gradient[linearIndex(j,i)] += ((i==j) ? 1 : 2) * gradientSpec[a]*eigenVectors[a][i]*eigenVectors[a][j];
return gradient;
}
friend auto derivative(const DamagedElasticEnergyDensityGradient& gradient) {
return DamagedElasticEnergyDensityHessian(*(gradient.f_));
}
......@@ -218,10 +201,14 @@ public:
+ hessianSpecM * (1.0 + f__->k_);
// psi,{eps eps} in physical space
Range result = spectralToPhysical(eigen, eigenV, damagedSigSpec, hessianSpec);
Range result;
SpectralToPhysicalTransform<RangeField,dim>::hessian(eigen, eigenV,
damagedSigSpec, hessianSpec,
result, 0);
// psi,{eps d}
auto undamagedGradientP = Gradient::spectralToPhysical(eigenV,sigP);
FieldVector<RangeField,size> undamagedGradientP;
SpectralToPhysicalTransform<RangeField,dim>::gradient(eigenV,sigP,undamagedGradientP,0);
for (std::size_t i=0;i < size-1; ++i)
result[size-1][i] = result[i][size-1] = -2.0*(1 - state[size-1])*undamagedGradientP[i];
......@@ -254,7 +241,9 @@ public:
FieldMatrix<double, dim,dim> hessianSpec = hessianSpecP * f__->degradation(jet)
+ hessianSpecM * (1.0 + f__->k_);
return spectralToPhysical(eigen, eigenV, damagedSigSpec, hessianSpec);
SpectralToPhysicalTransform<RangeField,dim>::hessian(eigen, eigenV,
damagedSigSpec, hessianSpec,
result, 0);
}
// psi,{eps d}
......@@ -264,7 +253,8 @@ public:
auto [sigP,sigM] = f__->mieheSplitting_.gradient(eigen);
auto undamagedGradientP = Gradient::spectralToPhysical(eigenV,sigP);
FieldVector<RangeField,size> undamagedGradientP;
SpectralToPhysicalTransform<RangeField,dim>::gradient(eigenV,sigP,undamagedGradientP,0);
for (std::size_t i=0;i < size-1; ++i)
result[i][size-1] = -2.0*(1 - jet[size-1])*undamagedGradientP[i];
}
......@@ -275,7 +265,8 @@ public:
auto [sigP,sigM] = f__->mieheSplitting_.gradient(eigen);
auto undamagedGradientP = Gradient::spectralToPhysical(eigenV,sigP);
FieldVector<RangeField,size> undamagedGradientP;
SpectralToPhysicalTransform<RangeField,dim>::gradient(eigenV,sigP,undamagedGradientP,0);
for (std::size_t i=0;i < size-1; ++i)
result[size-1][i] = -2.0*(1 - jet[size-1])*undamagedGradientP[i];
}
......@@ -307,73 +298,6 @@ public:
}
private:
/** \brief Transform Hesse matrix in spectral space to Hesse matrix in physical space
*
* 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.
*/
static Range spectralToPhysical(const FieldVector<RangeField,dim>& eigenValues,
const FieldMatrix<RangeField,dim,dim>& eigenVectors,
const FieldVector<RangeField,dim>& gradientSpectral,
const FieldMatrix<RangeField,dim,dim>& hessianSpectral)
{
Range result(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;
typename Range::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
result[linearIndex(i,j)][linearIndex(k,l)] += value;
if (linearIndex(i,j) !=linearIndex(k,l)) // Off-diagonal entry?
result[linearIndex(k,l)][linearIndex(i,j)] += value;
}
return result;
}
const DamagedElasticEnergyDensity* f__;
};
......@@ -416,7 +340,7 @@ public:
const auto gradientSpectral = derivative(*this).damagedSigSpec(state,eigenValues,sigP,sigM);
// transform sigma in physical space
gradient = DamagedElasticEnergyDensityGradient::spectralToPhysical(eigenVectors, gradientSpectral);
SpectralToPhysicalTransform<double,dim>::gradient(eigenVectors, gradientSpectral, gradient, 0);
// Damage block
const auto [psiP, psiM] = mieheSplitting_(eigenValues);
......@@ -437,10 +361,13 @@ public:
FieldMatrix<double, dim,dim> hessianSpec = hessianSpecP * degradation(state)
+ hessianSpecM * (1.0 + k_);
hessian = DamagedElasticEnergyDensityHessian::spectralToPhysical(eigenValues, eigenVectors, gradientSpectral, hessianSpec);
SpectralToPhysicalTransform<RangeField,dim>::hessian(eigenValues, eigenVectors,
gradientSpectral, hessianSpec,
hessian, 0);
// Displacement-damage and damage-displacement blocks
auto undamagedGradientP = DamagedElasticEnergyDensityGradient::spectralToPhysical(eigenVectors,sigP);
typename DamagedElasticEnergyDensityGradient::Range undamagedGradientP;
SpectralToPhysicalTransform<double,dim>::gradient(eigenVectors, sigP, undamagedGradientP, 0);
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];
......@@ -477,7 +404,7 @@ public:
const auto gradientSpectral = derivative(*this).damagedSigSpec(state,eigenValues,sigP,sigM);
// transform sigma in physical space
gradient = DamagedElasticEnergyDensityGradient::spectralToPhysical(eigenVectors, gradientSpectral);
SpectralToPhysicalTransform<RangeField,dim>::gradient(eigenVectors, gradientSpectral, gradient, 0);
// Compute Hesse matrix
// Positive and negative parts of the Hesse matrix in spectral space
......@@ -487,7 +414,9 @@ public:
FieldMatrix<double, dim,dim> hessianSpec = hessianSpecP * degradation(state)
+ hessianSpecM * (1.0 + k_);
hessian = DamagedElasticEnergyDensityHessian::spectralToPhysical(eigenValues, eigenVectors, gradientSpectral, hessianSpec);
SpectralToPhysicalTransform<RangeField,dim>::hessian(eigenValues, eigenVectors,
gradientSpectral, hessianSpec,
hessian, 0);
}
// Damage block
......
// -*- 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
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