Line data Source code
1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/CubicCrystal.hpp" 5 : 6 : #include <array> 7 : #include <pup.h> // IWYU pragma: keep 8 : 9 : #include "DataStructures/DataVector.hpp" 10 : #include "DataStructures/Tensor/Tensor.hpp" // IWYU pragma: keep 11 : #include "DataStructures/Tensor/TypeAliases.hpp" 12 : #include "Utilities/ErrorHandling/Assert.hpp" 13 : #include "Utilities/MakeWithValue.hpp" 14 : 15 : // IWYU pragma: no_forward_declare Tensor 16 : 17 : namespace Elasticity::ConstitutiveRelations { 18 : 19 : CubicCrystal::CubicCrystal(const double c_11, const double c_12, 20 : const double c_44) noexcept 21 : : c_11_(c_11), c_12_(c_12), c_44_(c_44) { 22 : ASSERT( 23 : c_11_ >= c_12_, 24 : "c_11 must be bigger than c_12, but are c_11=" 25 : << c_11 << " and c_12=" << c_12 26 : << ". This is because the youngs_modulus " 27 : "must be positive and the poisson ratio smaller or equal to 0.5."); 28 : } 29 : 30 : void CubicCrystal::stress(const gsl::not_null<tnsr::II<DataVector, 3>*> stress, 31 : const tnsr::ii<DataVector, 3>& strain, 32 : const tnsr::I<DataVector, 3>& /*x*/) const noexcept { 33 : get<2, 2>(*stress) = -get<0, 0>(strain); 34 : for (size_t i = 1; i < 3; ++i) { 35 : stress->get(2, 2) -= strain.get(i, i); 36 : } 37 : stress->get(2, 2) *= c_12_; 38 : for (size_t i = 0; i < 3; ++i) { 39 : stress->get(i, i) = stress->get(2, 2) - (c_11_ - c_12_) * strain.get(i, i); 40 : for (size_t j = 0; j < i; ++j) { 41 : stress->get(i, j) = -2. * c_44_ * strain.get(i, j); 42 : } 43 : } 44 : } 45 : 46 : double CubicCrystal::c_11() const noexcept { return c_11_; } 47 : 48 : double CubicCrystal::c_12() const noexcept { return c_12_; } 49 : 50 : double CubicCrystal::c_44() const noexcept { return c_44_; } 51 : 52 : // through \lambda = c_{12} = \frac{E\nu}{(1+\nu)(1-2\nu)} 53 : double CubicCrystal::youngs_modulus() const noexcept { 54 : return (c_11_ + 2. * c_12_) * (c_11_ - c_12_) / (c_11_ + c_12_); 55 : } 56 : 57 : double CubicCrystal::poisson_ratio() const noexcept { 58 : return 1. / (1. + (c_11_ / c_12_)); 59 : } 60 : 61 : /// \cond 62 : PUP::able::PUP_ID CubicCrystal::my_PUP_ID = 0; 63 : /// \endcond 64 : 65 : void CubicCrystal::pup(PUP::er& p) noexcept { 66 : p | c_11_; 67 : p | c_12_; 68 : p | c_44_; 69 : } 70 : 71 0 : bool operator==(const CubicCrystal& lhs, const CubicCrystal& rhs) noexcept { 72 : return lhs.c_11() == rhs.c_11() and lhs.c_12() == rhs.c_12() and 73 : lhs.c_44() == rhs.c_44(); 74 : } 75 : 76 0 : bool operator!=(const CubicCrystal& lhs, const CubicCrystal& rhs) noexcept { 77 : return not(lhs == rhs); 78 : } 79 : 80 : } // namespace Elasticity::ConstitutiveRelations