]> https://gitweb.dealii.org/ - dealii.git/commitdiff
generalize the portable MF methods to fe_degree < n_q_points_1d cases
authorYimin Jin <jinyimin14@mails.ucas.ac.cn>
Wed, 12 Feb 2025 08:24:21 +0000 (16:24 +0800)
committerYimin Jin <jinyimin14@mails.ucas.ac.cn>
Sat, 1 Mar 2025 15:44:22 +0000 (23:44 +0800)
include/deal.II/matrix_free/portable_evaluation_kernels.h [new file with mode: 0644]
include/deal.II/matrix_free/portable_fe_evaluation.h
include/deal.II/matrix_free/portable_matrix_free.h
include/deal.II/matrix_free/portable_matrix_free.templates.h
include/deal.II/matrix_free/portable_tensor_product_kernels.h
tests/matrix_free_kokkos/evaluate_1d_shape.cc
tests/matrix_free_kokkos/evaluate_1d_shape.output
tests/matrix_free_kokkos/evaluate_2d_shape.cc
tests/matrix_free_kokkos/evaluate_2d_shape.output
tests/tests.h

diff --git a/include/deal.II/matrix_free/portable_evaluation_kernels.h b/include/deal.II/matrix_free/portable_evaluation_kernels.h
new file mode 100644 (file)
index 0000000..42a477b
--- /dev/null
@@ -0,0 +1,620 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2017 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+#ifndef dealii__evaluation_kernels_h
+#define dealii__evaluation_kernels_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/portable_matrix_free.h>
+#include <deal.II/matrix_free/portable_tensor_product_kernels.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace Portable
+{
+  namespace internal
+  {
+    /**
+     * Helper function to specify whether a transformation to collocation should
+     * be used: It should give correct results (first condition), we need to be
+     * able to initialize the fields in shape_info.templates.h from the
+     * polynomials (second condition), and it should be the most efficient
+     * choice in terms of operation counts (third condition).
+     */
+    constexpr bool
+    use_collocation_evaluation(const unsigned int fe_degree,
+                               const unsigned int n_q_points_1d)
+    {
+      // TODO: are the conditions suit for GPU parallelization?
+      return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
+             (n_q_points_1d <= 3 * fe_degree / 2 + 1);
+    }
+
+
+
+    /**
+     * This struct performs the evaluation of function values and gradients for
+     * tensor-product finite elements. There are two specialized implementation
+     * classes FEEvaluationImplCollocation (for Gauss-Lobatto elements where the
+     * nodal points and the quadrature points coincide and the 'values'
+     * operation is identity) and FEEvaluationImplTransformToCollocation (which
+     * can be transformed to a collocation space and can then use the identity
+     * in these spaces), which both allow for shorter code.
+     */
+    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+    struct FEEvaluationImpl
+    {
+      using TeamHandle = Kokkos::TeamPolicy<
+        MemorySpace::Default::kokkos_space::execution_space>::member_type;
+      using SharedView = Kokkos::View<Number *,
+                                      MemorySpace::Default::kokkos_space::
+                                        execution_space::scratch_memory_space,
+                                      Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
+
+      DEAL_II_HOST_DEVICE static void
+      evaluate(const unsigned int                            n_components,
+               const EvaluationFlags::EvaluationFlags        evaluation_flag,
+               const typename MatrixFree<dim, Number>::Data *data,
+               SharedData<dim, Number>                      *shared_data)
+      {
+        if (evaluation_flag == EvaluationFlags::nothing)
+          return;
+
+        // the evaluator does not need temporary storage since no in-place
+        // operation takes place in this function
+        auto scratch_for_eval =
+          Kokkos::subview(shared_data->scratch_pad, Kokkos::make_pair(0, 0));
+        EvaluatorTensorProduct<EvaluatorVariant::evaluate_general,
+                               dim,
+                               fe_degree + 1,
+                               n_q_points_1d,
+                               Number>
+          eval(shared_data->team_member,
+               data->shape_values,
+               data->shape_gradients,
+               data->co_shape_gradients,
+               scratch_for_eval);
+
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            auto u      = Kokkos::subview(shared_data->values, Kokkos::ALL, c);
+            auto grad_u = Kokkos::subview(shared_data->gradients,
+                                          Kokkos::ALL,
+                                          Kokkos::ALL,
+                                          c);
+
+            if constexpr (dim == 1)
+              {
+                auto temp =
+                  Kokkos::subview(shared_data->scratch_pad,
+                                  Kokkos::make_pair(0, n_q_points_1d));
+
+                if (evaluation_flag & EvaluationFlags::gradients)
+                  eval.template gradients<0, true, false, false>(
+                    u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
+                if (evaluation_flag & EvaluationFlags::values)
+                  {
+                    eval.template values<0, true, false, false>(u, temp);
+                    populate_view<false>(shared_data->team_member,
+                                         u,
+                                         temp,
+                                         n_q_points_1d);
+                  }
+              }
+            else if constexpr (dim == 2)
+              {
+                constexpr int temp_size = (fe_degree + 1) * n_q_points_1d;
+                auto          temp = Kokkos::subview(shared_data->scratch_pad,
+                                            Kokkos::make_pair(0, temp_size));
+
+                // grad x
+                if (evaluation_flag & EvaluationFlags::gradients)
+                  {
+                    eval.template gradients<0, true, false, false>(u, temp);
+                    eval.template values<1, true, false, false>(
+                      temp, Kokkos::subview(grad_u, Kokkos::ALL, 0));
+                  }
+
+                // grad y
+                eval.template values<0, true, false, false>(u, temp);
+                if (evaluation_flag & EvaluationFlags::gradients)
+                  eval.template gradients<1, true, false, false>(
+                    temp, Kokkos::subview(grad_u, Kokkos::ALL, 1));
+
+                // val: can use values applied in x
+                if (evaluation_flag & EvaluationFlags::values)
+                  eval.template values<1, true, false, false>(temp, u);
+              }
+            else if constexpr (dim == 3)
+              {
+                constexpr int temp1_size = Utilities::pow(fe_degree + 1, 2) *
+                                           n_q_points_1d,
+                              temp2_size = Utilities::pow(n_q_points_1d, 2) *
+                                           (fe_degree + 1);
+
+                auto temp1 = Kokkos::subview(shared_data->scratch_pad,
+                                             Kokkos::make_pair(0, temp1_size));
+                auto temp2 =
+                  Kokkos::subview(shared_data->scratch_pad,
+                                  Kokkos::make_pair(temp1_size,
+                                                    temp1_size + temp2_size));
+
+                if (evaluation_flag & EvaluationFlags::gradients)
+                  {
+                    // grad x
+                    eval.template gradients<0, true, false, false>(u, temp1);
+                    eval.template values<1, true, false, false>(temp1, temp2);
+                    eval.template values<2, true, false, false>(
+                      temp2, Kokkos::subview(grad_u, Kokkos::ALL, 0));
+                  }
+
+                // grad y
+                eval.template values<0, true, false, false>(u, temp1);
+                if (evaluation_flag & EvaluationFlags::gradients)
+                  {
+                    eval.template gradients<1, true, false, false>(temp1,
+                                                                   temp2);
+                    eval.template values<2, true, false, false>(
+                      temp2, Kokkos::subview(grad_u, Kokkos::ALL, 1));
+                  }
+
+                // grad z: can use the values applied in x direction stored
+                // in temp1
+                eval.template values<1, true, false, false>(temp1, temp2);
+                if (evaluation_flag & EvaluationFlags::gradients)
+                  eval.template gradients<2, true, false, false>(
+                    temp2, Kokkos::subview(grad_u, Kokkos::ALL, 2));
+
+                // val: can use the values applied in x & y direction stored
+                // in temp2
+                if (evaluation_flag & EvaluationFlags::values)
+                  eval.template values<2, true, false, false>(temp2, u);
+              }
+            else
+              Assert(false, ExcMessage("dim must not exceed 3!"));
+          }
+      }
+
+
+
+      DEAL_II_HOST_DEVICE static void
+      integrate(const unsigned int                            n_components,
+                const EvaluationFlags::EvaluationFlags        integration_flag,
+                const typename MatrixFree<dim, Number>::Data *data,
+                SharedData<dim, Number>                      *shared_data)
+      {
+        if (integration_flag == EvaluationFlags::nothing)
+          return;
+
+        // the evaluator does not need temporary storage since no in-place
+        // operation takes place in this function
+        auto scratch_for_eval =
+          Kokkos::subview(shared_data->scratch_pad, Kokkos::make_pair(0, 0));
+        EvaluatorTensorProduct<EvaluatorVariant::evaluate_general,
+                               dim,
+                               fe_degree + 1,
+                               n_q_points_1d,
+                               Number>
+          eval(shared_data->team_member,
+               data->shape_values,
+               data->shape_gradients,
+               data->co_shape_gradients,
+               scratch_for_eval);
+
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            auto u      = Kokkos::subview(shared_data->values, Kokkos::ALL, c);
+            auto grad_u = Kokkos::subview(shared_data->gradients,
+                                          Kokkos::ALL,
+                                          Kokkos::ALL,
+                                          c);
+
+            if constexpr (dim == 1)
+              {
+                auto temp =
+                  Kokkos::subview(shared_data->scratch_pad,
+                                  Kokkos::make_pair(0, fe_degree + 1));
+
+                if ((integration_flag & EvaluationFlags::values) &&
+                    !(integration_flag & EvaluationFlags::gradients))
+                  {
+                    eval.template values<0, false, false, false>(u, temp);
+                    populate_view<false>(shared_data->team_member,
+                                         u,
+                                         temp,
+                                         fe_degree + 1);
+                  }
+                if (integration_flag & EvaluationFlags::gradients)
+                  {
+                    if (integration_flag & EvaluationFlags::values)
+                      {
+                        eval.template values<0, false, false, false>(u, temp);
+                        eval.template gradients<0, false, true, false>(
+                          Kokkos::subview(grad_u, Kokkos::ALL, 0), temp);
+                        populate_view<false>(shared_data->team_member,
+                                             u,
+                                             temp,
+                                             fe_degree + 1);
+                      }
+                    else
+                      eval.template gradients<0, false, false, false>(
+                        Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
+                  }
+              }
+            else if constexpr (dim == 2)
+              {
+                constexpr int temp_size = (fe_degree + 1) * n_q_points_1d;
+                auto          temp = Kokkos::subview(shared_data->scratch_pad,
+                                            Kokkos::make_pair(0, temp_size));
+
+                if ((integration_flag & EvaluationFlags::values) &&
+                    !(integration_flag & EvaluationFlags::gradients))
+                  {
+                    eval.template values<1, false, false, false>(u, temp);
+                    eval.template values<0, false, false, false>(temp, u);
+                  }
+                if (integration_flag & EvaluationFlags::gradients)
+                  {
+                    eval.template gradients<1, false, false, false>(
+                      Kokkos::subview(grad_u, Kokkos::ALL, 1), temp);
+                    if (integration_flag & EvaluationFlags::values)
+                      eval.template values<1, false, true, false>(u, temp);
+                    eval.template values<0, false, false, false>(temp, u);
+                    eval.template values<1, false, false, false>(
+                      Kokkos::subview(grad_u, Kokkos::ALL, 0), temp);
+                    eval.template gradients<0, false, true, false>(temp, u);
+                  }
+              }
+            else if constexpr (dim == 3)
+              {
+                constexpr int temp1_size = Utilities::pow(n_q_points_1d, 2) *
+                                           (fe_degree + 1),
+                              temp2_size = Utilities::pow(fe_degree + 1, 2) *
+                                           n_q_points_1d;
+
+                auto temp1 = Kokkos::subview(shared_data->scratch_pad,
+                                             Kokkos::make_pair(0, temp1_size));
+                auto temp2 =
+                  Kokkos::subview(shared_data->scratch_pad,
+                                  Kokkos::make_pair(temp1_size,
+                                                    temp1_size + temp2_size));
+
+                if ((integration_flag & EvaluationFlags::values) &&
+                    !(integration_flag & EvaluationFlags::gradients))
+                  {
+                    eval.template values<2, false, false, false>(u, temp1);
+                    eval.template values<1, false, false, false>(temp1, temp2);
+                    eval.template values<0, false, false, false>(temp2, u);
+                  }
+                if (integration_flag & EvaluationFlags::gradients)
+                  {
+                    eval.template gradients<2, false, false, false>(
+                      Kokkos::subview(grad_u, Kokkos::ALL, 2), temp1);
+                    if (integration_flag & EvaluationFlags::values)
+                      eval.template values<2, false, true, false>(u, temp1);
+                    eval.template values<1, false, false, false>(temp1, temp2);
+                    eval.template values<2, false, false, false>(
+                      Kokkos::subview(grad_u, Kokkos::ALL, 1), temp1);
+                    eval.template gradients<1, false, true, false>(temp1,
+                                                                   temp2);
+                    eval.template values<0, false, false, false>(temp2, u);
+                    eval.template values<2, false, false, false>(
+                      Kokkos::subview(grad_u, Kokkos::ALL, 0), temp1);
+                    eval.template values<1, false, false, false>(temp1, temp2);
+                    eval.template gradients<0, false, true, false>(temp2, u);
+                  }
+              }
+            else
+              Assert(false, ExcMessage("dim must not exceed 3!"));
+          }
+      }
+    };
+
+
+
+    /**
+     * This struct performs the evaluation of function values and gradients for
+     * tensor-product finite elements. This is a specialization for elements
+     * where the nodal points coincide with the quadrature points like FE_Q
+     * shape functions on Gauss-Lobatto elements integrated with Gauss-Lobatto
+     * quadrature. The assumption of this class is that the shape 'values'
+     * operation is identity, which allows us to write shorter code.
+     *
+     * In literature, this form of evaluation is often called spectral
+     * evaluation, spectral collocation or simply collocation, meaning the same
+     * location for shape functions and evaluation space (quadrature points).
+     */
+    template <int dim, int fe_degree, typename Number>
+    struct FEEvaluationImplCollocation
+    {
+      DEAL_II_HOST_DEVICE static void
+      evaluate(const unsigned int                            n_components,
+               const EvaluationFlags::EvaluationFlags        evaluation_flag,
+               const typename MatrixFree<dim, Number>::Data *data,
+               SharedData<dim, Number>                      *shared_data)
+      {
+        // since the dof values have already been stored in
+        // shared_data->values, there is nothing to do if the gradients are
+        // not required
+        if (!(evaluation_flag & EvaluationFlags::gradients))
+          return;
+
+        constexpr int n_points = Utilities::pow(fe_degree + 1, dim);
+        auto scratch_for_eval  = Kokkos::subview(shared_data->scratch_pad,
+                                                Kokkos::make_pair(0, n_points));
+
+        EvaluatorTensorProduct<EvaluatorVariant::evaluate_general,
+                               dim,
+                               fe_degree + 1,
+                               fe_degree + 1,
+                               Number>
+          eval(shared_data->team_member,
+               data->shape_values,
+               data->shape_gradients,
+               data->co_shape_gradients,
+               scratch_for_eval);
+
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            auto u      = Kokkos::subview(shared_data->values, Kokkos::ALL, c);
+            auto grad_u = Kokkos::subview(shared_data->gradients,
+                                          Kokkos::ALL,
+                                          Kokkos::ALL,
+                                          c);
+
+            eval.template co_gradients<0, true, false, false>(
+              u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
+            if constexpr (dim > 1)
+              eval.template co_gradients<1, true, false, false>(
+                u, Kokkos::subview(grad_u, Kokkos::ALL, 1));
+            if constexpr (dim > 2)
+              eval.template co_gradients<2, true, false, false>(
+                u, Kokkos::subview(grad_u, Kokkos::ALL, 2));
+          }
+      }
+
+
+      DEAL_II_HOST_DEVICE static void
+      integrate(const unsigned int                            n_components,
+                const EvaluationFlags::EvaluationFlags        integration_flag,
+                const typename MatrixFree<dim, Number>::Data *data,
+                SharedData<dim, Number>                      *shared_data)
+      {
+        // since the quad values have already been stored in
+        // shared_data->values, there is nothing to do if the gradients are
+        // not required
+        if (!(integration_flag & EvaluationFlags::gradients))
+          return;
+
+        constexpr int n_points = Utilities::pow(fe_degree + 1, dim);
+        auto scratch_for_eval  = Kokkos::subview(shared_data->scratch_pad,
+                                                Kokkos::make_pair(0, n_points));
+
+        EvaluatorTensorProduct<EvaluatorVariant::evaluate_general,
+                               dim,
+                               fe_degree + 1,
+                               fe_degree + 1,
+                               Number>
+          eval(shared_data->team_member,
+               data->shape_values,
+               data->shape_gradients,
+               data->co_shape_gradients,
+               scratch_for_eval);
+
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            auto u      = Kokkos::subview(shared_data->values, Kokkos::ALL, c);
+            auto grad_u = Kokkos::subview(shared_data->gradients,
+                                          Kokkos::ALL,
+                                          Kokkos::ALL,
+                                          c);
+
+            if constexpr (dim == 1)
+              {
+                if (integration_flag & EvaluationFlags::values)
+                  eval.template co_gradients<0, false, true, false>(
+                    Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
+                else
+                  eval.template co_gradients<2, false, false, false>(
+                    Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
+              }
+            else if constexpr (dim == 2)
+              {
+                if (integration_flag & EvaluationFlags::values)
+                  eval.template co_gradients<1, false, true, false>(
+                    Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
+                else
+                  eval.template co_gradients<1, false, false, false>(
+                    Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
+                eval.template co_gradients<0, false, true, false>(
+                  Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
+              }
+            else if constexpr (dim == 3)
+              {
+                if (integration_flag & EvaluationFlags::values)
+                  eval.template co_gradients<2, false, true, false>(
+                    Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
+                else
+                  eval.template co_gradients<2, false, false, false>(
+                    Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
+                eval.template co_gradients<1, false, true, false>(
+                  Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
+                eval.template co_gradients<0, false, true, false>(
+                  Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
+              }
+            else
+              Assert(false, ExcMessage("dim must not exceed 3!"));
+          }
+      }
+    };
+
+
+
+    /**
+     * This struct performs the evaluation of function values and gradients for
+     * tensor-product finite elements. This is a specialization for symmetric
+     * basis functions about the mid point 0.5 of the unit interval with the
+     * same number of quadrature points as degrees of freedom. In that case, we
+     * can first transform the basis to one that has the nodal points in the
+     * quadrature points (i.e., the collocation space) and then perform the
+     * evaluation of the first and second derivatives in this transformed space,
+     * using the identity operation for the shape values.
+     */
+    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+    struct FEEvaluationImplTransformToCollocation
+    {
+      DEAL_II_HOST_DEVICE static void
+      evaluate(const unsigned int                            n_components,
+               const EvaluationFlags::EvaluationFlags        evaluation_flag,
+               const typename MatrixFree<dim, Number>::Data *data,
+               SharedData<dim, Number>                      *shared_data)
+      {
+        constexpr int scratch_size = Utilities::pow(n_q_points_1d, dim);
+        auto          scratch_for_eval =
+          Kokkos::subview(shared_data->scratch_pad,
+                          Kokkos::make_pair(0, scratch_size));
+
+        EvaluatorTensorProduct<EvaluatorVariant::evaluate_general,
+                               dim,
+                               fe_degree + 1,
+                               n_q_points_1d,
+                               Number>
+          eval(shared_data->team_member,
+               data->shape_values,
+               data->shape_gradients,
+               data->co_shape_gradients,
+               scratch_for_eval);
+
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            auto u      = Kokkos::subview(shared_data->values, Kokkos::ALL, c);
+            auto grad_u = Kokkos::subview(shared_data->gradients,
+                                          Kokkos::ALL,
+                                          Kokkos::ALL,
+                                          c);
+
+            eval.template values<0, true, false, true>(u, u);
+            if constexpr (dim > 1)
+              eval.template values<1, true, false, true>(u, u);
+            if constexpr (dim > 2)
+              eval.template values<2, true, false, true>(u, u);
+
+            if (evaluation_flag & EvaluationFlags::gradients)
+              {
+                eval.template co_gradients<0, true, false, false>(
+                  u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
+                if constexpr (dim > 1)
+                  eval.template co_gradients<1, true, false, false>(
+                    u, Kokkos::subview(grad_u, Kokkos::ALL, 1));
+                if constexpr (dim > 2)
+                  eval.template co_gradients<2, true, false, false>(
+                    u, Kokkos::subview(grad_u, Kokkos::ALL, 2));
+              }
+          }
+      }
+
+
+      DEAL_II_HOST_DEVICE static void
+      integrate(const unsigned int                            n_components,
+                const EvaluationFlags::EvaluationFlags        integration_flag,
+                const typename MatrixFree<dim, Number>::Data *data,
+                const SharedData<dim, Number>                *shared_data)
+      {
+        constexpr int scratch_size = Utilities::pow(n_q_points_1d, dim);
+        auto          scratch_for_eval =
+          Kokkos::subview(shared_data->scratch_pad,
+                          Kokkos::make_pair(0, scratch_size));
+
+        EvaluatorTensorProduct<EvaluatorVariant::evaluate_general,
+                               dim,
+                               fe_degree + 1,
+                               n_q_points_1d,
+                               Number>
+          eval(shared_data->team_member,
+               data->shape_values,
+               data->shape_gradients,
+               data->co_shape_gradients,
+               scratch_for_eval);
+
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            auto u      = Kokkos::subview(shared_data->values, Kokkos::ALL, c);
+            auto grad_u = Kokkos::subview(shared_data->gradients,
+                                          Kokkos::ALL,
+                                          Kokkos::ALL,
+                                          c);
+
+            // apply derivatives in collocation space
+            if (integration_flag & EvaluationFlags::gradients)
+              {
+                if constexpr (dim == 1)
+                  {
+                    if (integration_flag & EvaluationFlags::values)
+                      eval.template co_gradients<0, false, true, false>(
+                        Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
+                    else
+                      eval.template co_gradients<2, false, false, false>(
+                        Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
+                  }
+                else if constexpr (dim == 2)
+                  {
+                    if (integration_flag & EvaluationFlags::values)
+                      eval.template co_gradients<1, false, true, false>(
+                        Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
+                    else
+                      eval.template co_gradients<1, false, false, false>(
+                        Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
+                    eval.template co_gradients<0, false, true, false>(
+                      Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
+                  }
+                else if constexpr (dim == 3)
+                  {
+                    if (integration_flag & EvaluationFlags::values)
+                      eval.template co_gradients<2, false, true, false>(
+                        Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
+                    else
+                      eval.template co_gradients<2, false, false, false>(
+                        Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
+                    eval.template co_gradients<1, false, true, false>(
+                      Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
+                    eval.template co_gradients<0, false, true, false>(
+                      Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
+                  }
+                else
+                  Assert(false, ExcMessage("dim must not exceed 3!"));
+              }
+
+            // transform back to the original space
+            if constexpr (dim > 2)
+              eval.template values<2, false, false, true>(u, u);
+            if constexpr (dim > 1)
+              eval.template values<1, false, false, true>(u, u);
+            eval.template values<0, false, false, true>(u, u);
+          }
+      }
+    };
+  } // end of namespace internal
+} // end of namespace Portable
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index dfa5a46ff6480a3fd4605591377b9cb5772c51c4..ae23f70eb8078b45a26a57dee67b43fd222a5368 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/portable_evaluation_kernels.h>
 #include <deal.II/matrix_free/portable_hanging_nodes_internal.h>
 #include <deal.II/matrix_free/portable_matrix_free.h>
 #include <deal.II/matrix_free/portable_matrix_free.templates.h>
@@ -106,10 +107,18 @@ namespace Portable
       Utilities::pow(n_q_points_1d, dim);
 
     /**
-     * Number of tensor degrees of freedoms per cell.
+     * Number of tensor degrees of freedom of a scalar component determined
+     * from the given template argument `fe_degree`.
+     */
+    static constexpr unsigned int tensor_dofs_per_component =
+      Utilities::pow(fe_degree + 1, dim);
+
+    /**
+     * Number of tensor degrees of freedom of all component determined from
+     * the given template argument `fe_degree`.
      */
     static constexpr unsigned int tensor_dofs_per_cell =
-      Utilities::pow(fe_degree + 1, dim) * n_components;
+      tensor_dofs_per_component * n_components;
 
     /**
      * Constructor.
@@ -271,13 +280,14 @@ namespace Portable
     read_dof_values(const Number *src)
   {
     // Populate the scratch memory
-    Kokkos::parallel_for(
-      Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
-      [&](const int &i) {
-        for (unsigned int c = 0; c < n_components_; ++c)
-          shared_data->values(i, c) = src[data->local_to_global(
-            cell_id, i + (tensor_dofs_per_cell / n_components) * c)];
-      });
+    Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member,
+                                                 tensor_dofs_per_component),
+                         [&](const int &i) {
+                           for (unsigned int c = 0; c < n_components_; ++c)
+                             shared_data->values(i, c) =
+                               src[data->local_to_global(
+                                 cell_id, i + tensor_dofs_per_component * c)];
+                         });
     shared_data->team_member.team_barrier();
 
     for (unsigned int c = 0; c < n_components_; ++c)
@@ -313,24 +323,25 @@ namespace Portable
     if (data->use_coloring)
       {
         Kokkos::parallel_for(
-          Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
+          Kokkos::TeamThreadRange(shared_data->team_member,
+                                  tensor_dofs_per_component),
           [&](const int &i) {
             for (unsigned int c = 0; c < n_components_; ++c)
-              dst[data->local_to_global(
-                cell_id, i + (tensor_dofs_per_cell / n_components) * c)] +=
+              dst[data->local_to_global(cell_id,
+                                        i + tensor_dofs_per_component * c)] +=
                 shared_data->values(i, c);
           });
       }
     else
       {
         Kokkos::parallel_for(
-          Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
+          Kokkos::TeamThreadRange(shared_data->team_member,
+                                  tensor_dofs_per_component),
           [&](const int &i) {
             for (unsigned int c = 0; c < n_components_; ++c)
-              Kokkos::atomic_add(
-                &dst[data->local_to_global(
-                  cell_id, i + (tensor_dofs_per_cell / n_components) * c)],
-                shared_data->values(i, c));
+              Kokkos::atomic_add(&dst[data->local_to_global(
+                                   cell_id, i + (tensor_dofs_per_component)*c)],
+                                 shared_data->values(i, c));
           });
       }
   }
@@ -340,50 +351,42 @@ namespace Portable
   template <int dim,
             int fe_degree,
             int n_q_points_1d,
-            int n_components_,
+            int n_components,
             typename Number>
   DEAL_II_HOST_DEVICE void
-  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::evaluate(
-    const EvaluationFlags::EvaluationFlags evaluate_flag)
+  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>::evaluate(
+    const EvaluationFlags::EvaluationFlags evaluation_flag)
   {
-    // First evaluate the gradients because it requires values that will be
-    // changed if evaluate_val is true
-    internal::EvaluatorTensorProduct<
-      internal::EvaluatorVariant::evaluate_general,
-      dim,
-      fe_degree,
-      n_q_points_1d,
-      Number>
-      evaluator_tensor_product(shared_data->team_member,
-                               data->shape_values,
-                               data->shape_gradients,
-                               data->co_shape_gradients);
+    using ElementType = ::dealii::internal::MatrixFreeFunctions::ElementType;
 
-    for (unsigned int c = 0; c < n_components_; ++c)
+    if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
+        data->element_type == ElementType::tensor_symmetric_collocation)
       {
-        if ((evaluate_flag & EvaluationFlags::values) &&
-            (evaluate_flag & EvaluationFlags::gradients))
-          {
-            evaluator_tensor_product.evaluate_values_and_gradients(
-              Kokkos::subview(shared_data->values, Kokkos::ALL, c),
-              Kokkos::subview(
-                shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
-            shared_data->team_member.team_barrier();
-          }
-        else if (evaluate_flag & EvaluationFlags::gradients)
-          {
-            evaluator_tensor_product.evaluate_gradients(
-              Kokkos::subview(shared_data->values, Kokkos::ALL, c),
-              Kokkos::subview(
-                shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
-            shared_data->team_member.team_barrier();
-          }
-        else if (evaluate_flag & EvaluationFlags::values)
-          {
-            evaluator_tensor_product.evaluate_values(
-              Kokkos::subview(shared_data->values, Kokkos::ALL, c));
-            shared_data->team_member.team_barrier();
-          }
+        internal::FEEvaluationImplCollocation<dim, fe_degree, Number>::evaluate(
+          n_components, evaluation_flag, data, shared_data);
+      }
+    // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
+    // shape_info.h for more details
+    else if (fe_degree >= 0 &&
+             internal::use_collocation_evaluation(fe_degree, n_q_points_1d) &&
+             data->element_type <= ElementType::tensor_symmetric)
+      {
+        internal::FEEvaluationImplTransformToCollocation<
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          Number>::evaluate(n_components, evaluation_flag, data, shared_data);
+      }
+    else if (fe_degree >= 0 &&
+             data->element_type <= ElementType::tensor_symmetric_no_collocation)
+      {
+        internal::FEEvaluationImpl<dim, fe_degree, n_q_points_1d, Number>::
+          evaluate(n_components, evaluation_flag, data, shared_data);
+      }
+    else
+      {
+        Kokkos::abort("The element type is not yet supported by the portable "
+                      "matrix-free module.");
       }
   }
 
@@ -415,42 +418,36 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::integrate(
     const EvaluationFlags::EvaluationFlags integration_flag)
   {
-    internal::EvaluatorTensorProduct<
-      internal::EvaluatorVariant::evaluate_general,
-      dim,
-      fe_degree,
-      n_q_points_1d,
-      Number>
-      evaluator_tensor_product(shared_data->team_member,
-                               data->shape_values,
-                               data->shape_gradients,
-                               data->co_shape_gradients);
-
+    using ElementType = ::dealii::internal::MatrixFreeFunctions::ElementType;
 
-    for (unsigned int c = 0; c < n_components_; ++c)
+    if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
+        data->element_type == ElementType::tensor_symmetric_collocation)
       {
-        if ((integration_flag & EvaluationFlags::values) &&
-            (integration_flag & EvaluationFlags::gradients))
-          {
-            evaluator_tensor_product.integrate_values_and_gradients(
-              Kokkos::subview(shared_data->values, Kokkos::ALL, c),
-              Kokkos::subview(
-                shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
-          }
-        else if (integration_flag & EvaluationFlags::values)
-          {
-            evaluator_tensor_product.integrate_values(
-              Kokkos::subview(shared_data->values, Kokkos::ALL, c));
-            shared_data->team_member.team_barrier();
-          }
-        else if (integration_flag & EvaluationFlags::gradients)
-          {
-            evaluator_tensor_product.template integrate_gradients<false>(
-              Kokkos::subview(shared_data->values, Kokkos::ALL, c),
-              Kokkos::subview(
-                shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
-            shared_data->team_member.team_barrier();
-          }
+        internal::FEEvaluationImplCollocation<dim, fe_degree, Number>::
+          integrate(n_components, integration_flag, data, shared_data);
+      }
+    // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
+    // shape_info.h for more details
+    else if (fe_degree >= 0 &&
+             internal::use_collocation_evaluation(fe_degree, n_q_points_1d) &&
+             data->element_type <= ElementType::tensor_symmetric)
+      {
+        internal::FEEvaluationImplTransformToCollocation<
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          Number>::integrate(n_components, integration_flag, data, shared_data);
+      }
+    else if (fe_degree >= 0 &&
+             data->element_type <= ElementType::tensor_symmetric_no_collocation)
+      {
+        internal::FEEvaluationImpl<dim, fe_degree, n_q_points_1d, Number>::
+          integrate(n_components, integration_flag, data, shared_data);
+      }
+    else
+      {
+        Kokkos::abort("The element type is not yet supported by the portable "
+                      "matrix-free module.");
       }
   }
 
@@ -486,6 +483,7 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::get_value(
     int q_point) const
   {
+    Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
     if constexpr (n_components_ == 1)
       {
         return shared_data->values(q_point, 0);
@@ -512,17 +510,18 @@ namespace Portable
                                             n_components_,
                                             Number>::value_type
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
-    get_dof_value(int q_point) const
+    get_dof_value(int dof) const
   {
+    Assert(dof >= 0 && dof < tensor_dofs_per_component, ExcInternalError());
     if constexpr (n_components_ == 1)
       {
-        return shared_data->values(q_point, 0);
+        return shared_data->values(dof, 0);
       }
     else
       {
         value_type result;
         for (unsigned int c = 0; c < n_components; ++c)
-          result[c] = shared_data->values(q_point, c);
+          result[c] = shared_data->values(dof, c);
         return result;
       }
   }
@@ -538,6 +537,7 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     submit_value(const value_type &val_in, int q_point)
   {
+    Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
     if constexpr (n_components_ == 1)
       {
         shared_data->values(q_point, 0) = val_in * data->JxW(cell_id, q_point);
@@ -559,16 +559,17 @@ namespace Portable
             typename Number>
   DEAL_II_HOST_DEVICE void
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
-    submit_dof_value(const value_type &val_in, int q_point)
+    submit_dof_value(const value_type &val_in, int dof)
   {
+    Assert(dof >= 0 && dof < tensor_dofs_per_component, ExcInternalError());
     if constexpr (n_components_ == 1)
       {
-        shared_data->values(q_point, 0) = val_in;
+        shared_data->values(dof, 0) = val_in;
       }
     else
       {
         for (unsigned int c = 0; c < n_components; ++c)
-          shared_data->values(q_point, c) = val_in[c];
+          shared_data->values(dof, c) = val_in[c];
       }
   }
 
@@ -587,6 +588,7 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     get_gradient(int q_point) const
   {
+    Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
     gradient_type grad;
 
     if constexpr (n_components_ == 1)
@@ -627,6 +629,7 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     submit_gradient(const gradient_type &grad_in, int q_point)
   {
+    Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
     if constexpr (n_components_ == 1)
       {
         for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
index dc1b9c22cf75d939696b2877967fa69ea545671b..dee00888322b7999d1eddd40a5f77bf06b12ff51 100644 (file)
@@ -35,6 +35,8 @@
 #include <deal.II/lac/affine_constraints.h>
 #include <deal.II/lac/la_parallel_vector.h>
 
+#include <deal.II/matrix_free/shape_info.h>
+
 #include <Kokkos_Core.hpp>
 
 
@@ -234,6 +236,18 @@ namespace Portable
        */
       bool use_coloring;
 
+      /**
+       * Encodes the type of element detected at construction. FEEvaluation
+       * will select the most efficient algorithm based on the given element
+       * type.
+       */
+      ::dealii::internal::MatrixFreeFunctions::ElementType element_type;
+
+      /**
+       * Size of the scratch pad for temporary storage in shared memory.
+       */
+      unsigned int scratch_pad_size;
+
       /**
        * Return the quadrature point index local. The index is
        * only unique for a given MPI process.
@@ -474,6 +488,18 @@ namespace Portable
      */
     types::global_dof_index n_dofs;
 
+    /**
+     * Encodes the type of element detected at construction. FEEvaluation
+     * will select the most efficient algorithm based on the given element
+     * type.
+     */
+    ::dealii::internal::MatrixFreeFunctions::ElementType element_type;
+
+    /**
+     * Size of the scratch pad for temporary storage in shared memory.
+     */
+    unsigned int scratch_pad_size;
+
     /**
      * Degree of the finite element used.
      */
@@ -627,14 +653,20 @@ namespace Portable
       Number ***,
       MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
       Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
+    using SharedViewScratchPad = Kokkos::View<
+      Number *,
+      MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
+      Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
 
     DEAL_II_HOST_DEVICE
-    SharedData(const TeamHandle          &team_member,
-               const SharedViewValues    &values,
-               const SharedViewGradients &gradients)
+    SharedData(const TeamHandle           &team_member,
+               const SharedViewValues     &values,
+               const SharedViewGradients  &gradients,
+               const SharedViewScratchPad &scratch_pad)
       : team_member(team_member)
       , values(values)
       , gradients(gradients)
+      , scratch_pad(scratch_pad)
     {}
 
     /**
@@ -651,6 +683,11 @@ namespace Portable
      * Memory for computed gradients in reference coordinate system.
      */
     SharedViewGradients gradients;
+
+    /**
+     * Memory for temporary arrays required by evaluation and integration.
+     */
+    SharedViewScratchPad scratch_pad;
   };
 
 
index e8c760923bba468ad85b539a08fce66b1cc1caad..3903861aa768013680a04bfc755ad0ee52b9f2fe 100644 (file)
@@ -345,6 +345,11 @@ namespace Portable
                      MemorySpace::Default::kokkos_space::execution_space::
                        scratch_memory_space,
                      Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
+      using SharedViewScratchPad =
+        Kokkos::View<Number *,
+                     MemorySpace::Default::kokkos_space::execution_space::
+                       scratch_memory_space,
+                     Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
 
       ApplyKernel(Functor                                      func,
                   const typename MatrixFree<dim, Number>::Data gpu_data,
@@ -367,13 +372,12 @@ namespace Portable
       size_t
       team_shmem_size(int /*team_size*/) const
       {
-        return SharedViewValues::shmem_size(Functor::n_local_dofs /
-                                              gpu_data.n_components,
+        return SharedViewValues::shmem_size(Functor::n_q_points,
                                             gpu_data.n_components) +
-               SharedViewGradients::shmem_size(Functor::n_local_dofs /
-                                                 gpu_data.n_components,
+               SharedViewGradients::shmem_size(Functor::n_q_points,
                                                dim,
-                                               gpu_data.n_components);
+                                               gpu_data.n_components) +
+               SharedViewScratchPad::shmem_size(gpu_data.scratch_pad_size);
       }
 
 
@@ -382,16 +386,20 @@ namespace Portable
       operator()(const TeamHandle &team_member) const
       {
         // Get the scratch memory
-        SharedViewValues    values(team_member.team_shmem(),
-                                Functor::n_local_dofs / gpu_data.n_components,
+        SharedViewValues     values(team_member.team_shmem(),
+                                Functor::n_q_points,
                                 gpu_data.n_components);
-        SharedViewGradients gradients(team_member.team_shmem(),
-                                      Functor::n_local_dofs /
-                                        gpu_data.n_components,
+        SharedViewGradients  gradients(team_member.team_shmem(),
+                                      Functor::n_q_points,
                                       dim,
                                       gpu_data.n_components);
+        SharedViewScratchPad scratch_pad(team_member.team_shmem(),
+                                         gpu_data.scratch_pad_size);
 
-        SharedData<dim, Number> shared_data(team_member, values, gradients);
+        SharedData<dim, Number> shared_data(team_member,
+                                            values,
+                                            gradients,
+                                            scratch_pad);
         func(team_member.league_rank(), &gpu_data, &shared_data, src, dst);
       }
     };
@@ -500,6 +508,8 @@ namespace Portable
     data_copy.padding_length     = padding_length;
     data_copy.row_start          = row_start[color];
     data_copy.use_coloring       = use_coloring;
+    data_copy.element_type       = element_type;
+    data_copy.scratch_pad_size   = scratch_pad_size;
 
     return data_copy;
   }
@@ -660,6 +670,49 @@ namespace Portable
 
 
 
+  namespace
+  {
+    /**
+     * Helper function for determining the scratch pad size.
+     */
+    inline unsigned int
+    compute_scratch_pad_size(
+      const ::dealii::internal::MatrixFreeFunctions::ElementType element_type,
+      const int                                                  dim,
+      const int                                                  fe_degree,
+      const int                                                  n_q_points_1d)
+    {
+      using ElementType = ::dealii::internal::MatrixFreeFunctions::ElementType;
+
+      if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric)
+        {
+          // evaluate/integrate with FEEvaluationImplCollocation or
+          // FEEvaluationImplTransformToCollocation
+          return Utilities::pow(n_q_points_1d, dim);
+        }
+      else if (fe_degree >= 0 &&
+               element_type <= ElementType::tensor_symmetric_no_collocation)
+        {
+          // evaluate/integrate with FEEvaluationImpl
+          if (dim == 1)
+            return n_q_points_1d;
+          else if (dim == 2)
+            return (fe_degree + 1) * n_q_points_1d;
+          else if (dim == 3)
+            return (fe_degree + 1) * n_q_points_1d *
+                   (fe_degree + 1 + n_q_points_1d);
+          else
+            AssertThrow(false, ExcNotImplemented());
+        }
+      else
+        AssertThrow(false, ExcNotImplemented());
+
+      return numbers::invalid_unsigned_int;
+    }
+  } // namespace
+
+
+
   template <int dim, typename Number>
   template <typename IteratorFiltersType>
   void
@@ -691,13 +744,15 @@ namespace Portable
     const unsigned int n_dofs_1d     = fe_degree + 1;
     const unsigned int n_q_points_1d = quad.size();
 
-    Assert(n_dofs_1d == n_q_points_1d,
-           ExcMessage("n_q_points_1d must be equal to fe_degree+1."));
+    // TODO remove the limitation in the future
+    AssertThrow(n_dofs_1d <= n_q_points_1d,
+                ExcMessage("n_q_points_1d must be greater than or equal to "
+                           "fe_degree + 1."));
 
     // Set padding length to the closest power of two larger than or equal to
     // the number of threads.
-    padding_length = 1 << static_cast<unsigned int>(
-                       std::ceil(dim * std::log2(fe_degree + 1.)));
+    padding_length = 1 << static_cast<unsigned int>(std::ceil(
+                       dim * std::log2(static_cast<float>(n_q_points_1d))));
 
     dofs_per_cell        = fe.n_dofs_per_cell();
     n_components         = fe.n_components();
@@ -707,6 +762,12 @@ namespace Portable
     ::dealii::internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info(quad,
                                                                           fe);
 
+    this->element_type     = shape_info.element_type;
+    this->scratch_pad_size = compute_scratch_pad_size(shape_info.element_type,
+                                                      dim,
+                                                      fe_degree,
+                                                      n_q_points_1d);
+
     unsigned int size_shape_values = n_dofs_1d * n_q_points_1d;
 
     shape_values = Kokkos::View<Number *, MemorySpace::Default::kokkos_space>(
index 6e29a3a5932855eba50939f550bdacd3cb2cad42..e168512b18e7fb57f0af30c75e0eeae82e2ff59a 100644 (file)
@@ -42,16 +42,46 @@ namespace Portable
 
 
 
+    /**
+     * Helper function that copies or adds the first N entries of src to
+     * dst, depending on the template argument "add".
+     */
+    template <bool add, typename ViewTypeIn, typename ViewTypeOut>
+    DEAL_II_HOST_DEVICE void
+    populate_view(
+      const Kokkos::TeamPolicy<
+        MemorySpace::Default::kokkos_space::execution_space>::member_type
+                      &team_member,
+      ViewTypeOut      dst,
+      const ViewTypeIn src,
+      const int        N)
+    {
+      Assert(dst.size() >= N, ExcInternalError());
+      Assert(src.size() >= N, ExcInternalError());
+
+      Kokkos::parallel_for(Kokkos::TeamVectorRange(team_member, N),
+                           [&](const int i) {
+                             if constexpr (add)
+                               Kokkos::atomic_add(&dst(i), src(i));
+                             else
+                               dst(i) = src(i);
+                           });
+
+      team_member.team_barrier();
+    }
+
+
+
 #if KOKKOS_VERSION >= 40000
     /**
-     * Helper function for values() and gradients() in 1D
+     * Helper function for apply() in 1D
      */
-    template <int n_q_points_1d,
+    template <int n_rows,
+              int n_columns,
+              int direction,
               typename Number,
-              int  direction,
-              bool dof_to_quad,
+              bool contract_over_rows,
               bool add,
-              bool in_place,
               typename ViewTypeIn,
               typename ViewTypeOut>
     DEAL_II_HOST_DEVICE void
@@ -63,49 +93,44 @@ namespace Portable
              const ViewTypeIn in,
              ViewTypeOut      out)
     {
-      Number t[n_q_points_1d];
-      Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, n_q_points_1d),
-                           [&](const int &q) {
-                             t[q] = 0;
-                             // This loop simply multiplies the shape function
-                             // at the quadrature point by the value finite
-                             // element coefficient.
-                             // FIXME check why using parallel_reduce
-                             // ThreadVector is slower
-                             for (int k = 0; k < n_q_points_1d; ++k)
+      constexpr int Nk = (contract_over_rows ? n_rows : n_columns),
+                    Nq = (contract_over_rows ? n_columns : n_rows);
+
+      Assert(shape_data.size() == n_rows * n_columns, ExcInternalError());
+      Assert(in.size() >= Nk, ExcInternalError());
+      Assert(out.size() >= Nq, ExcInternalError());
+
+      Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, Nq),
+                           [&](const int q) {
+                             Number sum = 0;
+                             for (int k = 0; k < Nk; ++k)
                                {
-                                 const unsigned int shape_idx =
-                                   dof_to_quad ? (q + k * n_q_points_1d) :
-                                                 (k + q * n_q_points_1d);
-                                 const unsigned int source_idx = k;
-                                 t[q] += shape_data[shape_idx] * in(source_idx);
+                                 const int shape_idx =
+                                   (contract_over_rows ? q + k * Nq :
+                                                         k + q * Nk);
+                                 sum += shape_data(shape_idx) * in(k);
                                }
-                           });
-
-      if constexpr (in_place)
-        team_member.team_barrier();
 
-      Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, n_q_points_1d),
-                           [&](const int &q) {
-                             const unsigned int destination_idx = q;
                              if constexpr (add)
-                               Kokkos::atomic_add(&out(destination_idx), t[q]);
+                               Kokkos::atomic_add(&out(q), sum);
                              else
-                               out(destination_idx) = t[q];
+                               out(q) = sum;
                            });
+
+      team_member.team_barrier();
     }
 
 
 
     /**
-     * Helper function for values() and gradients() in 2D
+     * Helper function for apply() in 2D
      */
-    template <int n_q_points_1d,
+    template <int n_rows,
+              int n_columns,
+              int direction,
               typename Number,
-              int  direction,
-              bool dof_to_quad,
+              bool contract_over_rows,
               bool add,
-              bool in_place,
               typename ViewTypeIn,
               typename ViewTypeOut>
     DEAL_II_HOST_DEVICE void
@@ -119,58 +144,70 @@ namespace Portable
     {
       using TeamType = Kokkos::TeamPolicy<
         MemorySpace::Default::kokkos_space::execution_space>::member_type;
-      constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 2);
 
-      Number t[n_q_points];
-      auto   thread_policy =
+      // Sizes of the input and output vectors:
+      // -----------------------------------------------------------
+      //   direction  |  contract_over_rows  |  !contract_over_rows
+      // -----------------------------------------------------------
+      //       0      |    m x m -> n x m    |     n x m -> m x m
+      // -----------------------------------------------------------
+      //       1      |    n x m -> n x n    |     n x n -> n x m
+      // -----------------------------------------------------------
+      //
+      // Directions of the cycle indices:
+      // -----------------------------
+      //   direction  |   j   |  q/k
+      // -----------------------------
+      //       0      |   1   |   0
+      // -----------------------------
+      //       1      |   0   |   1
+      // -----------------------------
+      constexpr int Nj = (direction < 1 ? n_rows : n_columns),
+                    Nk = (contract_over_rows ? n_rows : n_columns),
+                    Nq = (contract_over_rows ? n_columns : n_rows);
+
+      Assert(shape_data.size() == n_rows * n_columns, ExcInternalError());
+      Assert(in.size() >= Nj * Nk, ExcInternalError());
+      Assert(out.size() >= Nj * Nq, ExcInternalError());
+
+      auto thread_policy =
         Kokkos::TeamThreadMDRange<Kokkos::Rank<2>, TeamType>(team_member,
-                                                             n_q_points_1d,
-                                                             n_q_points_1d);
-      Kokkos::parallel_for(thread_policy, [&](const int i, const int j) {
-        int q_point = i + j * n_q_points_1d;
-
-        // This loop simply multiplies the shape function at the quadrature
-        // point by the value finite element coefficient.
-        // FIXME check why using parallel_reduce ThreadVector is slower
-        const int base_shape   = dof_to_quad ? j : j * n_q_points_1d;
-        const int stride_shape = dof_to_quad ? n_q_points_1d : 1;
-        const int base_in      = (direction == 0) ? (n_q_points_1d * i) : i;
-        const int stride_in    = Utilities::pow(n_q_points_1d, direction);
-        Number    sum          = shape_data[base_shape] * in(base_in);
-        for (int k = 1; k < n_q_points_1d; ++k)
-          {
-            sum += shape_data[base_shape + k * stride_shape] *
-                   in(base_in + k * stride_in);
-          }
-        t[q_point] = sum;
-      });
+                                                             Nj,
+                                                             Nq);
+      Kokkos::parallel_for(thread_policy, [&](const int j, const int q) {
+        const int base_shape   = contract_over_rows ? q : q * n_columns;
+        const int stride_shape = contract_over_rows ? n_columns : 1;
 
-      if (in_place)
-        team_member.team_barrier();
+        const int base_in   = (direction == 0 ? j * Nk : j);
+        const int stride_in = Utilities::pow(n_columns, direction);
 
-      Kokkos::parallel_for(thread_policy, [&](const int i, const int j) {
-        const int          q_point = i + j * n_q_points_1d;
-        const unsigned int destination_idx =
-          (direction == 0) ? (j + n_q_points_1d * i) : (i + n_q_points_1d * j);
+        Number sum = shape_data(base_shape) * in(base_in);
+        for (int k = 1; k < Nk; ++k)
+          sum += shape_data(base_shape + k * stride_shape) *
+                 in(base_in + k * stride_in);
 
-        if (add)
-          Kokkos::atomic_add(&out(destination_idx), t[q_point]);
+        const int index_out = (direction == 0 ? j * Nq + q : j + q * Nj);
+
+        if constexpr (add)
+          Kokkos::atomic_add(&out(index_out), sum);
         else
-          out(destination_idx) = t[q_point];
+          out(index_out) = sum;
       });
+
+      team_member.team_barrier();
     }
 
 
 
     /**
-     * Helper function for values() and gradients() in 3D
+     * Helper function for apply() in 3D
      */
-    template <int n_q_points_1d,
+    template <int n_rows,
+              int n_columns,
+              int direction,
               typename Number,
-              int  direction,
-              bool dof_to_quad,
+              bool contract_over_rows,
               bool add,
-              bool in_place,
               typename ViewTypeIn,
               typename ViewTypeOut>
     DEAL_II_HOST_DEVICE void
@@ -184,68 +221,78 @@ namespace Portable
     {
       using TeamType = Kokkos::TeamPolicy<
         MemorySpace::Default::kokkos_space::execution_space>::member_type;
-      constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 3);
 
-      Number t[n_q_points];
+      // Sizes of the input and output vectors:
+      // ------------------------------------------------------------------
+      //   direction  |    contract_over_rows    |   !contract_over_rows
+      // ------------------------------------------------------------------
+      //       0      |  m x m x m -> n x m x m  |  n x m x m -> m x m x m
+      // ------------------------------------------------------------------
+      //       1      |  n x m x m -> n x n x m  |  n x n x m -> n x m x m
+      // ------------------------------------------------------------------
+      //       2      |  n x n x m -> n x n x n  |  n x n x n -> n x n x m
+      // ------------------------------------------------------------------
+      //
+      // Directions of the cycle indices:
+      // -------------------------------------
+      //   direction  |   i   |   j   |  q/k
+      // -------------------------------------
+      //       0      |   2   |   1   |   0
+      // -------------------------------------
+      //       1      |   0   |   2   |   1
+      // -------------------------------------
+      //       2      |   1   |   0   |   2
+      // -------------------------------------
+      constexpr int Ni = (direction < 1 ? n_rows : n_columns),
+                    Nj = (direction < 2 ? n_rows : n_columns),
+                    Nk = (contract_over_rows ? n_rows : n_columns),
+                    Nq = (contract_over_rows ? n_columns : n_rows);
+
+      Assert(shape_data.size() == n_rows * n_columns, ExcInternalError());
+      Assert(in.size() >= Ni * Nj * Nk, ExcInternalError());
+      Assert(out.size() >= Ni * Nj * Nq, ExcInternalError());
+
       auto thread_policy = Kokkos::TeamThreadMDRange<Kokkos::Rank<3>, TeamType>(
-        team_member, n_q_points_1d, n_q_points_1d, n_q_points_1d);
+        team_member, Ni, Nj, Nq);
       Kokkos::parallel_for(
         thread_policy, [&](const int i, const int j, const int q) {
-          const int q_point =
-            i + j * n_q_points_1d + q * n_q_points_1d * n_q_points_1d;
-
-          // This loop simply multiplies the shape function at the quadrature
-          // point by the value finite element coefficient.
-          // FIXME check why using parallel_reduce ThreadVector is slower
-          const int base_shape   = dof_to_quad ? q : q * n_q_points_1d;
-          const int stride_shape = dof_to_quad ? n_q_points_1d : 1;
+          const int base_shape   = contract_over_rows ? q : q * n_columns;
+          const int stride_shape = contract_over_rows ? n_columns : 1;
+
           const int base_in =
-            (direction == 0 ?
-               (n_q_points_1d * (i + n_q_points_1d * j)) :
-               (direction == 1 ? (i + n_q_points_1d * n_q_points_1d * j) :
-                                 (i + n_q_points_1d * j)));
-          const int stride_in = Utilities::pow(n_q_points_1d, direction);
-          Number    sum       = shape_data[base_shape] * in(base_in);
-          for (int k = 1; k < n_q_points_1d; ++k)
-            {
-              sum += shape_data[base_shape + k * stride_shape] *
-                     in(base_in + k * stride_in);
-            }
-          t[q_point] = sum;
-        });
+            (direction == 0 ? (i * Nj + j) * Nk :
+                              (direction == 1 ? i + j * Ni * Nk : i * Nj + j));
+          const int stride_in = Utilities::pow(n_columns, direction);
 
-      if (in_place)
-        team_member.team_barrier();
+          Number sum = shape_data(base_shape) * in(base_in);
+          for (int k = 1; k < Nk; ++k)
+            sum += shape_data(base_shape + k * stride_shape) *
+                   in(base_in + k * stride_in);
 
-      Kokkos::parallel_for(
-        thread_policy, [&](const int i, const int j, const int q) {
-          const int q_point =
-            i + j * n_q_points_1d + q * n_q_points_1d * n_q_points_1d;
-          const unsigned int destination_idx =
-            (direction == 0) ? (q + n_q_points_1d * (i + n_q_points_1d * j)) :
-            (direction == 1) ? (i + n_q_points_1d * (q + n_q_points_1d * j)) :
-                               (i + n_q_points_1d * (j + n_q_points_1d * q));
-
-          if (add)
-            Kokkos::atomic_add(&out(destination_idx), t[q_point]);
+          const int index_out =
+            (direction == 0 ? (i * Nj + j) * Nq + q :
+                              (direction == 1 ? i + (j * Nq + q) * Ni :
+                                                (i + q * Ni) * Nj + j));
+
+          if constexpr (add)
+            Kokkos::atomic_add(&out(index_out), sum);
           else
-            out(destination_idx) = t[q_point];
+            out(index_out) = sum;
         });
+
+      team_member.team_barrier();
     }
 #endif
 
 
 
-    /**
-     * Helper function for values() and gradients().
-     */
     template <int dim,
-              int n_q_points_1d,
+              int n_rows,
+              int n_columns,
               typename Number,
               int  direction,
-              bool dof_to_quad,
+              bool contract_over_rows,
               bool add,
-              bool in_place,
               typename ViewTypeIn,
               typename ViewTypeOut>
     DEAL_II_HOST_DEVICE void
@@ -259,84 +306,68 @@ namespace Portable
     {
 #if KOKKOS_VERSION >= 40000
       if constexpr (dim == 1)
-        apply_1d<n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
+        apply_1d<n_rows, n_columns, direction, Number, contract_over_rows, add>(
           team_member, shape_data, in, out);
       if constexpr (dim == 2)
-        apply_2d<n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
+        apply_2d<n_rows, n_columns, direction, Number, contract_over_rows, add>(
           team_member, shape_data, in, out);
       if constexpr (dim == 3)
-        apply_3d<n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
+        apply_3d<n_rows, n_columns, direction, Number, contract_over_rows, add>(
           team_member, shape_data, in, out);
 #else
-      constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
+      // I: [0, m^{dim - direction - 1})
+      // J: [0, n^direction)
+      constexpr int NI = Utilities::pow(n_rows, dim - direction - 1);
+      constexpr int NJ = Utilities::pow(n_columns, direction);
+
+      constexpr int Nk = contract_over_rows ? n_rows : n_columns;
+      constexpr int Nq = contract_over_rows ? n_columns : n_rows;
+
+      Assert(shape_data.size() == n_rows * n_columns, ExcInternalError());
+      Assert(in.size() >= NI * NJ * Nk, ExcInternalError());
+      Assert(out.size() >= NI * NJ * Nq, ExcInternalError());
+
+      constexpr int N      = NI * NJ * Nq;
+      constexpr int stride = Utilities::pow(n_columns, direction);
 
-      Number t[n_q_points];
       Kokkos::parallel_for(
-        Kokkos::TeamThreadRange(team_member, n_q_points),
-        [&](const int &q_point) {
-          const unsigned int i = (dim == 1) ? 0 : q_point % n_q_points_1d;
-          const unsigned int j =
-            (dim == 3) ? (q_point / n_q_points_1d) % n_q_points_1d : 0;
-          const unsigned int q =
-            (dim == 1) ? q_point :
-            (dim == 2) ? (q_point / n_q_points_1d) % n_q_points_1d :
-                         q_point / (n_q_points_1d * n_q_points_1d);
-
-          // This loop simply multiplies the shape function at the quadrature
-          // point by the value finite element coefficient.
-          const int stride_shape = dof_to_quad ? n_q_points_1d : 1;
-          const int stride       = Utilities::pow(n_q_points_1d, direction);
-          const int base_shape   = dof_to_quad ? q : (q * n_q_points_1d);
-          const int base =
-            (direction == 0) ? (n_q_points_1d * (i + n_q_points_1d * j)) :
-            (direction == 1) ? (i + n_q_points_1d * (n_q_points_1d * j)) :
-                               (i + n_q_points_1d * j);
-          Number sum =
-            shape_data[base_shape] * (in_place ? out(base) : in(base));
-          for (int k = 1; k < n_q_points_1d; ++k)
+        Kokkos::TeamThreadRange(team_member, N), [&](const int index_out) {
+          // index_in  = (I Nk + k) n^direction + J
+          // index_out = (I Nq + q) n^direction + J
+          const int q = (index_out / stride) % Nq;
+          const int I = (index_out / stride) / Nq;
+          const int J = index_out % stride;
+
+          const int base_shape   = contract_over_rows ? q : q * n_columns;
+          const int stride_shape = contract_over_rows ? n_columns : 1;
+          const int base_in      = I * Nk * stride + J;
+
+          Number sum = shape_data(base_shape) * in(base_in);
+          for (int k = 1; k < Nk; ++k)
             {
-              sum +=
-                shape_data[base_shape + k * stride_shape] *
-                (in_place ? out(base + k * stride) : in(base + k * stride));
+              const int index_in = (I * Nk + k) * stride + J;
+              sum += shape_data(base_shape + k * stride_shape) * in(index_in);
             }
-          t[q_point] = sum;
-        });
-
-      if (in_place)
-        team_member.team_barrier();
 
-      Kokkos::parallel_for(
-        Kokkos::TeamThreadRange(team_member, n_q_points),
-        [&](const int &q_point) {
-          const unsigned int i = (dim == 1) ? 0 : q_point % n_q_points_1d;
-          const unsigned int j =
-            (dim == 3) ? (q_point / n_q_points_1d) % n_q_points_1d : 0;
-          const unsigned int q =
-            (dim == 1) ? q_point :
-            (dim == 2) ? (q_point / n_q_points_1d) % n_q_points_1d :
-                         q_point / (n_q_points_1d * n_q_points_1d);
-
-          const unsigned int destination_idx =
-            (direction == 0) ? (q + n_q_points_1d * (i + n_q_points_1d * j)) :
-            (direction == 1) ? (i + n_q_points_1d * (q + n_q_points_1d * j)) :
-                               (i + n_q_points_1d * (j + n_q_points_1d * q));
-
-          if (add)
-            Kokkos::atomic_add(&out(destination_idx), t[q_point]);
+          if constexpr (add)
+            Kokkos::atomic_add(&out(index_out), sum);
           else
-            out(destination_idx) = t[q_point];
+            out(index_out) = sum;
         });
+
+      team_member.team_barrier();
 #endif
     }
 
 
+
     /**
      * Generic evaluator framework.
      */
     template <EvaluatorVariant variant,
               int              dim,
-              int              fe_degree,
-              int              n_q_points_1d,
+              int              n_rows,
+              int              n_columns,
               typename Number>
     struct EvaluatorTensorProduct
     {};
@@ -347,17 +378,22 @@ namespace Portable
      * Internal evaluator for 1d-3d shape function using the tensor product form
      * of the basis functions.
      */
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+    template <int dim, int n_rows, int n_columns, typename Number>
     struct EvaluatorTensorProduct<evaluate_general,
                                   dim,
-                                  fe_degree,
-                                  n_q_points_1d,
+                                  n_rows,
+                                  n_columns,
                                   Number>
     {
     public:
       using TeamHandle = Kokkos::TeamPolicy<
         MemorySpace::Default::kokkos_space::execution_space>::member_type;
 
+      using SharedView = Kokkos::View<Number *,
+                                      MemorySpace::Default::kokkos_space::
+                                        execution_space::scratch_memory_space,
+                                      Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
+
       DEAL_II_HOST_DEVICE
       EvaluatorTensorProduct(
         const TeamHandle                                          &team_member,
@@ -365,53 +401,8 @@ namespace Portable
         Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
           shape_gradients,
         Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
-          co_shape_gradients);
-
-      /**
-       * Evaluate the finite element function at the quadrature points.
-       */
-      template <typename ViewType>
-      DEAL_II_HOST_DEVICE void
-      evaluate_values(ViewType u);
-
-      /**
-       * Evaluate the gradients of the finite element function at the quadrature
-       * points.
-       */
-      template <typename ViewTypeIn, typename ViewTypeOut>
-      DEAL_II_HOST_DEVICE void
-      evaluate_gradients(const ViewTypeIn u, ViewTypeOut grad_u);
-
-      /**
-       * Evaluate the values and the gradients of the finite element function at
-       * the quadrature points.
-       */
-      template <typename ViewType1, typename ViewType2>
-      DEAL_II_HOST_DEVICE void
-      evaluate_values_and_gradients(ViewType1 u, ViewType2 grad_u);
-
-      /**
-       * Helper function for integrate(). Integrate the finite element function.
-       */
-      template <typename ViewType>
-      DEAL_II_HOST_DEVICE void
-      integrate_values(ViewType u);
-
-      /**
-       * Helper function for integrate(). Integrate the gradients of the finite
-       * element function.
-       */
-      template <bool add, typename ViewType1, typename ViewType2>
-      DEAL_II_HOST_DEVICE void
-      integrate_gradients(ViewType1 u, ViewType2 grad_u);
-
-      /**
-       * Helper function for integrate(). Integrate the values and the gradients
-       * of the finite element function.
-       */
-      template <typename ViewType1, typename ViewType2>
-      DEAL_II_HOST_DEVICE void
-      integrate_values_and_gradients(ViewType1 u, ViewType2 grad_u);
+                   co_shape_gradients,
+        SharedView temp);
 
       /**
        * Evaluate/integrate the values of a finite element function at the
@@ -439,7 +430,6 @@ namespace Portable
       DEAL_II_HOST_DEVICE void
       gradients(const ViewTypeIn in, ViewTypeOut out) const;
 
-    public:
       /**
        * Evaluate the gradient of a finite element function at the quadrature
        * points for a given @p direction for collocation methods.
@@ -474,33 +464,36 @@ namespace Portable
        */
       Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
         co_shape_gradients;
+
+      /**
+       * Temporary storage for in-place evaluations.
+       */
+      SharedView temp;
     };
 
 
 
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+    template <int dim, int n_rows, int n_columns, typename Number>
     DEAL_II_HOST_DEVICE
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::
+    EvaluatorTensorProduct<evaluate_general, dim, n_rows, n_columns, Number>::
       EvaluatorTensorProduct(
         const TeamHandle                                          &team_member,
         Kokkos::View<Number *, MemorySpace::Default::kokkos_space> shape_values,
         Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
           shape_gradients,
         Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
-          co_shape_gradients)
+                   co_shape_gradients,
+        SharedView temp)
       : team_member(team_member)
       , shape_values(shape_values)
       , shape_gradients(shape_gradients)
       , co_shape_gradients(co_shape_gradients)
+      , temp(temp)
     {}
 
 
 
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+    template <int dim, int n_rows, int n_columns, typename Number>
     template <int  direction,
               bool dof_to_quad,
               bool add,
@@ -508,20 +501,24 @@ namespace Portable
               typename ViewTypeIn,
               typename ViewTypeOut>
     DEAL_II_HOST_DEVICE void
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::values(const ViewTypeIn in,
-                                           ViewTypeOut      out) const
+    EvaluatorTensorProduct<evaluate_general, dim, n_rows, n_columns, Number>::
+      values(const ViewTypeIn in, ViewTypeOut out) const
     {
-      apply<dim, n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
-        team_member, shape_values, in, out);
+      if constexpr (in_place)
+        {
+          apply<dim, n_rows, n_columns, Number, direction, dof_to_quad, false>(
+            team_member, shape_values, in, temp);
+
+          populate_view<add>(team_member, out, temp, out.extent(0));
+        }
+      else
+        apply<dim, n_rows, n_columns, Number, direction, dof_to_quad, add>(
+          team_member, shape_values, in, out);
     }
 
 
 
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+    template <int dim, int n_rows, int n_columns, typename Number>
     template <int  direction,
               bool dof_to_quad,
               bool add,
@@ -529,20 +526,24 @@ namespace Portable
               typename ViewTypeIn,
               typename ViewTypeOut>
     DEAL_II_HOST_DEVICE void
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::gradients(const ViewTypeIn in,
-                                              ViewTypeOut      out) const
+    EvaluatorTensorProduct<evaluate_general, dim, n_rows, n_columns, Number>::
+      gradients(const ViewTypeIn in, ViewTypeOut out) const
     {
-      apply<dim, n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
-        team_member, shape_gradients, in, out);
+      if constexpr (in_place)
+        {
+          apply<dim, n_rows, n_columns, Number, direction, dof_to_quad, false>(
+            team_member, shape_gradients, in, temp);
+
+          populate_view<add>(team_member, out, temp, out.extent(0));
+        }
+      else
+        apply<dim, n_rows, n_columns, Number, direction, dof_to_quad, add>(
+          team_member, shape_gradients, in, out);
     }
 
 
 
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+    template <int dim, int n_rows, int n_columns, typename Number>
     template <int  direction,
               bool dof_to_quad,
               bool add,
@@ -550,331 +551,24 @@ namespace Portable
               typename ViewTypeIn,
               typename ViewTypeOut>
     DEAL_II_HOST_DEVICE void
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::co_gradients(const ViewTypeIn in,
-                                                 ViewTypeOut      out) const
-    {
-      apply<dim, n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
-        team_member, co_shape_gradients, in, out);
-    }
-
-
-
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    template <typename ViewType>
-    DEAL_II_HOST_DEVICE inline void
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::evaluate_values(ViewType u)
-    {
-      if constexpr (dim == 1)
-        values<0, true, false, true>(u, u);
-      else if constexpr (dim == 2)
-        {
-          values<0, true, false, true>(u, u);
-          team_member.team_barrier();
-          values<1, true, false, true>(u, u);
-        }
-      else if constexpr (dim == 3)
-        {
-          values<0, true, false, true>(u, u);
-          team_member.team_barrier();
-          values<1, true, false, true>(u, u);
-          team_member.team_barrier();
-          values<2, true, false, true>(u, u);
-        }
-      else
-        Kokkos::abort("dim must not exceed 3!");
-    }
-
-
-
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    template <typename ViewType>
-    DEAL_II_HOST_DEVICE inline void
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::integrate_values(ViewType u)
+    EvaluatorTensorProduct<evaluate_general, dim, n_rows, n_columns, Number>::
+      co_gradients(const ViewTypeIn in, ViewTypeOut out) const
     {
-      if constexpr (dim == 1)
-        values<0, false, false, true>(u, u);
-      else if constexpr (dim == 2)
-        {
-          values<0, false, false, true>(u, u);
-          team_member.team_barrier();
-          values<1, false, false, true>(u, u);
-        }
-      else if constexpr (dim == 3)
-        {
-          values<0, false, false, true>(u, u);
-          team_member.team_barrier();
-          values<1, false, false, true>(u, u);
-          team_member.team_barrier();
-          values<2, false, false, true>(u, u);
-        }
-      else
-        Kokkos::abort("dim must not exceed 3!");
-    }
-
-
-
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    template <typename ViewTypeIn, typename ViewTypeOut>
-    DEAL_II_HOST_DEVICE inline void
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::evaluate_gradients(const ViewTypeIn u,
-                                                       ViewTypeOut      grad_u)
-    {
-      if constexpr (dim == 1)
-        {
-          gradients<0, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
-        }
-      else if constexpr (dim == 2)
-        {
-          gradients<0, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
-          values<0, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 1));
-
-          team_member.team_barrier();
-
-          values<1, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 0),
-                                       Kokkos::subview(grad_u, Kokkos::ALL, 0));
-          gradients<1, true, false, true>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 1),
-            Kokkos::subview(grad_u, Kokkos::ALL, 1));
-        }
-      else if constexpr (dim == 3)
-        {
-          gradients<0, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
-          values<0, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 1));
-          values<0, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 2));
-
-          team_member.team_barrier();
-
-          values<1, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 0),
-                                       Kokkos::subview(grad_u, Kokkos::ALL, 0));
-          gradients<1, true, false, true>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 1),
-            Kokkos::subview(grad_u, Kokkos::ALL, 1));
-          values<1, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 2),
-                                       Kokkos::subview(grad_u, Kokkos::ALL, 2));
-
-          team_member.team_barrier();
-
-          values<2, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 0),
-                                       Kokkos::subview(grad_u, Kokkos::ALL, 0));
-          values<2, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 1),
-                                       Kokkos::subview(grad_u, Kokkos::ALL, 1));
-          gradients<2, true, false, true>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 2),
-            Kokkos::subview(grad_u, Kokkos::ALL, 2));
-        }
-      else
-        Kokkos::abort("dim must not exceed 3!");
-    }
-
-
-
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    template <typename ViewType1, typename ViewType2>
-    DEAL_II_HOST_DEVICE inline void
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::evaluate_values_and_gradients(ViewType1 u,
-                                                                  ViewType2
-                                                                    grad_u)
-    {
-      if constexpr (dim == 1)
-        {
-          values<0, true, false, true>(u, u);
-          team_member.team_barrier();
-
-          co_gradients<0, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
-        }
-      else if constexpr (dim == 2)
-        {
-          values<0, true, false, true>(u, u);
-          team_member.team_barrier();
-          values<1, true, false, true>(u, u);
-          team_member.team_barrier();
-
-          co_gradients<0, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
-          co_gradients<1, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 1));
-        }
-      else if constexpr (dim == 3)
-        {
-          values<0, true, false, true>(u, u);
-          team_member.team_barrier();
-          values<1, true, false, true>(u, u);
-          team_member.team_barrier();
-          values<2, true, false, true>(u, u);
-          team_member.team_barrier();
-
-          co_gradients<0, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
-          co_gradients<1, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 1));
-          co_gradients<2, true, false, false>(
-            u, Kokkos::subview(grad_u, Kokkos::ALL, 2));
-        }
-      else
-        Kokkos::abort("dim must not exceed 3!");
-    }
-
-
-
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    template <bool add, typename ViewType1, typename ViewType2>
-    DEAL_II_HOST_DEVICE inline void
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::integrate_gradients(ViewType1 u,
-                                                        ViewType2 grad_u)
-    {
-      if constexpr (dim == 1)
-        {
-          gradients<0, false, add, false>(
-            Kokkos::subview(grad_u, Kokkos::ALL, dim), u);
-        }
-      else if constexpr (dim == 2)
-        {
-          gradients<0, false, false, true>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 0),
-            Kokkos::subview(grad_u, Kokkos::ALL, 0));
-          values<0, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 1),
-                                        Kokkos::subview(grad_u,
-                                                        Kokkos::ALL,
-                                                        1));
-
-          team_member.team_barrier();
-
-          values<1, false, add, false>(Kokkos::subview(grad_u, Kokkos::ALL, 0),
-                                       u);
-          team_member.team_barrier();
-          gradients<1, false, true, false>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
-        }
-      else if constexpr (dim == 3)
-        {
-          gradients<0, false, false, true>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 0),
-            Kokkos::subview(grad_u, Kokkos::ALL, 0));
-          values<0, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 1),
-                                        Kokkos::subview(grad_u,
-                                                        Kokkos::ALL,
-                                                        1));
-          values<0, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 2),
-                                        Kokkos::subview(grad_u,
-                                                        Kokkos::ALL,
-                                                        2));
-
-          team_member.team_barrier();
-
-          values<1, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 0),
-                                        Kokkos::subview(grad_u,
-                                                        Kokkos::ALL,
-                                                        0));
-          gradients<1, false, false, true>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 1),
-            Kokkos::subview(grad_u, Kokkos::ALL, 1));
-          values<1, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 2),
-                                        Kokkos::subview(grad_u,
-                                                        Kokkos::ALL,
-                                                        2));
-
-          team_member.team_barrier();
-
-          values<2, false, add, false>(Kokkos::subview(grad_u, Kokkos::ALL, 0),
-                                       u);
-          team_member.team_barrier();
-          values<2, false, true, false>(Kokkos::subview(grad_u, Kokkos::ALL, 1),
-                                        u);
-          team_member.team_barrier();
-          gradients<2, false, true, false>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
-        }
-      else
-        Kokkos::abort("dim must not exceed 3!");
-    }
-
-
-
-    template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    template <typename ViewType1, typename ViewType2>
-    DEAL_II_HOST_DEVICE inline void
-    EvaluatorTensorProduct<evaluate_general,
-                           dim,
-                           fe_degree,
-                           n_q_points_1d,
-                           Number>::integrate_values_and_gradients(ViewType1 u,
-                                                                   ViewType2
-                                                                     grad_u)
-    {
-      if constexpr (dim == 1)
-        {
-          co_gradients<0, false, true, false>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
-          team_member.team_barrier();
-
-          values<0, false, false, true>(u, u);
-        }
-      else if constexpr (dim == 2)
-        {
-          co_gradients<1, false, true, false>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
-          team_member.team_barrier();
-          co_gradients<0, false, true, false>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
-          team_member.team_barrier();
-
-          values<1, false, false, true>(u, u);
-          team_member.team_barrier();
-          values<0, false, false, true>(u, u);
-          team_member.team_barrier();
-        }
-      else if constexpr (dim == 3)
+      if constexpr (in_place)
         {
-          co_gradients<2, false, true, false>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
-          team_member.team_barrier();
-          co_gradients<1, false, true, false>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
-          team_member.team_barrier();
-          co_gradients<0, false, true, false>(
-            Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
-          team_member.team_barrier();
-
-          values<2, false, false, true>(u, u);
-          team_member.team_barrier();
-          values<1, false, false, true>(u, u);
-          team_member.team_barrier();
-          values<0, false, false, true>(u, u);
-          team_member.team_barrier();
+          apply<dim,
+                n_columns,
+                n_columns,
+                Number,
+                direction,
+                dof_to_quad,
+                false>(team_member, co_shape_gradients, in, temp);
+
+          populate_view<add>(team_member, out, temp, out.extent(0));
         }
       else
-        Kokkos::abort("dim must not exceed 3!");
+        apply<dim, n_columns, n_columns, Number, direction, dof_to_quad, add>(
+          team_member, co_shape_gradients, in, out);
     }
   } // namespace internal
 } // namespace Portable
