]> https://gitweb.dealii.org/ - dealii.git/commitdiff
FEFaceEvaluation: enable FE_DGP 18175/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 28 Feb 2025 08:27:07 +0000 (09:27 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 3 Mar 2025 18:54:43 +0000 (19:54 +0100)
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/evaluation_kernels_common.h [new file with mode: 0644]
include/deal.II/matrix_free/evaluation_kernels_face.h
include/deal.II/matrix_free/fe_evaluation_data.h
tests/matrix_free/matrix_vector_faces_36.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_faces_36.output [new file with mode: 0644]

index 236c157b1c0d01dbd673363eddbdc8e47498c6b4..4817976afe99826f984cb9cc47371a37bcfabcee 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/base/vectorization.h>
 
 #include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/evaluation_kernels_common.h>
 #include <deal.II/matrix_free/fe_evaluation_data.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
@@ -242,36 +243,19 @@ namespace internal
       (type == MatrixFreeFunctions::truncated_tensor) ?
         Utilities::pow(shape_data.front().fe_degree + 1, dim) :
         fe_eval.get_shape_info().dofs_per_component_on_cell;
-    const Number *values_dofs = values_dofs_actual;
+    const Number *values_dofs =
+      (type == MatrixFreeFunctions::truncated_tensor) ?
+        temp1 + 2 * (std::max<std::size_t>(
+                      fe_eval.get_shape_info().dofs_per_component_on_cell,
+                      n_q_points)) :
+        values_dofs_actual;
+
     if (type == MatrixFreeFunctions::truncated_tensor)
-      {
-        const std::size_t n_dofs_per_comp =
-          fe_eval.get_shape_info().dofs_per_component_on_cell;
-        Number *values_dofs_tmp =
-          temp1 + 2 * (std::max(n_dofs_per_comp, n_q_points));
-        const int degree =
-          fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
-        for (unsigned int c = 0; c < n_components; ++c)
-          for (int i = 0, count_p = 0, count_q = 0;
-               i < (dim > 2 ? degree + 1 : 1);
-               ++i)
-            {
-              for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
-                {
-                  for (int k = 0; k < degree + 1 - j - i;
-                       ++k, ++count_p, ++count_q)
-                    values_dofs_tmp[c * dofs_per_comp + count_q] =
-                      values_dofs_actual[c * n_dofs_per_comp + count_p];
-                  for (int k = degree + 1 - j - i; k < degree + 1;
-                       ++k, ++count_q)
-                    values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
-                }
-              for (int j = degree + 1 - i; j < degree + 1; ++j)
-                for (int k = 0; k < degree + 1; ++k, ++count_q)
-                  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
-            }
-        values_dofs = values_dofs_tmp;
-      }
+      embed_truncated_into_full_tensor_product<dim, fe_degree>(
+        n_components,
+        const_cast<Number *>(values_dofs),
+        values_dofs_actual,
+        fe_eval);
 
     Number *values_quad    = fe_eval.begin_values();
     Number *gradients_quad = fe_eval.begin_gradients();
@@ -436,6 +420,7 @@ namespace internal
       (type == MatrixFreeFunctions::truncated_tensor) ?
         Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
         fe_eval.get_shape_info().dofs_per_component_on_cell;
+
     // expand dof_values to tensor product for truncated tensor products
     Number *values_dofs =
       (type == MatrixFreeFunctions::truncated_tensor) ?
@@ -582,28 +567,11 @@ namespace internal
       }
 
     if (type == MatrixFreeFunctions::truncated_tensor)
