Merge pull request #16844 from kronbichler/eval_interface

Use ArrayView for interface to tensor-product point evaluator

include/deal.II/fe/mapping_q_internal.h

@@ -477,7 +477,7 @@ namespace internal
     do_transform_real_to_unit_cell_internal(
       const Point<spacedim, Number>                      &p,
       const Point<dim, Number>                           &initial_p_unit,
-      const std::vector<Point<spacedim>>                 &points,
+      const ArrayView<const Point<spacedim>>             &points,
       const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
       const std::vector<unsigned int>                    &renumber,
       const bool print_iterations_to_deallog = false)
@@ -734,7 +734,7 @@ namespace internal
     do_transform_real_to_unit_cell_internal_codim1(
       const Point<dim + 1>                               &p,
       const Point<dim>                                   &initial_p_unit,
-      const std::vector<Point<dim + 1>>                  &points,
+      const ArrayView<const Point<dim + 1>>              &points,
       const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
       const std::vector<unsigned int>                    &renumber)
     {
@@ -1245,9 +1245,9 @@ namespace internal
       std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
       std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians)
     {
-      const UpdateFlags                   update_flags = data.update_each;
-      const std::vector<Point<spacedim>> &support_points =
-        data.mapping_support_points;
+      const UpdateFlags                      update_flags = data.update_each;
+      const ArrayView<const Point<spacedim>> support_points(
+        data.mapping_support_points);
 
       const unsigned int n_points = unit_points.size();
       const unsigned int n_lanes  = VectorizedArray<double>::size();
@@ -1406,8 +1406,8 @@ namespace internal
     {
       if (data.update_each & update_jacobian_grads)
         {
-          const std::vector<Point<spacedim>> &support_points =
-            data.mapping_support_points;
+          const ArrayView<const Point<spacedim>> support_points(
+            data.mapping_support_points);
           const unsigned int n_q_points = jacobian_grads.size();
 
           if (cell_similarity != CellSimilarity::translation)
@@ -1448,8 +1448,8 @@ namespace internal
     {
       if (data.update_each & update_jacobian_pushed_forward_grads)
         {
-          const std::vector<Point<spacedim>> &support_points =
-            data.mapping_support_points;
+          const ArrayView<const Point<spacedim>> support_points(
+            data.mapping_support_points);
           const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
 
           if (cell_similarity != CellSimilarity::translation)
@@ -1563,8 +1563,8 @@ namespace internal
     {
       if (data.update_each & update_jacobian_2nd_derivatives)
         {
-          const std::vector<Point<spacedim>> &support_points =
-            data.mapping_support_points;
+          const ArrayView<const Point<spacedim>> support_points(
+            data.mapping_support_points);
           const unsigned int n_q_points = jacobian_2nd_derivatives.size();
 
           if (cell_similarity != CellSimilarity::translation)
@@ -1603,8 +1603,8 @@ namespace internal
     {
       if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives)
         {
-          const std::vector<Point<spacedim>> &support_points =
-            data.mapping_support_points;
+          const ArrayView<const Point<spacedim>> support_points(
+            data.mapping_support_points);
           const unsigned int n_q_points =
             jacobian_pushed_forward_2nd_derivatives.size();
 
@@ -1750,8 +1750,8 @@ namespace internal
     {
       if (data.update_each & update_jacobian_3rd_derivatives)
         {
-          const std::vector<Point<spacedim>> &support_points =
-            data.mapping_support_points;
+          const ArrayView<const Point<spacedim>> support_points(
+            data.mapping_support_points);
           const unsigned int n_q_points = jacobian_3rd_derivatives.size();
 
           if (cell_similarity != CellSimilarity::translation)
@@ -1790,8 +1790,8 @@ namespace internal
     {
       if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives)
         {
-          const std::vector<Point<spacedim>> &support_points =
-            data.mapping_support_points;
+          const ArrayView<const Point<spacedim>> support_points(
+            data.mapping_support_points);
           const unsigned int n_q_points =
             jacobian_pushed_forward_3rd_derivatives.size();
 

include/deal.II/matrix_free/tensor_product_point_kernels.h

@@ -461,7 +461,7 @@ namespace internal
     Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
   evaluate_tensor_product_value_and_gradient(
     const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
+    const ArrayView<const Number>                      &values,
     const Point<dim, Number2>                          &p,
     const bool                                          d_linear = false,
     const std::vector<unsigned int>                    &renumber = {})
@@ -682,7 +682,7 @@ namespace internal
   inline typename ProductTypeNoPoint<Number, Number2>::type
   evaluate_tensor_product_value(
     const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
+    const ArrayView<const Number>                      &values,
     const Point<dim, Number2>                          &p,
     const bool                                          d_linear = false,
     const std::vector<unsigned int>                    &renumber = {})
@@ -718,7 +718,7 @@ namespace internal
   inline Tensor<1, 1, typename ProductTypeNoPoint<Number, Number2>::type>
   evaluate_tensor_product_higher_derivatives(
     const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
+    const ArrayView<const Number>                      &values,
     const Point<1, Number2>                            &p,
     const std::vector<unsigned int>                    &renumber = {})
   {
@@ -758,7 +758,7 @@ namespace internal
                 typename ProductTypeNoPoint<Number, Number2>::type>
   evaluate_tensor_product_higher_derivatives(
     const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
+    const ArrayView<const Number>                      &values,
     const Point<2, Number2>                            &p,
     const std::vector<unsigned int>                    &renumber = {})
   {
@@ -810,7 +810,7 @@ namespace internal
                 typename ProductTypeNoPoint<Number, Number2>::type>
   evaluate_tensor_product_higher_derivatives(
     const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
+    const ArrayView<const Number>                      &values,
     const Point<3, Number2>                            &p,
     const std::vector<unsigned int>                    &renumber = {})
   {
@@ -874,7 +874,7 @@ namespace internal
   SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
   evaluate_tensor_product_hessian(
     const std::vector<Polynomials::Polynomial<double>> &poly,
-    const std::vector<Number>                          &values,
+    const ArrayView<const Number>                      &values,
     const Point<dim, Number2>                          &p,
     const std::vector<unsigned int>                    &renumber = {})
   {

source/fe/mapping_q.cc

@@ -292,7 +292,7 @@ MappingQ<dim, spacedim>::transform_unit_to_real_cell(
 {
   return Point<spacedim>(internal::evaluate_tensor_product_value(
     polynomials_1d,
-    this->compute_mapping_support_points(cell),
+    make_const_array_view(this->compute_mapping_support_points(cell)),
     p,
     polynomials_1d.size() == 2,
     renumber_lexicographic_to_hierarchic));
@@ -345,7 +345,7 @@ MappingQ<1, 1>::transform_real_to_unit_cell_internal(
     do_transform_real_to_unit_cell_internal<1>(
       p,
       initial_p_unit,
-      this->compute_mapping_support_points(cell),
+      make_const_array_view(this->compute_mapping_support_points(cell)),
       polynomials_1d,
       renumber_lexicographic_to_hierarchic);
 }
@@ -363,7 +363,7 @@ MappingQ<2, 2>::transform_real_to_unit_cell_internal(
     do_transform_real_to_unit_cell_internal<2>(
       p,
       initial_p_unit,
-      this->compute_mapping_support_points(cell),
+      make_const_array_view(this->compute_mapping_support_points(cell)),
       polynomials_1d,
       renumber_lexicographic_to_hierarchic);
 }
@@ -381,7 +381,7 @@ MappingQ<3, 3>::transform_real_to_unit_cell_internal(
     do_transform_real_to_unit_cell_internal<3>(
       p,
       initial_p_unit,
-      this->compute_mapping_support_points(cell),
+      make_const_array_view(this->compute_mapping_support_points(cell)),
       polynomials_1d,
       renumber_lexicographic_to_hierarchic);
 }
@@ -414,7 +414,7 @@ MappingQ<1, 2>::transform_real_to_unit_cell_internal(
     do_transform_real_to_unit_cell_internal_codim1<1>(
       p,
       initial_p_unit,
-      mdata->mapping_support_points,
+      make_const_array_view(mdata->mapping_support_points),
       polynomials_1d,
       renumber_lexicographic_to_hierarchic);
 }
@@ -447,7 +447,7 @@ MappingQ<2, 3>::transform_real_to_unit_cell_internal(
     do_transform_real_to_unit_cell_internal_codim1<2>(
       p,
       initial_p_unit,
-      mdata->mapping_support_points,
+      make_const_array_view(mdata->mapping_support_points),
       polynomials_1d,
       renumber_lexicographic_to_hierarchic);
 }

source/non_matching/quadrature_generator.cc

@@ -1631,12 +1631,12 @@ namespace NonMatching
             return this->is_fe_q_iso_q1() ?
                      dealii::internal::evaluate_tensor_product_value(
                        poly,
-                       local_dof_values_subcell,
+                       make_array_view(local_dof_values_subcell),
                        subcell_box.real_to_unit(point),
                        polynomials_are_hat_functions) :
                      dealii::internal::evaluate_tensor_product_value(
                        poly,
-                       local_dof_values,
+                       make_array_view(local_dof_values),
                        point,
                        polynomials_are_hat_functions,
                        renumber);
@@ -1670,13 +1670,13 @@ namespace NonMatching
                       dealii::internal::
                         evaluate_tensor_product_value_and_gradient(
                           poly,
-                          local_dof_values_subcell,
+                          make_array_view(local_dof_values_subcell),
                           subcell_box.real_to_unit(point),
                           polynomials_are_hat_functions) :
                       dealii::internal::
                         evaluate_tensor_product_value_and_gradient(
                           poly,
-                          local_dof_values,
+                          make_array_view(local_dof_values),
                           point,
                           polynomials_are_hat_functions,
                           renumber))
@@ -1710,10 +1710,13 @@ namespace NonMatching
             return this->is_fe_q_iso_q1() ?
                      dealii::internal::evaluate_tensor_product_hessian(
                        poly,
-                       local_dof_values_subcell,
+                       make_array_view(local_dof_values_subcell),
                        subcell_box.real_to_unit(point)) :
                      dealii::internal::evaluate_tensor_product_hessian(
-                       poly, local_dof_values, point, renumber);
+                       poly,
+                       make_array_view(local_dof_values),
+                       point,
+                       renumber);
           }
         else
           {

tests/mappings/mapping_q_real_to_unit_internal.cc

@@ -74,7 +74,7 @@ print_result(const unsigned int                  mapping_degree,
           do_transform_real_to_unit_cell_internal(
             p,
             cell->real_to_unit_cell_affine_approximation(p),
-            fe_values.get_quadrature_points(),
+            make_array_view(fe_values.get_quadrature_points()),
             polynomials,
             renumber,
             /* print_iterations = */ true);

tests/matrix_free/tensor_product_evaluate_01.cc

@@ -71,7 +71,11 @@ test(const unsigned int degree)
   for (const auto &p : evaluation_points)
     {
       const auto val = internal::evaluate_tensor_product_value_and_gradient(
-        polynomials, coefficients, p, false, renumbering);
+        polynomials,
+        make_const_array_view(coefficients),
+        p,
+        false,
+        renumbering);
       deallog << "Value " << val.first << " vs "
               << transform * (matrix * p + offset) << " ; gradient "
               << val.second << " vs " << transform * matrix << std::endl;
@@ -83,7 +87,11 @@ test(const unsigned int degree)
       for (const auto &p : evaluation_points)
         {
           const auto val = internal::evaluate_tensor_product_value_and_gradient(
-            polynomials, coefficients, p, true, renumbering);
+            polynomials,
+            make_const_array_view(coefficients),
+            p,
+            true,
+            renumbering);
           deallog << "Value " << val.first << " vs "
                   << transform * (matrix * p + offset) << " ; gradient "
                   << val.second << " vs " << transform * matrix << std::endl;

tests/matrix_free/tensor_product_evaluate_02.cc

@@ -66,7 +66,11 @@ test(const unsigned int degree)
   for (const auto &p : evaluation_points)
     {
       const auto val = internal::evaluate_tensor_product_value_and_gradient(
-        polynomials, coefficients, p, false, renumbering);
+        polynomials,
+        make_const_array_view(coefficients),
+        p,
+        false,
+        renumbering);
       deallog << "Value " << val.first << " vs " << matrix * p + offset
               << std::endl;
       deallog << "Gradient " << val.second << " vs " << transpose(matrix)
@@ -79,7 +83,11 @@ test(const unsigned int degree)
       for (const auto &p : evaluation_points)
         {
           const auto val = internal::evaluate_tensor_product_value_and_gradient(
-            polynomials, coefficients, p, true, renumbering);
+            polynomials,
+            make_const_array_view(coefficients),
+            p,
+            true,
+            renumbering);
           deallog << "Value " << val.first << " vs " << matrix * p + offset
                   << std::endl;
           deallog << "Gradient " << val.second << " vs " << transpose(matrix)

tests/matrix_free/tensor_product_evaluate_03.cc

@@ -72,7 +72,11 @@ test(const unsigned int degree)
           p_vec[d][v] = p[d] + 0.01 * v;
 
       const auto val = internal::evaluate_tensor_product_value_and_gradient(
-        polynomials, coefficients, p_vec, false, renumbering);
+        polynomials,
+        make_const_array_view(coefficients),
+        p_vec,
+        false,
+        renumbering);
 
       const auto error_vec = val.first - matrix * p_vec;
       double     error     = 0;
@@ -100,7 +104,11 @@ test(const unsigned int degree)
               p_vec[d][v] = p[d] + 0.01 * v;
 
           const auto val = internal::evaluate_tensor_product_value_and_gradient(
-            polynomials, coefficients, p_vec, true, renumbering);
+            polynomials,
+            make_const_array_view(coefficients),
+            p_vec,
+            true,
+            renumbering);
 
           const auto error_vec = val.first - matrix * p_vec;
           double     error     = 0;

tests/matrix_free/tensor_product_evaluate_04.cc

@@ -69,7 +69,7 @@ test(const unsigned int degree)
           p_vec[d][v] = p[d] + 0.01 * v;
 
       const auto val = internal::evaluate_tensor_product_value_and_gradient(
-        polynomials, coefficients, p_vec, false);
+        polynomials, make_const_array_view(coefficients), p_vec, false);
 
       const auto error_vec = val.first - matrix * p_vec;
       double     error     = 0;
@@ -97,7 +97,7 @@ test(const unsigned int degree)
               p_vec[d][v] = p[d] + 0.01 * v;
 
           const auto val = internal::evaluate_tensor_product_value_and_gradient(
-            polynomials, coefficients, p_vec, true);
+            polynomials, make_const_array_view(coefficients), p_vec, true);
 
           const auto error_vec = val.first - matrix * p_vec;
           double     error     = 0;

tests/matrix_free/tensor_product_evaluate_05.cc

@@ -69,10 +69,8 @@ test(const unsigned int degree)
         for (unsigned int d = 0; d < dim; ++d)
           p_vec[d][v] = p[d] + 0.01 * v;
 
-      const auto hess = internal::evaluate_tensor_product_hessian(polynomials,
-                                                                  coefficients,
-                                                                  p_vec,
-                                                                  renumbering);
+      const auto hess = internal::evaluate_tensor_product_hessian(
+        polynomials, make_const_array_view(coefficients), p_vec, renumbering);
 
       double error = 0;
       for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)

tests/matrix_free/tensor_product_evaluate_06.cc

@@ -65,9 +65,8 @@ test(const unsigned int degree)
         for (unsigned int d = 0; d < dim; ++d)
           p_vec[d][v] = p[d] + 0.01 * v;
 
-      const auto hess = internal::evaluate_tensor_product_hessian(polynomials,
-                                                                  coefficients,
-                                                                  p_vec);
+      const auto hess = internal::evaluate_tensor_product_hessian(
+        polynomials, make_const_array_view(coefficients), p_vec);
 
       std::cout << hess << "    " << matrix << std::endl;
 

tests/matrix_free/tensor_product_evaluate_07.cc

@@ -75,7 +75,7 @@ test(const unsigned int degree)
       deallog << "]: ";
       const auto derivative =
         internal::evaluate_tensor_product_higher_derivatives<derivative_order>(
-          polynomials, function_values, p, renumbering);
+          polynomials, make_const_array_view(function_values), p, renumbering);
 
       for (unsigned int d = 0; d < derivative.dimension; ++d)
         deallog << (std::abs(derivative[d]) < 1e-11 ? 0. : derivative[d])

tests/matrix_free/tensor_product_evaluate_08.cc

@@ -72,7 +72,7 @@ test(const unsigned int degree)
       deallog << "]: ";
       const auto derivative =
         internal::evaluate_tensor_product_higher_derivatives<derivative_order>(
-          polynomials, function_values, p);
+          polynomials, make_const_array_view(function_values), p);
 
       for (unsigned int d = 0; d < derivative.dimension; ++d)
         deallog << (std::abs(derivative[d]) < 1e-11 ? 0. : derivative[d])