index 7daa538697ddc412a3e362a190217b2252a34c83..037bd421093ea58c32dbf9c2f5fb37de78d9b715 100644 (file)
@@ -44,13 +44,23 @@ evaluate_tensor_product(
   Kokkos::View<double *, MemorySpace::Default::kokkos_space> dst,
   Kokkos::View<double *, MemorySpace::Default::kokkos_space> src)
 {
+  Kokkos::View<
+    double *,
+    MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
+    Kokkos::MemoryTraits<Kokkos::Unmanaged>>
+    dummy_scratch(team_member.team_shmem(), 0);
+
   Portable::internal::EvaluatorTensorProduct<
     Portable::internal::evaluate_general,
     1,
-    M - 1,
+    M,
     N,
     double>
-    evaluator(team_member, shape_values, shape_gradients, co_shape_gradients);
+    evaluator(team_member,
+              shape_values,
+              shape_gradients,
+              co_shape_gradients,
+              dummy_scratch);
 
   if (type == 0)
     evaluator.template values<0, dof_to_quad, add, false>(src, dst);
@@ -197,19 +207,19 @@ main()
   deallog.push("values");
   test<4, 4, 0, false>();
   test<3, 3, 0, false>();
+  test<3, 4, 0, false>();
+  test<3, 5, 0, false>();
   // Not supported right now
   // test<4,3,0,false>();
-  // test<3,4,0,false>();
-  // test<3,5,0,false>();
   deallog.pop();
 
   deallog.push("gradients");
   test<4, 4, 1, false>();
   test<3, 3, 1, false>();
+  test<3, 4, 1, false>();
+  test<3, 5, 1, false>();
   // Not supported right now
   // test<4,3,1,false>();
-  // test<3,4,1,false>();
-  // test<3,5,1,false>();
   deallog.pop();
 
   deallog.push("add");
@@ -217,19 +227,19 @@ main()
   deallog.push("values");
   test<4, 4, 0, true>();
   test<3, 3, 0, true>();
+  test<3, 4, 0, true>();
+  test<3, 5, 0, true>();
   // Not supported right now
   // test<4,3,0,true>();
-  // test<3,4,0,true>();
-  // test<3,5,0,true>();
   deallog.pop();
 
   deallog.push("gradients");
   test<4, 4, 1, true>();
   test<3, 3, 1, true>();
+  test<3, 4, 1, true>();
+  test<3, 5, 1, true>();
   // Not supported right now
   // test<4,3,1,true>();
-  // test<3,4,1,true>();
-  // test<3,5,1,true>();
   deallog.pop();
 
   deallog.pop();
index 5aa83953a84a83b1cdade8c0e44cd47f3a63e918..b18c6ceb08b87b63c55a874bf5fe9ea883315a37 100644 (file)
@@ -5,21 +5,45 @@ DEAL:values::Errors transpose:    0 0 0 0
 DEAL:values::Test 3 x 3
 DEAL:values::Errors no transpose: 0 0 0 
 DEAL:values::Errors transpose:    0 0 0 
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 0 
 DEAL:gradients::Test 4 x 4
 DEAL:gradients::Errors no transpose: 0 0 0 0 
 DEAL:gradients::Errors transpose:    0 0 0 0 
 DEAL:gradients::Test 3 x 3
 DEAL:gradients::Errors no transpose: 0 0 0 
 DEAL:gradients::Errors transpose:    0 0 0 
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 0 
 DEAL:add:values::Test 4 x 4
 DEAL:add:values::Errors no transpose: 0 0 0 0 
 DEAL:add:values::Errors transpose:    0 0 0 0 
 DEAL:add:values::Test 3 x 3
 DEAL:add:values::Errors no transpose: 0 0 0 
 DEAL:add:values::Errors transpose:    0 0 0 
+DEAL:add:values::Test 3 x 4
+DEAL:add:values::Errors no transpose: 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 
+DEAL:add:values::Test 3 x 5
+DEAL:add:values::Errors no transpose: 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 0 
 DEAL:add:gradients::Test 4 x 4
 DEAL:add:gradients::Errors no transpose: 0 0 0 0 
 DEAL:add:gradients::Errors transpose:    0 0 0 0 
 DEAL:add:gradients::Test 3 x 3
 DEAL:add:gradients::Errors no transpose: 0 0 0 
 DEAL:add:gradients::Errors transpose:    0 0 0 
+DEAL:add:gradients::Test 3 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 
+DEAL:add:gradients::Test 3 x 5
+DEAL:add:gradients::Errors no transpose: 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 0 
index a2f51058e46ac7cec971d8ddf41efe48676718a7..e82a4fea4f4f3d821e7398fd606a48de8d95ba0b 100644 (file)
@@ -41,27 +41,41 @@ evaluate_tensor_product(
   Kokkos::View<double *, MemorySpace::Default::kokkos_space> shape_gradients,
   Kokkos::View<double *, MemorySpace::Default::kokkos_space> co_shape_gradients,
   Kokkos::View<double *, MemorySpace::Default::kokkos_space> dst,
-  Kokkos::View<double *, MemorySpace::Default::kokkos_space> src)
+  Kokkos::View<double *, MemorySpace::Default::kokkos_space> src,
+  Kokkos::View<double *, MemorySpace::Default::kokkos_space> tmp)
 {
+  Kokkos::View<
+    double *,
+    MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
+    Kokkos::MemoryTraits<Kokkos::Unmanaged>>
+    dummy_scratch(team_member.team_shmem(), 0);
+
   Portable::internal::EvaluatorTensorProduct<
     Portable::internal::evaluate_general,
     2,
-    M - 1,
+    M,
     N,
     double>
-    evaluator(team_member, shape_values, shape_gradients, co_shape_gradients);
+    evaluator(team_member,
+              shape_values,
+              shape_gradients,
+              co_shape_gradients,
+              dummy_scratch);
+
+  constexpr int d0 = dof_to_quad ? 0 : 1;
+  constexpr int d1 = dof_to_quad ? 1 : 0;
 
   if (type == 0)
     {
-      evaluator.template values<0, dof_to_quad, false, false>(src, src);
+      evaluator.template values<d0, dof_to_quad, false, false>(src, tmp);
       team_member.team_barrier();
-      evaluator.template values<1, dof_to_quad, add, false>(src, dst);
+      evaluator.template values<d1, dof_to_quad, add, false>(tmp, dst);
     }
   if (type == 1)
     {
-      evaluator.template gradients<0, dof_to_quad, false, false>(src, src);
+      evaluator.template gradients<d0, dof_to_quad, false, false>(src, tmp);
       team_member.team_barrier();
-      evaluator.template gradients<1, dof_to_quad, add, false>(src, dst);
+      evaluator.template gradients<d1, dof_to_quad, add, false>(tmp, dst);
     }
 }
 
