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