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