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(const 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(const 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             :   template <typename T>
      64           0 :   tnsr::I<T, Dim, TargetFrame> operator()(
      65             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
      66             :       const double time = std::numeric_limits<double>::signaling_NaN(),
      67             :       const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
      68             :     auto block_source_point =
      69             :         apply_affine_transformation_to_point(source_point);
      70             :     return block_map_->operator()(std::move(block_source_point), time,
      71             :                                   functions_of_time);
      72             :   }
      73             : 
      74             :   template <typename T>
      75           0 :   tnsr::I<T, Dim, Frame::ElementLogical> inverse(
      76             :       tnsr::I<T, Dim, TargetFrame> target_point,
      77             :       const double time = std::numeric_limits<double>::signaling_NaN(),
      78             :       const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
      79             :     auto block_source_point{
      80             :         block_map_->inverse(std::move(target_point), time, functions_of_time)
      81             :             .value()};
      82             :     // Apply the affine map to the points
      83             :     tnsr::I<T, Dim, Frame::ElementLogical> source_point;
      84             :     for (size_t d = 0; d < Dim; ++d) {
      85             :       const double inv_jac = 1.0 / gsl::at(map_slope_, d);
      86             :       const double temp = gsl::at(map_offset_, d) * inv_jac;
      87             :       source_point.get(d) = inv_jac * block_source_point.get(d) - temp;
      88             :     }
      89             :     return source_point;
      90             :   }
      91             : 
      92             :   template <typename T>
      93           0 :   InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jacobian(
      94             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
      95             :       const double time = std::numeric_limits<double>::signaling_NaN(),
      96             :       const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
      97             :     auto block_source_point =
      98             :         apply_affine_transformation_to_point(source_point);
      99             :     auto block_inv_jac = block_map_->inv_jacobian(std::move(block_source_point),
     100             :                                                   time, functions_of_time);
     101             :     InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jac;
     102             :     for (size_t d = 0; d < Dim; ++d) {
     103             :       const double inv_jac_in_dir = 1.0 / gsl::at(map_slope_, d);
     104             :       for (size_t i = 0; i < Dim; ++i) {
     105             :         inv_jac.get(d, i) = inv_jac_in_dir * block_inv_jac.get(d, i);
     106             :       }
     107             :     }
     108             :     return inv_jac;
     109             :   }
     110             : 
     111             :   template <typename T>
     112           0 :   Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jacobian(
     113             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
     114             :       const double time = std::numeric_limits<double>::signaling_NaN(),
     115             :       const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
     116             :     auto block_source_point =
     117             :         apply_affine_transformation_to_point(source_point);
     118             :     auto block_jac = block_map_->jacobian(std::move(block_source_point), time,
     119             :                                           functions_of_time);
     120             :     Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jac;
     121             :     for (size_t d = 0; d < Dim; ++d) {
     122             :       for (size_t i = 0; i < Dim; ++i) {
     123             :         jac.get(i, d) = block_jac.get(i, d) * gsl::at(map_slope_, d);
     124             :       }
     125             :     }
     126             :     return jac;
     127             :   }
     128             : 
     129             :   // NOLINTNEXTLINE(google-runtime-references)
     130           0 :   void pup(PUP::er& p);
     131             : 
     132             :  private:
     133             :   template <typename T>
     134           0 :   tnsr::I<T, Dim, Frame::BlockLogical> apply_affine_transformation_to_point(
     135             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
     136             :     tnsr::I<T, Dim, Frame::BlockLogical> block_source_point;
     137             :     for (size_t d = 0; d < Dim; ++d) {
     138             :       block_source_point.get(d) = source_point.get(d) * gsl::at(map_slope_, d) +
     139             :                                   gsl::at(map_offset_, d);
     140             :     }
     141             :     return block_source_point;
     142             :   }
     143             : 
     144             :   std::unique_ptr<
     145             :       domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>>
     146           0 :       block_map_{nullptr};
     147             :   // map_slope_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) -
     148             :   //                        segment_ids[i].endpoint(Side::Lower))
     149             :   // Note: the map_slope_ is the jacobian_
     150           0 :   std::array<double, Dim> map_slope_{
     151             :       make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
     152             :   // map_offset_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) +
     153             :   //                         segment_ids[i].endpoint(Side::Lower))
     154           0 :   std::array<double, Dim> map_offset_{
     155             :       make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
     156             : };
     157             : 
     158             : template <size_t Dim, typename TargetFrame>
     159           0 : ElementMap(ElementId<Dim> element_id,
     160             :            std::unique_ptr<
     161             :                domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>>
     162             :                block_map) -> ElementMap<Dim, TargetFrame>;