From c00cd7802b52445eed3a83c12ac3b9a5b5f37876 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <martin.kronbichler@rub.de>
Date: Thu, 28 Mar 2024 15:57:15 +0100
Subject: [PATCH] MappingQ: Do not assume contiguous storage of Tensor<2, dim>

---
 include/deal.II/fe/mapping_q_internal.h | 47 +++++--------------------
 1 file changed, 8 insertions(+), 39 deletions(-)

diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h
index 5bc48229fa..90ba787dbc 100644
--- a/include/deal.II/fe/mapping_q_internal.h
+++ b/include/deal.II/fe/mapping_q_internal.h
@@ -1233,40 +1233,6 @@ namespace internal
 
 
 
-    template <int dim, int spacedim>
-    inline DEAL_II_ALWAYS_INLINE void
-    store_vectorized_tensor(
-      const unsigned int n_points,
-      const unsigned int cur_index,
-      const DerivativeForm<1, dim, spacedim, VectorizedArray<double>>
-                                                    &derivative,
-      std::vector<DerivativeForm<1, dim, spacedim>> &result_array)
-    {
-      AssertDimension(result_array.size(), n_points);
-      constexpr unsigned int n_lanes = VectorizedArray<double>::size();
-      if (cur_index + n_lanes <= n_points)
-        {
-          std::array<unsigned int, n_lanes> indices;
-          for (unsigned int j = 0; j < n_lanes; ++j)
-            indices[j] = j * dim * spacedim;
-          const unsigned int even       = (dim * spacedim) / 4 * 4;
-          double            *result_ptr = &result_array[cur_index][0][0];
-          const VectorizedArray<double> *derivative_ptr = &derivative[0][0];
-          vectorized_transpose_and_store(
-            false, even, derivative_ptr, indices.data(), result_ptr);
-          for (unsigned int d = even; d < dim * spacedim; ++d)
-            for (unsigned int j = 0; j < n_lanes; ++j)
-              result_ptr[j * dim * spacedim + d] = derivative_ptr[d][j];
-        }
-      else
-        for (unsigned int j = 0; j < n_lanes && cur_index + j < n_points; ++j)
-          for (unsigned int d = 0; d < spacedim; ++d)
-            for (unsigned int e = 0; e < dim; ++e)
-              result_array[cur_index + j][d][e] = derivative[d][e][j];
-    }
-
-
-
     template <int dim, int spacedim>
     inline void
     maybe_update_q_points_Jacobians_generic(
@@ -1354,7 +1320,10 @@ namespace internal
               continue;
 
             if (update_flags & update_contravariant_transformation)
-              store_vectorized_tensor(n_points, i, derivative, jacobians);
+              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                for (unsigned int d = 0; d < spacedim; ++d)
+                  for (unsigned int e = 0; e < dim; ++e)
+                    jacobians[i + j][d][e] = derivative[d][e][j];
 
             if (update_flags & update_volume_elements)
               {
@@ -1366,10 +1335,10 @@ namespace internal
             if (update_flags & update_covariant_transformation)
               {
                 const auto covariant = derivative.covariant_form();
-                store_vectorized_tensor(n_points,
-                                        i,
-                                        covariant.transpose(),
-                                        inverse_jacobians);
+                for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                  for (unsigned int d = 0; d < dim; ++d)
+                    for (unsigned int e = 0; e < spacedim; ++e)
+                      inverse_jacobians[i + j][d][e] = covariant[e][d][j];
               }
           }
         else
-- 
2.39.5