@@ -93,11 +107,14 @@ test()
 
   constexpr int                          M_2d = M * M;
   constexpr int                          N_2d = N * N;
+  constexpr int                          MN   = M * N;
   LinearAlgebra::ReadWriteVector<double> x_ref(N_2d), y_ref(M_2d);
   Kokkos::View<double[N_2d], MemorySpace::Default::kokkos_space> x_dev(
     Kokkos::view_alloc("x_dev", Kokkos::WithoutInitializing));
   Kokkos::View<double[M_2d], MemorySpace::Default::kokkos_space> y_dev(
     Kokkos::view_alloc("y_dev", Kokkos::WithoutInitializing));
+  Kokkos::View<double[MN], MemorySpace::Default::kokkos_space> tmp_dev(
+    Kokkos::view_alloc("tmp_dev", Kokkos::WithoutInitializing));
   auto x_host = Kokkos::create_mirror_view(x_dev);
   auto y_host = Kokkos::create_mirror_view(y_dev);
 
@@ -160,7 +177,8 @@ test()
                                                       shape_gradients,
                                                       co_shape_gradients,
                                                       y_dev,
-                                                      x_dev);
+                                                      x_dev,
+                                                      tmp_dev);
     });
 
   // Check the results on the host
@@ -197,7 +215,8 @@ test()
                                                      shape_gradients,
                                                      co_shape_gradients,
                                                      x_dev,
