From: Bruno Turcksin <bruno.turcksin@gmail.com>
Date: Thu, 16 Sep 2021 20:53:31 +0000 (+0000)
Subject: Refactor hanging_nodes_internals.h
X-Git-Tag: v9.4.0-rc1~983^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=19105e31ecb29751496ec7c59da0580fdc1ba5d0;p=dealii.git

Refactor hanging_nodes_internals.h
---

diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h
index 965ff1c27e..d913082f36 100644
--- a/include/deal.II/matrix_free/cuda_fe_evaluation.h
+++ b/include/deal.II/matrix_free/cuda_fe_evaluation.h
@@ -26,10 +26,10 @@
 #  include <deal.II/lac/cuda_atomic.h>
 #  include <deal.II/lac/cuda_vector.h>
 
+#  include <deal.II/matrix_free/cuda_hanging_nodes_internal.h>
 #  include <deal.II/matrix_free/cuda_matrix_free.h>
 #  include <deal.II/matrix_free/cuda_matrix_free.templates.h>
 #  include <deal.II/matrix_free/cuda_tensor_product_kernels.h>
-#  include <deal.II/matrix_free/hanging_nodes_internal.h>
 
 DEAL_II_NAMESPACE_OPEN
 
diff --git a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h
new file mode 100644
index 0000000000..d76f56fa25
--- /dev/null
+++ b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h
@@ -0,0 +1,349 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_cuda_hanging_nodes_internal_h
+#define dealii_cuda_hanging_nodes_internal_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_COMPILER_CUDA_AWARE
+
+#  include <deal.II/base/cuda_size.h>
+
+#  include <deal.II/matrix_free/hanging_nodes_internal.h>
+
+DEAL_II_NAMESPACE_OPEN
+namespace CUDAWrappers
+{
+  namespace internal
+  {
+    __constant__ double
+      constraint_weights[(CUDAWrappers::mf_max_elem_degree + 1) *
+                         (CUDAWrappers::mf_max_elem_degree + 1)];
+
+    //------------------------------------------------------------------------//
+    // Functions for resolving the hanging node constraints on the GPU        //
+    //------------------------------------------------------------------------//
+    template <unsigned int size>
+    __device__ inline unsigned int
+    index2(unsigned int i, unsigned int j)
+    {
+      return i + size * j;
+    }
+
+
+
+    template <unsigned int size>
+    __device__ inline unsigned int
+    index3(unsigned int i, unsigned int j, unsigned int k)
+    {
+      return i + size * j + size * size * k;
+    }
+
+
+
+    template <unsigned int fe_degree,
+              unsigned int direction,
+              bool         transpose,
+              typename Number>
+    __device__ inline void
+    interpolate_boundary_2d(
+      const dealii::internal::MatrixFreeFunctions::ConstraintKinds
+              constraint_mask,
+      Number *values)
+    {
+      const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
+      const unsigned int y_idx = threadIdx.y;
+
+      const auto this_type =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x :
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y;
+
+      const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
+
+      Number t = 0;
+      // Flag is true if dof is constrained for the given direction and the
+      // given face.
+      const bool constrained_face =
+        (constraint_mask &
+         (((direction == 0) ?
+             dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y :
+             dealii::internal::MatrixFreeFunctions::ConstraintKinds::
+               unconstrained) |
+          ((direction == 1) ?
+             dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
+             dealii::internal::MatrixFreeFunctions::ConstraintKinds::
+               unconstrained))) !=
+        dealii::internal::MatrixFreeFunctions::ConstraintKinds::unconstrained;
+
+      // Flag is true if for the given direction, the dof is constrained with
+      // the right type and is on the correct side (left (= 0) or right (=
+      // fe_degree))
+      const bool constrained_dof =
+        ((direction == 0) &&
+         (((constraint_mask &
+            dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y) !=
+           dealii::internal::MatrixFreeFunctions::ConstraintKinds::
+             unconstrained) ?
+            (y_idx == 0) :
+            (y_idx == fe_degree))) ||
+        ((direction == 1) &&
+         (((constraint_mask &
+            dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x) !=
+           dealii::internal::MatrixFreeFunctions::ConstraintKinds::
+             unconstrained) ?
+            (x_idx == 0) :
+            (x_idx == fe_degree)));
+
+      if (constrained_face && constrained_dof)
+        {
+          const bool type = (constraint_mask & this_type) !=
+                            dealii::internal::MatrixFreeFunctions::
+                              ConstraintKinds::unconstrained;
+
+          if (type)
+            {
+              for (unsigned int i = 0; i <= fe_degree; ++i)
+                {
+                  const unsigned int real_idx =
+                    (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
+                                       index2<fe_degree + 1>(x_idx, i);
+
+                  const Number w =
+                    transpose ?
+                      constraint_weights[i * (fe_degree + 1) + interp_idx] :
+                      constraint_weights[interp_idx * (fe_degree + 1) + i];
+                  t += w * values[real_idx];
+                }
+            }
+          else
+            {
+              for (unsigned int i = 0; i <= fe_degree; ++i)
+                {
+                  const unsigned int real_idx =
+                    (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
+                                       index2<fe_degree + 1>(x_idx, i);
+
+                  const Number w =
+                    transpose ?
+                      constraint_weights[(fe_degree - i) * (fe_degree + 1) +
+                                         fe_degree - interp_idx] :
+                      constraint_weights[(fe_degree - interp_idx) *
+                                           (fe_degree + 1) +
+                                         fe_degree - i];
+                  t += w * values[real_idx];
+                }
+            }
+        }
+
+      // The synchronization is done for all the threads in one block with
+      // each block being assigned to one element.
+      __syncthreads();
+      if (constrained_face && constrained_dof)
+        values[index2<fe_degree + 1>(x_idx, y_idx)] = t;
+
+      __syncthreads();
+    }
+
+
+
+    template <unsigned int fe_degree,
+              unsigned int direction,
+              bool         transpose,
+              typename Number>
+    __device__ inline void
+    interpolate_boundary_3d(
+      const dealii::internal::MatrixFreeFunctions::ConstraintKinds
+              constraint_mask,
+      Number *values)
+    {
+      const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
+      const unsigned int y_idx = threadIdx.y;
+      const unsigned int z_idx = threadIdx.z;
+
+      const auto this_type =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x :
+        (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y :
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z;
+      const auto face1_type =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y :
+        (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z :
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x;
+      const auto face2_type =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z :
+        (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x :
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y;
+
+      // If computing in x-direction, need to match against face_y or
+      // face_z
+      const auto face1 =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y :
+        (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x;
+      const auto face2 =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
+        (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y;
+      const auto edge =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_yz :
+        (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_zx :
+          dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_xy;
+      const auto constrained_face = constraint_mask & (face1 | face2 | edge);
+
+      const unsigned int interp_idx = (direction == 0) ? x_idx :
+                                      (direction == 1) ? y_idx :
+                                                         z_idx;
+      const unsigned int face1_idx  = (direction == 0) ? y_idx :
+                                      (direction == 1) ? z_idx :
+                                                         x_idx;
+      const unsigned int face2_idx  = (direction == 0) ? z_idx :
+                                      (direction == 1) ? x_idx :
+                                                         y_idx;
+
+      Number     t        = 0;
+      const bool on_face1 = ((constraint_mask & face1_type) !=
+                             dealii::internal::MatrixFreeFunctions::
+                               ConstraintKinds::unconstrained) ?
+                              (face1_idx == 0) :
+                              (face1_idx == fe_degree);
+      const bool on_face2 = ((constraint_mask & face2_type) !=
+                             dealii::internal::MatrixFreeFunctions::
+                               ConstraintKinds::unconstrained) ?
+                              (face2_idx == 0) :
+                              (face2_idx == fe_degree);
+      const bool constrained_dof =
+        ((((constraint_mask & face1) != dealii::internal::MatrixFreeFunctions::
+                                          ConstraintKinds::unconstrained) &&
+          on_face1) ||
+         (((constraint_mask & face2) != dealii::internal::MatrixFreeFunctions::
+                                          ConstraintKinds::unconstrained) &&
+          on_face2) ||
+         (((constraint_mask & edge) != dealii::internal::MatrixFreeFunctions::
+                                         ConstraintKinds::unconstrained) &&
+          on_face1 && on_face2));
+
+      if ((constrained_face != dealii::internal::MatrixFreeFunctions::
+                                 ConstraintKinds::unconstrained) &&
+          constrained_dof)
+        {
+          const bool type = (constraint_mask & this_type) !=
+                            dealii::internal::MatrixFreeFunctions::
+                              ConstraintKinds::unconstrained;
+          if (type)
+            {
+              for (unsigned int i = 0; i <= fe_degree; ++i)
+                {
+                  const unsigned int real_idx =
+                    (direction == 0) ? index3<fe_degree + 1>(i, y_idx, z_idx) :
+                    (direction == 1) ? index3<fe_degree + 1>(x_idx, i, z_idx) :
+                                       index3<fe_degree + 1>(x_idx, y_idx, i);
+
+                  const Number w =
+                    transpose ?
+                      constraint_weights[i * (fe_degree + 1) + interp_idx] :
+                      constraint_weights[interp_idx * (fe_degree + 1) + i];
+                  t += w * values[real_idx];
+                }
+            }
+          else
+            {
+              for (unsigned int i = 0; i <= fe_degree; ++i)
+                {
+                  const unsigned int real_idx =
+                    (direction == 0) ? index3<fe_degree + 1>(i, y_idx, z_idx) :
+                    (direction == 1) ? index3<fe_degree + 1>(x_idx, i, z_idx) :
+                                       index3<fe_degree + 1>(x_idx, y_idx, i);
+
+                  const Number w =
+                    transpose ?
+                      constraint_weights[(fe_degree - i) * (fe_degree + 1) +
+                                         fe_degree - interp_idx] :
+                      constraint_weights[(fe_degree - interp_idx) *
+                                           (fe_degree + 1) +
+                                         fe_degree - i];
+                  t += w * values[real_idx];
+                }
+            }
+        }
+
+      // The synchronization is done for all the threads in one block with
+      // each block being assigned to one element.
+      __syncthreads();
+
+      if ((constrained_face != dealii::internal::MatrixFreeFunctions::
+                                 ConstraintKinds::unconstrained) &&
+          constrained_dof)
+        values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = t;
+
+      __syncthreads();
+    }
+
+
+
+    /**
+     * This function resolves the hanging nodes using tensor product.
+     *
+     * The implementation of this class is explained in Section 3 of
+     * @cite ljungkvist2017matrix and in Section 3.4 of
+     * @cite kronbichler2019multigrid.
+     */
+    template <int dim, int fe_degree, bool transpose, typename Number>
+    __device__ void
+    resolve_hanging_nodes(
+      const dealii::internal::MatrixFreeFunctions::ConstraintKinds
+              constraint_mask,
+      Number *values)
+    {
+      if (dim == 2)
+        {
+          interpolate_boundary_2d<fe_degree, 0, transpose>(constraint_mask,
+                                                           values);
+
+          interpolate_boundary_2d<fe_degree, 1, transpose>(constraint_mask,
+                                                           values);
+        }
+      else if (dim == 3)
+        {
+          // Interpolate y and z faces (x-direction)
+          interpolate_boundary_3d<fe_degree, 0, transpose>(constraint_mask,
+                                                           values);
+          // Interpolate x and z faces (y-direction)
+          interpolate_boundary_3d<fe_degree, 1, transpose>(constraint_mask,
+                                                           values);
+          // Interpolate x and y faces (z-direction)
+          interpolate_boundary_3d<fe_degree, 2, transpose>(constraint_mask,
+                                                           values);
+        }
+    }
+  } // namespace internal
+} // namespace CUDAWrappers
+
+DEAL_II_NAMESPACE_CLOSE
+#endif
+
+#endif
diff --git a/include/deal.II/matrix_free/cuda_matrix_free.h b/include/deal.II/matrix_free/cuda_matrix_free.h
index 2fa34f40c1..31f322fa18 100644
--- a/include/deal.II/matrix_free/cuda_matrix_free.h
+++ b/include/deal.II/matrix_free/cuda_matrix_free.h
@@ -38,11 +38,19 @@
 #  include <deal.II/lac/cuda_vector.h>
 #  include <deal.II/lac/la_parallel_vector.h>
 
-#  include <deal.II/matrix_free/hanging_nodes_internal.h>
 
 
 DEAL_II_NAMESPACE_OPEN
 
+// Forward declaration
+namespace internal
+{
+  namespace MatrixFreeFunctions
+  {
+    enum class ConstraintKinds : std::uint16_t;
+  }
+} // namespace internal
+
 namespace CUDAWrappers
 {
   // forward declaration
diff --git a/include/deal.II/matrix_free/cuda_matrix_free.templates.h b/include/deal.II/matrix_free/cuda_matrix_free.templates.h
index 1bd8d75285..5039981471 100644
--- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h
+++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h
@@ -23,6 +23,7 @@
 
 #ifdef DEAL_II_COMPILER_CUDA_AWARE
 
+#  include <deal.II/base/cuda.h>
 #  include <deal.II/base/cuda_size.h>
 #  include <deal.II/base/graph_coloring.h>
 
@@ -31,6 +32,7 @@
 #  include <deal.II/fe/fe_dgq.h>
 #  include <deal.II/fe/fe_values.h>
 
+#  include <deal.II/matrix_free/cuda_hanging_nodes_internal.h>
 #  include <deal.II/matrix_free/shape_info.h>
 
 #  include <cuda_runtime_api.h>
diff --git a/include/deal.II/matrix_free/hanging_nodes_internal.h b/include/deal.II/matrix_free/hanging_nodes_internal.h
index f13aee5beb..b92fcf2316 100644
--- a/include/deal.II/matrix_free/hanging_nodes_internal.h
+++ b/include/deal.II/matrix_free/hanging_nodes_internal.h
@@ -18,7 +18,6 @@
 
 #include <deal.II/base/config.h>
 
-#include <deal.II/base/cuda_size.h>
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/dofs/dof_accessor.h>
@@ -855,327 +854,6 @@ namespace internal
 } // namespace internal
 
 
-
-#ifdef DEAL_II_COMPILER_CUDA_AWARE
-namespace CUDAWrappers
-{
-  namespace internal
-  {
-    __constant__ double
-      constraint_weights[(CUDAWrappers::mf_max_elem_degree + 1) *
-                         (CUDAWrappers::mf_max_elem_degree + 1)];
-
-    //------------------------------------------------------------------------//
-    // Functions for resolving the hanging node constraints on the GPU        //
-    //------------------------------------------------------------------------//
-    template <unsigned int size>
-    __device__ inline unsigned int
-    index2(unsigned int i, unsigned int j)
-    {
-      return i + size * j;
-    }
-
-
-
-    template <unsigned int size>
-    __device__ inline unsigned int
-    index3(unsigned int i, unsigned int j, unsigned int k)
-    {
-      return i + size * j + size * size * k;
-    }
-
-
-
-    template <unsigned int fe_degree,
-              unsigned int direction,
-              bool         transpose,
-              typename Number>
-    __device__ inline void
-    interpolate_boundary_2d(
-      const dealii::internal::MatrixFreeFunctions::ConstraintKinds
-              constraint_mask,
-      Number *values)
-    {
-      const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
-      const unsigned int y_idx = threadIdx.y;
-
-      const auto this_type =
-        (direction == 0) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x :
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y;
-
-      const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
-
-      Number t = 0;
-      // Flag is true if dof is constrained for the given direction and the
-      // given face.
-      const bool constrained_face =
-        (constraint_mask &
-         (((direction == 0) ?
-             dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y :
-             dealii::internal::MatrixFreeFunctions::ConstraintKinds::
-               unconstrained) |
-          ((direction == 1) ?
-             dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
-             dealii::internal::MatrixFreeFunctions::ConstraintKinds::
-               unconstrained))) !=
-        dealii::internal::MatrixFreeFunctions::ConstraintKinds::unconstrained;
-
-      // Flag is true if for the given direction, the dof is constrained with
-      // the right type and is on the correct side (left (= 0) or right (=
-      // fe_degree))
-      const bool constrained_dof =
-        ((direction == 0) &&
-         (((constraint_mask &
-            dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y) !=
-           dealii::internal::MatrixFreeFunctions::ConstraintKinds::
-             unconstrained) ?
-            (y_idx == 0) :
-            (y_idx == fe_degree))) ||
-        ((direction == 1) &&
-         (((constraint_mask &
-            dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x) !=
-           dealii::internal::MatrixFreeFunctions::ConstraintKinds::
-             unconstrained) ?
-            (x_idx == 0) :
-            (x_idx == fe_degree)));
-
-      if (constrained_face && constrained_dof)
-        {
-          const bool type = (constraint_mask & this_type) !=
-                            dealii::internal::MatrixFreeFunctions::
-                              ConstraintKinds::unconstrained;
-
-          if (type)
-            {
-              for (unsigned int i = 0; i <= fe_degree; ++i)
-                {
-                  const unsigned int real_idx =
-                    (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
-                                       index2<fe_degree + 1>(x_idx, i);
-
-                  const Number w =
-                    transpose ?
-                      constraint_weights[i * (fe_degree + 1) + interp_idx] :
-                      constraint_weights[interp_idx * (fe_degree + 1) + i];
-                  t += w * values[real_idx];
-                }
-            }
-          else
-            {
-              for (unsigned int i = 0; i <= fe_degree; ++i)
-                {
-                  const unsigned int real_idx =
-                    (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
-                                       index2<fe_degree + 1>(x_idx, i);
-
-                  const Number w =
-                    transpose ?
-                      constraint_weights[(fe_degree - i) * (fe_degree + 1) +
-                                         fe_degree - interp_idx] :
-                      constraint_weights[(fe_degree - interp_idx) *
-                                           (fe_degree + 1) +
-                                         fe_degree - i];
-                  t += w * values[real_idx];
-                }
-            }
-        }
-
-      // The synchronization is done for all the threads in one block with
-      // each block being assigned to one element.
-      __syncthreads();
-      if (constrained_face && constrained_dof)
-        values[index2<fe_degree + 1>(x_idx, y_idx)] = t;
-
-      __syncthreads();
-    }
-
-
-
-    template <unsigned int fe_degree,
-              unsigned int direction,
-              bool         transpose,
-              typename Number>
-    __device__ inline void
-    interpolate_boundary_3d(
-      const dealii::internal::MatrixFreeFunctions::ConstraintKinds
-              constraint_mask,
-      Number *values)
-    {
-      const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
-      const unsigned int y_idx = threadIdx.y;
-      const unsigned int z_idx = threadIdx.z;
-
-      const auto this_type =
-        (direction == 0) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x :
-        (direction == 1) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y :
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z;
-      const auto face1_type =
-        (direction == 0) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y :
-        (direction == 1) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z :
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x;
-      const auto face2_type =
-        (direction == 0) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z :
-        (direction == 1) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x :
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y;
-
-      // If computing in x-direction, need to match against face_y or
-      // face_z
-      const auto face1 =
-        (direction == 0) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y :
-        (direction == 1) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x;
-      const auto face2 =
-        (direction == 0) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
-        (direction == 1) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y;
-      const auto edge =
-        (direction == 0) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_yz :
-        (direction == 1) ?
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_zx :
-          dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_xy;
-      const auto constrained_face = constraint_mask & (face1 | face2 | edge);
-
-      const unsigned int interp_idx = (direction == 0) ? x_idx :
-                                      (direction == 1) ? y_idx :
-                                                         z_idx;
-      const unsigned int face1_idx  = (direction == 0) ? y_idx :
-                                      (direction == 1) ? z_idx :
-                                                         x_idx;
-      const unsigned int face2_idx  = (direction == 0) ? z_idx :
-                                      (direction == 1) ? x_idx :
-                                                         y_idx;
-
-      Number     t        = 0;
-      const bool on_face1 = ((constraint_mask & face1_type) !=
-                             dealii::internal::MatrixFreeFunctions::
-                               ConstraintKinds::unconstrained) ?
-                              (face1_idx == 0) :
-                              (face1_idx == fe_degree);
-      const bool on_face2 = ((constraint_mask & face2_type) !=
-                             dealii::internal::MatrixFreeFunctions::
-                               ConstraintKinds::unconstrained) ?
-                              (face2_idx == 0) :
-                              (face2_idx == fe_degree);
-      const bool constrained_dof =
-        ((((constraint_mask & face1) != dealii::internal::MatrixFreeFunctions::
-                                          ConstraintKinds::unconstrained) &&
-          on_face1) ||
-         (((constraint_mask & face2) != dealii::internal::MatrixFreeFunctions::
-                                          ConstraintKinds::unconstrained) &&
-          on_face2) ||
-         (((constraint_mask & edge) != dealii::internal::MatrixFreeFunctions::
-                                         ConstraintKinds::unconstrained) &&
-          on_face1 && on_face2));
-
-      if ((constrained_face != dealii::internal::MatrixFreeFunctions::
-                                 ConstraintKinds::unconstrained) &&
-          constrained_dof)
-        {
-          const bool type = (constraint_mask & this_type) !=
-                            dealii::internal::MatrixFreeFunctions::
-                              ConstraintKinds::unconstrained;
-          if (type)
-            {
-              for (unsigned int i = 0; i <= fe_degree; ++i)
-                {
-                  const unsigned int real_idx =
-                    (direction == 0) ? index3<fe_degree + 1>(i, y_idx, z_idx) :
-                    (direction == 1) ? index3<fe_degree + 1>(x_idx, i, z_idx) :
-                                       index3<fe_degree + 1>(x_idx, y_idx, i);
-
-                  const Number w =
-                    transpose ?
-                      constraint_weights[i * (fe_degree + 1) + interp_idx] :
-                      constraint_weights[interp_idx * (fe_degree + 1) + i];
-                  t += w * values[real_idx];
-                }
-            }
-          else
-            {
-              for (unsigned int i = 0; i <= fe_degree; ++i)
-                {
-                  const unsigned int real_idx =
-                    (direction == 0) ? index3<fe_degree + 1>(i, y_idx, z_idx) :
-                    (direction == 1) ? index3<fe_degree + 1>(x_idx, i, z_idx) :
-                                       index3<fe_degree + 1>(x_idx, y_idx, i);
-
-                  const Number w =
-                    transpose ?
-                      constraint_weights[(fe_degree - i) * (fe_degree + 1) +
-                                         fe_degree - interp_idx] :
-                      constraint_weights[(fe_degree - interp_idx) *
-                                           (fe_degree + 1) +
-                                         fe_degree - i];
-                  t += w * values[real_idx];
-                }
-            }
-        }
-
-      // The synchronization is done for all the threads in one block with
-      // each block being assigned to one element.
-      __syncthreads();
-
-      if ((constrained_face != dealii::internal::MatrixFreeFunctions::
-                                 ConstraintKinds::unconstrained) &&
-          constrained_dof)
-        values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = t;
-
-      __syncthreads();
-    }
-
-
-
-    /**
-     * This function resolves the hanging nodes using tensor product.
-     *
-     * The implementation of this class is explained in Section 3 of
-     * @cite ljungkvist2017matrix and in Section 3.4 of
-     * @cite kronbichler2019multigrid.
-     */
-    template <int dim, int fe_degree, bool transpose, typename Number>
-    __device__ void
-    resolve_hanging_nodes(
-      const dealii::internal::MatrixFreeFunctions::ConstraintKinds
-              constraint_mask,
-      Number *values)
-    {
-      if (dim == 2)
-        {
-          interpolate_boundary_2d<fe_degree, 0, transpose>(constraint_mask,
-                                                           values);
-
-          interpolate_boundary_2d<fe_degree, 1, transpose>(constraint_mask,
-                                                           values);
-        }
-      else if (dim == 3)
-        {
-          // Interpolate y and z faces (x-direction)
-          interpolate_boundary_3d<fe_degree, 0, transpose>(constraint_mask,
-                                                           values);
-          // Interpolate x and z faces (y-direction)
-          interpolate_boundary_3d<fe_degree, 1, transpose>(constraint_mask,
-                                                           values);
-          // Interpolate x and y faces (z-direction)
-          interpolate_boundary_3d<fe_degree, 2, transpose>(constraint_mask,
-                                                           values);
-        }
-    }
-  } // namespace internal
-} // namespace CUDAWrappers
-#endif
-
 DEAL_II_NAMESPACE_CLOSE
 
 #endif