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