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 <limits>
       9             : #include <memory>
      10             : 
      11             : #include "DataStructures/Tensor/Tensor.hpp"
      12             : #include "Domain/Block.hpp"
      13             : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
      14             : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
      15             : #include "Domain/Structure/ElementId.hpp"
      16             : #include "Utilities/Gsl.hpp"
      17             : #include "Utilities/MakeArray.hpp"
      18             : 
      19             : namespace PUP {
      20             : class er;
      21             : }  // namespace PUP
      22             : 
      23             : /*!
      24             :  * \ingroup ComputationalDomainGroup
      25             :  * \brief The CoordinateMap for the Element from the Logical frame to the
      26             :  * `TargetFrame`
      27             :  *
      28             :  * An ElementMap takes a CoordinateMap for a Block and an ElementId as input,
      29             :  * and then "prepends" the correct affine map to the CoordinateMap so that the
      30             :  * map corresponds to the coordinate map for the Element rather than the Block.
      31             :  * This allows DomainCreators to only specify the maps for the Blocks without
      32             :  * worrying about how the domain may be decomposed beyond that.
      33             :  */
      34             : template <size_t Dim, typename TargetFrame>
      35           1 : class ElementMap {
      36             :  public:
      37           0 :   static constexpr size_t dim = Dim;
      38           0 :   using source_frame = Frame::ElementLogical;
      39           0 :   using target_frame = TargetFrame;
      40             : 
      41             :   /// \cond HIDDEN_SYMBOLS
      42             :   ElementMap() = default;
      43             :   /// \endcond
      44             : 
      45           0 :   ElementMap(ElementId<Dim> element_id,
      46             :              std::unique_ptr<domain::CoordinateMapBase<Frame::BlockLogical,
      47             :                                                        TargetFrame, Dim>>
      48             :                  block_map);
      49             : 
      50             :   /// Construct from an `element_id` within the `block`. The (affine)
      51             :   /// ElementLogical to BlockLogical map is determined by the `element_id`. The
      52             :   /// BlockLogical to TargetFrame map is determined by the `block`:
      53             :   /// - If the block is time-independent: the `block.stationary_map()` is used.
      54             :   /// - If the block is time-dependent: The `block.moving_mesh_*_map()` maps
      55             :   ///   are used. Which maps are used depends on the TargetFrame.
      56           1 :   ElementMap(ElementId<Dim> element_id, const Block<Dim>& block);
      57             : 
      58             :   const domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>&
      59           0 :   block_map() const {
      60             :     return *block_map_;
      61             :   }
      62             : 
      63           0 :   const ElementId<Dim>& element_id() const { return element_id_; }
      64             : 
      65             :   template <typename T>
      66           0 :   tnsr::I<T, Dim, TargetFrame> operator()(
      67             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
      68             :       const double time = std::numeric_limits<double>::signaling_NaN(),
      69             :       const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
      70             :     auto block_source_point =
      71             :         apply_affine_transformation_to_point(source_point);
      72             :     return block_map_->operator()(std::move(block_source_point), time,
      73             :                                   functions_of_time);
      74             :   }
      75             : 
      76             :   template <typename T>
      77           0 :   tnsr::I<T, Dim, Frame::ElementLogical> inverse(
      78             :       tnsr::I<T, Dim, TargetFrame> target_point,
      79             :       const double time = std::numeric_limits<double>::signaling_NaN(),
      80             :       const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
      81             :     auto block_source_point{
      82             :         block_map_->inverse(std::move(target_point), time, functions_of_time)
      83             :             .value()};
      84             :     // Apply the affine map to the points
      85             :     tnsr::I<T, Dim, Frame::ElementLogical> source_point;
      86             :     for (size_t d = 0; d < Dim; ++d) {
      87             :       source_point.get(d) =
      88             :           block_source_point.get(d) * gsl::at(map_inverse_slope_, d) +
      89             :           gsl::at(map_inverse_offset_, d);
      90             :     }
      91             :     return source_point;
      92             :   }
      93             : 
      94             :   template <typename T>
      95           0 :   InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jacobian(
      96             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
      97             :       const double time = std::numeric_limits<double>::signaling_NaN(),
      98             :       const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
      99             :     auto block_source_point =
     100             :         apply_affine_transformation_to_point(source_point);
     101             :     auto block_inv_jac = block_map_->inv_jacobian(std::move(block_source_point),
     102             :                                                   time, functions_of_time);
     103             :     InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jac;
     104             :     for (size_t d = 0; d < Dim; ++d) {
     105             :       for (size_t i = 0; i < Dim; ++i) {
     106             :         inv_jac.get(d, i) =
     107             :             block_inv_jac.get(d, i) * gsl::at(inverse_jacobian_, d);
     108             :       }
     109             :     }
     110             :     return inv_jac;
     111             :   }
     112             : 
     113             :   template <typename T>
     114           0 :   Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jacobian(
     115             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
     116             :       const double time = std::numeric_limits<double>::signaling_NaN(),
     117             :       const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
     118             :     auto block_source_point =
     119             :         apply_affine_transformation_to_point(source_point);
     120             :     auto block_jac = block_map_->jacobian(std::move(block_source_point), time,
     121             :                                           functions_of_time);
     122             :     Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jac;
     123             :     for (size_t d = 0; d < Dim; ++d) {
     124             :       for (size_t i = 0; i < Dim; ++i) {
     125             :         jac.get(i, d) = block_jac.get(i, d) * gsl::at(jacobian_, d);
     126             :       }
     127             :     }
     128             :     return jac;
     129             :   }
     130             : 
     131             :   // NOLINTNEXTLINE(google-runtime-references)
     132           0 :   void pup(PUP::er& p);
     133             : 
     134             :  private:
     135             :   template <typename T>
     136           0 :   tnsr::I<T, Dim, Frame::BlockLogical> apply_affine_transformation_to_point(
     137             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
     138             :     tnsr::I<T, Dim, Frame::BlockLogical> block_source_point;
     139             :     for (size_t d = 0; d < Dim; ++d) {
     140             :       block_source_point.get(d) = source_point.get(d) * gsl::at(map_slope_, d) +
     141             :                                   gsl::at(map_offset_, d);
     142             :     }
     143             :     return block_source_point;
     144             :   }
     145             : 
     146             :   std::unique_ptr<
     147             :       domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>>
     148           0 :       block_map_{nullptr};
     149           0 :   ElementId<Dim> element_id_{};
     150             :   // map_slope_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) -
     151             :   //                        segment_ids[i].endpoint(Side::Lower))
     152           0 :   std::array<double, Dim> map_slope_{
     153             :       make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
     154             :   // map_offset_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) +
     155             :   //                         segment_ids[i].endpoint(Side::Lower))
     156           0 :   std::array<double, Dim> map_offset_{
     157             :       make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
     158             :   // map_inverse_slope_[i] = 1.0 / map_slope_[i]
     159           0 :   std::array<double, Dim> map_inverse_slope_{
     160             :       make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
     161             :   // map_inverse_offset_[i] = -map_offset_[i] / map_slope_[i]
     162           0 :   std::array<double, Dim> map_inverse_offset_{
     163             :       make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
     164             :   // Note: The Jacobian is diagonal
     165           0 :   std::array<double, Dim> jacobian_{map_slope_};
     166           0 :   std::array<double, Dim> inverse_jacobian_{map_inverse_slope_};
     167             : };
     168             : 
     169             : template <size_t Dim, typename TargetFrame>
     170           0 : ElementMap(ElementId<Dim> element_id,
     171             :            std::unique_ptr<
     172             :                domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>>
     173             :                block_map) -> ElementMap<Dim, TargetFrame>;