Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <memory>
       8             : #include <pup.h>
       9             : #include <tuple>
      10             : #include <type_traits>
      11             : #include <unordered_map>
      12             : #include <unordered_set>
      13             : #include <utility>
      14             : #include <vector>
      15             : 
      16             : #include "DataStructures/Tensor/Tensor.hpp"
      17             : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
      18             : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
      19             : #include "Utilities/Serialization/CharmPupable.hpp"
      20             : #include "Utilities/TMPL.hpp"
      21             : 
      22             : namespace domain::CoordinateMaps {
      23             : 
      24             : /*!
      25             :  * \brief A composition of coordinate maps at runtime
      26             :  *
      27             :  * Composes a sequence of `domain::CoordinateMapBase` that step through the
      28             :  * `Frames`. The result is another `domain::CoordinateMapBase`. This is
      29             :  * different to `domain::CoordinateMap`, which does the composition at compile
      30             :  * time. The use cases are different:
      31             :  *
      32             :  * - Use `domain::CoordinateMap` to compose maps at compile time to go from one
      33             :  *   (named) frame to another using any number of coordinate transformation. The
      34             :  *   coordinate transformations are concatenated to effectively form a single
      35             :  *   map, and intermediate frames have no meaning. This has the performance
      36             :  *   benefit that looking up pointers and calling into virtual member functions
      37             :  *   of intermediate maps is avoided, and it has the semantic benefit that
      38             :  *   intermediate frames without meaning are not named or even accessible.
      39             :  *   **Example:** A static BlockLogical -> Grid map that deforms the logical
      40             :  *   cube to a wedge, applies a radial redistribution of grid points, and
      41             :  *   translates + rotates the wedge.
      42             :  * - Use `domain::CoordinateMaps::Composition` (this class) to compose maps at
      43             :  *   runtime to step through a sequence of (named) frames.
      44             :  *   **Example:** A time-dependent ElementLogical -> BlockLogical -> Grid ->
      45             :  *   Inertial map that applies an affine transformation between the
      46             :  *   ElementLogical and BlockLogical frames (see
      47             :  *   domain::element_to_block_logical_map), then deforms the logical cube to a
      48             :  *   wedge using the map described above (Grid frame), and then rotates the grid
      49             :  *   with a time-dependent rotation map (Inertial frame).
      50             :  *
      51             :  * \warning Think about performance implications before using this Composition
      52             :  * class. In an evolution it's usually advantageous to keep at least some of the
      53             :  * maps separate, e.g. to avoid reevaluating the static maps in every time step.
      54             :  * You can also access individual components of the composition in this class.
      55             :  *
      56             :  * \tparam Frames The list of frames in the composition, as a tmpl::list<>.
      57             :  * The first entry in the list is the source frame, and the last entry is the
      58             :  * target frame of the composition. Maps in the composition step through the
      59             :  * `Frames`. For example, if `Frames = tmpl::list<Frame::ElementLogical,
      60             :  * Frame::BlockLogical, Frame::Inertial>`, then the composition has two maps:
      61             :  * ElementLogical -> BlockLogical and BlockLogical -> Inertial.
      62             :  */
      63             : template <typename Frames, size_t Dim,
      64             :           typename Is = std::make_index_sequence<tmpl::size<Frames>::value - 1>>
      65           1 : struct Composition;
      66             : 
      67             : template <typename Frames, size_t Dim, size_t... Is>
      68             : struct Composition<Frames, Dim, std::index_sequence<Is...>>
      69             :     : public CoordinateMapBase<tmpl::front<Frames>, tmpl::back<Frames>, Dim> {
      70             :   using frames = Frames;
      71             :   using SourceFrame = tmpl::front<Frames>;
      72             :   using TargetFrame = tmpl::back<Frames>;
      73             :   static constexpr size_t num_frames = tmpl::size<Frames>::value;
      74             :   using Base = CoordinateMapBase<SourceFrame, TargetFrame, Dim>;
      75             :   using FuncOfTimeMap = std::unordered_map<
      76             :       std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>;
      77             : 
      78             :   Composition() = default;
      79             :   Composition(const Composition& rhs) { *this = rhs; }
      80             :   Composition& operator=(const Composition& rhs);
      81             :   Composition(Composition&& /*rhs*/) = default;
      82             :   Composition& operator=(Composition&& /*rhs*/) = default;
      83             :   ~Composition() override = default;
      84             : 
      85             :   std::unique_ptr<Base> get_clone() const override {
      86             :     return std::make_unique<Composition>(*this);
      87             :   }
      88             : 
      89             :   /// \cond
      90             :   explicit Composition(CkMigrateMessage* /*m*/) {}
      91             :   using PUP::able::register_constructor;
      92             :   WRAPPED_PUPable_decl_template(Composition);  // NOLINT
      93             :   /// \endcond
      94             : 
      95             :   Composition(std::unique_ptr<CoordinateMapBase<
      96             :                   tmpl::at<frames, tmpl::size_t<Is>>,
      97             :                   tmpl::at<frames, tmpl::size_t<Is + 1>>, Dim>>... maps);
      98             : 
      99             :   const std::tuple<std::unique_ptr<
     100             :       CoordinateMapBase<tmpl::at<frames, tmpl::size_t<Is>>,
     101             :                         tmpl::at<frames, tmpl::size_t<Is + 1>>, Dim>>...>&
     102             :   maps() const {
     103             :     return maps_;
     104             :   }
     105             : 
     106             :   template <typename LocalSourceFrame, typename LocalTargetFrame>
     107             :   const CoordinateMapBase<LocalSourceFrame, LocalTargetFrame, Dim>&
     108             :   get_component() const {
     109             :     return *get<std::unique_ptr<
     110             :         CoordinateMapBase<LocalSourceFrame, LocalTargetFrame, Dim>>>(maps_);
     111             :   }
     112             : 
     113             :   bool is_identity() const override;
     114             : 
     115             :   bool inv_jacobian_is_time_dependent() const override;
     116             : 
     117             :   bool jacobian_is_time_dependent() const override;
     118             : 
     119             :   const std::unordered_set<std::string>& function_of_time_names()
     120             :       const override {
     121             :     return function_of_time_names_;
     122             :   }
     123             : 
     124             :   tnsr::I<double, Dim, TargetFrame> operator()(
     125             :       tnsr::I<double, Dim, SourceFrame> source_point,
     126             :       double time = std::numeric_limits<double>::signaling_NaN(),
     127             :       const FuncOfTimeMap& functions_of_time = {}) const override;
     128             : 
     129             :   tnsr::I<DataVector, Dim, TargetFrame> operator()(
     130             :       tnsr::I<DataVector, Dim, SourceFrame> source_point,
     131             :       double time = std::numeric_limits<double>::signaling_NaN(),
     132             :       const FuncOfTimeMap& functions_of_time = {}) const override;
     133             : 
     134             :   std::optional<tnsr::I<double, Dim, SourceFrame>> inverse(
     135             :       tnsr::I<double, Dim, TargetFrame> target_point,
     136             :       double time = std::numeric_limits<double>::signaling_NaN(),
     137             :       const FuncOfTimeMap& functions_of_time = {}) const override;
     138             : 
     139             :   InverseJacobian<double, Dim, SourceFrame, TargetFrame> inv_jacobian(
     140             :       tnsr::I<double, Dim, SourceFrame> source_point,
     141             :       double time = std::numeric_limits<double>::signaling_NaN(),
     142             :       const FuncOfTimeMap& functions_of_time = {}) const override;
     143             : 
     144             :   InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame> inv_jacobian(
     145             :       tnsr::I<DataVector, Dim, SourceFrame> source_point,
     146             :       double time = std::numeric_limits<double>::signaling_NaN(),
     147             :       const FuncOfTimeMap& functions_of_time = {}) const override;
     148             : 
     149             :   Jacobian<double, Dim, SourceFrame, TargetFrame> jacobian(
     150             :       tnsr::I<double, Dim, SourceFrame> source_point,
     151             :       double time = std::numeric_limits<double>::signaling_NaN(),
     152             :       const FuncOfTimeMap& functions_of_time = {}) const override;
     153             : 
     154             :   Jacobian<DataVector, Dim, SourceFrame, TargetFrame> jacobian(
     155             :       tnsr::I<DataVector, Dim, SourceFrame> source_point,
     156             :       double time = std::numeric_limits<double>::signaling_NaN(),
     157             :       const FuncOfTimeMap& functions_of_time = {}) const override;
     158             : 
     159             :   [[noreturn]] std::tuple<
     160             :       tnsr::I<double, Dim, TargetFrame>,
     161             :       InverseJacobian<double, Dim, SourceFrame, TargetFrame>,
     162             :       Jacobian<double, Dim, SourceFrame, TargetFrame>,
     163             :       tnsr::I<double, Dim, TargetFrame>>
     164             :   coords_frame_velocity_jacobians(
     165             :       tnsr::I<double, Dim, SourceFrame> source_point,
     166             :       double time = std::numeric_limits<double>::signaling_NaN(),
     167             :       const FuncOfTimeMap& functions_of_time = {}) const override;
     168             : 
     169             :   [[noreturn]] std::tuple<
     170             :       tnsr::I<DataVector, Dim, TargetFrame>,
     171             :       InverseJacobian<DataVector, Dim, SourceFrame, TargetFrame>,
     172             :       Jacobian<DataVector, Dim, SourceFrame, TargetFrame>,
     173             :       tnsr::I<DataVector, Dim, TargetFrame>>
     174             :   coords_frame_velocity_jacobians(
     175             :       tnsr::I<DataVector, Dim, SourceFrame> source_point,
     176             :       double time = std::numeric_limits<double>::signaling_NaN(),
     177             :       const FuncOfTimeMap& functions_of_time = {}) const override;
     178             : 
     179             :   [[noreturn]] std::unique_ptr<CoordinateMapBase<SourceFrame, Frame::Grid, Dim>>
     180             :   get_to_grid_frame() const override;
     181             : 
     182             :   // NOLINTNEXTLINE(google-runtime-references)
     183             :   void pup(PUP::er& p) override;
     184             : 
     185             :  private:
     186             :   template <typename DataType>
     187             :   tnsr::I<DataType, Dim, TargetFrame> call_impl(
     188             :       tnsr::I<DataType, Dim, SourceFrame> source_point,
     189             :       const double time = std::numeric_limits<double>::signaling_NaN(),
     190             :       const FuncOfTimeMap& functions_of_time = {}) const;
     191             : 
     192             :   template <typename DataType>
     193             :   std::optional<tnsr::I<DataType, Dim, SourceFrame>> inverse_impl(
     194             :       tnsr::I<DataType, Dim, TargetFrame> target_point,
     195             :       const double time = std::numeric_limits<double>::signaling_NaN(),
     196             :       const FuncOfTimeMap& functions_of_time = {}) const;
     197             : 
     198             :   template <typename DataType>
     199             :   InverseJacobian<DataType, Dim, SourceFrame, TargetFrame> inv_jacobian_impl(
     200             :       tnsr::I<DataType, Dim, SourceFrame> source_point,
     201             :       const double time = std::numeric_limits<double>::signaling_NaN(),
     202             :       const FuncOfTimeMap& functions_of_time = {}) const;
     203             : 
     204             :   template <typename DataType>
     205             :   Jacobian<DataType, Dim, SourceFrame, TargetFrame> jacobian_impl(
     206             :       tnsr::I<DataType, Dim, SourceFrame> source_point,
     207             :       const double time = std::numeric_limits<double>::signaling_NaN(),
     208             :       const FuncOfTimeMap& functions_of_time = {}) const;
     209             : 
     210             :   bool is_equal_to(const CoordinateMapBase<SourceFrame, TargetFrame, Dim>&
     211             :                        other) const override;
     212             : 
     213             :   std::tuple<std::unique_ptr<
     214             :       CoordinateMapBase<tmpl::at<frames, tmpl::size_t<Is>>,
     215             :                         tmpl::at<frames, tmpl::size_t<Is + 1>>, Dim>>...>
     216             :       maps_;
     217             :   std::unordered_set<std::string> function_of_time_names_;
     218             : };
     219             : 
     220             : // Template deduction guide
     221             : template <typename FirstMap, typename... Maps>
     222           0 : Composition(std::unique_ptr<FirstMap>, std::unique_ptr<Maps>... maps)
     223             :     -> Composition<tmpl::list<typename FirstMap::source_frame,
     224             :                               typename FirstMap::target_frame,
     225             :                               typename Maps::target_frame...>,
     226             :                    FirstMap::dim>;
     227             : 
     228             : /// \cond
     229             : template <typename Frames, size_t Dim, size_t... Is>
     230             : PUP::able::PUP_ID
     231             :     Composition<Frames, Dim, std::index_sequence<Is...>>::my_PUP_ID = 0;
     232             : /// \endcond
     233             : 
     234             : }  // namespace domain::CoordinateMaps