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