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 <iosfwd>
       9             : 
      10             : #include "DataStructures/Tensor/TypeAliases.hpp"
      11             : #include "Domain/Structure/Direction.hpp"
      12             : #include "Domain/Structure/SegmentId.hpp"
      13             : #include "Domain/Structure/Side.hpp"
      14             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      15             : #include "Utilities/Gsl.hpp"
      16             : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
      17             : 
      18             : namespace PUP {
      19             : class er;
      20             : }  // namespace PUP
      21             : 
      22             : /*!
      23             :  * \ingroup ComputationalDomainGroup
      24             :  * \brief A mapping of the logical coordinate axes of a host to the logical
      25             :  * coordinate axes of a neighbor of the host.
      26             :  *
      27             :  * Given a `size_t dimension`, a `Direction`, a `SegmentId`, or a `Mesh` of the
      28             :  * host, an `OrientationMap` will give the corresponding value in the neighbor.
      29             :  *
      30             :  * \tparam VolumeDim the dimension of the blocks.
      31             :  *
      32             :  * See the [tutorial](@ref tutorial_orientations) for information on how
      33             :  * OrientationMaps are used and constructed.
      34             :  *
      35             :  * \note If there is no discrete rotation between logical coordinates (e.g. the
      36             :  * angular coordinates of a spherical shell abutting a wedge of a cubed sphere)
      37             :  * specify Direction<VolumeDim>::self() as the mapped direction
      38             :  *
      39             :  */
      40             : template <size_t VolumeDim>
      41           1 : class OrientationMap {
      42             :  public:
      43           0 :   static constexpr uint16_t aligned_mask = 0b1000000000000000;
      44           0 :   static constexpr uint16_t version_mask = 0b0111000000000000;
      45             : 
      46             :   /// \brief Creates an OrientationMap in an uninitialized state.
      47             :   ///
      48             :   /// This can be helpful for debugging code. If you would like the identity
      49             :   /// map, please use `create_aligned()`.
      50           1 :   OrientationMap();
      51             :   /// Mapped directions relative to the positive (`Side::Upper`) direction in
      52             :   /// each logical direction.
      53           1 :   explicit OrientationMap(
      54             :       std::array<Direction<VolumeDim>, VolumeDim> mapped_directions);
      55           0 :   OrientationMap(
      56             :       const std::array<Direction<VolumeDim>, VolumeDim>& directions_in_host,
      57             :       const std::array<Direction<VolumeDim>, VolumeDim>&
      58             :           directions_in_neighbor);
      59           0 :   ~OrientationMap() = default;
      60           0 :   OrientationMap(const OrientationMap&) = default;
      61           0 :   OrientationMap& operator=(const OrientationMap&) = default;
      62           0 :   OrientationMap(OrientationMap&& /*rhs*/) = default;
      63           0 :   OrientationMap& operator=(OrientationMap&& /*rhs*/) = default;
      64             : 
      65             :   /// Creates an OrientationMap that is the identity map on directions.
      66             :   /// `is_aligned()` is `true` in this case.
      67           1 :   static OrientationMap<VolumeDim> create_aligned();
      68             : 
      69             :   /// True when mapped(Direction) == Direction
      70           1 :   bool is_aligned() const {
      71             :     ASSERT(bit_field_ != static_cast<uint16_t>(0b1 << 15),
      72             :            "Cannot use a default-constructed OrientationMap");
      73             :     return (bit_field_ bitand aligned_mask) == aligned_mask;
      74             :   }
      75             : 
      76             :   /// The corresponding dimension in the neighbor.
      77           1 :   size_t operator()(const size_t dim) const {
      78             :     ASSERT(bit_field_ != static_cast<uint16_t>(0b1 << 15),
      79             :            "Cannot use a default-constructed OrientationMap");
      80             :     const auto neighbor_direction = get_direction(dim);
      81             :     ASSERT(neighbor_direction.side() != Side::Self,
      82             :            "There is no corresponding dimension");
      83             :     return neighbor_direction.dimension();
      84             :   }
      85             : 
      86             :   /// The corresponding direction in the neighbor.
      87           1 :   Direction<VolumeDim> operator()(const Direction<VolumeDim>& direction) const {
      88             :     ASSERT(bit_field_ != static_cast<uint16_t>(0b1 << 15),
      89             :            "Cannot use a default-constructed OrientationMap");
      90             :     return direction.side() == Side::Upper
      91             :                ? get_direction(direction.dimension())
      92             :                : get_direction(direction.dimension()).opposite();
      93             :   }
      94             : 
      95             :   /// The corresponding SegmentIds in the neighbor.
      96           1 :   std::array<SegmentId, VolumeDim> operator()(
      97             :       const std::array<SegmentId, VolumeDim>& segmentIds) const;
      98             : 
      99             :   /// The corresponding Mesh in the neighbor
     100           1 :   Mesh<VolumeDim> operator()(const Mesh<VolumeDim>& mesh) const;
     101             : 
     102             :   /// An array whose elements are permuted such that
     103             :   /// `result[this->operator()(d)] = array_to_permute[d]`.
     104             :   ///
     105             :   /// \note the permutation depends only on how the dimension is mapped
     106             :   /// and ignores the side of the mapped direction.
     107             :   template <typename T>
     108           1 :   std::array<T, VolumeDim> permute_to_neighbor(
     109             :       const std::array<T, VolumeDim>& array_to_permute) const;
     110             : 
     111             :   /// An array whose elements are permuted such that
     112             :   /// `result[d] = array_in_neighbor[this->operator()(d)]`
     113             :   ///
     114             :   /// \note the permutation depends only on how the dimension is mapped
     115             :   /// and ignores the side of the mapped direction.
     116             :   template <typename T>
     117           1 :   std::array<T, VolumeDim> permute_from_neighbor(
     118             :       const std::array<T, VolumeDim>& array_in_neighbor) const;
     119             : 
     120             :   /// The corresponding Orientation of the host in the frame of the neighbor.
     121           1 :   OrientationMap<VolumeDim> inverse_map() const;
     122             : 
     123             :   /// Serialization for Charm++
     124             :   // NOLINTNEXTLINE(google-runtime-references)
     125           1 :   void pup(PUP::er& p);
     126             : 
     127             :  private:
     128           0 :   friend bool operator==(const OrientationMap& lhs, const OrientationMap& rhs) {
     129             :     return lhs.bit_field_ == rhs.bit_field_;
     130             :   }
     131             : 
     132           0 :   Direction<VolumeDim> get_direction(size_t dim) const;
     133           0 :   void set_direction(size_t dim, const Direction<VolumeDim>& direction);
     134           0 :   void set_aligned(bool is_aligned);
     135           0 :   std::set<size_t> set_of_dimensions() const;
     136             : 
     137           0 :   uint16_t bit_field_{0b1 << 15};
     138             : };
     139             : 
     140             : /// Output operator for OrientationMap.
     141             : template <size_t VolumeDim>
     142           1 : std::ostream& operator<<(std::ostream& os,
     143             :                          const OrientationMap<VolumeDim>& orientation);
     144             : 
     145             : template <size_t VolumeDim>
     146           0 : bool operator!=(const OrientationMap<VolumeDim>& lhs,
     147             :                 const OrientationMap<VolumeDim>& rhs) {
     148             :   return not(lhs == rhs);
     149             : }
     150             : 
     151             : template <size_t VolumeDim>
     152             : template <typename T>
     153             : std::array<T, VolumeDim> OrientationMap<VolumeDim>::permute_to_neighbor(
     154             :     const std::array<T, VolumeDim>& array_to_permute) const {
     155             :   std::array<T, VolumeDim> array_in_neighbor = array_to_permute;
     156             :   if (is_aligned() or VolumeDim <= 1) {
     157             :     return array_in_neighbor;
     158             :   }
     159             :   for (size_t i = 0; i < VolumeDim; i++) {
     160             :     gsl::at(array_in_neighbor, this->operator()(i)) =
     161             :         gsl::at(array_to_permute, i);
     162             :   }
     163             :   return array_in_neighbor;
     164             : }
     165             : 
     166             : template <size_t VolumeDim>
     167             : template <typename T>
     168             : std::array<T, VolumeDim> OrientationMap<VolumeDim>::permute_from_neighbor(
     169             :     const std::array<T, VolumeDim>& array_in_neighbor) const {
     170             :   std::array<T, VolumeDim> result = array_in_neighbor;
     171             :   if (not is_aligned() and VolumeDim > 1) {
     172             :     for (size_t i = 0; i < VolumeDim; i++) {
     173             :       gsl::at(result, i) = gsl::at(array_in_neighbor, this->operator()(i));
     174             :     }
     175             :   }
     176             :   return result;
     177             : }
     178             : 
     179             : /// \ingroup ComputationalDomainGroup
     180             : /// `OrientationMap`s define an active rotation of the logical axes that bring
     181             : /// the axes of a host block into alignment with the logical axes of the
     182             : /// neighbor block. `discrete_rotation` applies this active rotation on the
     183             : /// coordinates as opposed to the axes.
     184             : /// For a two-dimensional example, consider a host block and a neighbor block,
     185             : /// where the OrientationMap between them is \f$\{-\eta,+\xi\}\f$. A quarter-
     186             : /// turn counterclockwise of the host block's logical axes would bring them into
     187             : /// alignment with those of the neighbor. That is, after this active rotation,
     188             : /// the blocks would be Aligned. Now consider a point A with coordinates
     189             : /// (+1.0,-0.5). An active quarter-turn rotation counter-clockwise about the
     190             : /// origin, keeping the axes fixed, brings point A into the coordinates
     191             : /// (+0.5,+1.0). This is how `discrete_rotation` interprets the
     192             : /// `OrientationMap` passed to it.
     193             : template <size_t VolumeDim, typename T>
     194           1 : std::array<tt::remove_cvref_wrap_t<T>, VolumeDim> discrete_rotation(
     195             :     const OrientationMap<VolumeDim>& rotation,
     196             :     std::array<T, VolumeDim> source_coords);
     197             : 
     198             : /*!
     199             :  * \ingroup ComputationalDomainGroup
     200             :  * \brief Computes the Jacobian of the transformation that is computed by
     201             :  * `discrete_rotation()`
     202             :  *
     203             :  * \note This always returns a `double` because the Jacobian is spatially
     204             :  * constant.
     205             :  */
     206             : template <size_t VolumeDim>
     207           1 : tnsr::Ij<double, VolumeDim, Frame::NoFrame> discrete_rotation_jacobian(
     208             :     const OrientationMap<VolumeDim>& orientation);
     209             : 
     210             : /*!
     211             :  * \ingroup ComputationalDomainGroup
     212             :  * \brief Computes the inverse Jacobian of the transformation that is computed
     213             :  * by `discrete_rotation()`
     214             :  */
     215             : template <size_t VolumeDim>
     216           1 : tnsr::Ij<double, VolumeDim, Frame::NoFrame> discrete_rotation_inverse_jacobian(
     217             :     const OrientationMap<VolumeDim>& orientation);