Commit a7e5a266 authored by Patrick Jaap's avatar Patrick Jaap
Browse files

Replace local energies by local densities

 - wrap the densities by a generic Integral energy
 - this preserves all material properties
 - reduces redundant code
parent 93edeb51
# Master (will become release 2.8)
- Introduce class `LocalDensity` for material-specific implementations
- Introduce class `LocalIntegralEnergy` to work with the densities
## Deprecations and removals
- Redundant implementations of `LocalEnergy` classes are now replaced by combining `LocalDensity` and `LocalIntegralEnergy`
install(FILES
adolcmaterial.hh
adolcneohookeanmaterial.hh
exphenckydensity.hh
exphenckyenergy.hh
geomexactstvenantkirchhoffmaterial.hh
henckydensity.hh
henckyenergy.hh
localdensity.hh
localintegralenergy.hh
mooneyrivlindensity.hh
mooneyrivlinenergy.hh
neohookeanmaterial.hh
neohookedensity.hh
neohookeenergy.hh
neumannenergy.hh
ogdenmaterial.hh
stvenantkirchhoffdensity.hh
stvenantkirchhoffenergy.hh
sumenergy.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/elasticity/materials)
#ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYDENSITY_HH
#define DUNE_ELASTICITY_MATERIALS_EXPHENCKYDENSITY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/elasticity/materials/localdensity.hh>
namespace Dune::Elasticity {
template<int dim, class field_type = double>
class ExpHenckyDensity final
: public LocalDensity<dim,field_type>
{
public:
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
ExpHenckyDensity(const Dune::ParameterTree& parameters)
{
// Lame constants
mu_ = parameters.get<double>("mu");
lambda_ = parameters.get<double>("lambda");
kappa_ = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
// Exponential Hencky constants with naming convention from Neff, Ghiba and Lankeit,
// "The exponentiated Hencky-logarithmic strain energy. Part I: Constitutive issues
// and rank one convexity"
k_ = parameters.get<double>("k");
khat_ = parameters.get<double>("khat");
// weights for distortion and dilatation parts of the energy
alpha_ = mu_ / k_;
beta_ = kappa_ / (2.0 * khat_);
}
/** \brief Evaluation with the deformation gradient
*
* \param x Position of the gradient
* \param gradient Deformation gradient
*/
field_type operator() (const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
{
/////////////////////////////////////////////////////////
// compute C = F^TF
/////////////////////////////////////////////////////////
Dune::FieldMatrix<field_type,dim,dim> C(0);
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
for (int k=0; k<dim; k++)
C[i][j] += gradient[k][i] * gradient[k][j];
//////////////////////////////////////////////////////////
// Eigenvalues of FTF
//////////////////////////////////////////////////////////
// compute eigenvalues of C, i.e., squared singular values \sigma_i of F
Dune::FieldVector<field_type, dim> sigmaSquared;
FMatrixHelp::eigenValues(C, sigmaSquared);
// 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]);
field_type trLogU = 0;
for (int i = 0; i < dim; i++)
trLogU += lnSigma[i];
field_type normDevLogUSquared = 0;
for (int i = 0; i < dim; i++)
{
// the deviator shifts the
field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
}
// Add the local energy density
auto distortionEnergy = alpha_ * exp(khat_ * normDevLogUSquared);
auto dilatationEnergy = beta_ * exp(k_ * (trLogU * trLogU));
return distortionEnergy + dilatationEnergy;
}
/** \brief Lame constants */
field_type mu_, lambda_, kappa_, k_, khat_, alpha_, beta_;
};
} // namespace Dune::Elasticity
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGYDENSITY_HH
#ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localenergy.hh>
#include <dune/elasticity/materials/exphenckydensity.hh>
#include <dune/elasticity/materials/localintegralenergy.hh>
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
class ExpHenckyEnergy
: public Elasticity::LocalEnergy<GridView,
LocalFiniteElement,
std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::ExpHenckyDensity")
ExpHenckyEnergy
: public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
{
// grid types
typedef typename GridView::Grid::ctype DT;
typedef typename GridView::template Codim<0>::Entity Entity;
// some other sizes
enum {gridDim=GridView::dimension};
enum {dim=GridView::dimension};
using Base = Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>;
public:
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
// backwards compatibility: forward directly to integral energy
ExpHenckyEnergy(const Dune::ParameterTree& parameters)
{
// Lame constants
mu_ = parameters.get<double>("mu");
lambda_ = parameters.get<double>("lambda");
kappa_ = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
// Exponential Hencky constants with naming convention from Neff, Ghiba and Lankeit,
// "The exponentiated Hencky-logarithmic strain energy. Part I: Constitutive issues
// and rank one convexity"
k_ = parameters.get<double>("k");
khat_ = parameters.get<double>("khat");
// weights for distortion and dilatation parts of the energy
alpha_ = mu_ / k_;
beta_ = kappa_ / (2.0 * khat_);
}
/** \brief Assemble the energy for a single element */
field_type energy (const Entity& e,
const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
/** \brief Lame constants */
field_type mu_, lambda_, kappa_, k_, khat_, alpha_, beta_;
: Base(std::make_shared<Elasticity::ExpHenckyDensity<GridView::dimension,field_type>>(parameters))
{}
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
ExpHenckyEnergy<GridView, LocalFiniteElement, field_type>::
energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
{
assert(element.type() == localFiniteElement.type());
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
// store gradients of shape functions and base functions
std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim;
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
const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
const DT integrationElement = element.geometry().integrationElement(quadPos);
const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
DT weight = quad[pt].weight() * integrationElement;
// get gradients of shape functions
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
// compute gradients of base functions
for (size_t i=0; i<gradients.size(); ++i)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
for (size_t i=0; i<gradients.size(); i++)
for (int j=0; j<gridDim; j++)
derivative[j].axpy(localConfiguration[i][j], gradients[i]);
/////////////////////////////////////////////////////////
// compute C = F^T F
/////////////////////////////////////////////////////////
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++)
C[i][j] += derivative[k][i] * derivative[k][j];
//////////////////////////////////////////////////////////
// Eigenvalues of FTF
//////////////////////////////////////////////////////////
// compute eigenvalues of C, i.e., squared singular values \sigma_i of F
Dune::FieldVector<field_type, dim> sigmaSquared;
FMatrixHelp::eigenValues(C, sigmaSquared);
// 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]);
field_type trLogU = 0;
for (int i = 0; i < dim; i++)
trLogU += lnSigma[i];
field_type normDevLogUSquared = 0;
for (int i = 0; i < dim; i++)
{
// the deviator shifts the
field_type lnDevSigma_i = lnSigma[i] - (1.0 / double(dim)) * trLogU;
normDevLogUSquared += lnDevSigma_i * lnDevSigma_i;
}
// Add the local energy density
distortionEnergy += weight * alpha_ * exp(khat_ * normDevLogUSquared);
dilatationEnergy += weight * beta_ * exp(k_ * (trLogU * trLogU));
}
return distortionEnergy + dilatationEnergy;
}
} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_EXPHENCKYENERGY_HH
......
#ifndef DUNE_ELASTICITY_MATERIALS_HENCKYDENSITY_HH
#define DUNE_ELASTICITY_MATERIALS_HENCKYDENSITY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/elasticity/materials/localdensity.hh>
namespace Dune::Elasticity {
template<int dim, class field_type = double>
class HenckyDensity final
: public LocalDensity<dim,field_type>
{
public:
/** \brief Constructor with a set of material parameters
*
* \param parameters The material parameters
*/
HenckyDensity(const Dune::ParameterTree& parameters)
{
// Lame constants
mu_ = parameters.template get<double>("mu");
lambda_ = parameters.template get<double>("lambda");
kappa_ = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
}
/** \brief Evaluation with the deformation gradient
*
* \param x Position of the gradient
* \param gradient Deformation gradient
*/
field_type operator() (const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const
{
/////////////////////////////////////////////////////////
// compute C = F^TF
/////////////////////////////////////////////////////////
Dune::FieldMatrix<field_type,dim,dim> C(0);
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
for (int k=0; k<dim; k++)
C[i][j] += gradient[k][i] * gradient[k][j];
//////////////////////////////////////////////////////////
// Eigenvalues of C = F^TF
//////////////////////////////////////////////////////////
Dune::FieldVector<field_type,dim> sigmaSquared;
FMatrixHelp::eigenValues(C, sigmaSquared);
// 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]);
field_type trLogU = 0;
for (int i = 0; i < dim; i++)
trLogU += lnSigma[i];
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
auto distortionEnergy = mu_ * normDevLogUSquared;
auto dilatationEnergy = 0.5 * kappa_ * (trLogU * trLogU);
return dilatationEnergy + distortionEnergy;
}
/** \brief Lame constants */
double mu_, lambda_, kappa_;
};
} // namespace Dune::Elasticity
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_HENCKYDENSITY_HH
#ifndef DUNE_ELASTICITY_MATERIALS_HENCKYENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_HENCKYENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localenergy.hh>
#include <dune/elasticity/materials/henckydensity.hh>
#include <dune/elasticity/materials/localintegralenergy.hh>
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
class HenckyEnergy
: public Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension> > >
class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::HenckyDensity")
HenckyEnergy
: public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
{
// grid types
typedef typename GridView::Grid::ctype DT;
typedef typename GridView::template Codim<0>::Entity Entity;
// some other sizes
enum {gridDim=GridView::dimension};
enum {dim=GridView::dimension};
public: // for testing
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
HenckyEnergy(const Dune::ParameterTree& parameters)
{
// Lame constants
mu_ = parameters.template get<double>("mu");
lambda_ = parameters.template get<double>("lambda");
kappa_ = (2.0 * mu_ + 3.0 * lambda_) / 3.0;
}
using Base = Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>;
public:
/** \brief Assemble the energy for a single element */
field_type energy (const Entity& e,
const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
/** \brief Lame constants */
double mu_, lambda_, kappa_;
// backwards compatibility: forward directly to integral energy
HenckyEnergy(const Dune::ParameterTree& parameters)
: Base(std::make_shared<Elasticity::HenckyDensity<GridView::dimension,field_type>>(parameters))
{}
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
HenckyEnergy<GridView,LocalFiniteElement,field_type>::
energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const
{
assert(element.type() == localFiniteElement.type());
typedef typename GridView::template Codim<0>::Entity::Geometry Geometry;
field_type energy = 0;
// store gradients of shape functions and base functions
std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim;
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
const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
const DT integrationElement = element.geometry().integrationElement(quadPos);
const typename Geometry::JacobianInverseTransposed& jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
DT weight = quad[pt].weight() * integrationElement;
// get gradients of shape functions
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
// compute gradients of base functions
for (size_t i=0; i<gradients.size(); ++i)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
for (size_t i=0; i<gradients.size(); i++)
for (int j=0; j<gridDim; j++)
derivative[j].axpy(localConfiguration[i][j], gradients[i]);
/////////////////////////////////////////////////////////
// compute C = F^TF
/////////////////////////////////////////////////////////
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++)
C[i][j] += derivative[k][i] * derivative[k][j];
//////////////////////////////////////////////////////////
// Eigenvalues of C = F^TF
//////////////////////////////////////////////////////////
Dune::FieldVector<field_type,dim> sigmaSquared;
FMatrixHelp::eigenValues(C, sigmaSquared);
// 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]);
field_type trLogU = 0;
for (int i = 0; i < dim; i++)
trLogU += lnSigma[i];
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 distortionEnergy + dilatationEnergy;
}
} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_HENCKYENERGY_HH
......
#ifndef DUNE_ELASTICITY_MATERIALS_LOCALDENSITY_HH
#define DUNE_ELASTICITY_MATERIALS_LOCALDENSITY_HH
#include <dune/common/fmatrix.hh>
namespace Dune::Elasticity {
/** \brief A base class for energy densities to be evaluated in an integral energy */
template<int dim, class field_type = double>
class LocalDensity
{
public:
/** \brief Evaluation with the deformation gradient
*
* \param x Position of the gradient
* \param gradient Deformation gradient
*/
virtual field_type operator() (const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& gradient) const = 0;
};
} // namespace Dune::Elasticity
#endif // DUNE_ELASTICITY_MATERIALS_LOCALDENSITY_HH
#ifndef DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localenergy.hh>
#include <dune/elasticity/materials/localdensity.hh>
namespace Dune::Elasticity {
template<class GridView, class LocalFiniteElement, class field_type=double>
class LocalIntegralEnergy
: public Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,GridView::dimension>>>
{
// grid types
using DT = typename GridView::Grid::ctype;
using Entity = typename GridView::template Codim<0>::Entity;
// some other sizes
enum {gridDim=GridView::dimension};
enum {dim=GridView::dimension};
public:
/** \brief Constructor with a local energy density
*/
LocalIntegralEnergy(const std::shared_ptr<LocalDensity<dim,field_type>>& ld)
: localDensity_(ld)
{}
/** \brief Assemble the energy for a single element */
field_type energy (const Entity& e,
const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
protected:
const std::shared_ptr<LocalDensity<dim,field_type>> localDensity_ = nullptr;