Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <limits>
       7             : #include <pup.h>
       8             : 
       9             : #include "DataStructures/CachedTempBuffer.hpp"
      10             : #include "DataStructures/DataBox/Prefixes.hpp"
      11             : #include "DataStructures/DataVector.hpp"
      12             : #include "DataStructures/TaggedTuple.hpp"
      13             : #include "DataStructures/Tensor/Tensor.hpp"
      14             : #include "Elliptic/Systems/Elasticity/Tags.hpp"
      15             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      16             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      17             : #include "Options/String.hpp"
      18             : #include "PointwiseFunctions/Elasticity/ConstitutiveRelations/IsotropicHomogeneous.hpp"
      19             : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
      20             : #include "Utilities/ContainerHelpers.hpp"
      21             : #include "Utilities/Serialization/CharmPupable.hpp"
      22             : #include "Utilities/TMPL.hpp"
      23             : 
      24           0 : namespace Elasticity::Solutions {
      25             : 
      26             : namespace detail {
      27             : template <typename DataType>
      28             : struct BentBeamVariables {
      29             :   using Cache = CachedTempBuffer<
      30             :       Tags::Displacement<2>,
      31             :       ::Tags::deriv<Tags::Displacement<2>, tmpl::size_t<2>, Frame::Inertial>,
      32             :       Tags::Strain<2>, Tags::MinusStress<2>, Tags::PotentialEnergyDensity<2>,
      33             :       ::Tags::FixedSource<Tags::Displacement<2>>>;
      34             : 
      35             :   const tnsr::I<DataType, 2>& x;
      36             :   const double length;
      37             :   const double height;
      38             :   const double bending_moment;
      39             :   const ConstitutiveRelations::IsotropicHomogeneous<2>& constitutive_relation;
      40             : 
      41             :   void operator()(gsl::not_null<tnsr::I<DataType, 2>*> displacement,
      42             :                   gsl::not_null<Cache*> cache,
      43             :                   Tags::Displacement<2> /*meta*/) const;
      44             :   void operator()(gsl::not_null<tnsr::iJ<DataType, 2>*> deriv_displacement,
      45             :                   gsl::not_null<Cache*> cache,
      46             :                   ::Tags::deriv<Tags::Displacement<2>, tmpl::size_t<2>,
      47             :                                 Frame::Inertial> /*meta*/) const;
      48             :   void operator()(gsl::not_null<tnsr::ii<DataType, 2>*> strain,
      49             :                   gsl::not_null<Cache*> cache, Tags::Strain<2> /*meta*/) const;
      50             :   void operator()(gsl::not_null<tnsr::II<DataType, 2>*> minus_stress,
      51             :                   gsl::not_null<Cache*> cache,
      52             :                   Tags::MinusStress<2> /*meta*/) const;
      53             :   void operator()(gsl::not_null<Scalar<DataType>*> potential_energy_density,
      54             :                   gsl::not_null<Cache*> cache,
      55             :                   Tags::PotentialEnergyDensity<2> /*meta*/) const;
      56             :   void operator()(
      57             :       gsl::not_null<tnsr::I<DataType, 2>*> fixed_source_for_displacement,
      58             :       gsl::not_null<Cache*> cache,
      59             :       ::Tags::FixedSource<Tags::Displacement<2>> /*meta*/) const;
      60             : };
      61             : }  // namespace detail
      62             : 
      63             : /*!
      64             :  * \brief A state of pure bending of an elastic beam in 2D
      65             :  *
      66             :  * \details This solution describes a 2D slice through an elastic beam of length
      67             :  * \f$L\f$ and height \f$H\f$, centered around (0, 0), that is subject to a
      68             :  * bending moment \f$M=\int T^{xx}y\mathrm{d}y\f$ (see e.g.
      69             :  * \cite ThorneBlandford2017, Eq. 11.41c for a bending moment in 1D). The beam
      70             :  * material is characterized by an isotropic and homogeneous constitutive
      71             :  * relation \f$Y^{ijkl}\f$ in the plane-stress approximation (see
      72             :  * `Elasticity::ConstitutiveRelations::IsotropicHomogeneous`). In this scenario,
      73             :  * no body-forces \f$f_\mathrm{ext}^j\f$ act on the material, so the
      74             :  * \ref Elasticity equations reduce to \f$\nabla_i T^{ij}=0\f$, but the bending
      75             :  * moment \f$M\f$ generates the stress
      76             :  *
      77             :  * \f{align}
      78             :  * T^{xx} &= \frac{12 M}{H^3} y \\
      79             :  * T^{xy} &= 0 = T^{yy} \text{.}
      80             :  * \f}
      81             :  *
      82             :  * By fixing the rigid-body motions to
      83             :  *
      84             :  * \f[
      85             :  * \xi^x(0,y)=0 \quad \text{and} \quad \xi^y\left(\pm \frac{L}{2},0\right)=0
      86             :  * \f]
      87             :  *
      88             :  * we find that this stress is produced by the displacement field
      89             :  *
      90             :  * \f{align}
      91             :  * \xi^x&=-\frac{12 M}{EH^3}xy \\
      92             :  * \xi^y&=\frac{6 M}{EH^3}\left(x^2+\nu y^2-\frac{L^2}{4}\right)
      93             :  * \f}
      94             :  *
      95             :  * in terms of the Young's modulus \f$E\f$ and the Poisson ration \f$\nu\f$ of
      96             :  * the material. The corresponding strain \f$S_{ij}=\partial_{(i}\xi_{j)}\f$ is
      97             :  *
      98             :  * \f{align}
      99             :  * S_{xx} &= -\frac{12 M}{EH^3} y \\
     100             :  * S_{yy} &= \frac{12 M}{EH^3} \nu y \\
     101             :  * S_{xy} &= S_{yx} = 0
     102             :  * \f}
     103             :  *
     104             :  * and the potential energy stored in the entire infinitesimal slice is
     105             :  *
     106             :  * \f[
     107             :  * \int_{-L/2}^{L/2} \int_{-H/2}^{H/2} U dy\,dx = \frac{6M^2}{EH^3}L \text{.}
     108             :  * \f]
     109             :  */
     110           1 : class BentBeam : public elliptic::analytic_data::AnalyticSolution {
     111             :  public:
     112           0 :   using constitutive_relation_type =
     113             :       Elasticity::ConstitutiveRelations::IsotropicHomogeneous<2>;
     114             : 
     115           0 :   struct Length {
     116           0 :     using type = double;
     117           0 :     static constexpr Options::String help{"The beam length"};
     118           0 :     static type lower_bound() { return 0.0; }
     119             :   };
     120           0 :   struct Height {
     121           0 :     using type = double;
     122           0 :     static constexpr Options::String help{"The beam height"};
     123           0 :     static type lower_bound() { return 0.0; }
     124             :   };
     125           0 :   struct BendingMoment {
     126           0 :     using type = double;
     127           0 :     static constexpr Options::String help{
     128             :         "The bending moment applied to the beam"};
     129           0 :     static type lower_bound() { return 0.0; }
     130             :   };
     131           0 :   struct Material {
     132           0 :     using type = constitutive_relation_type;
     133           0 :     static constexpr Options::String help{
     134             :         "The material properties of the beam"};
     135             :   };
     136             : 
     137           0 :   using options = tmpl::list<Length, Height, BendingMoment, Material>;
     138           0 :   static constexpr Options::String help{
     139             :       "A 2D slice through an elastic beam which is subject to a bending "
     140             :       "moment. The bending moment is applied along the length of the beam, "
     141             :       "i.e. the x-axis, so that the beam's left and right ends are bent "
     142             :       "towards the positive y-axis. It is measured in units of force."};
     143             : 
     144           0 :   BentBeam() = default;
     145           0 :   BentBeam(const BentBeam&) = default;
     146           0 :   BentBeam& operator=(const BentBeam&) = default;
     147           0 :   BentBeam(BentBeam&&) = default;
     148           0 :   BentBeam& operator=(BentBeam&&) = default;
     149           0 :   ~BentBeam() override = default;
     150           0 :   std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
     151             :       const override {
     152             :     return std::make_unique<BentBeam>(*this);
     153             :   }
     154             : 
     155             :   /// \cond
     156             :   explicit BentBeam(CkMigrateMessage* m)
     157             :       : elliptic::analytic_data::AnalyticSolution(m) {}
     158             :   using PUP::able::register_constructor;
     159             :   WRAPPED_PUPable_decl_template(BentBeam);  // NOLINT
     160             :   /// \endcond
     161             : 
     162           0 :   BentBeam(double length, double height, double bending_moment,
     163             :            constitutive_relation_type constitutive_relation)
     164             :       : length_(length),
     165             :         height_(height),
     166             :         bending_moment_(bending_moment),
     167             :         constitutive_relation_(std::move(constitutive_relation)) {}
     168             : 
     169           0 :   double length() const { return length_; }
     170           0 :   double height() const { return height_; }
     171           0 :   double bending_moment() const { return bending_moment_; }
     172             : 
     173           0 :   const constitutive_relation_type& constitutive_relation() const {
     174             :     return constitutive_relation_;
     175             :   }
     176             : 
     177             :   /// Return potential energy integrated over the whole beam material
     178           1 :   double potential_energy() const {
     179             :     return 6. * length_ * square(bending_moment_) /
     180             :            (cube(height_) * constitutive_relation_.youngs_modulus());
     181             :   }
     182             : 
     183             :   template <typename DataType, typename... RequestedTags>
     184           0 :   tuples::TaggedTuple<RequestedTags...> variables(
     185             :       const tnsr::I<DataType, 2>& x,
     186             :       tmpl::list<RequestedTags...> /*meta*/) const {
     187             :     using VarsComputer = detail::BentBeamVariables<DataType>;
     188             :     typename VarsComputer::Cache cache{get_size(*x.begin())};
     189             :     const VarsComputer computer{x, length_, height_, bending_moment_,
     190             :                                 constitutive_relation_};
     191             :     return {cache.get_var(computer, RequestedTags{})...};
     192             :   }
     193             : 
     194             :   template <typename... RequestedTags>
     195           0 :   tuples::TaggedTuple<RequestedTags...> variables(
     196             :       const tnsr::I<DataVector, 2>& x, const Mesh<2>& /*mesh*/,
     197             :       const InverseJacobian<DataVector, 2, Frame::ElementLogical,
     198             :                             Frame::Inertial>& /*inv_jacobian*/,
     199             :       tmpl::list<RequestedTags...> meta) const {
     200             :     return variables(x, meta);
     201             :   }
     202             : 
     203             :   /// NOLINTNEXTLINE(google-runtime-references)
     204           1 :   void pup(PUP::er& p) override {
     205             :     elliptic::analytic_data::AnalyticSolution::pup(p);
     206             :     p | length_;
     207             :     p | height_;
     208             :     p | bending_moment_;
     209             :     p | constitutive_relation_;
     210             :   }
     211             : 
     212             :  private:
     213           0 :   double length_{std::numeric_limits<double>::signaling_NaN()};
     214           0 :   double height_{std::numeric_limits<double>::signaling_NaN()};
     215           0 :   double bending_moment_{std::numeric_limits<double>::signaling_NaN()};
     216           0 :   constitutive_relation_type constitutive_relation_{};
     217             : };
     218             : 
     219           0 : bool operator==(const BentBeam& lhs, const BentBeam& rhs);
     220           0 : bool operator!=(const BentBeam& lhs, const BentBeam& rhs);
     221             : 
     222             : }  // namespace Elasticity::Solutions