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 "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             : tnsr::II<DataVector, 3> CubicCrystal::stress(
      31             :     const tnsr::ii<DataVector, 3>& strain,
      32             :     const tnsr::I<DataVector, 3>& /*x*/) const noexcept {
      33             :   auto result = make_with_value<tnsr::II<DataVector, 3>>(strain, 0.);
      34             :   for (size_t i = 0; i < 3; i++) {
      35             :     result.get(2, 2) -= strain.get(i, i);
      36             :   }
      37             :   result.get(2, 2) *= c_12_;
      38             :   for (size_t i = 0; i < 3; i++) {
      39             :     result.get(i, i) = result.get(2, 2) - (c_11_ - c_12_) * strain.get(i, i);
      40             :     for (size_t j = 0; j < i; j++) {
      41             :       result.get(i, j) = -2. * c_44_ * strain.get(i, j);
      42             :     }
      43             :   }
      44             :   return result;
      45             : }
      46             : 
      47             : double CubicCrystal::c_11() const noexcept { return c_11_; }
      48             : 
      49             : double CubicCrystal::c_12() const noexcept { return c_12_; }
      50             : 
      51             : double CubicCrystal::c_44() const noexcept { return c_44_; }
      52             : 
      53             : // through \lambda = c_{12} = \frac{E\nu}{(1+\nu)(1-2\nu)}
      54             : double CubicCrystal::youngs_modulus() const noexcept {
      55             :   return (c_11_ + 2. * c_12_) * (c_11_ - c_12_) / (c_11_ + c_12_);
      56             : }
      57             : 
      58             : double CubicCrystal::poisson_ratio() const noexcept {
      59             :   return 1. / (1. + (c_11_ / c_12_));
      60             : }
      61             : 
      62             : /// \cond
      63             : PUP::able::PUP_ID CubicCrystal::my_PUP_ID = 0;
      64             : /// \endcond
      65             : 
      66             : void CubicCrystal::pup(PUP::er& p) noexcept {
      67             :   p | c_11_;
      68             :   p | c_12_;
      69             :   p | c_44_;
      70             : }
      71             : 
      72           0 : bool operator==(const CubicCrystal& lhs, const CubicCrystal& rhs) noexcept {
      73             :   return lhs.c_11() == rhs.c_11() and lhs.c_12() == rhs.c_12() and
      74             :          lhs.c_44() == rhs.c_44();
      75             : }
      76             : 
      77           0 : bool operator!=(const CubicCrystal& lhs, const CubicCrystal& rhs) noexcept {
      78             :   return not(lhs == rhs);
      79             : }
      80             : 
      81             : }  // namespace Elasticity::ConstitutiveRelations