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