Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <optional>
       8             : #include <unordered_map>
       9             : #include <vector>
      10             : 
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : 
      13             : /// \cond
      14             : namespace domain {
      15             : class BlockId;
      16             : }  // namespace domain
      17             : class DataVector;
      18             : template <size_t VolumeDim>
      19             : class ElementId;
      20             : template <typename IdType, typename DataType>
      21             : class IdPair;
      22             : /// \endcond
      23             : 
      24             : /*!
      25             :  * \brief Map block logical coordinates to element logical coordinates
      26             :  *
      27             :  * If the point is outside the element, returns std::nullopt. Otherwise, returns
      28             :  * the element logical coordinates of the point.
      29             :  *
      30             :  * Points on the element boundary are considered to be in the element and will
      31             :  * have element logical coordinates of -1 or 1. See the other function overload
      32             :  * for handling multiple points and disambiguating points on shared element
      33             :  * boundaries.
      34             :  */
      35             : template <size_t Dim>
      36             : std::optional<tnsr::I<double, Dim, Frame::ElementLogical>>
      37           1 : element_logical_coordinates(
      38             :     const tnsr::I<double, Dim, Frame::BlockLogical>& x_block_logical,
      39             :     const ElementId<Dim>& element_id);
      40             : 
      41             : /// \ingroup ComputationalDomainGroup
      42             : ///
      43             : /// Holds element logical coordinates of an arbitrary set of points on
      44             : /// a single `Element`.  The arbitrary set of points is assumed to be
      45             : /// a subset of a larger set of points spanning multiple `Element`s,
      46             : /// and this class holds `offsets` that index into that larger set of
      47             : /// points.
      48             : ///
      49             : /// \details `offsets.size()` is the same as the size of the `DataVector`
      50             : /// inside `element_logical_coords`.
      51             : ///
      52             : /// This is used during the process of interpolating volume quantities
      53             : /// on the `Elements` (e.g. the spatial metric) onto an arbitrary set
      54             : /// of points (e.g. the points on an apparent horizon or a
      55             : /// wave-extraction surface) expressed in some frame.  Here is an
      56             : /// outline of how this interpolation proceeds, and where
      57             : /// `element_logical_coordinates` and `block_logical_coordinates` fit
      58             : /// into the picture:
      59             : ///
      60             : /// Assume some component (e.g. HorizonA) has a `Tensor<DataVector>`
      61             : /// of target coordinate points in some coordinate frame.  The goal is
      62             : /// to determine the `Element` and logical coordinates of each point,
      63             : /// have each `Element` interpolate volume data onto the points
      64             : /// contained inside that `Element`, and send the interpolated data
      65             : /// back to the component.  The first step of this process is to
      66             : /// determine the block_id and block_logical_coordinates of each
      67             : /// point; this is done by the component (e.g. HorizonA), which calls
      68             : /// the function `block_logical_coordinates` on its full set of target
      69             : /// points.  The result of `block_logical_coordinates` is then
      70             : /// communicated to the members of a NodeGroup component
      71             : /// (e.g. HorizonManager).  Each node of the NodeGroup then calls
      72             : /// `element_logical_coordinates`, which returns a map of `ElementId`
      73             : /// to `ElementLogicalCoordHolder` for all the `Element`s on that node
      74             : /// that contain one or more of the target points. The NodeGroup
      75             : /// (which already has received the volume data from the `Elements` on
      76             : /// that node), interpolates the volume data to the element logical
      77             : /// coordinates for all of these `ElementId`s.  The `offsets` in the
      78             : /// `ElementLogicalCoordHolder` are the indices into the `DataVectors`
      79             : /// of the original target coordinates and will be used to assemble
      80             : /// the interpolated data into `Tensor<DataVector>`s that have the
      81             : /// same ordering as the original target coordinates. The NodeGroups
      82             : /// perform a reduction to get the data back to the original
      83             : /// component.
      84             : template <size_t Dim>
      85           1 : struct ElementLogicalCoordHolder {
      86           0 :   tnsr::I<DataVector, Dim, Frame::ElementLogical> element_logical_coords;
      87           0 :   std::vector<size_t> offsets;
      88             : };
      89             : 
      90             : /// \ingroup ComputationalDomainGroup
      91             : ///
      92             : /// Given a set of points in block logical coordinates and their
      93             : /// `BlockIds`, as returned from the function
      94             : /// `block_logical_coordinates`, determines which `Element`s in a list
      95             : /// of `ElementId`s contains each point, and determines the element
      96             : /// logical coordinates of each point.
      97             : ///
      98             : /// \details Returns a std::unordered_map from `ElementId`s to
      99             : /// `ElementLogicalCoordHolder`s.
     100             : /// It is expected that only a subset of the points will be found
     101             : /// in the given `Element`s.
     102             : /// Boundary points: If a point is on the boundary of an Element, it is
     103             : /// considered contained in that Element only if it is on the lower bound
     104             : /// of the Element, or if it is on the upper bound of the element and that
     105             : /// upper bound coincides with the upper bound of the containing Block.
     106             : /// This means that each boundary point is contained in one and only one
     107             : /// Element.  We assume that the input block_coord_holders associates
     108             : /// a point on a Block boundary with only a single Block, the one with
     109             : /// the smaller BlockId, which is always the lower-bounding Block.
     110             : ///
     111             : /// \code
     112             : ///  <---    Block 0   ---> <---   Block 1   --->
     113             : ///  |          |          |          |          |
     114             : /// P_0   E0   P_1   E1   P_2   E2   P_3   E3   P_4
     115             : ///  |          |          |          |          |
     116             : ///
     117             : /// For example, the above 1D diagram shows four Elements labeled E0
     118             : /// through E3, and five boundary points labeled P_0 through P_4 (where
     119             : /// P_0 and P_4 are external boundaries).  There are two Blocks.  This
     120             : /// algorithm assigns each boundary point to one and only one Element as
     121             : /// follows:
     122             : /// P_0 -> E0
     123             : /// P_1 -> E1
     124             : /// P_2 -> E1 (Note: block_coord_holders includes P_2 only in Block 0)
     125             : /// P_3 -> E3
     126             : /// P_4 -> E3
     127             : /// \endcode
     128             : template <size_t Dim>
     129           1 : auto element_logical_coordinates(
     130             :     const std::vector<ElementId<Dim>>& element_ids,
     131             :     const std::vector<std::optional<IdPair<
     132             :         domain::BlockId, tnsr::I<double, Dim, typename Frame::BlockLogical>>>>&
     133             :         block_coord_holders)
     134             :     -> std::unordered_map<ElementId<Dim>, ElementLogicalCoordHolder<Dim>>;