Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <array>
       7              : #include <cstddef>
       8              : #include <memory>
       9              : #include <pup.h>
      10              : 
      11              : #include "DataStructures/Tensor/TypeAliases.hpp"
      12              : #include "Options/Context.hpp"
      13              : #include "Options/String.hpp"
      14              : #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/ConstitutiveRelation.hpp"
      15              : #include "Utilities/Serialization/CharmPupable.hpp"
      16              : #include "Utilities/TMPL.hpp"
      17              : 
      18              : /// \cond
      19              : class DataVector;
      20              : /// \endcond
      21              : 
      22            1 : namespace Elasticity::ConstitutiveRelations {
      23              : 
      24              : /*!
      25              :  * \brief A general homogeneous anisotropic linear elastic material in 3D
      26              :  *
      27              :  * \details The material is specified by a symmetric 6x6 stiffness matrix in
      28              :  * Voigt notation. The stress-strain relation is
      29              :  *
      30              :  * \f[
      31              :  * T^{ij} = -Y^{ijkl} S_{kl}
      32              :  * \f]
      33              :  *
      34              :  * where \f$S_{ij}=\nabla_{(i}u_{j)}\f$ is the (symmetric) strain tensor.
      35              :  *
      36              :  * We interpret the Voigt strain vector with engineering shear strains:
      37              :  *
      38              :  * \f[
      39              :  * e = [S_{xx}, S_{yy}, S_{zz}, 2S_{yz}, 2S_{xz}, 2S_{xy}]
      40              :  * \f]
      41              :  *
      42              :  * and the Voigt stress vector as:
      43              :  *
      44              :  * \f[
      45              :  * s = [T_{xx}, T_{yy}, T_{zz}, T_{yz}, T_{xz}, T_{xy}]
      46              :  * \f]
      47              :  *
      48              :  * Then \f$ s = -C e \f$, where \f$C\f$ is the 6x6 stiffness matrix.
      49              :  *
      50              :  * The matrix must be symmetric and positive definite, so that the strain
      51              :  * energy density \f$e^T C e / 2\f$ is positive for every nonzero strain. The
      52              :  * positive-definiteness check is relative: the smallest eigenvalue of \f$C\f$
      53              :  * must exceed `min_relative_eigenvalue` times its largest eigenvalue, so the
      54              :  * check does not depend on the units of \f$C\f$ and numerically singular
      55              :  * matrices are rejected deterministically.
      56              :  */
      57            1 : class AnisotropicHomogeneous : public ConstitutiveRelation<3> {
      58              :  public:
      59            0 :   static constexpr size_t volume_dim = 3;
      60              : 
      61              :   /// The smallest eigenvalue of the Voigt matrix must exceed this fraction of
      62              :   /// the largest eigenvalue for the matrix to count as positive definite.
      63            1 :   static constexpr double min_relative_eigenvalue = 1.e-12;
      64              : 
      65            0 :   struct StiffnessMatrix {
      66            0 :     using type = std::array<std::array<double, 6>, 6>;
      67            0 :     static constexpr Options::String help = {
      68              :         "Symmetric, positive definite 6x6 stiffness matrix in Voigt "
      69              :         "notation (smallest eigenvalue > 1e-12 (`min_relative_eigenvalue`) "
      70              :         "times the largest). Interpreted with engineering shear strains: "
      71              :         "[Sxx,Syy,Szz,2Syz,2Sxz,2Sxy]. Entries are in the same units of "
      72              :         "stress as the other materials (e.g. GPa in the input files)."};
      73              :   };
      74              : 
      75            0 :   using options = tmpl::list<StiffnessMatrix>;
      76              : 
      77            0 :   static constexpr Options::String help = {
      78              :       "A constitutive relation that describes a homogeneous, anisotropic "
      79              :       "material in terms of a symmetric 6x6 stiffness matrix in Voigt "
      80              :       "notation. The matrix entries are in the same units of stress as the "
      81              :       "other materials (e.g. GPa in the input files)."};
      82              : 
      83            0 :   AnisotropicHomogeneous();
      84              : 
      85            0 :   AnisotropicHomogeneous(const AnisotropicHomogeneous&) = default;
      86            0 :   AnisotropicHomogeneous& operator=(const AnisotropicHomogeneous&) = default;
      87            0 :   AnisotropicHomogeneous(AnisotropicHomogeneous&&) = default;
      88            0 :   AnisotropicHomogeneous& operator=(AnisotropicHomogeneous&&) = default;
      89            0 :   ~AnisotropicHomogeneous() override = default;
      90              : 
      91              :   /// The stiffness matrix must be symmetric and positive definite (see the
      92              :   /// class documentation), which is checked in all build types (not only in
      93              :   /// Debug builds).
      94            1 :   explicit AnisotropicHomogeneous(
      95              :       std::array<std::array<double, 6>, 6> stiffness,
      96              :       const Options::Context& context = {});
      97              : 
      98            1 :   std::unique_ptr<ConstitutiveRelation<3>> get_clone() const override;
      99              : 
     100            0 :   void stress(gsl::not_null<tnsr::II<DataVector, 3>*> stress,
     101              :               const tnsr::ii<DataVector, 3>& strain,
     102              :               const tnsr::I<DataVector, 3>& x) const override;
     103              : 
     104            0 :   const std::array<std::array<double, 6>, 6>& stiffness_matrix() const;
     105              : 
     106              :   // clang-tidy: no runtime references
     107            0 :   void pup(PUP::er& /*p*/) override;  // NOLINT
     108              : 
     109            0 :   explicit AnisotropicHomogeneous(CkMigrateMessage* /*unused*/) {}
     110              : 
     111            0 :   WRAPPED_PUPable_decl_base_template(  // NOLINT
     112              :       SINGLE_ARG(ConstitutiveRelation<3>), AnisotropicHomogeneous);
     113              : 
     114              :  private:
     115            0 :   std::array<std::array<double, 6>, 6> stiffness_{};
     116              : };
     117              : 
     118            0 : bool operator==(const AnisotropicHomogeneous& lhs,
     119              :                 const AnisotropicHomogeneous& rhs);
     120            0 : bool operator!=(const AnisotropicHomogeneous& lhs,
     121              :                 const AnisotropicHomogeneous& rhs);
     122              : 
     123              : }  // namespace Elasticity::ConstitutiveRelations