-                                                     y_dev);
+                                                     y_dev,
+                                                     tmp_dev);
     });
 
   // Check the results on the host
@@ -219,11 +238,15 @@ main()
   deallog.push("values");
   test<4, 4, 0, false>();
   test<3, 3, 0, false>();
+  test<3, 4, 0, false>();
+  test<3, 5, 0, false>();
   deallog.pop();
 
   deallog.push("gradients");
   test<4, 4, 1, false>();
   test<3, 3, 1, false>();
+  test<3, 4, 1, false>();
+  test<3, 5, 1, false>();
   deallog.pop();
 
   deallog.push("add");
@@ -231,11 +254,15 @@ main()
   deallog.push("values");
   test<4, 4, 0, true>();
   test<3, 3, 0, true>();
+  test<3, 4, 0, true>();
+  test<3, 5, 0, true>();
   deallog.pop();
 
   deallog.push("gradients");
   test<4, 4, 1, true>();
   test<3, 3, 1, true>();
+  test<3, 4, 1, true>();
+  test<3, 5, 1, true>();
   deallog.pop();
 
   deallog.pop();
index 312a78a60e3c56469a75415736ea19fe1d6365f5..432869ab0c6593a9637518176c2026f20a5e05b2 100644 (file)
@@ -5,21 +5,45 @@ DEAL:values::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 DEAL:values::Test 3 x 3
 DEAL:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 
 DEAL:values::Errors transpose:    0 0 0 0 0 0 0 0 0 
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:values::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 DEAL:gradients::Test 4 x 4
 DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 DEAL:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 DEAL:gradients::Test 3 x 3
-DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0
-DEAL:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0
+DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 DEAL:add:values::Test 4 x 4
-DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-DEAL:add:values::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 DEAL:add:values::Test 3 x 3
-DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0
-DEAL:add:values::Errors transpose:    0 0 0 0 0 0 0 0 0
+DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 0 0 0 0 0 
+DEAL:add:values::Test 3 x 4
+DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:add:values::Test 3 x 5
+DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:add:values::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 DEAL:add:gradients::Test 4 x 4
-DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-DEAL:add:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 DEAL:add:gradients::Test 3 x 3
-DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0
-DEAL:add:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 
+DEAL:add:gradients::Test 3 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:add:gradients::Test 3 x 5
+DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 
+DEAL:add:gradients::Errors transpose:    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
index 49ff8f7c7c8223ff77bc51a9cb7917f63ab1dda3..b9b5fa3678d161044d8fbd66dd0552a7bb36db5a 100644 (file)
@@ -72,9 +72,13 @@ struct DisableWindowsDebugRuntimeDialog
 // Redefine Assert as AssertThrow to make sure that the code is tested similarly
 // in Release mode and in Debug mode. clang-format makes sure that this file is
 // included after all regular header files but before all the other local header
-// files.
-#undef Assert
-#define Assert AssertThrow
+// files. The redefinition is not applied when Kokkos GPU backend is enabled,
+// because AssertThrow is not defined for device codes.
+#if !defined(KOKKOS_ENABLE_CUDA) && !defined(KOKKOS_ENABLE_HIP) && \
+  !defined(KOKKOS_ENABLE_SYCL)
+#  undef Assert
+#  define Assert AssertThrow
+#endif
 
 // implicitly use the deal.II namespace everywhere, without us having to say
 // so in each and every testcase

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.