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