From 47ce740a0138701d92e87023a0055e8ad9d7443e Mon Sep 17 00:00:00 2001
From: Andreas Fischle <andreas.fischle@uni-due.de>
Date: Sun, 5 Apr 2015 16:30:28 +0200
Subject: [PATCH] Rewrite Hencky energy following Patrizio Neff's notation
---
dune/gfe/henckyenergy.hh | 47 +++++++++++++++++++++++-----------------
1 file changed, 27 insertions(+), 20 deletions(-)
diff --git a/dune/gfe/henckyenergy.hh b/dune/gfe/henckyenergy.hh
index 1b1ed5c7..51a6ca49 100644
--- a/dune/gfe/henckyenergy.hh
+++ b/dune/gfe/henckyenergy.hh
@@ -32,6 +32,7 @@ public: // for testing
// Lame constants
mu_ = parameters.template get<double>("mu");
lambda_ = parameters.template get<double>("lambda");
+ kappa_ = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
}
/** \brief Assemble the energy for a single element */
@@ -40,7 +41,7 @@ public: // for testing
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
/** \brief Lame constants */
- double mu_, lambda_;
+ double mu_, lambda_, kappa_;
};
template <class GridView, class LocalFiniteElement, class field_type>
@@ -65,6 +66,7 @@ energy(const Entity& element,
const Dune::QuadratureRule<DT, gridDim>& quad
= Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
+ field_type distortionEnergy = 0, dilatationEnergy = 0;
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
@@ -89,40 +91,45 @@ energy(const Entity& element,
derivative[j].axpy(localConfiguration[i][j], gradients[i]);
/////////////////////////////////////////////////////////
- // compute F^T F
+ // compute C = F^TF
/////////////////////////////////////////////////////////
- Dune::FieldMatrix<field_type,gridDim,gridDim> FTF(0);
+ Dune::FieldMatrix<field_type,gridDim,gridDim> C(0);
for (int i=0; i<gridDim; i++)
for (int j=0; j<gridDim; j++)
for (int k=0; k<gridDim; k++)
- FTF[i][j] += derivative[k][i] * derivative[k][j];
+ C[i][j] += derivative[k][i] * derivative[k][j];
//////////////////////////////////////////////////////////
- // Eigenvalues of FTF
+ // Eigenvalues of C = F^TF
//////////////////////////////////////////////////////////
- Dune::FieldVector<field_type,dim> lambda;
- FMatrixHelp::eigenValues(FTF, lambda);
+ Dune::FieldVector<field_type,dim> sigmaSquared;
+ FMatrixHelp::eigenValues(C, sigmaSquared);
- // logarithms of the eigenvalues
- std::array<field_type,dim> ln;
- for (int i=0; i<dim; i++)
- ln[i] = std::log(lambda[i]);
+ // logarithms of the singular values of F, i.e., eigenvalues of U
+ std::array<field_type, dim> lnSigma;
+ for (int i = 0; i < dim; i++)
+ lnSigma[i] = 0.5 * log(sigmaSquared[i]);
- // Add the local energy density
- for (int i=0; i<dim; i++)
- energy += weight * mu_ * ln[i]*ln[i];
-
- field_type trace = 0;
- for (int i=0; i<dim; i++)
- trace += ln[i];
+ field_type trLogU = 0;
+ for (int i = 0; i < dim; i++)
+ trLogU += lnSigma[i];
- energy += weight * 0.5 * lambda_ * trace * trace;
+ field_type normDevLogUSquared = 0;
+ for (int i = 0; i < dim; i++)
+ {
+ // the deviator shifts the spectrum by the trace
+ field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
+ normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
+ }
+ // Add the local energy density
+ distortionEnergy += weight * mu_ * normDevLogUSquared;
+ dilatationEnergy += weight * 0.5 * kappa_ * (trLogU * trLogU);
}
- return energy;
+ return distortionEnergy + dilatationEnergy;
}
} // namespace Dune
--
GitLab