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>>;
|
