diff --git a/docs/web/release-notes-5.1.md b/docs/web/release-notes-5.1.md index 0e4827130..12a26359c 100644 --- a/docs/web/release-notes-5.1.md +++ b/docs/web/release-notes-5.1.md @@ -76,6 +76,18 @@ $ tfel-config-5.1.0-release --python-module-suffix 5_1_0_release ~~~~ +# New `TFEL/Material` features + +## Homogenization + +### Ellipsoidal inclusion embedded in anisotropic matrix + +When \(\tenseur C_0\) is anisotropic, the Eshelby tensor can be computed +with `computeAnisotropicEshelbyTensor` in 3D and `compute2DAnisotropicEshelbyTensor` +in 2D. There are also `computeAnisotropicHillTensor`, `compute2DAnisotropicHillTensor`, +and also `computeAnisotropicLocalisationTensor` and `compute2DAnisotropicLocalisationTensor`. + + # Python bindings Python bindings are now generated using the diff --git a/docs/web/tfel-material.md b/docs/web/tfel-material.md index f522b8e3c..3e05aa0f0 100644 --- a/docs/web/tfel-material.md +++ b/docs/web/tfel-material.md @@ -384,7 +384,7 @@ A specialization for elasticity is defined: `tfel::material::homogenization::ela ## Eshelby tensors -The header `Eshelby.hxx` introduces +The header `IsotropicEshelbyTensor.hxx` introduces the function `computeEshelbyTensor` which computes the Eshelby tensor of an ellipsoid. If we consider a constant stress-free strain \(\tenseur \varepsilon^\mathrm{T}\) @@ -423,9 +423,15 @@ for `long double`. In the same way, the formulas for the axisymmetrical case are instable when the aspect ratio is near one, so a parameter allows to switch to the formula for a sphere. +When \(\tenseur C_0\) is anisotropic, the Eshelby tensor can be computed +with `computeAnisotropicEshelbyTensor` in 3D and `compute2DAnisotropicEshelbyTensor` +in 2D. There are also `computeAnisotropicHillTensor`, `compute2DAnisotropicHillTensor`, +and also `computeAnisotropicLocalisationTensor` and `compute2DAnisotropicLocalisationTensor`. +These functions are introduced by the header `AnisotropicEshelbyTensor.hxx`. + ## Strain localisation tensors -The header `Eshelby.hxx` also introduces +The header `IsotropicEshelbyTensor.hxx` also introduces three functions that compute the strain localisation tensor of an ellipsoid. If we consider an ellipsoid whose elasticity is \(\tenseur C_i\), embedded in an infinite homogeneous medium whose elasticity is \(\tenseur C_0\), @@ -446,9 +452,10 @@ The functions then return the localisation tensors taking into account the orien ## Homogenization schemes Different schemes are implemented and return the homogenized stiffness of the material. +These schemes are introduced by the header `LinearHomogenizationSchemes`. The scheme available are Mori-Tanaka scheme and dilute scheme. The available functions are `computeMoriTanakaScheme`, `computeDiluteScheme`, -computeSphereDiluteScheme, computeSphereMoriTanakaScheme. +`computeSphereDiluteScheme`, `computeSphereMoriTanakaScheme`. If a distribution of ellipsoids is considered, three types of distributions are considered : isotropic, transverse isotropic and with unique orientation. diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt index 97f1e4def..3a7d81189 100644 --- a/include/CMakeLists.txt +++ b/include/CMakeLists.txt @@ -500,10 +500,12 @@ install_header(TFEL/Material OrthotropicPlasticity.ixx) install_header(TFEL/Material Lame.hxx) install_header(TFEL/Material Hill.hxx) install_header(TFEL/Material Hill.ixx) -install_header(TFEL/Material Eshelby.hxx) -install_header(TFEL/Material Eshelby.ixx) -install_header(TFEL/Material EshelbyBasedHomogenization.ixx) -install_header(TFEL/Material EshelbyBasedHomogenization.hxx) +install_header(TFEL/Material IsotropicEshelbyTensor.hxx) +install_header(TFEL/Material IsotropicEshelbyTensor.ixx) +install_header(TFEL/Material AnisotropicEshelbyTensor.hxx) +install_header(TFEL/Material AnisotropicEshelbyTensor.ixx) +install_header(TFEL/Material LinearHomogenizationSchemes.ixx) +install_header(TFEL/Material LinearHomogenizationSchemes.hxx) install_header(TFEL/Material OrthotropicStressLinearTransformation.hxx) install_header(TFEL/Material OrthotropicStressLinearTransformation.ixx) install_header(TFEL/Material Hosford1972YieldCriterion.hxx) diff --git a/include/TFEL/Material/AnisotropicEshelbyTensor.hxx b/include/TFEL/Material/AnisotropicEshelbyTensor.hxx new file mode 100644 index 000000000..2a0451ef4 --- /dev/null +++ b/include/TFEL/Material/AnisotropicEshelbyTensor.hxx @@ -0,0 +1,192 @@ +/*! + * \file include/TFEL/Material/AnisotropicEshelbyTensor.hxx + * \author Antoine Martin + * \date 18 November 2024 + * \brief This file declares the Eshelby tensor for an ellipsoidal inclusion + * embedded in an anisotropic matrix. \copyright Copyright (C) 2006-2018 + * CEA/DEN, EDF R&D. All rights reserved. This project is publicly released + * under either the GNU GPL Licence or the CECILL-A licence. A copy of thoses + * licences are delivered with the sources of TFEL. CEA or EDF may also + * distribute this project under specific licensing conditions. + */ + +#ifndef LIB_TFEL_MATERIAL_ANISOTROPICESHELBYTENSOR_HXX +#define LIB_TFEL_MATERIAL_ANISOTROPICESHELBYTENSOR_HXX + +#include "TFEL/Math/st2tost2.hxx" +#include "TFEL/Material/StiffnessTensor.hxx" +#include + +namespace tfel::material::homogenization::elasticity { + + /*! + * This function builds the Hill tensor of a general 2d ellipse embedded + * in an ANISOTROPIC matrix. The + * ellipse also has a specific orientation given by the vector + * \f$n_a\f$. A LOCAL basis can then be defined. + * The function returns the Hill tensor in th GLOBAL basis. + * \return an object of type + * st2tost2<2u,tfel::math::invert_type> \tparam real: + * underlying type \tparam StressType: type of the elastic constants + * \param[in] C: stiffness tensor of the matrix + * \param [in] n_a: direction of the principal axis whose length is + * \f$a\f$ \param [in] a: length of semi-axis relative to the direction + * \f$n_a\f$ \param [in] b: length of semi-axis relative to the other + * direction \param[in] max_it: maximal number of iterations for integration + */ + template + TFEL_HOST_DEVICE static tfel::math::st2tost2< + 2u, + typename tfel::math::invert_type> + compute2DAnisotropicHillTensor(const tfel::math::st2tost2<2u, StressType>&, + const real&, + const std::size_t max_it = 12); + + /*! + * This function builds the Eshelby tensor of a general 2d ellipse + * embedded in an ANISOTROPIC matrix. The + * ellipse also has a specific orientation given by the vector + * \f$n_a\f$. A LOCAL basis can then be defined. + * The function returns the Hill tensor in th GLOBAL basis. + * The function returns the Eshelby + * tensor in the global basis \return an object of + * type st2tost2<2u,real> \tparam real: underlying type \tparam + * StressType: type of the elastic constants \param[in] C: stiffness + * tensor of the matrix + * * \param [in] n_a: direction of the principal axis whose length is + * \f$a\f$ \param [in] a: length of semi-axis relative to the direction + * \f$n_a\f$ \param [in] b: length of semi-axis relative to the other + * direction \param[in] max_it: maximal number of iterations for integration + */ + template + TFEL_HOST_DEVICE static tfel::math::st2tost2<2u, real> + compute2DAnisotropicEshelbyTensor(const tfel::math::st2tost2<2u, StressType>&, + const real&, + const std::size_t max_it = 12); + + /*! + * This function builds the Hill tensor of a general ellipsoid embedded in + * an ANISOTROPIC matrix. The function returns the Hill tensor in the + * global basis \return an object + * of type st2tost2<3u,tfel::math::invert_type> \tparam real: + * underlying type \tparam LengthType: type of the dimensions of the + * ellipsoid \tparam StressType: type of the elastic constants \param[in] + * C: stiffness tensor of the matrix in the global basis + * \param [in] n_a: direction of the principal axis whose length is + * \f$a\f$ \param [in] a: length of semi-axis relative to the direction + * \f$n_a\f$ \param [in] n_b: direction of the principal axis whose length + * is \f$b\f$ \param [in] b: length of semi-axis relative to the direction + * \f$n_b\f$ \param [in] c: length of the remaining semi-axis + * \param[in] max_it: maximal number of iterations + * for integration + * + */ + template + TFEL_HOST_DEVICE static tfel::math:: + st2tost2<3u, typename tfel::math::invert_type> + computeAnisotropicHillTensor(const tfel::math::st2tost2<3u, StressType>&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const LengthType&, + const std::size_t max_it = 12); + + /*! + * This function builds the Eshelby tensor of a general ellipsoid embedded + * in an ANISOTROPIC matrix. The function returns the Eshelby tensor in + * the global basi + * \return an object of type st2tost2<3u,real> + * \tparam real: underlying type + * \tparam LengthType: type of the dimensions of the ellipsoid + * \tparam StressType: type of the elastic constants + * \param[in] C: stiffness tensor of the + * matrix in the global basis + * \param [in] n_a: direction of the principal axis whose length is + * \f$a\f$ \param [in] a: length of semi-axis relative to the direction + * \f$n_a\f$ \param [in] n_b: direction of the principal axis whose length + * is \f$b\f$ \param [in] b: length of semi-axis relative to the direction + * \f$n_b\f$ \param [in] c: length of the remaining semi-axis + * \param[in] max_it: maximal number of iterations for integration + * + */ + template + TFEL_HOST_DEVICE static tfel::math::st2tost2<3u, real> + computeAnisotropicEshelbyTensor(const tfel::math::st2tost2<3u, StressType>&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const LengthType&, + const std::size_t max_it = 12); + + /*! + * This function builds the strain localisation tensor of a general + * ellipsoid embedded in an ANISOTROPIC matrix. The localisation tensor + * \f$A\f$ is defined as follows : \f[\epsilon = A:E_0\f] where \f$E_0\f$ + * is the uniform strain tensor imposed at infinity, and \f$\epsilon\f$ is + * the strain tensor solution of Eshelby problem for the ellipsoid. The + * ellipsoid also has a specific orientation given by the vectors + * \f$n_a\f$, \f$n_b\f$. A LOCAL basis can then be defined as n_a, n_b, + * n_c where n_c=cross_product(n_a,n_b) \return an object of type + * st2tost2<3u,real>, which is the fourth-order localisation tensor + * \f$A\f$ in the GLOBAL BASIS \tparam real: underlying type \tparam + * StressType: type of the elastic constants related to the matrix and the + * ellipsoid \tparam LengthType: type of the dimensions of the ellipsoid + * \param[in] C_0_glob: stiffness tensor of the matrix in the GLOBAL basis + * \param[in] C_i_loc: stiffness tensor of the inclusion in the LOCAL + * basis, which is the basis (n_a,n_b,n_c) + * \param [in] n_a: direction of the principal axis whose length is + * \f$a\f$ \param [in] a: length of semi-axis relative to the direction + * \f$n_a\f$ \param [in] n_b: direction of the principal axis whose length + * is \f$b\f$ \param [in] b: length of semi-axis relative to the direction + * \f$n_b\f$ \param [in] c: length of the remaining semi-axis + * \param[in] max_it: maximal number of iterations for integration + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> + computeAnisotropicLocalisationTensor( + const tfel::math::st2tost2<3u, StressType>&, + const tfel::math::st2tost2<3u, StressType>&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const LengthType&, + const std::size_t max_it = 12); + + /*! + * This function builds the strain localisation tensor of a general 2d ellipse + * embedded in an ANISOTROPIC matrix. The localisation tensor + * \f$A\f$ is defined as follows : \f[\epsilon = A:E_0\f] where \f$E_0\f$ + * is the uniform strain tensor imposed at infinity, and \f$\epsilon\f$ is + * the strain tensor solution of Eshelby problem for the ellipse. The + * ellipse also has a specific orientation given by the vector + * \f$n_a\f$. A LOCAL basis can then be defined. + * \return an object of type + * st2tost2<2u,real>, which is the fourth-order localisation tensor + * \f$A\f$ in the GLOBAL BASIS \tparam real: underlying type \tparam + * StressType: type of the elastic constants + * \param[in] C_0_glob: stiffness tensor of the matrix in the GLOBAL basis + * \param[in] C_i_loc: stiffness tensor of the inclusion in the LOCAL + * basis, which is the basis related to \f$n_a\f$ + * \param [in] n_a: direction of the principal axis whose length is + * \f$a\f$ \param [in] a: length of semi-axis relative to the direction + * \f$n_a\f$ \param [in] b: length of semi-axis relative to the other + * direction \param[in] max_it: maximal number of iterations for integration + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> + compute2DAnisotropicLocalisationTensor( + const tfel::math::st2tost2<2u, StressType>&, + const tfel::math::st2tost2<2u, StressType>&, + const tfel::math::tvector<2u, real>&, + const LengthType&, + const LengthType&, + const std::size_t max_it = 12); + +} // namespace tfel::material::homogenization::elasticity + +#include "TFEL/Material/AnisotropicEshelbyTensor.ixx" + +#endif /* LIB_TFEL_MATERIAL_ANISOTROPICESHELBYTENSOR_HXX */ diff --git a/include/TFEL/Material/AnisotropicEshelbyTensor.ixx b/include/TFEL/Material/AnisotropicEshelbyTensor.ixx new file mode 100644 index 000000000..5ad11a4eb --- /dev/null +++ b/include/TFEL/Material/AnisotropicEshelbyTensor.ixx @@ -0,0 +1,451 @@ +/*! + * \file include/TFEL/Material/AnisotropicEshelbyTensor.ixx + * \author Antoine Martin + * \date 18 November 2024 + * \brief This file defines the Eshelby tensor for an ellipsoidal inclusion + * embedded in an anisotropic matrix. \copyright Copyright (C) 2006-2018 + * CEA/DEN, EDF R&D. All rights reserved. This project is publicly released + * under either the GNU GPL Licence or the CECILL-A licence. A copy of thoses + * licences are delivered with the sources of TFEL. CEA or EDF may also + * distribute this project under specific licensing conditions. + */ + +#ifndef LIB_TFEL_MATERIAL_ANISOTROPICESHELBYTENSOR_IXX +#define LIB_TFEL_MATERIAL_ANISOTROPICESHELBYTENSOR_IXX + +#include +#include +#include + +namespace tfel::material::homogenization::elasticity { + + namespace internals { + + template + TFEL_HOST_DEVICE auto integrate1D( + const std::function f, + const real a, + const real b, + const std::size_t max_it = 12, + const real tol = std::sqrt(std::numeric_limits::epsilon())) { + auto h = (b - a) / 2; + auto I0 = (f(a) + f(b)) * h; + auto IL0 = (std::abs(f(a)) + std::abs(f(b))) * h; + auto I1 = I0 / 2 + f(a + h) * h; + auto IL1 = IL0 / 2 + std::abs(f(a + h)) * h; + + std::size_t k = 2; + std::size_t p = 2; + real error = std::abs(I0 - I1); + while (k < 4 || (k < max_it && error > tol * IL1)) { + I0 = I1; + IL0 = IL1; + + I1 = I0 / 2; + IL1 = IL0 / 2; + p *= 2; + h /= 2; + + for (std::size_t j = 1; j < p; j += 2) { + I1 += f(a + j * h) * h; + IL1 += std::abs(f(a + j * h)) * h; + } + + ++k; + error = std::abs(I0 - I1); + } + return I1; + } + + template + TFEL_HOST_DEVICE auto dblintegrate( + const std::function f, + const real a, + const real b, + const real c, + const real d, + const std::size_t max_it = 12) { + const auto f_ = [f, c, d, max_it](const real& x) { + return integrate1D([x, f](const real& y) { return f(x, y); }, c, + d, max_it); + }; + return integrate1D(f_, a, b, max_it); + } + + TFEL_HOST_DEVICE int vi(unsigned short int i, unsigned short int j) { + if ((i == 0) and (j == 0)) { + return 0; + } else if ((i == 1) and (j == 1)) { + return 1; + } else if ((i == 2) and (j == 2)) { + return 2; + } else if (((i == 0) and (j == 1)) || ((i == 1) and (j == 0))) { + return 3; + } else if (((i == 0) and (j == 2)) || ((i == 2) and (j == 0))) { + return 4; + } else if (((i == 1) and (j == 2)) || ((i == 2) and (j == 1))) { + return 5; + } + } + + template + TFEL_HOST_DEVICE Type + getST2toST2Component(const tfel::math::st2tost2& A, + unsigned short int i, + unsigned short int j, + unsigned short int k, + unsigned short int l) { + const int I = vi(i, j); + const int J = vi(k, l); + auto fac = real(1); + constexpr auto icste = tfel::math::Cste::isqrt2; + if (I > 2) { + fac *= icste; + }; + if (J > 2) { + fac *= icste; + }; + return fac * A(I, J); + } + + template + TFEL_HOST_DEVICE void setST2toST2Component(tfel::math::st2tost2& A, + unsigned short int i, + unsigned short int j, + unsigned short int k, + unsigned short int l, + Type Aijkl) { + const int I = vi(i, j); + const int J = vi(k, l); + auto fac = real(1); + constexpr auto cste = tfel::math::Cste::sqrt2; + if (I > 2) { + fac *= cste; + }; + if (J > 2) { + fac *= cste; + }; + A(I, J) = fac * Aijkl; + } + + template + TFEL_HOST_DEVICE void setStensorComponent(tfel::math::stensor& A, + unsigned short int i, + unsigned short int j, + Type Aij) { + const int I = vi(i, j); + auto fac = real(1); + constexpr auto cste = tfel::math::Cste::sqrt2; + if (I > 2) { + fac *= cste; + }; + A(I) = fac * Aij; + } + + template + TFEL_HOST_DEVICE Type + getStensorComponent(const tfel::math::stensor& A, + unsigned short int i, + unsigned short int j) { + const int I = vi(i, j); + auto fac = real(1); + constexpr auto icste = tfel::math::Cste::isqrt2; + if (I > 2) { + fac *= icste; + }; + return fac * A(I); + } + + template + TFEL_HOST_DEVICE tfel::math::tmatrix Acoustic( + const tfel::math::st2tost2& C, + const tfel::math::tvector& X) { + tfel::math::tmatrix A; + for (int i = 0; i < N; i++) + for (int k = i; k < N; k++) { + real A_ik = real(0); + for (int j = 0; j < N; j++) + for (int l = 0; l < N; l++) { + A_ik += real( + (getST2toST2Component(C, i, j, k, l)) / + StressType(1) * X[j] * X[l]); + } + A(i, k) = A_ik; + if (i != k) { + A(k, i) = A_ik; + } + } + return A; + } + + template + TFEL_HOST_DEVICE real + p_ijkl_2D(const tfel::math::st2tost2<2u, StressType>& C, + const LengthType& a, + const LengthType& b, + const real theta, + const int& i, + const int& j, + const int& k, + const int& l) { + const real pi = std::numbers::pi_v; + const tfel::math::tvector<2u, real> X = { + std::cos(theta) / a * LengthType(1), + std::sin(theta) / b * LengthType(1)}; + auto A_inv = Acoustic<2u, real, StressType>(C, X); + tfel::math::TinyMatrixInvert<2u, real>::exe(A_inv); + const auto Mijkl = + (A_inv(j, k) * X[i] * X[l] + A_inv(i, k) * X[j] * X[l] + + A_inv(j, l) * X[i] * X[k] + A_inv(i, l) * X[j] * X[k]) / + 4; + return Mijkl / 2 / pi; + } + + template + TFEL_HOST_DEVICE real p_ijkl(const tfel::math::st2tost2<3u, StressType>& C, + const real theta, + const real phi, + const LengthType& a, + const LengthType& b, + const LengthType& c, + const int& i, + const int& j, + const int& k, + const int& l) { + const real pi = std::numbers::pi_v; + const tfel::math::tvector<3u, real> X = { + std::sin(theta) * std::cos(phi) / a * LengthType(1), + std::sin(theta) * std::sin(phi) / b * LengthType(1), + std::cos(theta) / c * LengthType(1)}; + auto A_inv = Acoustic<3u, real, StressType>(C, X); + tfel::math::TinyMatrixInvert<3u, real>::exe(A_inv); + const auto Mijkl = + (A_inv(j, k) * X[i] * X[l] + A_inv(i, k) * X[j] * X[l] + + A_inv(j, l) * X[i] * X[k] + A_inv(i, l) * X[j] * X[k]) / + 4; + return Mijkl * std::sin(theta) / 4 / pi; + } + + } // end of namespace internals + + template + TFEL_HOST_DEVICE static tfel::math::st2tost2< + 2u, + typename tfel::math::invert_type> + compute2DAnisotropicHillTensor(const tfel::math::st2tost2<2u, StressType>& C, + const tfel::math::tvector<2u, real>& n_a, + const LengthType& a, + const LengthType& b, + const std::size_t max_it) { + if (not((a > LengthType{0}) and (b > LengthType{0}))) { + tfel::reportContractViolation("a<=0 or b<=0"); + }; + if (tfel::math::ieee754::fpclassify(norm(n_a)) == FP_ZERO) { + tfel::reportContractViolation("n_a is null"); + }; + using namespace tfel::math; + const auto n_a_ = n_a / norm(n_a); + tfel::math::tvector<2u, real> n_b_; + if (n_a_[1] != real(0)) { + n_b_ = {real(1), -n_a_[0] / n_a_[1]}; + } else { + n_b_ = {real(0), real(1)}; + } + const rotation_matrix r_glob_loc = {n_a_[0], n_b_[0], real(0), + n_a_[1], n_b_[1], real(0), + real(0), real(0), real(1)}; + const rotation_matrix r_loc_glob = transpose(r_glob_loc); + const auto C_loc = + StressType(1) * change_basis(C / StressType(1), r_glob_loc); + const real pi = std::numbers::pi_v; + const real zero = real(0); + using compliance = typename tfel::math::invert_type; + tfel::math::st2tost2<2u, compliance> P; + for (int i = 0; i < 2; i++) + for (int j = i; j < 2; j++) + for (int k = 0; k < 2; k++) + for (int l = k; l < 2; l++) { + const int I = internals::vi(i, j); + const int J = internals::vi(k, l); + const auto p_ = [C, a, b, i, j, k, l](const real& theta) { + return internals::p_ijkl_2D(C, a, b, theta, i, + j, k, l); + }; + const auto int_p = compliance( + internals::integrate1D(p_, zero, 2 * pi, max_it)); + internals::setST2toST2Component<2u, compliance, real>(P, i, j, k, l, + int_p); + } + return change_basis(P * StressType(1), r_loc_glob) / StressType(1); + } + + template + TFEL_HOST_DEVICE + tfel::math::st2tost2<3u, typename tfel::math::invert_type> + computeAnisotropicHillTensor( + const tfel::math::st2tost2<3u, StressType>& C, + const tfel::math::tvector<3u, real>& n_a, + const LengthType& a, + const tfel::math::tvector<3u, real>& n_b, + const LengthType& b, + const LengthType& c, + const std::size_t max_it) { + if (not((a > LengthType{0}) and (b > LengthType{0}) and + (c > LengthType{0}))) { + tfel::reportContractViolation("a<=0 or b<=0 or c<=0"); + }; + if (not(tfel::math::ieee754::fpclassify( + tfel::math::VectorVectorDotProduct::exe< + real, tfel::math::tvector<3u, real>, + tfel::math::tvector<3u, real>>(n_a, n_b)) == FP_ZERO)) { + tfel::reportContractViolation("n_a and n_b not normals"); + }; + if (tfel::math::ieee754::fpclassify(norm(n_a)) == FP_ZERO) { + tfel::reportContractViolation("n_a is null"); + }; + if (tfel::math::ieee754::fpclassify(norm(n_b)) == FP_ZERO) { + tfel::reportContractViolation("n_b is null"); + }; + using namespace tfel::math; + const auto n_a_ = n_a / norm(n_a); + const auto n_b_ = n_b / norm(n_b); + const auto n_c_ = cross_product(n_a_, n_b_); + + const rotation_matrix r_glob_loc = {n_a_[0], n_b_[0], n_c_[0], + n_a_[1], n_b_[1], n_c_[1], + n_a_[2], n_b_[2], n_c_[2]}; + const rotation_matrix r_loc_glob = transpose(r_glob_loc); + const auto C_loc = + StressType(1) * change_basis(C / StressType(1), r_glob_loc); + + const real pi = std::numbers::pi_v; + const real zero = real(0); + using compliance = typename tfel::math::invert_type; + tfel::math::st2tost2<3u, compliance> P; + for (int i = 0; i < 3; i++) + for (int j = i; j < 3; j++) + for (int k = 0; k < 3; k++) + for (int l = k; l < 3; l++) { + const int I = internals::vi(i, j); + const int J = internals::vi(k, l); + const auto p_ = [C_loc, a, b, c, i, j, k, l](const real& theta, + const real& phi) { + return internals::p_ijkl( + C_loc, theta, phi, a, b, c, i, j, k, l); + }; + const auto int_p = compliance(internals::dblintegrate( + p_, zero, pi, zero, 2 * pi, max_it)); + internals::setST2toST2Component<3u, compliance, real>(P, i, j, k, l, + int_p); + } + return change_basis(P * StressType(1), r_loc_glob) / StressType(1); + } + + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> + computeAnisotropicEshelbyTensor(const tfel::math::st2tost2<3u, StressType>& C, + const tfel::math::tvector<3u, real>& n_a, + const LengthType& a, + const tfel::math::tvector<3u, real>& n_b, + const LengthType& b, + const LengthType& c, + const std::size_t max_it) { + return computeAnisotropicHillTensor( + C, n_a, a, n_b, b, c, max_it) * + C; + } + + template + TFEL_HOST_DEVICE tfel::math::st2tost2<2u, real> + compute2DAnisotropicEshelbyTensor( + const tfel::math::st2tost2<2u, StressType>& C, + const tfel::math::tvector<2u, real>& n_a, + const LengthType& a, + const LengthType& b, + const std::size_t max_it) { + return compute2DAnisotropicHillTensor( + C, n_a, a, b, max_it) * + C; + } + + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> + computeAnisotropicLocalisationTensor( + const tfel::math::st2tost2<3u, StressType>& C_0_glob, + const tfel::math::st2tost2<3u, StressType>& C_i_loc, + const tfel::math::tvector<3u, real>& n_a, + const LengthType& a, + const tfel::math::tvector<3u, real>& n_b, + const LengthType& b, + const LengthType& c, + const std::size_t max_it) { + if (not(tfel::math::ieee754::fpclassify( + tfel::math::VectorVectorDotProduct::exe< + real, tfel::math::tvector<3u, real>, + tfel::math::tvector<3u, real>>(n_a, n_b)) == FP_ZERO)) { + tfel::reportContractViolation("n_a and n_b not normals"); + }; + if (tfel::math::ieee754::fpclassify(norm(n_a)) == FP_ZERO) { + tfel::reportContractViolation("n_a is null"); + }; + if (tfel::math::ieee754::fpclassify(norm(n_b)) == FP_ZERO) { + tfel::reportContractViolation("n_b is null"); + }; + const auto n_a_ = n_a / norm(n_a); + const auto n_b_ = n_b / norm(n_b); + const auto n_c_ = tfel::math::cross_product(n_a_, n_b_); + + using namespace tfel::math; + const rotation_matrix r_glob_loc = {n_a_[0], n_b_[0], n_c_[0], + n_a_[1], n_b_[1], n_c_[1], + n_a_[2], n_b_[2], n_c_[2]}; + const rotation_matrix r_loc_glob = transpose(r_glob_loc); + const auto P_0_glob = + computeAnisotropicHillTensor( + C_0_glob, n_a_, a, n_b_, b, c, max_it); + const auto C_i_glob = + change_basis(C_i_loc / StressType(1), r_loc_glob) * StressType(1); + const st2tost2<3u, StressType> C = C_i_glob - C_0_glob; + const auto Pr = P_0_glob * C; + const auto A = invert(st2tost2<3u, real>::Id() + Pr); + return A; + } + + template + TFEL_HOST_DEVICE tfel::math::st2tost2<2u, real> + compute2DAnisotropicLocalisationTensor( + const tfel::math::st2tost2<2u, StressType>& C_0_glob, + const tfel::math::st2tost2<2u, StressType>& C_i_loc, + const tfel::math::tvector<2u, real>& n_a, + const LengthType& a, + const LengthType& b, + const std::size_t max_it) { + if (tfel::math::ieee754::fpclassify(norm(n_a)) == FP_ZERO) { + tfel::reportContractViolation("n_a is null"); + }; + const auto n_a_ = n_a / norm(n_a); + tfel::math::tvector<2u, real> n_b_; + if (n_a_[1] != real(0)) { + n_b_ = {real(1), -n_a_[0] / n_a_[1]}; + } else { + n_b_ = {real(0), real(1)}; + } + using namespace tfel::math; + const rotation_matrix r_glob_loc = {n_a_[0], n_b_[0], real(0), + n_a_[1], n_b_[1], real(0), + real(0), real(0), real(1)}; + const rotation_matrix r_loc_glob = transpose(r_glob_loc); + const auto P_0_glob = + compute2DAnisotropicHillTensor( + C_0_glob, n_a_, a, b, max_it); + const auto C_i_glob = + change_basis(C_i_loc / StressType(1), r_loc_glob) * StressType(1); + const st2tost2<2u, StressType> C = C_i_glob - C_0_glob; + const auto Pr = P_0_glob * C; + const auto A = invert(st2tost2<2u, real>::Id() + Pr); + return A; + } + +} // namespace tfel::material::homogenization::elasticity + +#endif /* LIB_TFEL_MATERIAL_ANISOTROPICESHELBYTENSOR_IXX */ diff --git a/include/TFEL/Material/Eshelby.hxx b/include/TFEL/Material/Eshelby.hxx deleted file mode 100644 index fdd86a21a..000000000 --- a/include/TFEL/Material/Eshelby.hxx +++ /dev/null @@ -1,190 +0,0 @@ -/*! - * \file include/TFEL/Material/Eshelby.hxx - * \author Antoine Martin - * \date 15 October 2024 - * \brief This file declares the Eshelby tensor for an ellipsoidal inclusion - * embedded in an isotropic matrix. \copyright Copyright (C) 2006-2018 CEA/DEN, - * EDF R&D. All rights reserved. This project is publicly released under either - * the GNU GPL Licence or the CECILL-A licence. A copy of thoses licences are - * delivered with the sources of TFEL. CEA or EDF may also distribute this - * project under specific licensing conditions. - */ - -#ifndef LIB_TFEL_MATERIAL_ESHELBY_HXX -#define LIB_TFEL_MATERIAL_ESHELBY_HXX - -#include "TFEL/Math/st2tost2.hxx" -#include "TFEL/Material/StiffnessTensor.hxx" -#include - -namespace tfel::material { - - namespace homogenization { - namespace elasticity { - - /*! - * This function builds the Eshelby tensor of a circular cylinder embedded - * in an isotropic matrix, considering a PLANE STRAIN modelling hypothesis - * \return an object of type st2tost2<2u,real> - * \tparam real: underlying type - * \param[in] nu: Poisson's ratio of the matrix - */ - template - TFEL_HOST_DEVICE tfel::math::st2tost2<2u, real> - computeCircularCylinderEshelbyTensor(const real&); - - /*! - * This function builds the Eshelby tensor of an elliptic cylinder - * embedded in an isotropic matrix, considering a PLANE STRAIN modelling - * hypothesis The function returns the Eshelby tensor in the basis (e1,e2) - * where e1 corresponds to the biggest axis \return an object of type - * st2tost2<2u,real> \tparam real: underlying type \param[in] nu: - * Poisson's ratio of the matrix \param[in] e: aspect ratio of the - * elliptic basis - */ - template - TFEL_HOST_DEVICE tfel::math::st2tost2<2u, real> - computeEllipticCylinderEshelbyTensor(const real&, const real&); - - /*! - * This function builds the Eshelby tensor of a sphere embedded in an - * isotropic matrix. \return an object of type st2tost2<3u,real> \tparam - * real: underlying type \param[in] nu: Poisson's ratio of the matrix - */ - template - TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> - computeSphereEshelbyTensor(const real&); - - /*! - * This function builds the Eshelby tensor of an axisymmetrical ellipsoid - * embedded in an isotropic matrix. The function returns the Eshelby - * tensor in the basis (e1,e2,e3) where e1 corresponds to (one of) the - * biggest ax(es) \return an object of type st2tost2<3u,real> \tparam - * real: underlying type \param[in] nu: Poisson's ratio of the matrix - * \param[in] e: aspect ratio of the ellipsoid (e>1 : prolate, e<1 : - * oblate) \param[in] precf, precd, precld: default arguments which aim at - * preventing the numerical instability of the formula when the ellipsoid - * is almost a sphere. When the absolute value of (e-1) is below precf - * (resp. precd, precld) for real=float (resp. real= double, long double), - * the returned tensor is computeSphereEshelbyTensor(nu). - * - * The expressions can be found in Torquato, Random Heterogeneous - * Materials (2002). - */ - template - TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> - computeAxisymmetricalEshelbyTensor(const real&, - const real&, - const real = real{8e-3}, - const real = real{1.5e-4}, - const real = real{1e-5}); - - /*! - * This function builds the Eshelby tensor of a general ellipsoid embedded - * in an isotropic matrix. The function returns the Eshelby tensor in the - * basis (e1,e2,e3) where e1 (resp. e2, e3) is aligned with the axis with - * the first (resp. the second and third) biggest length \return an object - * of type st2tost2<3u,real> \tparam real: underlying type \tparam - * LengthType: type of the dimensions of the ellipsoid \param[in] nu: - * Poisson's ratio of the matrix \param[in] a: length of the first - * semi-axis \param[in] b: length of the second semi-axis \param[in] c: - * length of the third semi-axis \param[in] precf, precd, precld: default - * arguments which aim at preventing the numerical instability of the - * formula when the ellipsoid is almost axisymmetrical. When the absolute - * value of (a-b)/c (or (a-c)/b or (b-c)/a) is below precf (resp. precd, - * precld) for real=float (resp. real= double, long double), the returned - * tensor is computeAxisymmetricalEshelbyTensor. - * - * The expressions for the case of three distinct semi-axes can be found - * in Eshelby (1957). - */ - template - TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> computeEshelbyTensor( - const real&, - const LengthType&, - const LengthType&, - const LengthType&, - const real = real{8e-3}, - const real = real{1.5e-4}, - const real = real{1e-5}); - - /*! - * This function builds the strain localisation tensor of a general - * ellipsoid with a general elasticity, embedded in an isotropic matrix. - * The localisation tensor \f$A\f$ is defined as follows : \f[\epsilon = - * A:E_0\f] where \f$E_0\f$ is the uniform strain tensor imposed at - * infinity, and \f$\epsilon\f$ is the strain tensor solution of Eshelby - * problem for the ellipsoid. The ellipsoid also has a specific - * orientation given by the vectors \f$n_a\f$, \f$n_b\f$. \return an - * object of type st2tost2, which is the fourth-order localisation tensor - * \f$A\f$ \tparam real: underlying type \tparam StressType: type of the - * elastic constants related to the matrix and the ellipsoid \tparam - * LengthType: type of the dimensions of the ellipsoid \param [in] - * young,nu: Young modulus and Poisson's ratio of the matrix \param [in] - * young_i,nu_i: Young modulus and Poisson's ratio of the inclusions - * \param [in] n_a: direction of the principal axis whose length is - * \f$a\f$ \param [in] a: length of semi-axis relative to the direction - * \f$n_a\f$ \param [in] n_b: direction of the principal axis whose length - * is \f$b\f$ \param [in] b: length of semi-axis relative to the direction - * \f$n_b\f$ \param [in] c: length of the remaining semi-axis - */ - template - TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> - computeEllipsoidLocalisationTensor(const StressType&, - const real&, - const StressType&, - const real&, - const tfel::math::tvector<3u, real>&, - const LengthType&, - const tfel::math::tvector<3u, real>&, - const LengthType&, - const LengthType&); - - /*! - * This function builds the strain localisation tensor of an - * axisymmetrical ellipsoid with a general elasticity, embedded in an - * isotropic matrix. The ellipsoid also has a specific orientation given - * by the vector \f$n_a\f$, axis of the ellipsoid, whose semi-length is - * \f$a\f$. \return an object of type st2tost2 \tparam real: underlying - * type \tparam StressType: type of the elastic constants related to the - * matrix and the ellipsoid \param [in] young,nu: Young modulus and - * Poisson's ratio of the matrix \param [in] young_i,nu_i: Young modulus - * and Poisson's ratio of the inclusions \param [in] n_a: direction of the - * axis of the ellipsoid (whose semi-length is \f$a\f$) \param [in] a: - * length of semi-axis relative to the direction \f$n_a\f$ - */ - template - TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> - computeAxisymmetricalEllipsoidLocalisationTensor( - const StressType&, - const real&, - const StressType&, - const real&, - const tfel::math::tvector<3u, real>&, - const real&); - - /*! - * This function builds the strain localisation tensor of a sphere - * with a general elasticity, embedded in an isotropic matrix. - * \return an object of type st2tost2 - * \tparam real: underlying type - * \tparam StressType: type of the elastic constants related to the matrix - * and the ellipsoid \param [in] young,nu: Young modulus and Poisson's - * ratio of the matrix \param [in] young_i,nu_i: Young modulus and - * Poisson's ratio of the inclusions - */ - template - TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> - computeSphereLocalisationTensor(const StressType&, - const real&, - const StressType&, - const real&); - - } // end of namespace elasticity - } // end of namespace homogenization - -} // end of namespace tfel::material - -#include "TFEL/Material/Eshelby.ixx" - -#endif /* LIB_TFEL_MATERIAL_ESHELBY_HXX */ diff --git a/include/TFEL/Material/EshelbyBasedHomogenization.hxx b/include/TFEL/Material/EshelbyBasedHomogenization.hxx deleted file mode 100644 index b83fb5a10..000000000 --- a/include/TFEL/Material/EshelbyBasedHomogenization.hxx +++ /dev/null @@ -1,271 +0,0 @@ -/*! - * \file include/TFEL/Material/EshelbyBasedHomogenization.hxx - * \author Antoine Martin - * \date 24 October 2024 - * \brief This file declares some homogenization schemes based on solution of - * Eshelby problem. \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All - * rights reserved. This project is publicly released under either the GNU GPL - * Licence or the CECILL-A licence. A copy of thoses licences are delivered with - * the sources of TFEL. CEA or EDF may also distribute this project under - * specific licensing conditions. - */ - -#ifndef LIB_TFEL_MATERIAL_ESHELBYBASEDHOMOGENIZATION_HXX -#define LIB_TFEL_MATERIAL_ESHELBYBASEDHOMOGENIZATION_HXX - -#include "TFEL/Math/st2tost2.hxx" -#include "TFEL/Material/Eshelby.hxx" - -namespace tfel::material { - - namespace homogenization { - namespace elasticity { - /*! - * This function gives the homogenized stiffness for a dilute scheme, - * knowing the strain localisation tensor. \tparam real: underlying type - * \tparam StressType: type of the elastic constants related to the matrix - * and the inclusions \return an object of type st2tost2<3u,StressType> - * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix - * \param [in] f: volumic fraction of inclusions - * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the - * inclusions \param [in] A: mean strain localisation tensor of inclusions - */ - template - TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> - computeDiluteScheme(const StressType&, - const real&, - const real&, - const StressType&, - const real&, - const tfel::math::st2tost2<3u, real>&); - - /*! - * This function gives the homogenized stiffness for a Mori-Tanaka scheme - * knowing the strain localisation tensor. \tparam real: underlying type - * \tparam StressType: type of the elastic constants related to the matrix - * and the inclusions \return an object of type st2tost2<3u,StressType> - * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix - * \param [in] f: volumic fraction of inclusions - * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the - * inclusions \param [in] A: mean strain localisation tensor of inclusions - */ - template - TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> - computeMoriTanakaScheme(const StressType&, - const real&, - const real&, - const StressType&, - const real&, - const tfel::math::st2tost2<3u, real>&); - - /*! - * This function gives the homogenized moduli for a dilute scheme, for - * spheres \tparam real: underlying type \tparam StressType: type of the - * elastic constants related to the matrix and the spheres \return an - * object of type std::pair \param [in] young,nu: Young - * modulus and Poisson's ratio of the matrix \param [in] f: volumic - * fraction of spheres \param [in] young_i,nu_i: Young modulus and - * Poisson's ratio of the inclusions - */ - template - TFEL_HOST_DEVICE const std::pair - computeSphereDiluteScheme(const StressType&, - const real&, - const real&, - const StressType&, - const real&); - - /*! - * This function gives the homogenized moduli for a Mori-Tanaka scheme, - * for spheres \tparam real: underlying type \tparam StressType: type of - * the elastic constants related to the matrix and the spheres \return an - * object of type std::pair \param [in] young,nu: Young - * modulus and Poisson's ratio of the matrix \param [in] f: volumic - * fraction of spheres \param [in] young_i,nu_i: Young modulus and - * Poisson's ratio of the inclusions - */ - template - TFEL_HOST_DEVICE const std::pair - computeSphereMoriTanakaScheme(const StressType&, - const real&, - const real&, - const StressType&, - const real&); - - /*! - * This function gives the homogenized moduli for a dilute scheme, for an - * isotropic distribution of ellipsoids \tparam real: underlying type - * \tparam StressType: type of the elastic constants related to the matrix - * and the ellipsoids \tparam LengthType: type of the dimensions of the - * ellipsoids \return an object of type std::pair \param - * [in] young,nu: Young modulus and Poisson's ratio of the matrix \param - * [in] f: volumic fraction of ellipsoids \param [in] young_i,nu_i: Young - * modulus and Poisson's ratio of the inclusions \param[in] a: length of - * the first semi-axis \param[in] b: length of the second semi-axis - * \param[in] c: length of the third semi-axis - */ - template - TFEL_HOST_DEVICE const std::pair - computeIsotropicDiluteScheme(const StressType&, - const real&, - const real&, - const StressType&, - const real&, - const LengthType&, - const LengthType&, - const LengthType&); - - /*! - * This function gives the homogenized stiffness for a dilute scheme, for - * a transverse isotropic distribution of ellipsoids. One principal axis - * (the same principal axis for all ellipsoids) is oriented in a fixed - * direction n. The distribution of the other axes are in the plane normal - * to n, and the distribution of these axes inside this plane is - * isotropic. \tparam real: underlying type \tparam StressType: type of - * the elastic constants related to the matrix and the ellipsoids \tparam - * LengthType: type of the dimensions of the ellipsoids \return an object - * of type st2tost2<3u,StressType> \param [in] young,nu: Young modulus and - * Poisson's ratio of the matrix \param [in] f: volumic fraction of - * ellipsoids \param [in] young_i,nu_i: Young modulus and Poisson's ratio - * of the inclusions \param [in] n_a: direction of the principal axis - * which has a fixed orientation and a semi-length \f$a\f$ \param[in] a: - * length of semi-axis relative to the direction \f$n_a\f$ \param[in] b: - * length of the second semi-axis \param[in] c: length of the third - * semi-axis - */ - template - TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> - computeTransverseIsotropicDiluteScheme( - const StressType&, - const real&, - const real&, - const StressType&, - const real&, - const tfel::math::tvector<3u, real>&, - const LengthType&, - const LengthType&, - const LengthType&); - - /*! - * This function gives the homogenized stiffness for a dilute scheme, for - * a distribution of oriented ellipsoids : all principal axes have the - * same fixed orientation. \tparam real: underlying type \tparam - * StressType: type of the elastic constants related to the matrix and the - * ellipsoids \tparam LengthType: type of the dimensions of the ellipsoids - * \return an object of type st2tost2<3u,StressType> - * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix - * \param [in] f: volumic fraction of ellipsoids - * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the - * inclusions \param [in] n_a: direction of the principal axis whose - * semi-length is \f$a\f$ \param[in] a: length of semi-axis relative to - * the direction \f$n_a\f$ \param [in] n_b: direction of the principal - * axis whose semi-length is \f$b\f$ \param[in] b: length of the semi-axis - * relative to the direction \f$n_b\f$ \param[in] c: length of the third - * semi-axis - */ - template - TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> - computeOrientedDiluteScheme(const StressType&, - const real&, - const real&, - const StressType&, - const real&, - const tfel::math::tvector<3u, real>&, - const LengthType&, - const tfel::math::tvector<3u, real>&, - const LengthType&, - const LengthType&); - - /*! - * This function gives the homogenized moduli for a Mori-Tanaka scheme, - * for an isotropic distribution of ellipsoids \tparam real: underlying - * type \tparam StressType: type of the elastic constants related to the - * matrix and the ellipsoids \tparam LengthType: type of the dimensions of - * the ellipsoids \return an object of type std::pair - * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix - * \param [in] f: volumic fraction of ellipsoids - * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the - * inclusions \param[in] a: length of the first semi-axis \param[in] b: - * length of the second semi-axis \param[in] c: length of the third - * semi-axis - */ - template - TFEL_HOST_DEVICE const std::pair - computeIsotropicMoriTanakaScheme(const StressType&, - const real&, - const real&, - const StressType&, - const real&, - const LengthType&, - const LengthType&, - const LengthType&); - - /*! - * This function gives the homogenized stiffness for a Mori-Tanaka scheme, - * for a transverse isotropic distribution of ellipsoids. One principal - * axis (the same principal axis for all ellipsoids) is oriented in a - * fixed direction n. The distribution of the other axes are in the plane - * normal to n, and the distribution of these axes inside this plane is - * isotropic. \tparam real: underlying type \tparam StressType: type of - * the elastic constants related to the matrix and the ellipsoids \tparam - * LengthType: type of the dimensions of the ellipsoids \return an object - * of type st2tost2<3u,StressType> \param [in] young,nu: Young modulus and - * Poisson's ratio of the matrix \param [in] f: volumic fraction of - * ellipsoids \param [in] young_i,nu_i: Young modulus and Poisson's ratio - * of the inclusions \param [in] n_a: direction of the principal axis - * which has a fixed orientation and a semi-length \f$a\f$ \param[in] a: - * length of semi-axis relative to the direction \f$n_a\f$ \param[in] b: - * length of the second semi-axis \param[in] c: length of the third - * semi-axis - */ - template - TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> - computeTransverseIsotropicMoriTanakaScheme( - const StressType&, - const real&, - const real&, - const StressType&, - const real&, - const tfel::math::tvector<3u, real>&, - const LengthType&, - const LengthType&, - const LengthType&); - - /*! - * This function gives the homogenized stiffness for a Mori-Tanaka scheme, - * for a distribution of oriented ellipsoids : all principal axes have the - * same fixed orientation. \tparam real: underlying type \tparam - * StressType: type of the elastic constants related to the matrix and the - * ellipsoids \tparam LengthType: type of the dimensions of the ellipsoids - * \return an object of type st2tost2<3u,StressType> - * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix - * \param [in] f: volumic fraction of ellipsoids - * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the - * inclusions \param [in] n_a: direction of the principal axis whose - * semi-length is \f$a\f$ \param[in] a: length of semi-axis relative to - * the direction \f$n_a\f$ \param [in] n_b: direction of the principal - * axis whose semi-length is \f$b\f$ \param[in] b: length of the semi-axis - * relative to the direction \f$n_b\f$ \param[in] c: length of the third - * semi-axis - */ - template - TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> - computeOrientedMoriTanakaScheme(const StressType&, - const real&, - const real&, - const StressType&, - const real&, - const tfel::math::tvector<3u, real>&, - const LengthType&, - const tfel::math::tvector<3u, real>&, - const LengthType&, - const LengthType&); - - } // end of namespace elasticity - } // end of namespace homogenization - -} // end of namespace tfel::material - -#include "TFEL/Material/EshelbyBasedHomogenization.ixx" - -#endif /* LIB_TFEL_MATERIAL_ESHELBYBASEDHOMOGENIZATION_HXX */ diff --git a/include/TFEL/Material/EshelbyBasedHomogenization.ixx b/include/TFEL/Material/EshelbyBasedHomogenization.ixx deleted file mode 100644 index f1ba883f5..000000000 --- a/include/TFEL/Material/EshelbyBasedHomogenization.ixx +++ /dev/null @@ -1,555 +0,0 @@ -/*! - * \file include/TFEL/Material/EshelbyBasedHomogenization.ixx - * \author Antoine Martin - * \date 24 October 2024 - * \brief This file defines some homogenization schemes based on solution of - * Eshelby problem. \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All - * rights reserved. This project is publicly released under either the GNU GPL - * Licence or the CECILL-A licence. A copy of thoses licences are delivered with - * the sources of TFEL. CEA or EDF may also distribute this project under - * specific licensing conditions. - */ - -#ifndef LIB_TFEL_MATERIAL_ESHELBYBASEDHOMOGENIZATION_IXX -#define LIB_TFEL_MATERIAL_ESHELBYBASEDHOMOGENIZATION_IXX - -#include - -namespace tfel::material { - - namespace homogenization { - namespace elasticity { - - template - struct EllipsoidMeanLocalisator { - static constexpr auto eps = std::numeric_limits::epsilon(); - - TFEL_HOST_DEVICE static const std::pair Isotropic( - const StressType& young, - const real& nu, - const StressType& young_i, - const real& nu_i, - const LengthType& a, - const LengthType& b, - const LengthType& c) { - if ((nu > 0.5) || (nu < -1)) { - tfel::reportContractViolation("nu>0.5 or nu<-1"); - } - if (not(young > StressType{0})) { - tfel::reportContractViolation("E<=0"); - } - if (not((a > LengthType{0}) and (b > LengthType{0}) and - (c > LengthType{0}))) { - tfel::reportContractViolation("a<=0 or b<=0 or c<=0"); - } - const tfel::math::tvector<3u, real> n_1 = {1., 0., 0.}; - const tfel::math::tvector<3u, real> n_2 = {0., 1., 0.}; - real mu; - real ka; - if ((a == b) and (b == c)) { - const auto A = computeSphereLocalisationTensor( - young, nu, young_i, nu_i); - mu = A(3, 3) / 2; - ka = A(0, 0) - 4 * mu / 3; - } else if ((a == b) || (a == c) || (b == c)) { - tfel::math::st2tost2<3u, real> A_; - if (a == b) { - A_ = computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_1, c / a); - } else if (a == c) { - A_ = computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_1, b / a); - } else { - A_ = computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_1, a / b); - } - const auto A1111 = 8 * A_(1, 1) / 15 + A_(0, 0) / 5 + - 2 * (A_(2, 0) + A_(0, 2) + 2 * A_(4, 4)) / 15; - const auto A1122 = A_(1, 1) / 15 + A_(0, 0) / 15 + A_(1, 2) / 3 + - 4 * A_(0, 1) / 15 + 4 * A_(1, 0) / 15 - - 2 * A_(3, 3) / 15; - mu = (A1111 - A1122) / 2; - ka = (A1111 + 2 * A1122) / 3; - } else { - const auto A_ = computeEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_1, a, n_2, b, c); - const auto A1111 = A_(0, 0) / 5 + A_(1, 1) / 5 + A_(2, 2) / 5 + - (A_(0, 1) + A_(1, 0) + 2 * A_(3, 3)) / 15 + - (A_(0, 2) + A_(2, 0) + 2 * A_(4, 4)) / 15 + - (A_(1, 2) + A_(2, 1) + 2 * A_(5, 5)) / 15; - const auto A1122 = - A_(0, 0) / 15 + A_(1, 1) / 15 + A_(2, 2) / 15 + - 2 * A_(0, 1) / 15 - 2 * A_(3, 3) / 30 + 2 * A_(1, 0) / 15 + - 2 * A_(2, 1) / 15 - 2 * A_(5, 5) / 30 + 2 * A_(1, 2) / 15 + - 2 * A_(2, 0) / 15 - 2 * A_(4, 4) / 30 + 2 * A_(0, 2) / 15; - mu = (A1111 - A1122) / 2; - ka = (A1111 + 2 * A1122) / 3; - } - return {ka, mu}; - } // end of Isotropic - - TFEL_HOST_DEVICE static const tfel::math::st2tost2<3u, real> - TransverseIsotropic(const StressType& young, - const real& nu, - const StressType& young_i, - const real& nu_i, - const tfel::math::tvector<3u, real>& n_a, - const LengthType& a, - const LengthType& b, - const LengthType& c) { - if ((nu > 0.5) || (nu < -1)) { - tfel::reportContractViolation("nu>0.5 or nu<-1"); - } - if (not(young > StressType{0})) { - tfel::reportContractViolation("E<=0"); - } - if (not((a > LengthType{0}) and (b > LengthType{0}) and - (c > LengthType{0}))) { - tfel::reportContractViolation("a<=0 or b<=0 or c<=0"); - } - if (tfel::math::ieee754::fpclassify(norm(n_a)) == FP_ZERO) { - tfel::reportContractViolation("n_a is null"); - } - using namespace tfel::math; - tvector<3u, real> n_; - if ((n_a[1] != 0.) || (n_a[2] != 0.)) { - n_ = {1., 0., 0.}; - } else { - n_ = {0., 1., 0.}; - } - const auto n_2 = cross_product(n_a, n_); - const auto n_3 = cross_product(n_a, n_2); - const tvector<3u, real> n_x = {1., 0., 0.}; - const tvector<3u, real> n_y = {0., 1., 0.}; - const tvector<3u, real> n_z = {0., 0., 1.}; - st2tost2<3u, real> A; - if ((a == b) and (b == c)) { - A = computeSphereLocalisationTensor( - young, nu, young_i, nu_i); - } else if (b == c) { - const auto A_ = - computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_z, a / b); - const rotation_matrix r = {n_2[0], n_2[1], n_2[2], - n_3[0], n_3[1], n_3[2], - n_a[0], n_a[1], n_a[2]}; - A = change_basis(A_, r); - } else { - st2tost2<3u, real> A_; - if (a == b) { - A_ = computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_y, c / a); - } else if (a == c) { - A_ = computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_x, b / a); - } else { - A_ = computeEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_z, a, n_x, b, c); - } - const auto A11 = 3 * (A_(0, 0) + A_(1, 1)) / 8 + - (A_(0, 1) + A_(1, 0) + 2 * A_(3, 3)) / 8; - const auto A22 = A11; - const auto A33 = A_(2, 2); - const auto A12 = (A_(0, 0) + A_(1, 1)) / 8 + - 3 * (A_(0, 1) + A_(1, 0)) / 8 - 2 * A_(3, 3) / 8; - const auto A21 = A12; - const auto A13 = (A_(0, 2) + A_(1, 2)) / 2; - const auto A23 = A13; - const auto A31 = (A_(2, 1) + A_(2, 0)) / 2; - const auto A32 = A31; - const auto A44 = - 2 * ((A_(0, 0) + A_(1, 1) - A_(0, 1) - A_(1, 0)) / 8 + - A_(3, 3) / 4); - const auto A55 = 2 * (A_(4, 4) / 4 + A_(5, 5) / 4); - const auto A66 = A55; - const auto zero = real{0}; - const st2tost2<3u, real> A_moy = { - A11, A12, A13, zero, zero, zero, A21, A22, A23, - zero, zero, zero, A31, A32, A33, zero, zero, zero, - zero, zero, zero, A44, zero, zero, zero, zero, zero, - zero, A55, zero, zero, zero, zero, zero, zero, A66}; - const rotation_matrix r = {n_2[0], n_2[1], n_2[2], - n_3[0], n_3[1], n_3[2], - n_a[0], n_a[1], n_a[2]}; - A = change_basis(A_moy, r); - } - return A; - } // end of TransverseIsotropic - - TFEL_HOST_DEVICE static const tfel::math::st2tost2<3u, real> Oriented( - const StressType& young, - const real& nu, - const StressType& young_i, - const real& nu_i, - const tfel::math::tvector<3u, real>& n_a, - const LengthType& a, - const tfel::math::tvector<3u, real>& n_b, - const LengthType& b, - const LengthType& c) { - if ((nu > 0.5) || (nu < -1)) { - tfel::reportContractViolation("nu>0.5 or nu<-1"); - } - if (not(young > StressType{0})) { - tfel::reportContractViolation("E<=0"); - } - if (not((a > LengthType{0}) and (b > LengthType{0}) and - (c > LengthType{0}))) { - tfel::reportContractViolation("a<=0 or b<=0 or c<=0"); - } - if (not(tfel::math::ieee754::fpclassify( - tfel::math::VectorVectorDotProduct::exe< - real, tfel::math::tvector<3u, real>, - tfel::math::tvector<3u, real>>(n_a, n_b)) == - FP_ZERO)) { - tfel::reportContractViolation("n_a and n_b not normals"); - } - if (tfel::math::ieee754::fpclassify(norm(n_a)) == FP_ZERO) { - tfel::reportContractViolation("n_a is null"); - } - if (tfel::math::ieee754::fpclassify(norm(n_b)) == FP_ZERO) { - tfel::reportContractViolation("n_b is null"); - } - using namespace tfel::math; - st2tost2<3u, real> A; - if ((a == b) and (b == c)) { - A = computeSphereLocalisationTensor( - young, nu, young_i, nu_i); - } else if (a == b) { - const tvector<3u, real> n_1 = tfel::math::cross_product(n_a, n_b); - A = computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_1, c / a); - } else if (a == c) { - A = computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_b, b / a); - } else if (b == c) { - A = computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_a, a / b); - } else { - A = computeEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_a, a, n_b, b, c); - } - return A; - } // end of Oriented - - }; // end of struct EllipsoidMeanLocalisator ; - - template - TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> - computeDiluteScheme(const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i, - const tfel::math::st2tost2<3u, real>& A) { - if ((nu > 0.5) || (nu < -1)) { - tfel::reportContractViolation("nu>0.5 or nu<-1"); - } - if (not(young > StressType{0})) { - tfel::reportContractViolation("E<=0"); - } - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - if ((nu_i > 0.5) || (nu_i < -1)) { - tfel::reportContractViolation("nu_i>0.5 or nu_i<-1"); - } - if (not(young_i > StressType{0})) { - tfel::reportContractViolation("E_i<=0"); - } - using namespace tfel::math; - st2tost2<3u, StressType> C_0; - static constexpr auto value = - StiffnessTensorAlterationCharacteristic::UNALTERED; - computeIsotropicStiffnessTensorII<3u, value, StressType, real>( - C_0, young, nu); - st2tost2<3u, StressType> C_i; - computeIsotropicStiffnessTensorII<3u, value, StressType, real>( - C_i, young_i, nu_i); - const auto C = C_i - C_0; - const auto Pr = C * A; - return C_0 + f * Pr; - } // end of DiluteScheme - - template - TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> - computeMoriTanakaScheme(const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i, - const tfel::math::st2tost2<3u, real>& A) { - if ((nu > 0.5) || (nu < -1)) { - tfel::reportContractViolation("nu>0.5 or nu<-1"); - } - if (not(young > StressType{0})) { - tfel::reportContractViolation("E<=0"); - } - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - if ((nu_i > 0.5) || (nu_i < -1)) { - tfel::reportContractViolation("nu_i>0.5 or nu_i<-1"); - } - if (not(young_i > StressType{0})) { - tfel::reportContractViolation("E_i<=0"); - } - using namespace tfel::math; - st2tost2<3u, StressType> C_0; - static constexpr auto value = - StiffnessTensorAlterationCharacteristic::UNALTERED; - computeIsotropicStiffnessTensorII<3u, value, StressType, real>( - C_0, young, nu); - st2tost2<3u, StressType> C_i; - computeIsotropicStiffnessTensorII<3u, value, StressType, real>( - C_i, young_i, nu_i); - const auto C = C_i - C_0; - const auto Pr = C * A; - const auto inv = invert(f * A + (1 - f) * st2tost2<3u, real>::Id()); - const auto PPr = Pr * inv; - return C_0 + f * PPr; - } // end of MoriTanakaScheme - - template - TFEL_HOST_DEVICE const std::pair - computeSphereDiluteScheme(const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i) { - if ((nu > 0.5) || (nu < -1)) { - tfel::reportContractViolation("nu>0.5 or nu<-1"); - } - if (not(young > StressType{0})) { - tfel::reportContractViolation("E<=0"); - } - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - - const auto k0 = young / 3 / (1 - 2 * nu); - const auto mu0 = young / 2 / (1 + nu); - const auto k_i = young_i / 3 / (1 - 2 * nu_i); - const auto mu_i = young_i / 2 / (1 + nu_i); - const auto alpha0 = 3 * k0 / (3 * k0 + 4 * mu0); - const auto beta0 = 6 * (k0 + 2 * mu0) / 5 / (3 * k0 + 4 * mu0); - const auto muhom = - mu0 + f * (mu_i - mu0) / (1 + beta0 * (mu_i - mu0) / mu0); - const auto khom = k0 + f * (k_i - k0) / (1 + alpha0 * (k_i - k0) / k0); - const auto nuhom = (3 * khom - 2 * muhom) / (2 * muhom + 6 * khom); - const auto Ehom = 2 * muhom * (1 + nuhom); - return {Ehom, nuhom}; - }; - - template - TFEL_HOST_DEVICE const std::pair - computeSphereMoriTanakaScheme(const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i) { - if ((nu > 0.5) || (nu < -1)) { - tfel::reportContractViolation("nu>0.5 or nu<-1"); - } - if (not(young > StressType{0})) { - tfel::reportContractViolation("E<=0"); - } - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - const auto k0 = young / 3 / (1 - 2 * nu); - const auto mu0 = young / 2 / (1 + nu); - const auto k_i = young_i / 3 / (1 - 2 * nu_i); - const auto mu_i = young_i / 2 / (1 + nu_i); - const auto muhom = - mu0 + - f * (mu_i - mu0) / - (1 + (1 - f) * (mu_i - mu0) / - (mu0 + mu0 * (9 * k0 + 8 * mu0) / 6 / (k0 + 2 * mu0))); - const auto khom = - k0 + - f * (k_i - k0) / (1 + (1 - f) * (k_i - k0) / (k0 + 4 * mu0 / 3)); - const auto nuhom = (3 * khom - 2 * muhom) / (2 * muhom + 6 * khom); - const auto Ehom = 2 * muhom * (1 + nuhom); - return {Ehom, nuhom}; - }; - - template - const std::pair computeIsotropicDiluteScheme( - const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i, - const LengthType& a, - const LengthType& b, - const LengthType& c) { - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - - const auto pair = - EllipsoidMeanLocalisator<3u, real, StressType, - LengthType>::Isotropic(young, nu, young_i, - nu_i, a, b, c); - const auto ka = std::get<0>(pair); - const auto mu = std::get<1>(pair); - const auto k0 = young / 3 / (1 - 2 * nu); - const auto mu0 = young / 2 / (1 + nu); - const auto k_i = young_i / 3 / (1 - 2 * nu_i); - const auto mu_i = young_i / 2 / (1 + nu_i); - const auto muhom = mu0 + 2 * f * (mu_i - mu0) * mu; - const auto khom = k0 + 3 * f * (k_i - k0) * ka; - const auto nuhom = (3 * khom - 2 * muhom) / (2 * muhom + 6 * khom); - const auto Ehom = 2 * muhom * (1 + nuhom); - return {Ehom, nuhom}; - }; - - template - const tfel::math::st2tost2<3u, StressType> - computeTransverseIsotropicDiluteScheme( - const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i, - const tfel::math::tvector<3u, real>& n_a, - const LengthType& a, - const LengthType& b, - const LengthType& c) { - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - const auto A = - EllipsoidMeanLocalisator<3u, real, StressType, - LengthType>::TransverseIsotropic(young, nu, - young_i, - nu_i, n_a, - a, b, c); - return computeDiluteScheme(young, nu, f, young_i, - nu_i, A); - }; - - template - const tfel::math::st2tost2<3u, StressType> computeOrientedDiluteScheme( - const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i, - const tfel::math::tvector<3u, real>& n_a, - const LengthType& a, - const tfel::math::tvector<3u, real>& n_b, - const LengthType& b, - const LengthType& c) { - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - const auto A = - EllipsoidMeanLocalisator<3u, real, StressType, - LengthType>::Oriented(young, nu, young_i, - nu_i, n_a, a, n_b, b, - c); - return computeDiluteScheme(young, nu, f, young_i, - nu_i, A); - }; - - template - const std::pair computeIsotropicMoriTanakaScheme( - const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i, - const LengthType& a, - const LengthType& b, - const LengthType& c) { - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - const auto pair = - EllipsoidMeanLocalisator<3u, real, StressType, - LengthType>::Isotropic(young, nu, young_i, - nu_i, a, b, c); - const auto ka = std::get<0>(pair); - const auto mu = std::get<1>(pair); - const auto k0 = young / 3 / (1 - 2 * nu); - const auto mu0 = young / 2 / (1 + nu); - const auto k_i = young_i / 3 / (1 - 2 * nu_i); - const auto mu_i = young_i / 2 / (1 + nu_i); - const auto muhom = mu0 + f * (mu_i - mu0) * mu / (f * mu + (1 - f) / 2); - const auto khom = k0 + f * (k_i - k0) * ka / (ka * f + (1 - f) / 3); - const auto nuhom = (3 * khom - 2 * muhom) / (2 * muhom + 6 * khom); - const auto Ehom = 2 * muhom * (1 + nuhom); - return {Ehom, nuhom}; - }; - - template - const tfel::math::st2tost2<3u, StressType> - computeTransverseIsotropicMoriTanakaScheme( - const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i, - const tfel::math::tvector<3u, real>& n_a, - const LengthType& a, - const LengthType& b, - const LengthType& c) { - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - const auto A = - EllipsoidMeanLocalisator<3u, real, StressType, - LengthType>::TransverseIsotropic(young, nu, - young_i, - nu_i, n_a, - a, b, c); - return computeMoriTanakaScheme(young, nu, f, young_i, - nu_i, A); - }; - - template - const tfel::math::st2tost2<3u, StressType> - computeOrientedMoriTanakaScheme(const StressType& young, - const real& nu, - const real& f, - const StressType& young_i, - const real& nu_i, - const tfel::math::tvector<3u, real>& n_a, - const LengthType& a, - const tfel::math::tvector<3u, real>& n_b, - const LengthType& b, - const LengthType& c) { - if ((f < 0) || (f > 1)) { - tfel::reportContractViolation("f<0 or f>1"); - } - const auto A = - EllipsoidMeanLocalisator<3u, real, StressType, - LengthType>::Oriented(young, nu, young_i, - nu_i, n_a, a, n_b, b, - c); - return computeMoriTanakaScheme(young, nu, f, young_i, - nu_i, A); - }; - - } // end of namespace elasticity - } // end of namespace homogenization - -} // end of namespace tfel::material - -#endif /* LIB_TFEL_MATERIAL_ESHELBYBASEDHOMOGENIZATION_IXX */ diff --git a/include/TFEL/Material/IsotropicEshelbyTensor.hxx b/include/TFEL/Material/IsotropicEshelbyTensor.hxx new file mode 100644 index 000000000..7ae8107db --- /dev/null +++ b/include/TFEL/Material/IsotropicEshelbyTensor.hxx @@ -0,0 +1,184 @@ +/*! + * \file include/TFEL/Material/IsotropicEshelbyTensor.hxx + * \author Antoine Martin + * \date 15 October 2024 + * \brief This file declares the Eshelby tensor for an ellipsoidal inclusion + * embedded in an isotropic matrix. \copyright Copyright (C) 2006-2018 CEA/DEN, + * EDF R&D. All rights reserved. This project is publicly released under either + * the GNU GPL Licence or the CECILL-A licence. A copy of thoses licences are + * delivered with the sources of TFEL. CEA or EDF may also distribute this + * project under specific licensing conditions. + */ + +#ifndef LIB_TFEL_MATERIAL_ISOTROPICESHELBYTENSOR_HXX +#define LIB_TFEL_MATERIAL_ISOTROPICESHELBYTENSOR_HXX + +#include "TFEL/Math/st2tost2.hxx" +#include "TFEL/Material/StiffnessTensor.hxx" +#include + +namespace tfel::material::homogenization::elasticity { + + /*! + * This function builds the Eshelby tensor of a circular cylinder embedded + * in an isotropic matrix, considering a PLANE STRAIN modelling hypothesis + * \return an object of type st2tost2<2u,real> + * \tparam real: underlying type + * \param[in] nu: Poisson's ratio of the matrix + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<2u, real> + computeCircularCylinderEshelbyTensor(const real&); + + /*! + * This function builds the Eshelby tensor of an elliptic cylinder + * embedded in an isotropic matrix, considering a PLANE STRAIN modelling + * hypothesis The function returns the Eshelby tensor in the basis (e1,e2) + * where e1 corresponds to the biggest axis \return an object of type + * st2tost2<2u,real> \tparam real: underlying type \param[in] nu: + * Poisson's ratio of the matrix \param[in] e: aspect ratio of the + * elliptic basis + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<2u, real> + computeEllipticCylinderEshelbyTensor(const real&, const real&); + + /*! + * This function builds the Eshelby tensor of a sphere embedded in an + * isotropic matrix. \return an object of type st2tost2<3u,real> \tparam + * real: underlying type \param[in] nu: Poisson's ratio of the matrix + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> computeSphereEshelbyTensor( + const real&); + + /*! + * This function builds the Eshelby tensor of an axisymmetrical ellipsoid + * embedded in an isotropic matrix. The function returns the Eshelby + * tensor in the basis (e1,e2,e3) where e1 corresponds to (one of) the + * biggest ax(es) \return an object of type st2tost2<3u,real> \tparam + * real: underlying type \param[in] nu: Poisson's ratio of the matrix + * \param[in] e: aspect ratio of the ellipsoid (e>1 : prolate, e<1 : + * oblate) \param[in] precf, precd, precld: default arguments which aim at + * preventing the numerical instability of the formula when the ellipsoid + * is almost a sphere. When the absolute value of (e-1) is below precf + * (resp. precd, precld) for real=float (resp. real= double, long double), + * the returned tensor is computeSphereEshelbyTensor(nu). + * + * The expressions can be found in Torquato, Random Heterogeneous + * Materials (2002). + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> + computeAxisymmetricalEshelbyTensor(const real&, + const real&, + const real = real{8e-3}, + const real = real{1.5e-4}, + const real = real{1e-5}); + + /*! + * This function builds the Eshelby tensor of a general ellipsoid embedded + * in an isotropic matrix. The function returns the Eshelby tensor in the + * basis (e1,e2,e3) where e1 (resp. e2, e3) is aligned with the axis with + * the first (resp. the second and third) biggest length \return an object + * of type st2tost2<3u,real> \tparam real: underlying type \tparam + * LengthType: type of the dimensions of the ellipsoid \param[in] nu: + * Poisson's ratio of the matrix \param[in] a: length of the first + * semi-axis \param[in] b: length of the second semi-axis \param[in] c: + * length of the third semi-axis \param[in] precf, precd, precld: default + * arguments which aim at preventing the numerical instability of the + * formula when the ellipsoid is almost axisymmetrical. When the absolute + * value of (a-b)/c (or (a-c)/b or (b-c)/a) is below precf (resp. precd, + * precld) for real=float (resp. real= double, long double), the returned + * tensor is computeAxisymmetricalEshelbyTensor. + * + * The expressions for the case of three distinct semi-axes can be found + * in Eshelby (1957). + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> computeEshelbyTensor( + const real&, + const LengthType&, + const LengthType&, + const LengthType&, + const real = real{8e-3}, + const real = real{1.5e-4}, + const real = real{1e-5}); + + /*! + * This function builds the strain localisation tensor of a general + * ellipsoid with a general elasticity, embedded in an isotropic matrix. + * The localisation tensor \f$A\f$ is defined as follows : \f[\epsilon = + * A:E_0\f] where \f$E_0\f$ is the uniform strain tensor imposed at + * infinity, and \f$\epsilon\f$ is the strain tensor solution of Eshelby + * problem for the ellipsoid. The ellipsoid also has a specific + * orientation given by the vectors \f$n_a\f$, \f$n_b\f$. \return an + * object of type st2tost2, which is the fourth-order localisation tensor + * \f$A\f$ \tparam real: underlying type \tparam StressType: type of the + * elastic constants related to the matrix and the ellipsoid \tparam + * LengthType: type of the dimensions of the ellipsoid \param [in] + * young,nu: Young modulus and Poisson's ratio of the matrix \param [in] + * young_i,nu_i: Young modulus and Poisson's ratio of the inclusions + * \param [in] n_a: direction of the principal axis whose length is + * \f$a\f$ \param [in] a: length of semi-axis relative to the direction + * \f$n_a\f$ \param [in] n_b: direction of the principal axis whose length + * is \f$b\f$ \param [in] b: length of semi-axis relative to the direction + * \f$n_b\f$ \param [in] c: length of the remaining semi-axis + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> + computeEllipsoidLocalisationTensor(const StressType&, + const real&, + const StressType&, + const real&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const LengthType&); + + /*! + * This function builds the strain localisation tensor of an + * axisymmetrical ellipsoid with a general elasticity, embedded in an + * isotropic matrix. The ellipsoid also has a specific orientation given + * by the vector \f$n_a\f$, axis of the ellipsoid, whose semi-length is + * \f$a\f$. \return an object of type st2tost2 \tparam real: underlying + * type \tparam StressType: type of the elastic constants related to the + * matrix and the ellipsoid \param [in] young,nu: Young modulus and + * Poisson's ratio of the matrix \param [in] young_i,nu_i: Young modulus + * and Poisson's ratio of the inclusions \param [in] n_a: direction of the + * axis of the ellipsoid (whose semi-length is \f$a\f$) \param [in] a: + * length of semi-axis relative to the direction \f$n_a\f$ + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> + computeAxisymmetricalEllipsoidLocalisationTensor( + const StressType&, + const real&, + const StressType&, + const real&, + const tfel::math::tvector<3u, real>&, + const real&); + + /*! + * This function builds the strain localisation tensor of a sphere + * with a general elasticity, embedded in an isotropic matrix. + * \return an object of type st2tost2 + * \tparam real: underlying type + * \tparam StressType: type of the elastic constants related to the matrix + * and the ellipsoid \param [in] young,nu: Young modulus and Poisson's + * ratio of the matrix \param [in] young_i,nu_i: Young modulus and + * Poisson's ratio of the inclusions + */ + template + TFEL_HOST_DEVICE tfel::math::st2tost2<3u, real> + computeSphereLocalisationTensor(const StressType&, + const real&, + const StressType&, + const real&); + +} // end of namespace tfel::material::homogenization::elasticity + +#include "TFEL/Material/IsotropicEshelbyTensor.ixx" + +#endif /* LIB_TFEL_MATERIAL_ESHELBY_HXX */ diff --git a/include/TFEL/Material/Eshelby.ixx b/include/TFEL/Material/IsotropicEshelbyTensor.ixx similarity index 96% rename from include/TFEL/Material/Eshelby.ixx rename to include/TFEL/Material/IsotropicEshelbyTensor.ixx index ae8ba7e84..4beeca607 100644 --- a/include/TFEL/Material/Eshelby.ixx +++ b/include/TFEL/Material/IsotropicEshelbyTensor.ixx @@ -1,5 +1,5 @@ /*! - * \file include/TFEL/Material/Eshelby.ixx + * \file include/TFEL/Material/IsotropicEshelbyTensor.ixx * \author Antoine Martin * \date 15 October 2024 * \brief This file defines the Eshelby tensor for an ellipsoidal inclusion @@ -10,8 +10,8 @@ * project under specific licensing conditions. */ -#ifndef LIB_TFEL_MATERIAL_ESHELBY_IXX -#define LIB_TFEL_MATERIAL_ESHELBY_IXX +#ifndef LIB_TFEL_MATERIAL_ISOTROPICESHELBYTENSOR_IXX +#define LIB_TFEL_MATERIAL_ISOTROPICESHELBYTENSOR_IXX #include #include @@ -101,11 +101,14 @@ namespace tfel::material::homogenization::elasticity { tfel::reportContractViolation("e<=0"); } static constexpr auto eps = std::numeric_limits::epsilon(); + const auto precision = [precf, precd, precld]() { - if (eps == std::numeric_limits>::epsilon()) { - return precf; + if (typeid(real) == typeid(long double)) { + return precld; + } else if (typeid(real) == typeid(double)) { + return precd; } - return precld; + return precf; }(); if (std::abs(e - 1) < precision) { return computeSphereEshelbyTensor(nu); @@ -201,11 +204,14 @@ namespace tfel::material::homogenization::elasticity { tfel::reportContractViolation("a<=0 or b<=0 or c<=0"); } static constexpr auto eps = std::numeric_limits::epsilon(); + const auto precision = [precf, precd, precld]() { - if (eps == std::numeric_limits::epsilon()) { - return precf; + if (typeid(real) == typeid(long double)) { + return precld; + } else if (typeid(real) == typeid(double)) { + return precd; } - return precld; + return precf; }(); if (std::abs((b - a) / c) < precision || std::abs((a - c) / b) < precision || @@ -424,4 +430,4 @@ namespace tfel::material::homogenization::elasticity { } // end of namespace tfel::material::homogenization::elasticity -#endif /* LIB_TFEL_MATERIAL_ESHELBY_IXX */ +#endif /* LIB_TFEL_MATERIAL_ISOTROPICESHELBYTENSOR_IXX */ diff --git a/include/TFEL/Material/LinearHomogenizationSchemes.hxx b/include/TFEL/Material/LinearHomogenizationSchemes.hxx new file mode 100644 index 000000000..881e4a50e --- /dev/null +++ b/include/TFEL/Material/LinearHomogenizationSchemes.hxx @@ -0,0 +1,265 @@ +/*! + * \file include/TFEL/Material/LinearHomogenizationSchemes.hxx + * \author Antoine Martin + * \date 24 October 2024 + * \brief This file declares some homogenization schemes based on solution of + * Eshelby problem. \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All + * rights reserved. This project is publicly released under either the GNU GPL + * Licence or the CECILL-A licence. A copy of thoses licences are delivered with + * the sources of TFEL. CEA or EDF may also distribute this project under + * specific licensing conditions. + */ + +#ifndef LIB_TFEL_MATERIAL_LINEARHOMOGENIZATIONSCHEMES_HXX +#define LIB_TFEL_MATERIAL_LINEARHOMOGENIZATIONSCHEMES_HXX + +#include "TFEL/Math/st2tost2.hxx" +#include "TFEL/Material/IsotropicEshelbyTensor.hxx" + +namespace tfel::material::homogenization::elasticity { + + /*! + * This function gives the homogenized stiffness for a dilute scheme, + * knowing the strain localisation tensor. \tparam real: underlying type + * \tparam StressType: type of the elastic constants related to the matrix + * and the inclusions \return an object of type st2tost2<3u,StressType> + * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix + * \param [in] f: volumic fraction of inclusions + * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the + * inclusions \param [in] A: mean strain localisation tensor of inclusions + */ + template + TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> + computeDiluteScheme(const StressType&, + const real&, + const real&, + const StressType&, + const real&, + const tfel::math::st2tost2<3u, real>&); + + /*! + * This function gives the homogenized stiffness for a Mori-Tanaka scheme + * knowing the strain localisation tensor. \tparam real: underlying type + * \tparam StressType: type of the elastic constants related to the matrix + * and the inclusions \return an object of type st2tost2<3u,StressType> + * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix + * \param [in] f: volumic fraction of inclusions + * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the + * inclusions \param [in] A: mean strain localisation tensor of inclusions + */ + template + TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> + computeMoriTanakaScheme(const StressType&, + const real&, + const real&, + const StressType&, + const real&, + const tfel::math::st2tost2<3u, real>&); + + /*! + * This function gives the homogenized moduli for a dilute scheme, for + * spheres \tparam real: underlying type \tparam StressType: type of the + * elastic constants related to the matrix and the spheres \return an + * object of type std::pair \param [in] young,nu: Young + * modulus and Poisson's ratio of the matrix \param [in] f: volumic + * fraction of spheres \param [in] young_i,nu_i: Young modulus and + * Poisson's ratio of the inclusions + */ + template + TFEL_HOST_DEVICE const std::pair computeSphereDiluteScheme( + const StressType&, + const real&, + const real&, + const StressType&, + const real&); + + /*! + * This function gives the homogenized moduli for a Mori-Tanaka scheme, + * for spheres \tparam real: underlying type \tparam StressType: type of + * the elastic constants related to the matrix and the spheres \return an + * object of type std::pair \param [in] young,nu: Young + * modulus and Poisson's ratio of the matrix \param [in] f: volumic + * fraction of spheres \param [in] young_i,nu_i: Young modulus and + * Poisson's ratio of the inclusions + */ + template + TFEL_HOST_DEVICE const std::pair + computeSphereMoriTanakaScheme(const StressType&, + const real&, + const real&, + const StressType&, + const real&); + + /*! + * This function gives the homogenized moduli for a dilute scheme, for an + * isotropic distribution of ellipsoids \tparam real: underlying type + * \tparam StressType: type of the elastic constants related to the matrix + * and the ellipsoids \tparam LengthType: type of the dimensions of the + * ellipsoids \return an object of type std::pair \param + * [in] young,nu: Young modulus and Poisson's ratio of the matrix \param + * [in] f: volumic fraction of ellipsoids \param [in] young_i,nu_i: Young + * modulus and Poisson's ratio of the inclusions \param[in] a: length of + * the first semi-axis \param[in] b: length of the second semi-axis + * \param[in] c: length of the third semi-axis + */ + template + TFEL_HOST_DEVICE const std::pair + computeIsotropicDiluteScheme(const StressType&, + const real&, + const real&, + const StressType&, + const real&, + const LengthType&, + const LengthType&, + const LengthType&); + + /*! + * This function gives the homogenized stiffness for a dilute scheme, for + * a transverse isotropic distribution of ellipsoids. One principal axis + * (the same principal axis for all ellipsoids) is oriented in a fixed + * direction n. The distribution of the other axes are in the plane normal + * to n, and the distribution of these axes inside this plane is + * isotropic. \tparam real: underlying type \tparam StressType: type of + * the elastic constants related to the matrix and the ellipsoids \tparam + * LengthType: type of the dimensions of the ellipsoids \return an object + * of type st2tost2<3u,StressType> \param [in] young,nu: Young modulus and + * Poisson's ratio of the matrix \param [in] f: volumic fraction of + * ellipsoids \param [in] young_i,nu_i: Young modulus and Poisson's ratio + * of the inclusions \param [in] n_a: direction of the principal axis + * which has a fixed orientation and a semi-length \f$a\f$ \param[in] a: + * length of semi-axis relative to the direction \f$n_a\f$ \param[in] b: + * length of the second semi-axis \param[in] c: length of the third + * semi-axis + */ + template + TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> + computeTransverseIsotropicDiluteScheme(const StressType&, + const real&, + const real&, + const StressType&, + const real&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const LengthType&, + const LengthType&); + + /*! + * This function gives the homogenized stiffness for a dilute scheme, for + * a distribution of oriented ellipsoids : all principal axes have the + * same fixed orientation. \tparam real: underlying type \tparam + * StressType: type of the elastic constants related to the matrix and the + * ellipsoids \tparam LengthType: type of the dimensions of the ellipsoids + * \return an object of type st2tost2<3u,StressType> + * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix + * \param [in] f: volumic fraction of ellipsoids + * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the + * inclusions \param [in] n_a: direction of the principal axis whose + * semi-length is \f$a\f$ \param[in] a: length of semi-axis relative to + * the direction \f$n_a\f$ \param [in] n_b: direction of the principal + * axis whose semi-length is \f$b\f$ \param[in] b: length of the semi-axis + * relative to the direction \f$n_b\f$ \param[in] c: length of the third + * semi-axis + */ + template + TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> + computeOrientedDiluteScheme(const StressType&, + const real&, + const real&, + const StressType&, + const real&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const LengthType&); + + /*! + * This function gives the homogenized moduli for a Mori-Tanaka scheme, + * for an isotropic distribution of ellipsoids \tparam real: underlying + * type \tparam StressType: type of the elastic constants related to the + * matrix and the ellipsoids \tparam LengthType: type of the dimensions of + * the ellipsoids \return an object of type std::pair + * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix + * \param [in] f: volumic fraction of ellipsoids + * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the + * inclusions \param[in] a: length of the first semi-axis \param[in] b: + * length of the second semi-axis \param[in] c: length of the third + * semi-axis + */ + template + TFEL_HOST_DEVICE const std::pair + computeIsotropicMoriTanakaScheme(const StressType&, + const real&, + const real&, + const StressType&, + const real&, + const LengthType&, + const LengthType&, + const LengthType&); + + /*! + * This function gives the homogenized stiffness for a Mori-Tanaka scheme, + * for a transverse isotropic distribution of ellipsoids. One principal + * axis (the same principal axis for all ellipsoids) is oriented in a + * fixed direction n. The distribution of the other axes are in the plane + * normal to n, and the distribution of these axes inside this plane is + * isotropic. \tparam real: underlying type \tparam StressType: type of + * the elastic constants related to the matrix and the ellipsoids \tparam + * LengthType: type of the dimensions of the ellipsoids \return an object + * of type st2tost2<3u,StressType> \param [in] young,nu: Young modulus and + * Poisson's ratio of the matrix \param [in] f: volumic fraction of + * ellipsoids \param [in] young_i,nu_i: Young modulus and Poisson's ratio + * of the inclusions \param [in] n_a: direction of the principal axis + * which has a fixed orientation and a semi-length \f$a\f$ \param[in] a: + * length of semi-axis relative to the direction \f$n_a\f$ \param[in] b: + * length of the second semi-axis \param[in] c: length of the third + * semi-axis + */ + template + TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> + computeTransverseIsotropicMoriTanakaScheme( + const StressType&, + const real&, + const real&, + const StressType&, + const real&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const LengthType&, + const LengthType&); + + /*! + * This function gives the homogenized stiffness for a Mori-Tanaka scheme, + * for a distribution of oriented ellipsoids : all principal axes have the + * same fixed orientation. \tparam real: underlying type \tparam + * StressType: type of the elastic constants related to the matrix and the + * ellipsoids \tparam LengthType: type of the dimensions of the ellipsoids + * \return an object of type st2tost2<3u,StressType> + * \param [in] young,nu: Young modulus and Poisson's ratio of the matrix + * \param [in] f: volumic fraction of ellipsoids + * \param [in] young_i,nu_i: Young modulus and Poisson's ratio of the + * inclusions \param [in] n_a: direction of the principal axis whose + * semi-length is \f$a\f$ \param[in] a: length of semi-axis relative to + * the direction \f$n_a\f$ \param [in] n_b: direction of the principal + * axis whose semi-length is \f$b\f$ \param[in] b: length of the semi-axis + * relative to the direction \f$n_b\f$ \param[in] c: length of the third + * semi-axis + */ + template + TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> + computeOrientedMoriTanakaScheme(const StressType&, + const real&, + const real&, + const StressType&, + const real&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const tfel::math::tvector<3u, real>&, + const LengthType&, + const LengthType&); + +} // end of namespace tfel::material::homogenization::elasticity + +#include "TFEL/Material/LinearHomogenizationSchemes.ixx" + +#endif /* LIB_TFEL_MATERIAL_LINEARHOMOGENIZATIONSCHEMES_HXX */ diff --git a/include/TFEL/Material/LinearHomogenizationSchemes.ixx b/include/TFEL/Material/LinearHomogenizationSchemes.ixx new file mode 100644 index 000000000..9c23cf93b --- /dev/null +++ b/include/TFEL/Material/LinearHomogenizationSchemes.ixx @@ -0,0 +1,530 @@ +/*! + * \file include/TFEL/Material/LinearHomogenizationSchemes.ixx + * \author Antoine Martin + * \date 24 October 2024 + * \brief This file defines some homogenization schemes based on solution of + * Eshelby problem. \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All + * rights reserved. This project is publicly released under either the GNU GPL + * Licence or the CECILL-A licence. A copy of thoses licences are delivered with + * the sources of TFEL. CEA or EDF may also distribute this project under + * specific licensing conditions. + */ + +#ifndef LIB_TFEL_MATERIAL_LINEARHOMOGENIZATIONSCHEMES_IXX +#define LIB_TFEL_MATERIAL_LINEARHOMOGENIZATIONSCHEMES_IXX + +#include + +namespace tfel::material::homogenization::elasticity { + + template + struct EllipsoidMeanLocalisator { + static constexpr auto eps = std::numeric_limits::epsilon(); + + TFEL_HOST_DEVICE static const std::pair Isotropic( + const StressType& young, + const real& nu, + const StressType& young_i, + const real& nu_i, + const LengthType& a, + const LengthType& b, + const LengthType& c) { + if ((nu > 0.5) || (nu < -1)) { + tfel::reportContractViolation("nu>0.5 or nu<-1"); + } + if (not(young > StressType{0})) { + tfel::reportContractViolation("E<=0"); + } + if (not((a > LengthType{0}) and (b > LengthType{0}) and + (c > LengthType{0}))) { + tfel::reportContractViolation("a<=0 or b<=0 or c<=0"); + } + const tfel::math::tvector<3u, real> n_1 = {1., 0., 0.}; + const tfel::math::tvector<3u, real> n_2 = {0., 1., 0.}; + real mu; + real ka; + if ((a == b) and (b == c)) { + const auto A = computeSphereLocalisationTensor( + young, nu, young_i, nu_i); + mu = A(3, 3) / 2; + ka = A(0, 0) - 4 * mu / 3; + } else if ((a == b) || (a == c) || (b == c)) { + tfel::math::st2tost2<3u, real> A_; + if (a == b) { + A_ = computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_1, c / a); + } else if (a == c) { + A_ = computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_1, b / a); + } else { + A_ = computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_1, a / b); + } + const auto A1111 = 8 * A_(1, 1) / 15 + A_(0, 0) / 5 + + 2 * (A_(2, 0) + A_(0, 2) + 2 * A_(4, 4)) / 15; + const auto A1122 = A_(1, 1) / 15 + A_(0, 0) / 15 + A_(1, 2) / 3 + + 4 * A_(0, 1) / 15 + 4 * A_(1, 0) / 15 - + 2 * A_(3, 3) / 15; + mu = (A1111 - A1122) / 2; + ka = (A1111 + 2 * A1122) / 3; + } else { + const auto A_ = + computeEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_1, a, n_2, b, c); + const auto A1111 = A_(0, 0) / 5 + A_(1, 1) / 5 + A_(2, 2) / 5 + + (A_(0, 1) + A_(1, 0) + 2 * A_(3, 3)) / 15 + + (A_(0, 2) + A_(2, 0) + 2 * A_(4, 4)) / 15 + + (A_(1, 2) + A_(2, 1) + 2 * A_(5, 5)) / 15; + const auto A1122 = + A_(0, 0) / 15 + A_(1, 1) / 15 + A_(2, 2) / 15 + 2 * A_(0, 1) / 15 - + 2 * A_(3, 3) / 30 + 2 * A_(1, 0) / 15 + 2 * A_(2, 1) / 15 - + 2 * A_(5, 5) / 30 + 2 * A_(1, 2) / 15 + 2 * A_(2, 0) / 15 - + 2 * A_(4, 4) / 30 + 2 * A_(0, 2) / 15; + mu = (A1111 - A1122) / 2; + ka = (A1111 + 2 * A1122) / 3; + } + return {ka, mu}; + } // end of Isotropic + + TFEL_HOST_DEVICE static const tfel::math::st2tost2<3u, real> + TransverseIsotropic(const StressType& young, + const real& nu, + const StressType& young_i, + const real& nu_i, + const tfel::math::tvector<3u, real>& n_a, + const LengthType& a, + const LengthType& b, + const LengthType& c) { + if ((nu > 0.5) || (nu < -1)) { + tfel::reportContractViolation("nu>0.5 or nu<-1"); + } + if (not(young > StressType{0})) { + tfel::reportContractViolation("E<=0"); + } + if (not((a > LengthType{0}) and (b > LengthType{0}) and + (c > LengthType{0}))) { + tfel::reportContractViolation("a<=0 or b<=0 or c<=0"); + } + if (tfel::math::ieee754::fpclassify(norm(n_a)) == FP_ZERO) { + tfel::reportContractViolation("n_a is null"); + } + using namespace tfel::math; + tvector<3u, real> n_; + if ((n_a[1] != 0.) || (n_a[2] != 0.)) { + n_ = {1., 0., 0.}; + } else { + n_ = {0., 1., 0.}; + } + const auto n_2 = cross_product(n_a, n_); + const auto n_3 = cross_product(n_a, n_2); + const tvector<3u, real> n_x = {1., 0., 0.}; + const tvector<3u, real> n_y = {0., 1., 0.}; + const tvector<3u, real> n_z = {0., 0., 1.}; + st2tost2<3u, real> A; + if ((a == b) and (b == c)) { + A = computeSphereLocalisationTensor(young, nu, + young_i, nu_i); + } else if (b == c) { + const auto A_ = + computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_z, a / b); + const rotation_matrix r = {n_2[0], n_2[1], n_2[2], n_3[0], n_3[1], + n_3[2], n_a[0], n_a[1], n_a[2]}; + A = change_basis(A_, r); + } else { + st2tost2<3u, real> A_; + if (a == b) { + A_ = computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_y, c / a); + } else if (a == c) { + A_ = computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_x, b / a); + } else { + A_ = computeEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_z, a, n_x, b, c); + } + const auto A11 = 3 * (A_(0, 0) + A_(1, 1)) / 8 + + (A_(0, 1) + A_(1, 0) + 2 * A_(3, 3)) / 8; + const auto A22 = A11; + const auto A33 = A_(2, 2); + const auto A12 = (A_(0, 0) + A_(1, 1)) / 8 + + 3 * (A_(0, 1) + A_(1, 0)) / 8 - 2 * A_(3, 3) / 8; + const auto A21 = A12; + const auto A13 = (A_(0, 2) + A_(1, 2)) / 2; + const auto A23 = A13; + const auto A31 = (A_(2, 1) + A_(2, 0)) / 2; + const auto A32 = A31; + const auto A44 = 2 * ((A_(0, 0) + A_(1, 1) - A_(0, 1) - A_(1, 0)) / 8 + + A_(3, 3) / 4); + const auto A55 = 2 * (A_(4, 4) / 4 + A_(5, 5) / 4); + const auto A66 = A55; + const auto zero = real{0}; + const st2tost2<3u, real> A_moy = { + A11, A12, A13, zero, zero, zero, A21, A22, A23, + zero, zero, zero, A31, A32, A33, zero, zero, zero, + zero, zero, zero, A44, zero, zero, zero, zero, zero, + zero, A55, zero, zero, zero, zero, zero, zero, A66}; + const rotation_matrix r = {n_2[0], n_2[1], n_2[2], n_3[0], n_3[1], + n_3[2], n_a[0], n_a[1], n_a[2]}; + A = change_basis(A_moy, r); + } + return A; + } // end of TransverseIsotropic + + TFEL_HOST_DEVICE static const tfel::math::st2tost2<3u, real> Oriented( + const StressType& young, + const real& nu, + const StressType& young_i, + const real& nu_i, + const tfel::math::tvector<3u, real>& n_a, + const LengthType& a, + const tfel::math::tvector<3u, real>& n_b, + const LengthType& b, + const LengthType& c) { + if ((nu > 0.5) || (nu < -1)) { + tfel::reportContractViolation("nu>0.5 or nu<-1"); + } + if (not(young > StressType{0})) { + tfel::reportContractViolation("E<=0"); + } + if (not((a > LengthType{0}) and (b > LengthType{0}) and + (c > LengthType{0}))) { + tfel::reportContractViolation("a<=0 or b<=0 or c<=0"); + } + if (not(tfel::math::ieee754::fpclassify( + tfel::math::VectorVectorDotProduct::exe< + real, tfel::math::tvector<3u, real>, + tfel::math::tvector<3u, real>>(n_a, n_b)) == FP_ZERO)) { + tfel::reportContractViolation("n_a and n_b not normals"); + } + if (tfel::math::ieee754::fpclassify(norm(n_a)) == FP_ZERO) { + tfel::reportContractViolation("n_a is null"); + } + if (tfel::math::ieee754::fpclassify(norm(n_b)) == FP_ZERO) { + tfel::reportContractViolation("n_b is null"); + } + using namespace tfel::math; + st2tost2<3u, real> A; + if ((a == b) and (b == c)) { + A = computeSphereLocalisationTensor(young, nu, + young_i, nu_i); + } else if (a == b) { + const tvector<3u, real> n_1 = tfel::math::cross_product(n_a, n_b); + A = computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_1, c / a); + } else if (a == c) { + A = computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_b, b / a); + } else if (b == c) { + A = computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_a, a / b); + } else { + A = computeEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_a, a, n_b, b, c); + } + return A; + } // end of Oriented + + }; // end of struct EllipsoidMeanLocalisator + + template + TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> + computeDiluteScheme(const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i, + const tfel::math::st2tost2<3u, real>& A) { + if ((nu > 0.5) || (nu < -1)) { + tfel::reportContractViolation("nu>0.5 or nu<-1"); + } + if (not(young > StressType{0})) { + tfel::reportContractViolation("E<=0"); + } + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + if ((nu_i > 0.5) || (nu_i < -1)) { + tfel::reportContractViolation("nu_i>0.5 or nu_i<-1"); + } + if (not(young_i > StressType{0})) { + tfel::reportContractViolation("E_i<=0"); + } + using namespace tfel::math; + st2tost2<3u, StressType> C_0; + static constexpr auto value = + StiffnessTensorAlterationCharacteristic::UNALTERED; + computeIsotropicStiffnessTensorII<3u, value, StressType, real>(C_0, young, + nu); + st2tost2<3u, StressType> C_i; + computeIsotropicStiffnessTensorII<3u, value, StressType, real>(C_i, young_i, + nu_i); + const auto C = C_i - C_0; + const auto Pr = C * A; + return C_0 + f * Pr; + } // end of DiluteScheme + + template + TFEL_HOST_DEVICE const tfel::math::st2tost2<3u, StressType> + computeMoriTanakaScheme(const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i, + const tfel::math::st2tost2<3u, real>& A) { + if ((nu > 0.5) || (nu < -1)) { + tfel::reportContractViolation("nu>0.5 or nu<-1"); + } + if (not(young > StressType{0})) { + tfel::reportContractViolation("E<=0"); + } + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + if ((nu_i > 0.5) || (nu_i < -1)) { + tfel::reportContractViolation("nu_i>0.5 or nu_i<-1"); + } + if (not(young_i > StressType{0})) { + tfel::reportContractViolation("E_i<=0"); + } + using namespace tfel::math; + st2tost2<3u, StressType> C_0; + static constexpr auto value = + StiffnessTensorAlterationCharacteristic::UNALTERED; + computeIsotropicStiffnessTensorII<3u, value, StressType, real>(C_0, young, + nu); + st2tost2<3u, StressType> C_i; + computeIsotropicStiffnessTensorII<3u, value, StressType, real>(C_i, young_i, + nu_i); + const auto C = C_i - C_0; + const auto Pr = C * A; + const auto inv = invert(f * A + (1 - f) * st2tost2<3u, real>::Id()); + const auto PPr = Pr * inv; + return C_0 + f * PPr; + } // end of MoriTanakaScheme + + template + TFEL_HOST_DEVICE const std::pair computeSphereDiluteScheme( + const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i) { + if ((nu > 0.5) || (nu < -1)) { + tfel::reportContractViolation("nu>0.5 or nu<-1"); + } + if (not(young > StressType{0})) { + tfel::reportContractViolation("E<=0"); + } + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + + const auto k0 = young / 3 / (1 - 2 * nu); + const auto mu0 = young / 2 / (1 + nu); + const auto k_i = young_i / 3 / (1 - 2 * nu_i); + const auto mu_i = young_i / 2 / (1 + nu_i); + const auto alpha0 = 3 * k0 / (3 * k0 + 4 * mu0); + const auto beta0 = 6 * (k0 + 2 * mu0) / 5 / (3 * k0 + 4 * mu0); + const auto muhom = + mu0 + f * (mu_i - mu0) / (1 + beta0 * (mu_i - mu0) / mu0); + const auto khom = k0 + f * (k_i - k0) / (1 + alpha0 * (k_i - k0) / k0); + const auto nuhom = (3 * khom - 2 * muhom) / (2 * muhom + 6 * khom); + const auto Ehom = 2 * muhom * (1 + nuhom); + return {Ehom, nuhom}; + } + + template + TFEL_HOST_DEVICE const std::pair + computeSphereMoriTanakaScheme(const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i) { + if ((nu > 0.5) || (nu < -1)) { + tfel::reportContractViolation("nu>0.5 or nu<-1"); + } + if (not(young > StressType{0})) { + tfel::reportContractViolation("E<=0"); + } + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + const auto k0 = young / 3 / (1 - 2 * nu); + const auto mu0 = young / 2 / (1 + nu); + const auto k_i = young_i / 3 / (1 - 2 * nu_i); + const auto mu_i = young_i / 2 / (1 + nu_i); + const auto muhom = + mu0 + + f * (mu_i - mu0) / + (1 + (1 - f) * (mu_i - mu0) / + (mu0 + mu0 * (9 * k0 + 8 * mu0) / 6 / (k0 + 2 * mu0))); + const auto khom = + k0 + f * (k_i - k0) / (1 + (1 - f) * (k_i - k0) / (k0 + 4 * mu0 / 3)); + const auto nuhom = (3 * khom - 2 * muhom) / (2 * muhom + 6 * khom); + const auto Ehom = 2 * muhom * (1 + nuhom); + return {Ehom, nuhom}; + } + + template + const std::pair computeIsotropicDiluteScheme( + const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i, + const LengthType& a, + const LengthType& b, + const LengthType& c) { + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + + const auto pair = + EllipsoidMeanLocalisator<3u, real, StressType, LengthType>::Isotropic( + young, nu, young_i, nu_i, a, b, c); + const auto ka = std::get<0>(pair); + const auto mu = std::get<1>(pair); + const auto k0 = young / 3 / (1 - 2 * nu); + const auto mu0 = young / 2 / (1 + nu); + const auto k_i = young_i / 3 / (1 - 2 * nu_i); + const auto mu_i = young_i / 2 / (1 + nu_i); + const auto muhom = mu0 + 2 * f * (mu_i - mu0) * mu; + const auto khom = k0 + 3 * f * (k_i - k0) * ka; + const auto nuhom = (3 * khom - 2 * muhom) / (2 * muhom + 6 * khom); + const auto Ehom = 2 * muhom * (1 + nuhom); + return {Ehom, nuhom}; + } + + template + const tfel::math::st2tost2<3u, StressType> + computeTransverseIsotropicDiluteScheme( + const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i, + const tfel::math::tvector<3u, real>& n_a, + const LengthType& a, + const LengthType& b, + const LengthType& c) { + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + const auto A = + EllipsoidMeanLocalisator<3u, real, StressType, + LengthType>::TransverseIsotropic(young, nu, + young_i, nu_i, + n_a, a, b, c); + return computeDiluteScheme(young, nu, f, young_i, nu_i, + A); + } + + template + const tfel::math::st2tost2<3u, StressType> computeOrientedDiluteScheme( + const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i, + const tfel::math::tvector<3u, real>& n_a, + const LengthType& a, + const tfel::math::tvector<3u, real>& n_b, + const LengthType& b, + const LengthType& c) { + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + const auto A = + EllipsoidMeanLocalisator<3u, real, StressType, LengthType>::Oriented( + young, nu, young_i, nu_i, n_a, a, n_b, b, c); + return computeDiluteScheme(young, nu, f, young_i, nu_i, + A); + } + + template + const std::pair computeIsotropicMoriTanakaScheme( + const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i, + const LengthType& a, + const LengthType& b, + const LengthType& c) { + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + const auto pair = + EllipsoidMeanLocalisator<3u, real, StressType, LengthType>::Isotropic( + young, nu, young_i, nu_i, a, b, c); + const auto ka = std::get<0>(pair); + const auto mu = std::get<1>(pair); + const auto k0 = young / 3 / (1 - 2 * nu); + const auto mu0 = young / 2 / (1 + nu); + const auto k_i = young_i / 3 / (1 - 2 * nu_i); + const auto mu_i = young_i / 2 / (1 + nu_i); + const auto muhom = mu0 + f * (mu_i - mu0) * mu / (f * mu + (1 - f) / 2); + const auto khom = k0 + f * (k_i - k0) * ka / (ka * f + (1 - f) / 3); + const auto nuhom = (3 * khom - 2 * muhom) / (2 * muhom + 6 * khom); + const auto Ehom = 2 * muhom * (1 + nuhom); + return {Ehom, nuhom}; + } + + template + const tfel::math::st2tost2<3u, StressType> + computeTransverseIsotropicMoriTanakaScheme( + const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i, + const tfel::math::tvector<3u, real>& n_a, + const LengthType& a, + const LengthType& b, + const LengthType& c) { + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + const auto A = + EllipsoidMeanLocalisator<3u, real, StressType, + LengthType>::TransverseIsotropic(young, nu, + young_i, nu_i, + n_a, a, b, c); + return computeMoriTanakaScheme(young, nu, f, young_i, + nu_i, A); + } + + template + const tfel::math::st2tost2<3u, StressType> computeOrientedMoriTanakaScheme( + const StressType& young, + const real& nu, + const real& f, + const StressType& young_i, + const real& nu_i, + const tfel::math::tvector<3u, real>& n_a, + const LengthType& a, + const tfel::math::tvector<3u, real>& n_b, + const LengthType& b, + const LengthType& c) { + if ((f < 0) || (f > 1)) { + tfel::reportContractViolation("f<0 or f>1"); + } + const auto A = + EllipsoidMeanLocalisator<3u, real, StressType, LengthType>::Oriented( + young, nu, young_i, nu_i, n_a, a, n_b, b, c); + return computeMoriTanakaScheme(young, nu, f, young_i, + nu_i, A); + } + +} // end of namespace tfel::material::homogenization::elasticity + +#endif /* LIB_TFEL_MATERIAL_LINEARHOMOGENIZATIONSCHEMES_IXX */ diff --git a/tests/Material/AnisotropicEshelbyTensor.cxx b/tests/Material/AnisotropicEshelbyTensor.cxx new file mode 100644 index 000000000..7b35fedc1 --- /dev/null +++ b/tests/Material/AnisotropicEshelbyTensor.cxx @@ -0,0 +1,213 @@ +/*! + * \file tests/Material/AnisotropicEshelbyTensor.cxx + * \brief + * \author Thomas Helfer + * \date 22/06/2021 + * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights + * reserved. + * This project is publicly released under either the GNU GPL Licence with + * linking exception or the CECILL-A licence. A copy of thoses licences are + * delivered with the sources of TFEL. CEA or EDF may also distribute this + * project under specific licensing conditions. + */ + +#ifdef NDEBUG +#undef NDEBUG +#endif NDEBUG + +#include +#include +#include +#include "TFEL/Config/TFELTypes.hxx" +#include "TFEL/Math/qt.hxx" +#include "TFEL/Material/IsotropicEshelbyTensor.hxx" +#include "TFEL/Material/AnisotropicEshelbyTensor.hxx" +#include "TFEL/Tests/TestCase.hxx" +#include "TFEL/Tests/TestProxy.hxx" +#include "TFEL/Tests/TestManager.hxx" +#include "TFEL/Material/StiffnessTensor.hxx" + +template +static constexpr T my_abs(const T& v) noexcept { + return v < T(0) ? -v : v; +} + +struct AnisotropicEshelbyTensorTest final : public tfel::tests::TestCase { + AnisotropicEshelbyTensorTest() + : tfel::tests::TestCase("TFEL/Material", "AnisotropicEshelbyTensor") { + } // end of AnisotropicEshelbyTensorTest + + tfel::tests::TestResult execute() override { + this->template errors(); + this->template test_Eshelby(); + this->template test_Eshelby2D(); + this->template test_localisator(); + return this->result; + } + + // must return a warning + private: + template + void errors() { + using namespace tfel::material::homogenization::elasticity; + using lg = typename tfel::config::Types<1u, NumericType, use_qt>::length; + using real = NumericType; + using stress = typename tfel::config::Types<1u, real, use_qt>::stress; + const auto nu = real{0.3}; + const auto young = stress{150e9}; + tfel::math::st2tost2<3u, stress> C_0; + const tfel::math::tvector<3u, real> n_a = {0., 0., 1.}; + const tfel::math::tvector<3u, real> n_b = {1., 0., 0.}; + using namespace tfel::material; + static constexpr auto value = + StiffnessTensorAlterationCharacteristic::UNALTERED; + computeIsotropicStiffnessTensorII<3u, value, stress, real>(C_0, young, nu); + // const auto S1 = + // computeAnisotropicEshelbyTensor(C_0,n_a,lg{0},n_b,lg{30},lg{3}); + // const auto S1 = + // computeAnisotropicEshelbyTensor(C_0,n_a,lg{-2},n_b,lg{30},lg{3}); + } + + private: + template + void test_Eshelby() { + using namespace tfel::material::homogenization::elasticity; + using real = NumericType; + using lg = typename tfel::config::Types<1u, real, use_qt>::length; + using stress = typename tfel::config::Types<1u, real, use_qt>::stress; + + static constexpr auto eps = std::sqrt(std::numeric_limits::epsilon()); + const auto nu = real{0.3}; + const auto young = stress{150e9}; + const tfel::math::tvector<3u, real> n_a = {1., 0., 0.}; + const tfel::math::tvector<3u, real> n_b = {0., 1., 0.}; + tfel::math::st2tost2<3u, stress> C_0; + using namespace tfel::material; + static constexpr auto value = + StiffnessTensorAlterationCharacteristic::UNALTERED; + computeIsotropicStiffnessTensorII<3u, value, stress, real>(C_0, young, nu); + const auto S1 = computeAnisotropicEshelbyTensor( + C_0, n_a, lg{3}, n_b, lg{2}, lg{1}, 14); + const auto S2 = computeEshelbyTensor(nu, lg{1}, lg{2}, lg{3}); + const auto S3 = computeEshelbyTensor(nu, lg{2}, lg{3}, lg{1}); + for (int i : {0, 1, 2, 3, 4, 5}) { + for (int j : {0, 1, 2, 3, 4, 5}) { + TFEL_TESTS_ASSERT(std::abs(S1(i, j) - S2(i, j)) < eps); + TFEL_TESTS_ASSERT(std::abs(S1(i, j) - S3(i, j)) < eps); + // std::cout << S1(i,j)-S2(i,j)<< " "<< eps<<'\n'; + // std::cout << S1(i,j)-S3(i,j) << " "<< eps<<'\n'; + } + } + + const auto SAxi1 = computeAnisotropicEshelbyTensor( + C_0, n_a, lg{30}, n_b, lg{3}, lg{3}, 14); + const auto SAxi2 = computeEshelbyTensor(nu, lg{30}, lg{3}, lg{3}); + for (int i : {0, 1, 2, 3, 4, 5}) { + for (int j : {0, 1, 2, 3, 4, 5}) { + TFEL_TESTS_ASSERT(std::abs(SAxi1(i, j) - SAxi2(i, j)) < eps); + // std::cout << SAxi1(i,j)-SAxi2(i,j) << " "<< eps <<'\n'; + } + } + + const auto SSph1 = computeAnisotropicEshelbyTensor( + C_0, n_a, lg{3}, n_b, lg{3}, lg{3}, 14); + const auto SSph2 = computeEshelbyTensor(nu, lg{30}, lg{30}, lg{30}); + for (int i : {0, 1, 2, 3, 4, 5}) { + for (int j : {0, 1, 2, 3, 4, 5}) { + TFEL_TESTS_ASSERT(std::abs(SSph1(i, j) - SSph2(i, j)) < eps); + // std::cout << SSph1(i,j)-SSph2(i,j) << " "<< eps <<'\n'; + } + } + + } // end of test_Eshelby + + private: + template + void test_Eshelby2D() { + using namespace tfel::material::homogenization::elasticity; + using real = NumericType; + using lg = typename tfel::config::Types<1u, real, use_qt>::length; + using stress = typename tfel::config::Types<1u, real, use_qt>::stress; + + static constexpr auto eps = std::sqrt(std::numeric_limits::epsilon()); + const auto nu = real{0.3}; + const auto young = stress{150e9}; + const tfel::math::tvector<2u, real> n_a = {1., 0.}; + const tfel::math::tvector<2u, real> n_b = {0., 1.}; + tfel::math::st2tost2<2u, stress> C_0; + using namespace tfel::material; + static constexpr auto value = + StiffnessTensorAlterationCharacteristic::UNALTERED; + computeIsotropicStiffnessTensorII<2u, value, stress, real>(C_0, young, nu); + const auto S2D1 = compute2DAnisotropicEshelbyTensor( + C_0, n_a, lg{3}, lg{2}, 14); + const auto S2D2 = computeEllipticCylinderEshelbyTensor(nu, real{1.5}); + for (int i : {0, 1, 2, 3}) { + for (int j : {0, 1, 2, 3}) { + TFEL_TESTS_ASSERT(std::abs(S2D1(i, j) - S2D2(i, j)) < eps); + // std::cout << S2D1(i,j) <<" " <( + C_0, n_a, lg{3}, lg{3}, 14); + const auto S2DCir2 = computeCircularCylinderEshelbyTensor(nu); + for (int i : {0, 1, 2, 3}) { + for (int j : {0, 1, 2, 3}) { + TFEL_TESTS_ASSERT(std::abs(S2DCir1(i, j) - S2DCir2(i, j)) < eps); + // std::cout << i <<" " << j <<" " << S2DCir1(i,j) <<" " << S2DCir2(i,j) + // << " "<< eps <<'\n'; + } + } + + } // end of test_Eshelby2D + + private: + template + void test_localisator() { + using namespace tfel::material::homogenization::elasticity; + using real = NumericType; + using lg = typename tfel::config::Types<1u, real, use_qt>::length; + using stress = typename tfel::config::Types<1u, real, use_qt>::stress; + + static constexpr auto eps = std::sqrt(std::numeric_limits::epsilon()); + const auto nu = real{0.3}; + const auto young = stress{1e9}; + const auto young_i = stress{150e9}; + const auto nu_i = real{0.1}; + const tfel::math::tvector<3u, real> n_a = {0., 0., 1.}; + const tfel::math::tvector<3u, real> n_b = {1., 0., 0.}; + tfel::math::st2tost2<3u, stress> C_0; + tfel::math::st2tost2<3u, stress> C_i; + using namespace tfel::material; + static constexpr auto value = + StiffnessTensorAlterationCharacteristic::UNALTERED; + computeIsotropicStiffnessTensorII<3u, value, stress, real>(C_0, young, nu); + computeIsotropicStiffnessTensorII<3u, value, stress, real>(C_i, young_i, + nu_i); + + const auto AAxis_0 = computeAnisotropicLocalisationTensor( + C_0, C_i, n_a, lg{20}, n_b, lg{2}, lg{1}); + const auto AAxis_1 = computeEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_a, lg{20}, n_b, lg{2}, lg{1}); + for (int i : {0, 1, 2, 3, 4, 5}) { + for (int j : {0, 1, 2, 3, 4, 5}) { + TFEL_TESTS_ASSERT(std::abs(AAxis_0(i, j) - AAxis_1(i, j)) < eps); + // std::cout << AAxis_0(i,j)-AAxis_1(i,j) << " "<< eps<<'\n'; + } + } + + } // end of test_localisator + +}; // end of struct AnisotropicEshelbyTensorTest + +TFEL_TESTS_GENERATE_PROXY(AnisotropicEshelbyTensorTest, + "AnisotropicEshelbyTensor"); + +/* coverity [UNCAUGHT_EXCEPT]*/ +int main() { + auto& m = tfel::tests::TestManager::getTestManager(); + m.addTestOutput(std::cout); + m.addXMLTestOutput("AnisotropicEshelbyTensor.xml"); + return m.execute().success() ? EXIT_SUCCESS : EXIT_FAILURE; +} // end of main diff --git a/tests/Material/CMakeLists.txt b/tests/Material/CMakeLists.txt index 34ca7ab55..ea0d20141 100644 --- a/tests/Material/CMakeLists.txt +++ b/tests/Material/CMakeLists.txt @@ -11,9 +11,10 @@ macro(tests_material test_arg) TFELUtilities TFELException TFELTests) endmacro(tests_material) -tests_material(EshelbyBasedHomogenization) -tests_material(Eshelby) +tests_material(LinearHomogenizationSchemes) +tests_material(IsotropicEshelbyTensor) tests_material(Lame) +tests_material(AnisotropicEshelbyTensor) tests_material(boundsCheck) tests_material(OrthotropicAxesConventionTest) tests_material(AbaqusTangentOperator) diff --git a/tests/Material/Eshelby.cxx b/tests/Material/IsotropicEshelbyTensor.cxx similarity index 86% rename from tests/Material/Eshelby.cxx rename to tests/Material/IsotropicEshelbyTensor.cxx index 1c4c98862..cf27fa905 100644 --- a/tests/Material/Eshelby.cxx +++ b/tests/Material/IsotropicEshelbyTensor.cxx @@ -1,5 +1,5 @@ /*! - * \file tests/Material/Eshelby.cxx + * \file tests/Material/IsotropicEshelbyTensor.cxx * \brief * \author Thomas Helfer * \date 22/06/2021 @@ -18,9 +18,10 @@ #include #include #include +#include #include "TFEL/Config/TFELTypes.hxx" #include "TFEL/Math/qt.hxx" -#include "TFEL/Material/Eshelby.hxx" +#include "TFEL/Material/IsotropicEshelbyTensor.hxx" #include "TFEL/Tests/TestCase.hxx" #include "TFEL/Tests/TestProxy.hxx" #include "TFEL/Tests/TestManager.hxx" @@ -31,10 +32,10 @@ static constexpr T my_abs(const T& v) noexcept { return v < T(0) ? -v : v; } -struct EshelbyTest final : public tfel::tests::TestCase { - EshelbyTest() - : tfel::tests::TestCase("TFEL/Material", "Eshelby") { - } // end of EshelbyTest +struct IsotropicEshelbyTensorTest final : public tfel::tests::TestCase { + IsotropicEshelbyTensorTest() + : tfel::tests::TestCase("TFEL/Material", "IsotropicEshelbyTensor") { + } // end of IsotropicEshelbyTensorTest tfel::tests::TestResult execute() override { this->template test_Eshelby_2d(); @@ -71,10 +72,10 @@ struct EshelbyTest final : public tfel::tests::TestCase { for (int j : {0, 1, 2, 3}) { TFEL_TESTS_ASSERT(std::abs(S2d_1(i, j) - S2d_2(i, j)) < eps); // std::cout << S1(i,j) << " " << S2(i,j) << " " << value << '\n'; - }; - }; + } + } } - }; + } // tests compilation of the Eshelby tensor functions private: @@ -98,7 +99,7 @@ struct EshelbyTest final : public tfel::tests::TestCase { const auto S2 = computeEshelbyTensor(nu, lg{1}, lg{1}, lg{2}); const auto S3 = computeEshelbyTensor(nu, lg{1}, lg{1}, lg{1}); } - }; + } // must return an error private: @@ -110,12 +111,13 @@ struct EshelbyTest final : public tfel::tests::TestCase { using namespace tfel::material::homogenization::elasticity; { // const auto S1 = computeSphereEshelbyTensor(real{1}); + // const auto S2 = computeAxisymmetricalEshelbyTensor(nu,real{0}); // const auto S3 = computeAxisymmetricalEshelbyTensor(nu,real{-10}); // const auto S4 = computeEshelbyTensor(nu,lg{-1},lg{3},lg{2}); // const auto S5 = computeEshelbyTensor(nu,lg{-1},lg{-1},lg{1}); } - }; + } private: template @@ -138,8 +140,8 @@ struct EshelbyTest final : public tfel::tests::TestCase { // std::cout << S1(i,j) << " " << S2(i,j) << " " << value << '\n'; TFEL_TESTS_ASSERT(std::abs(SSphere_1(i, j) - SSphere_3(i, j)) < eps); // std::cout << "13" << value << '\n'; - }; - }; + } + } } // AxisymmetricalEshelbyTensor with e=10 is equal to EshelbyTensor(30,3,3) // or EshelbyTensor(3,3,30) @@ -152,10 +154,10 @@ struct EshelbyTest final : public tfel::tests::TestCase { for (int j : {0, 1, 2, 3, 4, 5}) { TFEL_TESTS_ASSERT(std::abs(SAxis_1(i, j) - SAxis_2(i, j)) < eps); TFEL_TESTS_ASSERT(std::abs(SAxis_1(i, j) - SAxis_3(i, j)) < eps); - }; - }; + } + } } - }; + } private: template @@ -184,8 +186,8 @@ struct EshelbyTest final : public tfel::tests::TestCase { TFEL_TESTS_ASSERT( std::abs(SSphere_lim_1(i, j) - SSphere_lim_2(i, j)) < eps); // std::cout << S1(i,j) << " " << S2(i,j) << " " << value << '\n'; - }; - }; + } + } } // EshelbyTensor when a is near b must be near AxisymmetricalEshelbyTensor @@ -206,10 +208,10 @@ struct EshelbyTest final : public tfel::tests::TestCase { TFEL_TESTS_ASSERT(std::abs(SAxis_lim_1(i, j) - SAxis_lim_2(i, j)) < eps); // std::cout << S1(i,j) << " " << S2(i,j) << " " << value << '\n'; - }; - }; + } + } } - }; + } // tests compilation of localisation tensor functions private: @@ -240,7 +242,7 @@ struct EshelbyTest final : public tfel::tests::TestCase { const auto A3 = computeEllipsoidLocalisationTensor( young, nu, young_i, nu_i, n_a, a, n_b, b, c); } - }; + } // must return a warning private: @@ -269,6 +271,7 @@ struct EshelbyTest final : public tfel::tests::TestCase { using namespace tfel::material::homogenization::elasticity; { // const auto A1 = // computeSphereLocalisationTensor(young,real{1},young_i,nu_i); + // const auto A2 = // computeAxisymmetricalEllipsoidLocalisationTensor(young,nu,young_i,nu_i,n_a,real{0}); // const auto A3 = @@ -282,7 +285,7 @@ struct EshelbyTest final : public tfel::tests::TestCase { // const auto A7 = // computeEllipsoidLocalisationTensor(young,nu,young_i,nu_i,n_a,a,n__,b,c); } - }; + } // These functions must return the same thing private: @@ -307,7 +310,7 @@ struct EshelbyTest final : public tfel::tests::TestCase { young, nu, young_i, nu_i); const auto ASphere_2 = computeAxisymmetricalEllipsoidLocalisationTensor( - young, nu, young_i, nu_i, n_a, 1); + young, nu, young_i, nu_i, n_a, real{1}); const auto ASphere_3 = computeEllipsoidLocalisationTensor( young, nu, young_i, nu_i, n_a, lg{2}, n_b, lg{2.00000001}, @@ -316,8 +319,8 @@ struct EshelbyTest final : public tfel::tests::TestCase { for (int j : {0, 1, 2, 3, 4, 5}) { TFEL_TESTS_ASSERT(std::abs(ASphere_1(i, j) - ASphere_2(i, j)) < eps); TFEL_TESTS_ASSERT(std::abs(ASphere_1(i, j) - ASphere_3(i, j)) < eps); - }; - }; + } + } } { const auto AAxis_1 = @@ -331,19 +334,44 @@ struct EshelbyTest final : public tfel::tests::TestCase { TFEL_TESTS_ASSERT(std::abs(AAxis_1(i, j) - AAxis_2(i, j)) < eps); // std::cout << A1(i,j) <<" " << i<< " " << j<< " " << A2(i,j) << // value << '\n'; - }; - }; + } + } } - }; -}; // end of struct EshelbyTest + // this test checks that the rotation of the principal axis of the + // axisymmetrical ellipsoid in the direction n_aa is ok : + // A_nnnn when n=n_aa is much bigger than A_tttt where t=n_bb + // + { + const tfel::math::tvector<3u, real> n_aa = {std::sqrt(2) / 2, 0., + std::sqrt(2) / 2}; + const tfel::math::tvector<3u, real> n_bb = {-std::sqrt(2) / 2, 0., + std::sqrt(2) / 2}; + const auto dem = real{1} / 2; + const tfel::math::stensor<3u, real> n_aa_n_aa = {dem, 0., dem, + 0., dem, 0.}; + const tfel::math::stensor<3u, real> n_bb_n_bb = {dem, 0., dem, + 0., -dem, 0.}; + + const auto Arot1 = + computeAxisymmetricalEllipsoidLocalisationTensor( + young, nu, young_i, nu_i, n_aa, real{20}); + const auto Annnn = n_aa_n_aa * (Arot1 * n_aa_n_aa); + const auto Atttt = n_bb_n_bb * (Arot1 * n_bb_n_bb); + const auto Attnn = n_bb_n_bb * (Arot1 * n_aa_n_aa); + // std::cout << Atttt(0) << Annnn(0) << Attnn(0) << std::endl; + TFEL_TESTS_ASSERT(Annnn(0) - 5 * Atttt(0) > 0); + } + } + +}; // end of struct IsotropicEshelbyTensorTest -TFEL_TESTS_GENERATE_PROXY(EshelbyTest, "Eshelby"); +TFEL_TESTS_GENERATE_PROXY(IsotropicEshelbyTensorTest, "IsotropicEshelbyTensor"); /* coverity [UNCAUGHT_EXCEPT]*/ int main() { auto& m = tfel::tests::TestManager::getTestManager(); m.addTestOutput(std::cout); - m.addXMLTestOutput("Eshelby.xml"); + m.addXMLTestOutput("IsotropicEshelbyTensor.xml"); return m.execute().success() ? EXIT_SUCCESS : EXIT_FAILURE; } // end of main diff --git a/tests/Material/EshelbyBasedHomogenization.cxx b/tests/Material/LinearHomogenizationSchemes.cxx similarity index 96% rename from tests/Material/EshelbyBasedHomogenization.cxx rename to tests/Material/LinearHomogenizationSchemes.cxx index 5da36092c..0e18f1478 100644 --- a/tests/Material/EshelbyBasedHomogenization.cxx +++ b/tests/Material/LinearHomogenizationSchemes.cxx @@ -1,5 +1,5 @@ /*! - * \file tests/Material/EshelbyBasedHomogenization.cxx + * \file tests/Material/LinearHomogenizationSchemes.cxx * \brief * \author Antoine Martin * \date 25/10/2024 @@ -22,7 +22,7 @@ #include "TFEL/Config/TFELTypes.hxx" #include "TFEL/Math/qt.hxx" #include "TFEL/Math/General/ConstExprMathFunctions.hxx" -#include "TFEL/Material/EshelbyBasedHomogenization.hxx" +#include "TFEL/Material/LinearHomogenizationSchemes.hxx" #include "TFEL/Tests/TestCase.hxx" #include "TFEL/Tests/TestProxy.hxx" #include "TFEL/Tests/TestManager.hxx" @@ -33,15 +33,16 @@ static constexpr T my_abs(const T& v) noexcept { return v < T(0) ? -v : v; } -struct EshelbyBasedHomogenizationTest final : public tfel::tests::TestCase { - EshelbyBasedHomogenizationTest() - : tfel::tests::TestCase("TFEL/Material", "EshelbyBasedHomogenization") { - } // end of EshelbyBasedHomogenizationTest +struct LinearHomogenizationSchemesTest final : public tfel::tests::TestCase { + LinearHomogenizationSchemesTest() + : tfel::tests::TestCase("TFEL/Material", "LinearHomogenizationSchemes") { + } // end of LinearHomogenizationSchemesTest tfel::tests::TestResult execute() override { using real = double; constexpr bool qt = true; using stress = typename tfel::config::Types<1u, real, qt>::stress; using length = typename tfel::config::Types<1u, real, qt>::length; + this->template test1(); this->template errors(); this->template test3(); @@ -397,15 +398,15 @@ struct EshelbyBasedHomogenizationTest final : public tfel::tests::TestCase { } } -}; // end of struct EshelbyBasedHomogenizationTest +}; // end of struct LinearHomogenizationSchemesTest -TFEL_TESTS_GENERATE_PROXY(EshelbyBasedHomogenizationTest, - "EshelbyBasedHomogenization"); +TFEL_TESTS_GENERATE_PROXY(LinearHomogenizationSchemesTest, + "LinearHomogenizationSchemes"); /* coverity [UNCAUGHT_EXCEPT]*/ int main() { auto& m = tfel::tests::TestManager::getTestManager(); m.addTestOutput(std::cout); - m.addXMLTestOutput("EshelbyBasedHomogenization.xml"); + m.addXMLTestOutput("LinearHomogenizationSchemes.xml"); return m.execute().success() ? EXIT_SUCCESS : EXIT_FAILURE; } // end of main