Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <limits>
9 : #include <memory>
10 :
11 : #include "DataStructures/Tensor/Tensor.hpp"
12 : #include "Domain/Block.hpp"
13 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
14 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
15 : #include "Domain/Structure/ElementId.hpp"
16 : #include "Utilities/Gsl.hpp"
17 : #include "Utilities/MakeArray.hpp"
18 :
19 : namespace PUP {
20 : class er;
21 : } // namespace PUP
22 :
23 : /*!
24 : * \ingroup ComputationalDomainGroup
25 : * \brief The CoordinateMap for the Element from the Logical frame to the
26 : * `TargetFrame`
27 : *
28 : * An ElementMap takes a CoordinateMap for a Block and an ElementId as input,
29 : * and then "prepends" the correct affine map to the CoordinateMap so that the
30 : * map corresponds to the coordinate map for the Element rather than the Block.
31 : * This allows DomainCreators to only specify the maps for the Blocks without
32 : * worrying about how the domain may be decomposed beyond that.
33 : */
34 : template <size_t Dim, typename TargetFrame>
35 1 : class ElementMap {
36 : public:
37 0 : static constexpr size_t dim = Dim;
38 0 : using source_frame = Frame::ElementLogical;
39 0 : using target_frame = TargetFrame;
40 :
41 : /// \cond HIDDEN_SYMBOLS
42 : ElementMap() = default;
43 : /// \endcond
44 :
45 0 : ElementMap(const ElementId<Dim>& element_id,
46 : std::unique_ptr<domain::CoordinateMapBase<Frame::BlockLogical,
47 : TargetFrame, Dim>>
48 : block_map);
49 :
50 : /// Construct from an `element_id` within the `block`. The (affine)
51 : /// ElementLogical to BlockLogical map is determined by the `element_id`. The
52 : /// BlockLogical to TargetFrame map is determined by the `block`:
53 : /// - If the block is time-independent: the `block.stationary_map()` is used.
54 : /// - If the block is time-dependent: The `block.moving_mesh_*_map()` maps
55 : /// are used. Which maps are used depends on the TargetFrame.
56 1 : ElementMap(const ElementId<Dim>& element_id, const Block<Dim>& block);
57 :
58 : const domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>&
59 0 : block_map() const {
60 : return *block_map_;
61 : }
62 :
63 : template <typename T>
64 0 : tnsr::I<T, Dim, TargetFrame> operator()(
65 : const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
66 : const double time = std::numeric_limits<double>::signaling_NaN(),
67 : const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
68 : auto block_source_point =
69 : apply_affine_transformation_to_point(source_point);
70 : return block_map_->operator()(std::move(block_source_point), time,
71 : functions_of_time);
72 : }
73 :
74 : template <typename T>
75 0 : tnsr::I<T, Dim, Frame::ElementLogical> inverse(
76 : tnsr::I<T, Dim, TargetFrame> target_point,
77 : const double time = std::numeric_limits<double>::signaling_NaN(),
78 : const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
79 : auto block_source_point{
80 : block_map_->inverse(std::move(target_point), time, functions_of_time)
81 : .value()};
82 : // Apply the affine map to the points
83 : tnsr::I<T, Dim, Frame::ElementLogical> source_point;
84 : for (size_t d = 0; d < Dim; ++d) {
85 : const double inv_jac = 1.0 / gsl::at(map_slope_, d);
86 : const double temp = gsl::at(map_offset_, d) * inv_jac;
87 : source_point.get(d) = inv_jac * block_source_point.get(d) - temp;
88 : }
89 : return source_point;
90 : }
91 :
92 : template <typename T>
93 0 : InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jacobian(
94 : const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
95 : const double time = std::numeric_limits<double>::signaling_NaN(),
96 : const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
97 : auto block_source_point =
98 : apply_affine_transformation_to_point(source_point);
99 : auto block_inv_jac = block_map_->inv_jacobian(std::move(block_source_point),
100 : time, functions_of_time);
101 : InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jac;
102 : for (size_t d = 0; d < Dim; ++d) {
103 : const double inv_jac_in_dir = 1.0 / gsl::at(map_slope_, d);
104 : for (size_t i = 0; i < Dim; ++i) {
105 : inv_jac.get(d, i) = inv_jac_in_dir * block_inv_jac.get(d, i);
106 : }
107 : }
108 : return inv_jac;
109 : }
110 :
111 : template <typename T>
112 0 : Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jacobian(
113 : const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
114 : const double time = std::numeric_limits<double>::signaling_NaN(),
115 : const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
116 : auto block_source_point =
117 : apply_affine_transformation_to_point(source_point);
118 : auto block_jac = block_map_->jacobian(std::move(block_source_point), time,
119 : functions_of_time);
120 : Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jac;
121 : for (size_t d = 0; d < Dim; ++d) {
122 : for (size_t i = 0; i < Dim; ++i) {
123 : jac.get(i, d) = block_jac.get(i, d) * gsl::at(map_slope_, d);
124 : }
125 : }
126 : return jac;
127 : }
128 :
129 : #ifdef SPECTRE_AUTODIFF
130 : /*!
131 : * \brief The inverse Hessian of the element map.
132 : *
133 : * Let \f$ \xi^i \f$ be the element logical coordinates, \f$ x^i \f$ be the
134 : * block logical coordinates, and \f$ y^i \f$ be the inertial coordinates.
135 : * The element map maps \f$ \xi^i \rightarrow x^i \rightarrow y^i \f$. Thus,
136 : * the inverse Hessian is
137 : * \f[
138 : * \frac{\partial^2\xi^i}{\partial y^j \partial y^k} =
139 : * \frac{\partial}{\partial y^j}\frac{\partial x^l}{\partial y^k}
140 : * \frac{\partial \xi^i}{\partial x^l}} =
141 : * \frac{\partial^2 x^l}{\partial y^j \partial y^k}
142 : * \frac{\partial \xi^i}{\partial x^l}
143 : * \f]
144 : * where we have used the fact
145 : * \f[
146 : * \frac{\partial^2 \xi^i}{\partial x^l \partial x^n} = 0.
147 : * \f]
148 : */
149 : template <typename T>
150 1 : InverseHessian<T, Dim, Frame::ElementLogical, TargetFrame> inv_hessian(
151 : const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
152 : const double time = std::numeric_limits<double>::signaling_NaN(),
153 : const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
154 : auto block_source_point =
155 : apply_affine_transformation_to_point(source_point);
156 : auto block_inv_jac =
157 : block_map_->inv_jacobian(block_source_point, time, functions_of_time);
158 : auto block_inv_hes = block_map_->inv_hessian(
159 : std::move(block_source_point), block_inv_jac, time, functions_of_time);
160 : InverseHessian<T, Dim, Frame::ElementLogical, TargetFrame> inv_hes;
161 : for (size_t i = 0; i < Dim; ++i) {
162 : for (size_t j = 0; j < Dim; ++j) {
163 : for (size_t k = j; k < Dim; ++k) {
164 : inv_hes.get(i, j, k) =
165 : block_inv_hes.get(i, j, k) / gsl::at(map_slope_, i);
166 : }
167 : }
168 : }
169 : return inv_hes;
170 : }
171 : #endif // SPECTRE_AUTODIFF
172 :
173 : // NOLINTNEXTLINE(google-runtime-references)
174 0 : void pup(PUP::er& p);
175 :
176 : private:
177 : template <typename T>
178 0 : tnsr::I<T, Dim, Frame::BlockLogical> apply_affine_transformation_to_point(
179 : const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
180 : tnsr::I<T, Dim, Frame::BlockLogical> block_source_point;
181 : for (size_t d = 0; d < Dim; ++d) {
182 : block_source_point.get(d) = source_point.get(d) * gsl::at(map_slope_, d) +
183 : gsl::at(map_offset_, d);
184 : }
185 : return block_source_point;
186 : }
187 :
188 : std::unique_ptr<
189 : domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>>
190 0 : block_map_{nullptr};
191 : // map_slope_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) -
192 : // segment_ids[i].endpoint(Side::Lower))
193 : // Note: the map_slope_ is the jacobian_
194 0 : std::array<double, Dim> map_slope_{
195 : make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
196 : // map_offset_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) +
197 : // segment_ids[i].endpoint(Side::Lower))
198 0 : std::array<double, Dim> map_offset_{
199 : make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
200 : };
201 :
202 : template <size_t Dim, typename TargetFrame>
203 0 : ElementMap(ElementId<Dim> element_id,
204 : std::unique_ptr<
205 : domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>>
206 : block_map) -> ElementMap<Dim, TargetFrame>;
|
