Commit bf351cb3 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 SolidDeltaResidual: incorporate some methods that were in hyperelasticity...

#820 SolidDeltaResidual: incorporate some methods that were in hyperelasticity policy in SecondPiolKirchhoffStressTensor.
parent e20b13f8
......@@ -103,6 +103,14 @@ namespace HappyHeart
using deriv_green_lagrange_type =
::HappyHeart::LocalVariationalOperatorNS::Private::DerivativeGreenLagrange<::HappyHeart::LocalVariationalOperatorNS::Private::GreenLagrangeOrEta::green_lagrange>;
//! Alias to invariant manager.
using invariant_holder_type = InvariantHolder;
//! Alias to Cauchy Green tensor.
using cauchy_green_tensor_type = CauchyGreenTensor;
public:
/// \name Special members.
......@@ -170,19 +178,13 @@ namespace HappyHeart
* \brief Part of ComputeEltArray once all the dimension depending operations have been performed.
*/
void ComputeAtQuadraturePoint(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const LocalMatrix& tangent_matrix,
const LocalVector& rhs_part,
elementary_data_type& elementary_data) const;
//! Compute the derivates of W.
template<unsigned int DimensionT>
void ComputeWDerivates(const InformationsAtQuadraturePoint& infos_at_quad_pt,
void ComputeWDerivates(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& ref_felt,
const LocalMatrix& gradient_displacement,
const CauchyGreenTensor& cauchy_green_tensor,
const LocalMatrix& De,
const LocalMatrix& transposed_De,
const cauchy_green_tensor_type& cauchy_green_tensor,
LocalVector& dW,
LocalMatrix& d2W);
......@@ -251,6 +253,36 @@ namespace HappyHeart
};
///@}
private:
//! Access to the invariant manager.
const invariant_holder_type& GetInvariantHolder() const noexcept;
//! Non constant access to the invariant manager.
invariant_holder_type& GetNonCstInvariantHolder() noexcept;
//! Helper matrix to store outer product.
LocalMatrix& GetNonCstWorkMatrixOuterProduct() noexcept;
//! Underlying hyperelasticity law.
const HyperelasticLawT& GetHyperelasticLaw() const noexcept;
private:
//! Invariant manager.
invariant_holder_type::unique_ptr invariant_holder_ = nullptr;
private:
//! Helper matrix to store outer product.
LocalMatrix work_matrix_outer_product_;
//! Underlying hyperelasticity law.
const HyperelasticLawT* const hyperelastic_law_;
};
......
......@@ -45,6 +45,7 @@ namespace HappyHeart
const TimeManager& time_manager,
const HyperelasticLawT* hyperelastic_law)
: parent(unknown_list, std::move(a_elementary_data)),
hyperelastic_law_(hyperelastic_law),
matrix_parent(),
vector_parent()
{
......@@ -101,6 +102,22 @@ namespace HappyHeart
}});
former_local_displacement_.resize(elementary_data.Ndof());
switch(mesh_dimension)
{
case 2u:
GetNonCstWorkMatrixOuterProduct().Resize(3, 3);
break;
case 3u:
GetNonCstWorkMatrixOuterProduct().Resize(6, 6);
break;
default:
assert(false);
break;
}
invariant_holder_ = std::make_unique<invariant_holder_type>(mesh_dimension,
InvariantHolderNS::Content::invariants_and_first_and_second_deriv);
}
......@@ -154,7 +171,6 @@ namespace HappyHeart
template <class HyperelasticLawT>
void SolidDeltaResidual<HyperelasticLawT>
::ComputeAtQuadraturePoint(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const LocalMatrix& tangent_matrix,
const LocalVector& rhs_part,
elementary_data_type& elementary_data) const
{
......@@ -194,29 +210,108 @@ namespace HappyHeart
template <class HyperelasticLawT>
template<unsigned int DimensionT>
void SolidDeltaResidual<HyperelasticLawT>
::ComputeWDerivates(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& ref_felt,
const LocalMatrix& gradient_displacement,
const CauchyGreenTensor& cauchy_green_tensor,
const LocalMatrix& De,
const LocalMatrix& transposed_De,
LocalVector& dW,
LocalMatrix& d2W)
void SolidDeltaResidual<HyperelasticLawT>::ComputeWDerivates(const QuadraturePoint& quad_pt,
const GeometricElt& geom_elt,
const cauchy_green_tensor_type& cauchy_green_tensor,
LocalVector& dW,
LocalMatrix& d2W)
{
static_cast<void>(gradient_displacement);
// static_cast<void>(pressure); // See #7.
auto& invariant_holder = this->GetNonCstInvariantHolder();
invariant_holder.template Update(cauchy_green_tensor);
const auto& law = GetHyperelasticLaw();
const auto& dI1dC = invariant_holder.template GetFirstDerivativeWrtCauchyGreen<1>();
const auto& dI2dC = invariant_holder.template GetFirstDerivativeWrtCauchyGreen<2>();
const auto& dI3dC = invariant_holder.template GetFirstDerivativeWrtCauchyGreen<3>();
const double dWdI1 = law.FirstDerivativeWFirstInvariant(invariant_holder, quad_pt, geom_elt);
const double dWdI2 = law.FirstDerivativeWSecondInvariant(invariant_holder, quad_pt, geom_elt);
const double dWdI3 = law.FirstDerivativeWThirdInvariant(invariant_holder, quad_pt, geom_elt);
{
const auto& d2I1dC = invariant_holder.template GetSecondDerivativeWrtCauchyGreen<1>();
const auto& d2I2dC = invariant_holder.template GetSecondDerivativeWrtCauchyGreen<2>();
const auto& d2I3dC = invariant_holder.template GetSecondDerivativeWrtCauchyGreen<3>();
{
assert(Wrappers::Seldon::IsZeroMatrix(d2W));
Seldon::Add(dWdI1, d2I1dC, d2W);
Seldon::Add(dWdI2, d2I2dC, d2W);
Seldon::Add(dWdI3, d2I3dC, d2W);
# ifndef NDEBUG
Wrappers::Seldon::ThrowIfNanInside(d2W, __FILE__, __LINE__);
# endif // NDEBUG
}
{
const double d2WdI1dI1 =
law.SecondDerivativeWFirstInvariant(invariant_holder, quad_pt, geom_elt);
const double d2WdI2dI2 =
law.SecondDerivativeWSecondInvariant(invariant_holder, quad_pt, geom_elt);
const double d2WdI3dI3 =
law.SecondDerivativeWThirdInvariant(invariant_holder, quad_pt, geom_elt);
const double d2WdI1dI2 =
law.SecondDerivativeWFirstAndSecondInvariant(invariant_holder, quad_pt, geom_elt);
const double d2WdI1dI3 =
law.SecondDerivativeWFirstAndThirdInvariant(invariant_holder, quad_pt, geom_elt);
const double d2WdI2dI3 =
law.SecondDerivativeWSecondAndThirdInvariant(invariant_holder, quad_pt, geom_elt);
using namespace Wrappers::Seldon;
auto& outer_prod = GetNonCstWorkMatrixOuterProduct();
OuterProd(dI1dC, dI1dC, outer_prod);
Seldon::Add(d2WdI1dI1, outer_prod, d2W);
OuterProd(dI1dC, dI2dC, outer_prod);
Seldon::Add(d2WdI1dI2, outer_prod, d2W);
OuterProd(dI1dC, dI3dC, outer_prod);
Seldon::Add(d2WdI1dI3, outer_prod, d2W);
OuterProd(dI2dC, dI1dC, outer_prod);
Seldon::Add(d2WdI1dI2, outer_prod, d2W);
OuterProd(dI2dC, dI2dC, outer_prod);
Seldon::Add(d2WdI2dI2, outer_prod, d2W);
OuterProd(dI2dC, dI3dC, outer_prod);
Seldon::Add(d2WdI2dI3, outer_prod, d2W);
OuterProd(dI3dC, dI1dC, outer_prod);
Seldon::Add(d2WdI1dI3, outer_prod, d2W);
OuterProd(dI3dC, dI2dC, outer_prod);
Seldon::Add(d2WdI2dI3, outer_prod, d2W);
OuterProd(dI3dC, dI3dC, outer_prod);
Seldon::Add(d2WdI3dI3, outer_prod, d2W);
}
Seldon::Mlt(4., d2W);
}
const auto& hyperelasticity_ptr = static_cast<HyperelasticLawT*>(this);
auto& hyperelasticity = *hyperelasticity_ptr;
hyperelasticity.ComputeWDerivates(infos_at_quad_pt.GetQuadraturePoint(),
geom_elt,
cauchy_green_tensor,
dW,
d2W);
{
assert(Wrappers::Seldon::IsZeroVector(dW));
Seldon::Add(dWdI1, dI1dC, dW);
Seldon::Add(dWdI2, dI2dC, dW);
Seldon::Add(dWdI3, dI3dC, dW);
Seldon::Mlt(2., dW);
}
#ifndef NDEBUG
invariant_holder.Reset();
#endif // NDEBUG
}
template <class HyperelasticLawT>
......@@ -330,17 +425,12 @@ namespace HappyHeart
dW.Zero();
d2W.Zero();
ComputeWDerivates<DimensionT>(infos_at_quad_pt,
geom_elt,
ref_felt,
gradient_displacement,
cauchy_green_tensor,
De,
transposed_De,
dW,
d2W);
ComputeWDerivates(infos_at_quad_pt.GetQuadraturePoint(),
geom_elt,
cauchy_green_tensor,
dW,
d2W);
# ifndef NDEBUG
Wrappers::Seldon::ThrowIfNanInside(d2W, __FILE__, __LINE__);
# endif // NDEBUG
......@@ -354,6 +444,43 @@ namespace HappyHeart
}
template <class HyperelasticLawT>
inline const typename SolidDeltaResidual<HyperelasticLawT>::invariant_holder_type&
SolidDeltaResidual<HyperelasticLawT>
::GetInvariantHolder() const noexcept
{
assert(!(!invariant_holder_));
return *invariant_holder_;
}
template <class HyperelasticLawT>
inline typename SolidDeltaResidual<HyperelasticLawT>::invariant_holder_type&
SolidDeltaResidual<HyperelasticLawT>
::GetNonCstInvariantHolder() noexcept
{
return const_cast<typename SolidDeltaResidual<HyperelasticLawT>::invariant_holder_type&>(GetInvariantHolder());
}
template <class HyperelasticLawT>
inline LocalMatrix& SolidDeltaResidual<HyperelasticLawT>::GetNonCstWorkMatrixOuterProduct() noexcept
{
return work_matrix_outer_product_;
}
template <class HyperelasticLawT>
inline const HyperelasticLawT& SolidDeltaResidual<HyperelasticLawT>::GetHyperelasticLaw() const noexcept
{
assert(!(!hyperelastic_law_));
return *hyperelastic_law_;
}
} // namespace LocalVariationalOperatorNS
......
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