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             : #ifdef SPECTRE_AUTODIFF
     130             :   /*!
     131             :    * \brief The inverse Hessian of the element map.
     132             :    *
     133             :    * Let \f$ \xi^i \f$ be the element logical coordinates, \f$ x^i \f$ be the
     134             :    * block logical coordinates, and \f$ y^i \f$ be the inertial coordinates.
     135             :    * The element map maps \f$ \xi^i \rightarrow x^i \rightarrow y^i \f$. Thus,
     136             :    * the inverse Hessian is
     137             :    * \f[
     138             :    * \frac{\partial^2\xi^i}{\partial y^j \partial y^k} =
     139             :    * \frac{\partial}{\partial y^j}\frac{\partial x^l}{\partial y^k}
     140             :    * \frac{\partial \xi^i}{\partial x^l}} =
     141             :    * \frac{\partial^2 x^l}{\partial y^j \partial y^k}
     142             :    * \frac{\partial \xi^i}{\partial x^l}
     143             :    * \f]
     144             :    * where we have used the fact
     145             :    * \f[
     146             :    * \frac{\partial^2 \xi^i}{\partial x^l \partial x^n} = 0.
     147             :    * \f]
     148             :    */
     149             :   template <typename T>
     150           1 :   InverseHessian<T, Dim, Frame::ElementLogical, TargetFrame> inv_hessian(
     151             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
     152             :       const double time = std::numeric_limits<double>::signaling_NaN(),
     153             :       const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
     154             :     auto block_source_point =
     155             :         apply_affine_transformation_to_point(source_point);
     156             :     auto block_inv_jac =
     157             :         block_map_->inv_jacobian(block_source_point, time, functions_of_time);
     158             :     auto block_inv_hes = block_map_->inv_hessian(
     159             :         std::move(block_source_point), block_inv_jac, time, functions_of_time);
     160             :     InverseHessian<T, Dim, Frame::ElementLogical, TargetFrame> inv_hes;
     161             :     for (size_t i = 0; i < Dim; ++i) {
     162             :       for (size_t j = 0; j < Dim; ++j) {
     163             :         for (size_t k = j; k < Dim; ++k) {
     164             :           inv_hes.get(i, j, k) =
     165             :               block_inv_hes.get(i, j, k) / gsl::at(map_slope_, i);
     166             :         }
     167             :       }
     168             :     }
     169             :     return inv_hes;
     170             :   }
     171             : #endif  // SPECTRE_AUTODIFF
     172             : 
     173             :   // NOLINTNEXTLINE(google-runtime-references)
     174           0 :   void pup(PUP::er& p);
     175             : 
     176             :  private:
     177             :   template <typename T>
     178           0 :   tnsr::I<T, Dim, Frame::BlockLogical> apply_affine_transformation_to_point(
     179             :       const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
     180             :     tnsr::I<T, Dim, Frame::BlockLogical> block_source_point;
     181             :     for (size_t d = 0; d < Dim; ++d) {
     182             :       block_source_point.get(d) = source_point.get(d) * gsl::at(map_slope_, d) +
     183             :                                   gsl::at(map_offset_, d);
     184             :     }
     185             :     return block_source_point;
     186             :   }
     187             : 
     188             :   std::unique_ptr<
     189             :       domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>>
     190           0 :       block_map_{nullptr};
     191             :   // map_slope_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) -
     192             :   //                        segment_ids[i].endpoint(Side::Lower))
     193             :   // Note: the map_slope_ is the jacobian_
     194           0 :   std::array<double, Dim> map_slope_{
     195             :       make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
     196             :   // map_offset_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) +
     197             :   //                         segment_ids[i].endpoint(Side::Lower))
     198           0 :   std::array<double, Dim> map_offset_{
     199             :       make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
     200             : };
     201             : 
     202             : template <size_t Dim, typename TargetFrame>
     203           0 : ElementMap(ElementId<Dim> element_id,
     204             :            std::unique_ptr<
     205             :                domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>>
     206             :                block_map) -> ElementMap<Dim, TargetFrame>;