9 #include "DataStructures/ApplyMatrices.hpp"
10 #include "DataStructures/DataVector.hpp"
12 #include "DataStructures/SliceIterator.hpp"
24 template <
typename TagsList>
25 struct NumberOfIndependentComponents;
27 template <
typename... Tags>
28 struct NumberOfIndependentComponents<tmpl::list<Tags...>> {
29 static constexpr
size_t value = (... + Tags::type::size());
52 template <
typename VolumeVarsTagsList,
typename FaceVarsTagsList,
size_t Dim>
54 const gsl::not_null<Variables<FaceVarsTagsList>*> face_fields,
55 const Variables<VolumeVarsTagsList>& volume_fields,
57 static_assert(tmpl::size<VolumeVarsTagsList>::value != 0,
58 "Must have non-zero number of volume fields");
59 static_assert(tmpl::size<FaceVarsTagsList>::value >=
60 tmpl::size<VolumeVarsTagsList>::value,
61 "There must not be more volume tags than there are face tags.");
63 tmpl::list_contains_v<FaceVarsTagsList, tmpl::front<VolumeVarsTagsList>>,
64 "The first tag of VolumeVarsTagsList is not in the face tags. The "
65 "VolumeVarsTagsList must be a subset of the FaceVarsTagsList");
67 tmpl::list_contains_v<FaceVarsTagsList, tmpl::back<VolumeVarsTagsList>>,
68 "The last tag of VolumeVarsTagsList is not in the face tags. The "
69 "VolumeVarsTagsList must be a subset of the FaceVarsTagsList");
70 using face_vars_excluding_extras_at_end = tmpl::front<
71 tmpl::split_at<FaceVarsTagsList,
72 tmpl::next<tmpl::index_of<
73 FaceVarsTagsList, tmpl::back<VolumeVarsTagsList>>>>>;
74 using front_face_vars_split = tmpl::split_at<
75 face_vars_excluding_extras_at_end,
76 tmpl::index_of<FaceVarsTagsList, tmpl::front<VolumeVarsTagsList>>>;
77 using volume_vars_face_subset_list = tmpl::back<front_face_vars_split>;
79 std::is_same_v<volume_vars_face_subset_list, VolumeVarsTagsList>,
80 "The VolumeVarsTagsList must be a subset of the FaceVarsTagsList.");
81 constexpr
const size_t number_of_independent_components =
82 Variables<VolumeVarsTagsList>::number_of_independent_components;
83 using first_volume_tag = tmpl::front<VolumeVarsTagsList>;
84 const size_t sliced_dim = direction.dimension();
85 if (volume_mesh.quadrature(sliced_dim) == Spectral::Quadrature::Gauss) {
87 auto interpolation_matrices = make_array<Dim>(std::cref(
identity));
88 const auto& matrix = Spectral::boundary_interpolation_matrices(
89 volume_mesh.slice_through(sliced_dim));
90 gsl::at(interpolation_matrices, sliced_dim) =
91 direction.side() == Side::Upper ? matrix.second : matrix.first;
93 auto& first_face_field = get<first_volume_tag>(*face_fields);
94 auto& first_volume_field = get<first_volume_tag>(volume_fields);
101 first_face_field[0].data(),
102 first_face_field[0].size() * number_of_independent_components};
106 DataVector{
const_cast<double*
>(first_volume_field[0].data()),
107 first_volume_field[0].size() *
108 number_of_independent_components},
109 volume_mesh.extents());
111 const size_t fixed_index = direction.side() == Side::Upper
112 ? volume_mesh.extents(sliced_dim) - 1
115 const size_t interface_grid_points =
116 volume_mesh.extents().slice_away(sliced_dim).product();
117 const size_t volume_grid_points = volume_mesh.number_of_grid_points();
119 const double* vars_data = volume_fields.data();
122 double* interface_vars_data = get<first_volume_tag>(*face_fields)[0].data();
132 for (
SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
134 for (
size_t i = 0; i < number_of_independent_components; ++i) {
136 interface_vars_data[si.slice_offset() +
137 i * interface_grid_points] =
138 vars_data[si.volume_offset() + i * volume_grid_points];
151 template <
typename TagsToProjectList,
typename VolumeVarsTagsList,
152 typename FaceVarsTagsList,
size_t Dim>
154 const gsl::not_null<Variables<FaceVarsTagsList>*> face_fields,
155 const Variables<VolumeVarsTagsList>& volume_fields,
157 static_assert(tmpl::size<VolumeVarsTagsList>::value != 0,
158 "Must have non-zero number of volume fields");
159 static_assert(tmpl::size<FaceVarsTagsList>::value >=
160 tmpl::size<VolumeVarsTagsList>::value,
161 "There must not be more volume tags than there are face tags.");
164 tmpl::list_difference<TagsToProjectList, FaceVarsTagsList>>::value ==
166 "All of the tags in TagsToProjectList must be in FaceVarsTagsList");
168 tmpl::size<tmpl::list_difference<TagsToProjectList,
169 VolumeVarsTagsList>>::value == 0,
170 "All of the tags in TagsToProjectList must be in VolumeVarsTagsList");
171 const size_t sliced_dim = direction.dimension();
172 if (volume_mesh.quadrature(sliced_dim) == Spectral::Quadrature::Gauss) {
174 auto interpolation_matrices = make_array<Dim>(std::cref(
identity));
176 Spectral::boundary_interpolation_matrices(
177 volume_mesh.slice_through(sliced_dim));
178 gsl::at(interpolation_matrices, sliced_dim) =
179 direction.side() == Side::Upper ? matrices.second : matrices.first;
180 tmpl::for_each<TagsToProjectList>([&face_fields, &interpolation_matrices,
182 &volume_mesh](
auto tag_v) noexcept {
183 using tag =
typename decltype(tag_v)::type;
184 auto& face_field = get<tag>(*face_fields);
185 const auto& volume_field = get<tag>(volume_fields);
187 face_field[0].size() * face_field.size()};
190 DataVector{
const_cast<double*
>(volume_field[0].data()),
191 volume_field[0].size() * volume_field.size()},
192 volume_mesh.extents());
195 const size_t fixed_index = direction.side() == Side::Upper
196 ? volume_mesh.extents(sliced_dim) - 1
199 const size_t interface_grid_points =
200 volume_mesh.extents().slice_away(sliced_dim).product();
201 const size_t volume_grid_points = volume_mesh.number_of_grid_points();
205 for (
SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
207 tmpl::for_each<TagsToProjectList>([&face_fields, interface_grid_points,
210 auto tag_v) noexcept {
211 using tag =
typename decltype(tag_v)::type;
213 const double* vars_data = get<tag>(volume_fields)[0].data();
214 double* interface_vars_data = get<tag>(*face_fields)[0].data();
215 static constexpr
size_t number_of_independent_components_in_tensor =
218 for (
size_t i = 0; i < number_of_independent_components_in_tensor;
221 interface_vars_data[si.slice_offset() +
222 i * interface_grid_points] =
223 vars_data[si.volume_offset() + i * volume_grid_points];
235 template <
typename Symm,
typename IndexList,
size_t Dim>
237 const gsl::not_null<Tensor<DataVector, Symm, IndexList>*> face_field,
238 const Tensor<DataVector, Symm, IndexList>& volume_field,
240 const size_t sliced_dim = direction.dimension();
241 if (volume_mesh.quadrature(sliced_dim) ==
242 Spectral::Quadrature::Gauss) {
244 auto interpolation_matrices = make_array<Dim>(std::cref(
identity));
246 Spectral::boundary_interpolation_matrices(
247 volume_mesh.slice_through(sliced_dim));
248 gsl::at(interpolation_matrices, sliced_dim) =
249 direction.side() == Side::Upper ? matrices.second : matrices.first;
250 for (
size_t tensor_storage_index = 0;
251 tensor_storage_index < volume_field.size(); ++tensor_storage_index) {
253 interpolation_matrices, volume_field[tensor_storage_index],
254 volume_mesh.extents());
257 const size_t fixed_index = direction.side() == Side::Upper
258 ? volume_mesh.extents(sliced_dim) - 1
263 for (
SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
265 for (
size_t tensor_storage_index = 0;
266 tensor_storage_index < volume_field.size(); ++tensor_storage_index) {
267 (*face_field)[tensor_storage_index][si.slice_offset()] =
268 volume_field[tensor_storage_index][si.volume_offset()];