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 <cstdint>
       9              : #include <optional>
      10              : #include <type_traits>
      11              : 
      12              : #include "DataStructures/DataBox/DataBox.hpp"
      13              : #include "DataStructures/DataVector.hpp"
      14              : #include "DataStructures/Tensor/Tensor.hpp"
      15              : #include "DataStructures/Transpose.hpp"
      16              : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      17              : #include "Evolution/DgSubcell/Tags/ReconstructionOrder.hpp"
      18              : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      19              : #include "Utilities/Gsl.hpp"
      20              : #include "Utilities/SetNumberOfGridPoints.hpp"
      21              : 
      22              : namespace evolution::dg::subcell {
      23              : /// Copies the `reconstruction_order` into the DataBox tag
      24              : /// `evolution::dg::subcell::Tags::ReconstructionOrder` if
      25              : /// `reconstruction_order` has a value.
      26              : template <size_t Dim, typename DbTagsList>
      27            1 : void store_reconstruction_order_in_databox(
      28              :     const gsl::not_null<db::DataBox<DbTagsList>*> box,
      29              :     const std::optional<std::array<gsl::span<std::uint8_t>, Dim>>&
      30              :         reconstruction_order) {
      31              :   if (reconstruction_order.has_value()) {
      32              :     db::mutate<evolution::dg::subcell::Tags::ReconstructionOrder<Dim>>(
      33              :         [&reconstruction_order](const auto recons_order_ptr,
      34              :                                 const ::Mesh<Dim>& subcell_mesh) {
      35              :           if (UNLIKELY(not recons_order_ptr->has_value())) {
      36              :             (*recons_order_ptr) = typename std::decay_t<
      37              :                 decltype(*recons_order_ptr)>::value_type{};
      38              :           }
      39              :           set_number_of_grid_points(
      40              :               make_not_null(&(recons_order_ptr->value())),
      41              :               subcell_mesh.number_of_grid_points());
      42              :           for (size_t d = Dim - 1; d < Dim; --d) {
      43              :             // The data is always in d-varies fastest ordering
      44              :             //
      45              :             // We always copy into the `d=0` buffer, then do transposes later.
      46              :             for (size_t j = 0;
      47              :                  j < subcell_mesh.slice_away(d).number_of_grid_points(); ++j) {
      48              :               for (size_t i = 0; i < subcell_mesh.extents(d); ++i) {
      49              :                 get<0>(recons_order_ptr
      50              :                            ->value())[j * subcell_mesh.extents(d) + i] =
      51              :                     gsl::at(reconstruction_order.value(),
      52              :                             d)[j * (subcell_mesh.extents(d) + 2) + i + 1];
      53              :               }
      54              :             }
      55              :             if (d == 1) {
      56              :               // Order is (y, z, x) -> (x, y, z)
      57              :               transpose(make_not_null(&recons_order_ptr->value().get(d)),
      58              :                         get<0>(recons_order_ptr->value()),
      59              :                         subcell_mesh.extents(1) *
      60              :                             (Dim > 2 ? subcell_mesh.extents(2) : 1),
      61              :                         subcell_mesh.number_of_grid_points() /
      62              :                             (subcell_mesh.extents(1) *
      63              :                              (Dim > 2 ? subcell_mesh.extents(2) : 1)));
      64              :             } else if (d == 2) {
      65              :               // Order is (z, x, y) -> (x, y, z)
      66              :               transpose(make_not_null(&recons_order_ptr->value().get(d)),
      67              :                         get<0>(recons_order_ptr->value()),
      68              :                         subcell_mesh.extents(2),
      69              :                         subcell_mesh.number_of_grid_points() /
      70              :                             subcell_mesh.extents(2));
      71              :             }
      72              :             // }
      73              :           }
      74              :         },
      75              :         box, db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(*box));
      76              :   }
      77              : }
      78              : }  // namespace evolution::dg::subcell