]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use ArrayView for interface to tensor-product point evaluator 16844/head
authorMartin Kronbichler <martin.kronbichler@rub.de>
Wed, 3 Apr 2024 11:39:35 +0000 (13:39 +0200)
committerMartin Kronbichler <martin.kronbichler@rub.de>
Wed, 3 Apr 2024 13:43:42 +0000 (15:43 +0200)
13 files changed:
include/deal.II/fe/mapping_q_internal.h
include/deal.II/matrix_free/tensor_product_point_kernels.h
source/fe/mapping_q.cc
source/non_matching/quadrature_generator.cc
tests/mappings/mapping_q_real_to_unit_internal.cc
tests/matrix_free/tensor_product_evaluate_01.cc
tests/matrix_free/tensor_product_evaluate_02.cc
tests/matrix_free/tensor_product_evaluate_03.cc
tests/matrix_free/tensor_product_evaluate_04.cc
tests/matrix_free/tensor_product_evaluate_05.cc
tests/matrix_free/tensor_product_evaluate_06.cc
tests/matrix_free/tensor_product_evaluate_07.cc
tests/matrix_free/tensor_product_evaluate_08.cc

index 90ba787dbc934d5cb63759df179d499d0924188c..d2e8cd9053f7398fc77bae492ffbcdcaa9b7eada 100644 (file)
@@ -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();
 
index d6dbfccaca1b85ed26a3f6468c41d562775d9461..0fed110392e9e3bba4b36b15f7341ac5875394bf 100644 (file)
@@ -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 = {})
   {
index f1631738add1e694b3a99e27ccc91fd1fd18fffe..38f5f2d7ae7bf16ab34657b7bf31584484348259 100644 (file)
@@ -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);
 }
index 3dc9d32c95470dc38bc00b327264d7293f99f550..05589996884947e09ca3f4e10f1e3a33ee48e280 100644 (file)
@@ -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
           {
index d1e568866245926194462c0fdb88b36eb49e0afb..5e02266e61230c7113fcc25b6c6ae2c95deeffb5 100644 (file)
@@ -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);
index 952f1e5fec11a25a8b00842a45e3ae7abefd195c..fcc477d74a22dd5511367924c234bf36d59f48c2 100644 (file)
@@ -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;
index eeb7f37f2cdd105510dbe44f48b6affe547b4b24..d41525d234e1583b2a009c0adbe12836ae0e9b98 100644 (file)
@@ -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)
index 2e97ca301de0882ce4f9e4797600f68c6649ab4c..7507516f063d0a502f73d23594addbfb3e59fd78 100644 (file)
@@ -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;
index b27fca1b43956b8db99c56c86260a5b3c53c3e7d..e2b4c19598af8f5bfb6171bb0bd0580fd885ddc1 100644 (file)
@@ -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;
index a3f784dd3a0ee174921cfffb63470d4d71413911..d1c4c8685a37c7116db4562a030c80d38d0e1c39 100644 (file)
@@ -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)
index ae8a9aa6fc80414f0004b70efd7896d3f7f8f3fe..754a53d486a862757c891fc855a60663d81fd3a1 100644 (file)
@@ -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;
 
index ab3e09b33d06c8657c11adbd351cf79ef44ddd7d..2dc825b1d04124af1e937ef1def27fcb3ce4d077 100644 (file)
@@ -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])
index 3b9fad9c218ac597253772479a21e34ba1f68e93..b561a04cbc4ff3dc75d8656ef548ae46d31d8e3c 100644 (file)
@@ -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])

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.