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