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