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