-      {
-        const std::size_t n_dofs_per_comp =
-          fe_eval.get_shape_info().dofs_per_component_on_cell;
-        values_dofs -= dofs_per_comp * n_components;
-        const int degree =
-          fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
-        for (unsigned int c = 0; c < n_components; ++c)
-          for (int i = 0, count_p = 0, count_q = 0;
-               i < (dim > 2 ? degree + 1 : 1);
-               ++i)
-            {
-              for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
-                {
-                  for (int k = 0; k < degree + 1 - j - i;
-                       ++k, ++count_p, ++count_q)
-                    values_dofs_actual[c * n_dofs_per_comp + count_p] =
-                      values_dofs[c * dofs_per_comp + count_q];
-                  count_q += j + i;
-                }
-              count_q += i * (degree + 1);
-            }
-      }
+      truncate_tensor_product_to_complete_degrees<dim, fe_degree>(
+        n_components,
+        values_dofs_actual,
+        values_dofs - dofs_per_comp * n_components,
+        fe_eval);
   }
 
 
diff --git a/include/deal.II/matrix_free/evaluation_kernels_common.h b/include/deal.II/matrix_free/evaluation_kernels_common.h
new file mode 100644 (file)
index 0000000..07a5d15
--- /dev/null
@@ -0,0 +1,101 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2025 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_matrix_free_evaluation_kernels_common_h
+#define dealii_matrix_free_evaluation_kernels_common_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/matrix_free/fe_evaluation_data.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+  template <int dim, int fe_degree, typename Number, bool is_face>
+  inline void
+  embed_truncated_into_full_tensor_product(
+    const unsigned int                      n_components,
+    Number                                 *values_dofs,
+    const Number                           *values_dofs_actual,
+    FEEvaluationData<dim, Number, is_face> &fe_eval)
+  {
+    const auto &shape_info = fe_eval.get_shape_info();
+    const auto &shape_data = shape_info.data.front();
+
+    const std::size_t dofs_per_comp =
+      Utilities::pow(shape_data.fe_degree + 1, dim);
+    const std::size_t n_dofs_per_comp = shape_info.dofs_per_component_on_cell;
+    const int degree = fe_degree != -1 ? fe_degree : shape_data.fe_degree;
+
+    for (unsigned int c = 0; c < n_components; ++c)
+      for (int i = 0, count_p = 0, count_q = 0; i < (dim > 2 ? degree + 1 : 1);
+           ++i)
+        {
+          for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
+            {
+              for (int k = 0; k < degree + 1 - j - i; ++k, ++count_p, ++count_q)
+                values_dofs[c * dofs_per_comp + count_q] =
+                  values_dofs_actual[c * n_dofs_per_comp + count_p];
+              for (int k = degree + 1 - j - i; k < degree + 1; ++k, ++count_q)
+                values_dofs[c * dofs_per_comp + count_q] = Number();
+            }
+          for (int j = degree + 1 - i; j < degree + 1; ++j)
+            for (int k = 0; k < degree + 1; ++k, ++count_q)
+              values_dofs[c * dofs_per_comp + count_q] = Number();
+        }
+  }
+
+
+
+  template <int dim, int fe_degree, typename Number, bool is_face>
+  inline void
+  truncate_tensor_product_to_complete_degrees(
+    const unsigned int                      n_components,
+    Number                                 *values_dofs_actual,
+    const Number                           *values_dofs,
+    FEEvaluationData<dim, Number, is_face> &fe_eval)
+  {
+    const auto &shape_info = fe_eval.get_shape_info();
+    const auto &shape_data = shape_info.data.front();
+
+    const unsigned int dofs_per_comp =
+      Utilities::fixed_power<dim>(shape_data.fe_degree + 1);
+    const std::size_t n_dofs_per_comp = shape_info.dofs_per_component_on_cell;
+    const int degree = fe_degree != -1 ? fe_degree : shape_data.fe_degree;
+    for (unsigned int c = 0; c < n_components; ++c)
+      for (int i = 0, count_p = 0, count_q = 0; i < (dim > 2 ? degree + 1 : 1);
+           ++i)
+        {
+          for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
+            {
+              for (int k = 0; k < degree + 1 - j - i; ++k, ++count_p, ++count_q)
+                values_dofs_actual[c * n_dofs_per_comp + count_p] =
+                  values_dofs[c * dofs_per_comp + count_q];
+              count_q += j + i;
+            }
+          count_q += i * (degree + 1);
+        }
+  }
+} // end of namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index e3f2e118bd9beb210511c2f3409f68eb38fe1e19..6214c6585cbf14f42b8a5b8a005e3b3f594e3c37 100644 (file)
@@ -25,6 +25,7 @@
 
 #include <deal.II/matrix_free/dof_info.h>
 #include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/evaluation_kernels_common.h>
 #include <deal.II/matrix_free/fe_evaluation_data.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
@@ -874,16 +875,20 @@ namespace internal
         interpolate_raviart_thomas<do_evaluate, add_into_output>(
           n_components, input, output, flags, face_no, shape_info);
       else
-        interpolate_generic<do_evaluate, add_into_output>(
-          n_components,
-          input,
-          output,
-          flags,
-          face_no,
-          shape_info.data.front().fe_degree + 1,
-          shape_info.data.front().shape_data_on_face,
-          shape_info.dofs_per_component_on_cell,
-          3 * shape_info.dofs_per_component_on_face);
+        {
+          const unsigned int fe_degree_ = shape_info.data.front().fe_degree;
+
+          interpolate_generic<do_evaluate, add_into_output>(
+            n_components,
+            input,
+            output,
+            flags,
+            face_no,
+            fe_degree_ + 1,
+            shape_info.data.front().shape_data_on_face,
+            Utilities::pow(fe_degree_ + 1, dim),
+            3 * Utilities::pow(fe_degree_ + 1, dim - 1));
+        }
     }
 
     /**
@@ -1751,7 +1756,7 @@ namespace internal
     static bool
     evaluate_tensor(const unsigned int                     n_components,
                     const EvaluationFlags::EvaluationFlags evaluation_flag,
-                    const Number                          *values_dofs,
+                    const Number                          *values_dofs_actual,
                     FEEvaluationData<dim, Number, true>   &fe_eval)
     {
       const auto &shape_info = fe_eval.get_shape_info();
@@ -1764,8 +1769,22 @@ namespace internal
 
       // Note: we always keep storage of values, 1st and 2nd derivatives in an
       // array, so reserve space for all three here
-      Number *temp         = fe_eval.get_scratch_data().begin();
-      Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+      Number *temp1 = fe_eval.get_scratch_data().begin();
+      Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face;
+
+      const Number *values_dofs =
+        (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ?
+          temp2 +
+            2 * (std::max(fe_eval.get_shape_info().dofs_per_component_on_cell,
+                          shape_info.n_q_points)) :
+          values_dofs_actual;
+
+      if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor)
+        embed_truncated_into_full_tensor_product<dim, fe_degree>(
+          n_components,
+          const_cast<Number *>(values_dofs),
+          values_dofs_actual,
+          fe_eval);
 
       bool use_vectorization = true;
       if (fe_eval.get_dof_access_index() ==
@@ -1781,15 +1800,15 @@ namespace internal
                                  values_dofs,
                                  fe_eval,
                                  use_vectorization,
-                                 temp,
-                                 scratch_data);
+                                 temp1,
+                                 temp2);
 
       evaluate_in_face<fe_degree, n_q_points_1d>(
-        n_components, evaluation_flag, fe_eval, temp, scratch_data);
+        n_components, evaluation_flag, fe_eval, temp1, temp2);
 
       if (dim == 3)
         adjust_quadrature_for_face_orientation(
-          n_components, evaluation_flag, fe_eval, use_vectorization, temp);
+          n_components, evaluation_flag, fe_eval, use_vectorization, temp1);
 
       return false;
     }
@@ -2290,7 +2309,7 @@ namespace internal
     static bool
     integrate_tensor(const unsigned int                     n_components,
                      const EvaluationFlags::EvaluationFlags integration_flag,
-                     Number                                *values_dofs,
+                     Number                                *values_dofs_actual,
                      FEEvaluationData<dim, Number, true>   &fe_eval,
                      const bool                             sum_into_values)
     {
@@ -2302,8 +2321,16 @@ namespace internal
           Utilities::pow(fe_degree + 1, dim - 1) :
           Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
 
-      Number *temp         = fe_eval.get_scratch_data().begin();
-      Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+      Number *temp1 = fe_eval.get_scratch_data().begin();
+      Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face;
+
+      // expand dof_values to tensor product for truncated tensor products
+      Number *values_dofs =
+        (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ?
+          temp2 + 2 * (std::max<std::size_t>(
+                        fe_eval.get_shape_info().dofs_per_component_on_cell,
+                        fe_eval.get_shape_info().n_q_points)) :
+          values_dofs_actual;
 
       bool use_vectorization = true;
 
@@ -2321,20 +2348,25 @@ namespace internal
 
       if (dim == 3)
         adjust_quadrature_for_face_orientation(
-          n_components, integration_flag, fe_eval, use_vectorization, temp);
+          n_components, integration_flag, fe_eval, use_vectorization, temp1);
 
       integrate_in_face<fe_degree, n_q_points_1d>(
-        n_components, integration_flag, fe_eval, temp, scratch_data);
+        n_components, integration_flag, fe_eval, temp1, temp2);
 
       collect_from_face<fe_degree>(n_components,
                                    integration_flag,
                                    values_dofs,
                                    fe_eval,
                                    use_vectorization,
-                                   temp,
-                                   scratch_data,
+                                   temp1,
+                                   temp2,
                                    sum_into_values);
 
+
+      if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor)
+        truncate_tensor_product_to_complete_degrees<dim, fe_degree>(
+          n_components, values_dofs_actual, values_dofs, fe_eval);
+
       return false;
     }
 
index 7a54ba8c560b82a73b34a52bbee6acaa20ea1432..517e6cd9a774350855439e70ea75521a99213a2d 100644 (file)
@@ -1201,7 +1201,7 @@ FEEvaluationData<dim, Number, is_face>::set_data_pointers(
 
   const unsigned int size_scratch_data =
     std::max(tensor_dofs_per_component + 1, dofs_per_component) * n_components *
-      3 +
+      4 +
     2 * n_quadrature_points;
   const unsigned int size_data_arrays =
     n_components * dofs_per_component +
diff --git a/tests/matrix_free/matrix_vector_faces_36.cc b/tests/matrix_free/matrix_vector_faces_36.cc
new file mode 100644 (file)
index 0000000..d562b65
--- /dev/null
@@ -0,0 +1,42 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2025 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.
+//
+// ------------------------------------------------------------------------
+
+
+
+// Similar to matrix_vector_faces_36 but with FE_DGP.
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/fe/fe_dgp.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_faces_common.h"
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(5 - dim);
+
+  FE_DGP<dim>     fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  do_test<dim, fe_degree, fe_degree + 1, double>(dof, constraints, true);
+}
diff --git a/tests/matrix_free/matrix_vector_faces_36.output b/tests/matrix_free/matrix_vector_faces_36.output
new file mode 100644 (file)
index 0000000..c4fe177
--- /dev/null
@@ -0,0 +1,17 @@
+
+DEAL:2d::Testing FE_DGP<2>(1)
+DEAL:2d::Norm of difference:          0
+DEAL:2d::Norm of difference parallel: 0
+DEAL:2d::
+DEAL:2d::Testing FE_DGP<2>(2)
+DEAL:2d::Norm of difference:          0
+DEAL:2d::Norm of difference parallel: 0
+DEAL:2d::
+DEAL:3d::Testing FE_DGP<3>(1)
+DEAL:3d::Norm of difference:          0
+DEAL:3d::Norm of difference parallel: 0
+DEAL:3d::
+DEAL:3d::Testing FE_DGP<3>(2)
+DEAL:3d::Norm of difference:          0
+DEAL:3d::Norm of difference parallel: 0
+DEAL:3d::

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.