]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Split up evaluation_template_factory 12778/head
authorPeter Munch <peterrmuench@gmail.com>
Tue, 28 Sep 2021 10:51:56 +0000 (12:51 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Tue, 28 Sep 2021 12:54:17 +0000 (14:54 +0200)
16 files changed:
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h [new file with mode: 0644]
include/deal.II/matrix_free/evaluation_selector.h
include/deal.II/matrix_free/evaluation_template_factory.templates.h
include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h [new file with mode: 0644]
include/deal.II/matrix_free/evaluation_template_factory_internal.h [new file with mode: 0644]
source/matrix_free/CMakeLists.txt
source/matrix_free/evaluation_template_factory.inst.in
source/matrix_free/evaluation_template_factory_hanging_nodes.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_factory_hanging_nodes.inst.in [new file with mode: 0644]
source/matrix_free/evaluation_template_factory_hanging_nodes_inst2.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_factory_hanging_nodes_inst3.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_factory_hanging_nodes_inst4.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_factory_hanging_nodes_inst5.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_factory_hanging_nodes_inst6.cc [new file with mode: 0644]
tests/matrix_free/hanging_node_kernels_01.cc

index 2fa379d2b9ad9ddb07b05eac377aaf5a91c8eb17..215a3ecd48266ceef43a0a27d99c8c4c5fbe3a2c 100644 (file)
@@ -4906,587 +4906,6 @@ namespace internal
     }
   };
 
-
-
-  template <int dim, typename Number, bool is_face>
-  struct FEEvaluationImplHangingNodes
-  {
-    template <int fe_degree, int n_q_points_1d>
-    static bool
-    run(const unsigned int                  n_desired_components,
-        const FEEvaluationBaseData<dim,
-                                   typename Number::value_type,
-                                   is_face,
-                                   Number> &fe_eval,
-        const bool                          transpose,
-        const std::array<MatrixFreeFunctions::ConstraintKinds, Number::size()>
-          &     c_mask,
-        Number *values)
-    {
-      if (transpose)
-        run_internal<fe_degree, true>(n_desired_components,
-                                      fe_eval,
-                                      c_mask,
-                                      values);
-      else
-        run_internal<fe_degree, false>(n_desired_components,
-                                       fe_eval,
-                                       c_mask,
-                                       values);
-
-      return false;
-    }
-
-  private:
-    template <int fe_degree, unsigned int side, bool transpose>
-    static void
-    interpolate_2D(const unsigned int given_degree,
-                   const unsigned int v,
-                   const Number *     weight,
-                   Number *           values)
-    {
-      typename Number::value_type temp[40];
-
-      const unsigned int points =
-        (fe_degree != -1 ? fe_degree : given_degree) + 1;
-
-      AssertIndexRange(points, 40);
-
-      const unsigned int d = side / 2; // direction
-      const unsigned int s = side % 2; // left or right surface
-
-      const unsigned int offset = dealii::Utilities::pow(points, d + 1);
-      const unsigned int stride =
-        (s == 0 ? 0 : (points - 1)) * dealii::Utilities::pow(points, d);
-
-      const unsigned int r1 = dealii::Utilities::pow(points, dim - d - 1);
-      const unsigned int r2 = dealii::Utilities::pow(points, d);
-
-      // copy result back
-      for (unsigned int i = 0, k = 0; i < r1; ++i)
-        for (unsigned int j = 0; j < r2; ++j, ++k)
-          temp[k] = values[i * offset + stride + j][v];
-
-      // perform interpolation point by point (note: r1 * r2 == points^(dim-1))
-      for (unsigned int i = 0, k = 0; i < r1; ++i)
-        for (unsigned int j = 0; j < r2; ++j, ++k)
-          {
-            typename Number::value_type sum = 0.0;
-            for (unsigned int h = 0; h < points; ++h)
-              sum += weight[(transpose ? 1 : points) * k +
-                            (transpose ? points : 1) * h][v] *
-                     temp[h];
-            values[i * offset + stride + j][v] = sum;
-          }
-    }
-
-    template <int          fe_degree,
-              unsigned int direction,
-              unsigned int d,
-              bool         transpose>
-    static void
-    interpolate_3D_face(const unsigned int dof_offset,
-                        const unsigned int given_degree,
-                        const unsigned int v,
-                        const Number *     weight,
-                        Number *           values)
-    {
-      typename Number::value_type temp[40];
-
-      const unsigned int points =
-        (fe_degree != -1 ? fe_degree : given_degree) + 1;
-
-      AssertIndexRange(points, 40);
-
-      const unsigned int stride = Utilities::pow(points, direction);
-
-      // direction   side0   side1   side2
-      // 0             -      p^2      p
-      // 1            p^2      -       1
-      // 2             p       -       1
-      const unsigned int stride2 =
-        ((direction == 0 && d == 1) || (direction == 1 && d == 0)) ?
-          (points * points) :
-          (((direction == 0 && d == 2) || (direction == 2 && d == 0)) ? points :
-                                                                        1);
-
-      for (unsigned int g = 1; g < points - 1; ++g)
-        {
-          // copy result back
-          for (unsigned int k = 0; k < points; ++k)
-            temp[k] = values[dof_offset + k * stride + stride2 * g][v];
-
-          // perform interpolation point by point
-          for (unsigned int k = 0; k < points; ++k)
-            {
-              typename Number::value_type sum = 0.0;
-              for (unsigned int h = 0; h < points; ++h)
-                sum += weight[(transpose ? 1 : points) * k +
-                              (transpose ? points : 1) * h][v] *
-                       temp[h];
-              values[dof_offset + k * stride + stride2 * g][v] = sum;
-            }
-        }
-    }
-
-    template <int fe_degree, unsigned int direction, bool transpose>
-    static void
-    interpolate_3D_edge(const unsigned int p,
-                        const unsigned int given_degree,
-                        const unsigned int v,
-                        const Number *     weight,
-                        Number *           values)
-    {
-      typename Number::value_type temp[40];
-
-      const unsigned int points =
-        (fe_degree != -1 ? fe_degree : given_degree) + 1;
-
-      AssertIndexRange(points, 40);
-
-      const unsigned int stride = Utilities::pow(points, direction);
-
-      // copy result back
-      for (unsigned int k = 0; k < points; ++k)
-        temp[k] = values[p + k * stride][v];
-
-      // perform interpolation point by point
-      for (unsigned int k = 0; k < points; ++k)
-        {
-          typename Number::value_type sum = 0.0;
-          for (unsigned int h = 0; h < points; ++h)
-            sum += weight[(transpose ? 1 : points) * k +
-                          (transpose ? points : 1) * h][v] *
-                   temp[h];
-          values[p + k * stride][v] = sum;
-        }
-    }
-
-    template <int fe_degree, bool transpose>
-    static void
-    run_internal(const unsigned int                  n_desired_components,
-                 const FEEvaluationBaseData<dim,
-                                            typename Number::value_type,
-                                            is_face,
-                                            Number> &fe_eval,
-                 const std::array<MatrixFreeFunctions::ConstraintKinds,
-                                  Number::size()> &  constraint_mask,
-                 Number *                            values)
-    {
-      const unsigned int given_degree =
-        fe_degree != -1 ? fe_degree :
-                          fe_eval.get_shape_info().data.front().fe_degree;
-
-      const auto &interpolation_matrices =
-        fe_eval.get_shape_info().data.front().subface_interpolation_matrices;
-
-      const auto is_set = [](const auto a, const auto b) -> bool {
-        return (a & b) == b;
-      };
-
-      const auto not_set = [](const auto a, const auto b) -> bool {
-        return (a & b) == MatrixFreeFunctions::ConstraintKinds::unconstrained;
-      };
-
-      const unsigned int points = given_degree + 1;
-
-      for (unsigned int c = 0; c < n_desired_components; ++c)
-        {
-          for (unsigned int v = 0; v < Number::size(); ++v)
-            {
-              const auto mask = constraint_mask[v];
-
-              if (mask == MatrixFreeFunctions::ConstraintKinds::unconstrained)
-                continue;
-
-              if (dim == 2) // 2D: only faces
-                {
-                  // direction 0:
-                  if ((mask & MatrixFreeFunctions::ConstraintKinds::face_y) !=
-                      MatrixFreeFunctions::ConstraintKinds::unconstrained)
-                    {
-                      const bool is_subface_0 =
-                        (mask & MatrixFreeFunctions::ConstraintKinds::type_x) ==
-                        MatrixFreeFunctions::ConstraintKinds::unconstrained;
-
-                      const Number *weights =
-                        interpolation_matrices[is_subface_0].data();
-
-                      if (is_set(mask,
-                                 MatrixFreeFunctions::ConstraintKinds::type_y))
-                        interpolate_2D<fe_degree, 2, transpose>(
-                          given_degree,
-                          v,
-                          weights,
-                          values); // face 2
-                      else
-                        interpolate_2D<fe_degree, 3, transpose>(
-                          given_degree,
-                          v,
-                          weights,
-                          values); // face 3
-                    }
-
-                  // direction 1:
-                  if ((mask & MatrixFreeFunctions::ConstraintKinds::face_x) !=
-                      MatrixFreeFunctions::ConstraintKinds::unconstrained)
-                    {
-                      const bool is_subface_0 =
-                        (mask & MatrixFreeFunctions::ConstraintKinds::type_y) ==
-                        MatrixFreeFunctions::ConstraintKinds::unconstrained;
-
-                      const Number *weights =
-                        interpolation_matrices[is_subface_0].data();
-
-                      if (is_set(mask,
-                                 MatrixFreeFunctions::ConstraintKinds::type_x))
-                        interpolate_2D<fe_degree, 0, transpose>(
-                          given_degree,
-                          v,
-                          weights,
-                          values); // face 0
-                      else
-                        interpolate_2D<fe_degree, 1, transpose>(
-                          given_degree,
-                          v,
-                          weights,
-                          values); // face 1
-                    }
-                }
-              else if (dim == 3) // 3D faces and edges
-                {
-                  const unsigned int p0 = 0;
-                  const unsigned int p1 = points - 1;
-                  const unsigned int p2 = points * points - points;
-                  const unsigned int p3 = points * points - 1;
-                  const unsigned int p4 =
-                    points * points * points - points * points;
-                  const unsigned int p5 =
-                    points * points * points - points * points + points - 1;
-                  const unsigned int p6 = points * points * points - points;
-
-                  const bool is_face_0 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::face_x) &&
-                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_x);
-                  const bool is_face_1 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::face_x) &&
-                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_x);
-                  const bool is_face_2 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::face_y) &&
-                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
-                  const bool is_face_3 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::face_y) &&
-                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
-                  const bool is_face_4 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::face_z) &&
-                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-                  const bool is_face_5 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::face_z) &&
-                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-
-                  const bool is_edge_2 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_yz) &&
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::type_y) &&
-                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-                  const bool is_edge_3 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_yz) &&
-                    not_set(mask,
-                            MatrixFreeFunctions::ConstraintKinds::type_y) &&
-                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-                  const bool is_edge_6 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_yz) &&
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::type_y) &&
-                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-                  const bool is_edge_7 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_yz) &&
-                    not_set(mask,
-                            MatrixFreeFunctions::ConstraintKinds::type_y) &&
-                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-
-                  const bool is_edge_0 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_zx) &&
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::type_x) &&
-                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-                  const bool is_edge_1 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_zx) &&
-                    not_set(mask,
-                            MatrixFreeFunctions::ConstraintKinds::type_x) &&
-                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-                  const bool is_edge_4 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_zx) &&
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::type_x) &&
-                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-                  const bool is_edge_5 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_zx) &&
-                    not_set(mask,
-                            MatrixFreeFunctions::ConstraintKinds::type_x) &&
-                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
-
-                  const bool is_edge_8 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_xy) &&
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::type_x) &&
-                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
-                  const bool is_edge_9 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_xy) &&
-                    not_set(mask,
-                            MatrixFreeFunctions::ConstraintKinds::type_x) &&
-                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
-                  const bool is_edge_10 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_xy) &&
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::type_x) &&
-                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
-                  const bool is_edge_11 =
-                    is_set(mask,
-                           MatrixFreeFunctions::ConstraintKinds::edge_xy) &&
-                    not_set(mask,
-                            MatrixFreeFunctions::ConstraintKinds::type_x) &&
-                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
-
-                  // direction 0:
-                  {
-                    const bool is_subface_0 =
-                      (mask & MatrixFreeFunctions::ConstraintKinds::type_x) ==
-                      MatrixFreeFunctions::ConstraintKinds::unconstrained;
-
-                    const Number *weights =
-                      interpolation_matrices[is_subface_0].data();
-
-                    // ... faces
-                    if (is_face_2)
-                      interpolate_3D_face<fe_degree, 0, 1, transpose>(
-                        p0,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 2
-                    else if (is_face_3)
-                      interpolate_3D_face<fe_degree, 0, 1, transpose>(
-                        p2,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 3
-                    if (is_face_4)
-                      interpolate_3D_face<fe_degree, 0, 2, transpose>(
-                        p0,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 4
-                    else if (is_face_5)
-                      interpolate_3D_face<fe_degree, 0, 2, transpose>(
-                        p4,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 5
-
-                    // ... edges
-                    if (is_face_2 || is_face_4 || is_edge_2)
-                      interpolate_3D_edge<fe_degree, 0, transpose>(
-                        p0,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 2
-                    if (is_face_3 || is_face_4 || is_edge_3)
-                      interpolate_3D_edge<fe_degree, 0, transpose>(
-                        p2,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 3
-                    if (is_face_2 || is_face_5 || is_edge_6)
-                      interpolate_3D_edge<fe_degree, 0, transpose>(
-                        p4,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 6
-                    if (is_face_3 || is_face_5 || is_edge_7)
-                      interpolate_3D_edge<fe_degree, 0, transpose>(
-                        p6,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 7
-                  }
-
-                  // direction 1:
-                  {
-                    const bool is_subface_0 =
-                      (mask & MatrixFreeFunctions::ConstraintKinds::type_y) ==
-                      MatrixFreeFunctions::ConstraintKinds::unconstrained;
-
-                    const Number *weights =
-                      interpolation_matrices[is_subface_0].data();
-
-                    // ... faces
-                    if (is_face_0)
-                      interpolate_3D_face<fe_degree, 1, 0, transpose>(
-                        p0,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 0
-                    else if (is_face_1)
-                      interpolate_3D_face<fe_degree, 1, 0, transpose>(
-                        p1,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 1
-                    if (is_face_4)
-                      interpolate_3D_face<fe_degree, 1, 2, transpose>(
-                        p0,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 4
-                    else if (is_face_5)
-                      interpolate_3D_face<fe_degree, 1, 2, transpose>(
-                        p4,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 5
-
-                    // ... edges
-                    if (is_face_0 || is_face_4 || is_edge_0)
-                      interpolate_3D_edge<fe_degree, 1, transpose>(
-                        p0,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 0
-                    if (is_face_1 || is_face_4 || is_edge_1)
-                      interpolate_3D_edge<fe_degree, 1, transpose>(
-                        p1,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 1
-                    if (is_face_0 || is_face_5 || is_edge_4)
-                      interpolate_3D_edge<fe_degree, 1, transpose>(
-                        p4,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 4
-                    if (is_face_1 || is_face_5 || is_edge_5)
-                      interpolate_3D_edge<fe_degree, 1, transpose>(
-                        p5,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 5
-                  }
-
-                  // direction 2:
-                  {
-                    const bool is_subface_0 =
-                      (mask & MatrixFreeFunctions::ConstraintKinds::type_z) ==
-                      MatrixFreeFunctions::ConstraintKinds::unconstrained;
-
-                    const Number *weights =
-                      interpolation_matrices[is_subface_0].data();
-
-                    // ... faces
-                    if (is_face_0)
-                      interpolate_3D_face<fe_degree, 2, 0, transpose>(
-                        p0,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 0
-                    else if (is_face_1)
-                      interpolate_3D_face<fe_degree, 2, 0, transpose>(
-                        p1,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 1
-                    if (is_face_2)
-                      interpolate_3D_face<fe_degree, 2, 1, transpose>(
-                        p0,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 2
-                    else if (is_face_3)
-                      interpolate_3D_face<fe_degree, 2, 1, transpose>(
-                        p2,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // face 3
-
-                    // ... edges
-                    if (is_face_0 || is_face_2 || is_edge_8)
-                      interpolate_3D_edge<fe_degree, 2, transpose>(
-                        p0,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 8
-                    if (is_face_1 || is_face_2 || is_edge_9)
-                      interpolate_3D_edge<fe_degree, 2, transpose>(
-                        p1,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 9
-                    if (is_face_0 || is_face_3 || is_edge_10)
-                      interpolate_3D_edge<fe_degree, 2, transpose>(
-                        p2,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 10
-                    if (is_face_1 || is_face_3 || is_edge_11)
-                      interpolate_3D_edge<fe_degree, 2, transpose>(
-                        p3,
-                        given_degree,
-                        v,
-                        weights,
-                        values); // edge 11
-                  }
-                }
-              else
-                {
-                  Assert(false, ExcNotImplemented());
-                }
-            }
-
-          values += fe_eval.get_shape_info().dofs_per_component_on_cell;
-        }
-    }
-  };
-
-
 } // end of namespace internal
 
 
diff --git a/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h b/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h
new file mode 100644 (file)
index 0000000..6b9deb6
--- /dev/null
@@ -0,0 +1,622 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_evaluation_kernels_hanging_nodes_h
+#define dealii_matrix_free_evaluation_kernels_hanging_nodes_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/matrix_free/hanging_nodes_internal.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+// forward declaration
+template <int, typename, bool, typename>
+class FEEvaluationBaseData;
+
+
+
+namespace internal
+{
+  template <int dim, typename Number, bool is_face>
+  struct FEEvaluationImplHangingNodes
+  {
+    template <int fe_degree, int n_q_points_1d>
+    static bool
+    run(const unsigned int                  n_desired_components,
+        const FEEvaluationBaseData<dim,
+                                   typename Number::value_type,
+                                   is_face,
+                                   Number> &fe_eval,
+        const bool                          transpose,
+        const std::array<MatrixFreeFunctions::ConstraintKinds, Number::size()>
+          &     c_mask,
+        Number *values)
+    {
+      if (transpose)
+        run_internal<fe_degree, true>(n_desired_components,
+                                      fe_eval,
+                                      c_mask,
+                                      values);
+      else
+        run_internal<fe_degree, false>(n_desired_components,
+                                       fe_eval,
+                                       c_mask,
+                                       values);
+
+      return false;
+    }
+
+  private:
+    template <int fe_degree, unsigned int side, bool transpose>
+    static void
+    interpolate_2D(const unsigned int given_degree,
+                   const unsigned int v,
+                   const Number *     weight,
+                   Number *           values)
+    {
+      typename Number::value_type temp[40];
+
+      const unsigned int points =
+        (fe_degree != -1 ? fe_degree : given_degree) + 1;
+
+      AssertIndexRange(points, 40);
+
+      const unsigned int d = side / 2; // direction
+      const unsigned int s = side % 2; // left or right surface
+
+      const unsigned int offset = dealii::Utilities::pow(points, d + 1);
+      const unsigned int stride =
+        (s == 0 ? 0 : (points - 1)) * dealii::Utilities::pow(points, d);
+
+      const unsigned int r1 = dealii::Utilities::pow(points, dim - d - 1);
+      const unsigned int r2 = dealii::Utilities::pow(points, d);
+
+      // copy result back
+      for (unsigned int i = 0, k = 0; i < r1; ++i)
+        for (unsigned int j = 0; j < r2; ++j, ++k)
+          temp[k] = values[i * offset + stride + j][v];
+
+      // perform interpolation point by point (note: r1 * r2 == points^(dim-1))
+      for (unsigned int i = 0, k = 0; i < r1; ++i)
+        for (unsigned int j = 0; j < r2; ++j, ++k)
+          {
+            typename Number::value_type sum = 0.0;
+            for (unsigned int h = 0; h < points; ++h)
+              sum += weight[(transpose ? 1 : points) * k +
+                            (transpose ? points : 1) * h][v] *
+                     temp[h];
+            values[i * offset + stride + j][v] = sum;
+          }
+    }
+
+    template <int          fe_degree,
+              unsigned int direction,
+              unsigned int d,
+              bool         transpose>
+    static void
+    interpolate_3D_face(const unsigned int dof_offset,
+                        const unsigned int given_degree,
+                        const unsigned int v,
+                        const Number *     weight,
+                        Number *           values)
+    {
+      typename Number::value_type temp[40];
+
+      const unsigned int points =
+        (fe_degree != -1 ? fe_degree : given_degree) + 1;
+
+      AssertIndexRange(points, 40);
+
+      const unsigned int stride = Utilities::pow(points, direction);
+
+      // direction   side0   side1   side2
+      // 0             -      p^2      p
+      // 1            p^2      -       1
+      // 2             p       -       1
+      const unsigned int stride2 =
+        ((direction == 0 && d == 1) || (direction == 1 && d == 0)) ?
+          (points * points) :
+          (((direction == 0 && d == 2) || (direction == 2 && d == 0)) ? points :
+                                                                        1);
+
+      for (unsigned int g = 1; g < points - 1; ++g)
+        {
+          // copy result back
+          for (unsigned int k = 0; k < points; ++k)
+            temp[k] = values[dof_offset + k * stride + stride2 * g][v];
+
+          // perform interpolation point by point
+          for (unsigned int k = 0; k < points; ++k)
+            {
+              typename Number::value_type sum = 0.0;
+              for (unsigned int h = 0; h < points; ++h)
+                sum += weight[(transpose ? 1 : points) * k +
+                              (transpose ? points : 1) * h][v] *
+                       temp[h];
+              values[dof_offset + k * stride + stride2 * g][v] = sum;
+            }
+        }
+    }
+
+    template <int fe_degree, unsigned int direction, bool transpose>
+    static void
+    interpolate_3D_edge(const unsigned int p,
+                        const unsigned int given_degree,
+                        const unsigned int v,
+                        const Number *     weight,
+                        Number *           values)
+    {
+      typename Number::value_type temp[40];
+
+      const unsigned int points =
+        (fe_degree != -1 ? fe_degree : given_degree) + 1;
+
+      AssertIndexRange(points, 40);
+
+      const unsigned int stride = Utilities::pow(points, direction);
+
+      // copy result back
+      for (unsigned int k = 0; k < points; ++k)
+        temp[k] = values[p + k * stride][v];
+
+      // perform interpolation point by point
+      for (unsigned int k = 0; k < points; ++k)
+        {
+          typename Number::value_type sum = 0.0;
+          for (unsigned int h = 0; h < points; ++h)
+            sum += weight[(transpose ? 1 : points) * k +
+                          (transpose ? points : 1) * h][v] *
+                   temp[h];
+          values[p + k * stride][v] = sum;
+        }
+    }
+
+    template <int fe_degree, bool transpose>
+    static void
+    run_internal(const unsigned int                  n_desired_components,
+                 const FEEvaluationBaseData<dim,
+                                            typename Number::value_type,
+                                            is_face,
+                                            Number> &fe_eval,
+                 const std::array<MatrixFreeFunctions::ConstraintKinds,
+                                  Number::size()> &  constraint_mask,
+                 Number *                            values)
+    {
+      const unsigned int given_degree =
+        fe_degree != -1 ? fe_degree :
+                          fe_eval.get_shape_info().data.front().fe_degree;
+
+      const auto &interpolation_matrices =
+        fe_eval.get_shape_info().data.front().subface_interpolation_matrices;
+
+      const auto is_set = [](const auto a, const auto b) -> bool {
+        return (a & b) == b;
+      };
+
+      const auto not_set = [](const auto a, const auto b) -> bool {
+        return (a & b) == MatrixFreeFunctions::ConstraintKinds::unconstrained;
+      };
+
+      const unsigned int points = given_degree + 1;
+
+      for (unsigned int c = 0; c < n_desired_components; ++c)
+        {
+          for (unsigned int v = 0; v < Number::size(); ++v)
+            {
+              const auto mask = constraint_mask[v];
+
+              if (mask == MatrixFreeFunctions::ConstraintKinds::unconstrained)
+                continue;
+
+              if (dim == 2) // 2D: only faces
+                {
+                  // direction 0:
+                  if ((mask & MatrixFreeFunctions::ConstraintKinds::face_y) !=
+                      MatrixFreeFunctions::ConstraintKinds::unconstrained)
+                    {
+                      const bool is_subface_0 =
+                        (mask & MatrixFreeFunctions::ConstraintKinds::type_x) ==
+                        MatrixFreeFunctions::ConstraintKinds::unconstrained;
+
+                      const Number *weights =
+                        interpolation_matrices[is_subface_0].data();
+
+                      if (is_set(mask,
+                                 MatrixFreeFunctions::ConstraintKinds::type_y))
+                        interpolate_2D<fe_degree, 2, transpose>(
+                          given_degree,
+                          v,
+                          weights,
+                          values); // face 2
+                      else
+                        interpolate_2D<fe_degree, 3, transpose>(
+                          given_degree,
+                          v,
+                          weights,
+                          values); // face 3
+                    }
+
+                  // direction 1:
+                  if ((mask & MatrixFreeFunctions::ConstraintKinds::face_x) !=
+                      MatrixFreeFunctions::ConstraintKinds::unconstrained)
+                    {
+                      const bool is_subface_0 =
+                        (mask & MatrixFreeFunctions::ConstraintKinds::type_y) ==
+                        MatrixFreeFunctions::ConstraintKinds::unconstrained;
+
+                      const Number *weights =
+                        interpolation_matrices[is_subface_0].data();
+
+                      if (is_set(mask,
+                                 MatrixFreeFunctions::ConstraintKinds::type_x))
+                        interpolate_2D<fe_degree, 0, transpose>(
+                          given_degree,
+                          v,
+                          weights,
+                          values); // face 0
+                      else
+                        interpolate_2D<fe_degree, 1, transpose>(
+                          given_degree,
+                          v,
+                          weights,
+                          values); // face 1
+                    }
+                }
+              else if (dim == 3) // 3D faces and edges
+                {
+                  const unsigned int p0 = 0;
+                  const unsigned int p1 = points - 1;
+                  const unsigned int p2 = points * points - points;
+                  const unsigned int p3 = points * points - 1;
+                  const unsigned int p4 =
+                    points * points * points - points * points;
+                  const unsigned int p5 =
+                    points * points * points - points * points + points - 1;
+                  const unsigned int p6 = points * points * points - points;
+
+                  const bool is_face_0 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::face_x) &&
+                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_x);
+                  const bool is_face_1 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::face_x) &&
+                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_x);
+                  const bool is_face_2 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::face_y) &&
+                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
+                  const bool is_face_3 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::face_y) &&
+                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
+                  const bool is_face_4 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::face_z) &&
+                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+                  const bool is_face_5 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::face_z) &&
+                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+
+                  const bool is_edge_2 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_yz) &&
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::type_y) &&
+                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+                  const bool is_edge_3 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_yz) &&
+                    not_set(mask,
+                            MatrixFreeFunctions::ConstraintKinds::type_y) &&
+                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+                  const bool is_edge_6 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_yz) &&
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::type_y) &&
+                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+                  const bool is_edge_7 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_yz) &&
+                    not_set(mask,
+                            MatrixFreeFunctions::ConstraintKinds::type_y) &&
+                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+
+                  const bool is_edge_0 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_zx) &&
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::type_x) &&
+                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+                  const bool is_edge_1 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_zx) &&
+                    not_set(mask,
+                            MatrixFreeFunctions::ConstraintKinds::type_x) &&
+                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+                  const bool is_edge_4 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_zx) &&
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::type_x) &&
+                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+                  const bool is_edge_5 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_zx) &&
+                    not_set(mask,
+                            MatrixFreeFunctions::ConstraintKinds::type_x) &&
+                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_z);
+
+                  const bool is_edge_8 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_xy) &&
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::type_x) &&
+                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
+                  const bool is_edge_9 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_xy) &&
+                    not_set(mask,
+                            MatrixFreeFunctions::ConstraintKinds::type_x) &&
+                    is_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
+                  const bool is_edge_10 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_xy) &&
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::type_x) &&
+                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
+                  const bool is_edge_11 =
+                    is_set(mask,
+                           MatrixFreeFunctions::ConstraintKinds::edge_xy) &&
+                    not_set(mask,
+                            MatrixFreeFunctions::ConstraintKinds::type_x) &&
+                    not_set(mask, MatrixFreeFunctions::ConstraintKinds::type_y);
+
+                  // direction 0:
+                  {
+                    const bool is_subface_0 =
+                      (mask & MatrixFreeFunctions::ConstraintKinds::type_x) ==
+                      MatrixFreeFunctions::ConstraintKinds::unconstrained;
+
+                    const Number *weights =
+                      interpolation_matrices[is_subface_0].data();
+
+                    // ... faces
+                    if (is_face_2)
+                      interpolate_3D_face<fe_degree, 0, 1, transpose>(
+                        p0,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 2
+                    else if (is_face_3)
+                      interpolate_3D_face<fe_degree, 0, 1, transpose>(
+                        p2,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 3
+                    if (is_face_4)
+                      interpolate_3D_face<fe_degree, 0, 2, transpose>(
+                        p0,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 4
+                    else if (is_face_5)
+                      interpolate_3D_face<fe_degree, 0, 2, transpose>(
+                        p4,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 5
+
+                    // ... edges
+                    if (is_face_2 || is_face_4 || is_edge_2)
+                      interpolate_3D_edge<fe_degree, 0, transpose>(
+                        p0,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 2
+                    if (is_face_3 || is_face_4 || is_edge_3)
+                      interpolate_3D_edge<fe_degree, 0, transpose>(
+                        p2,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 3
+                    if (is_face_2 || is_face_5 || is_edge_6)
+                      interpolate_3D_edge<fe_degree, 0, transpose>(
+                        p4,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 6
+                    if (is_face_3 || is_face_5 || is_edge_7)
+                      interpolate_3D_edge<fe_degree, 0, transpose>(
+                        p6,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 7
+                  }
+
+                  // direction 1:
+                  {
+                    const bool is_subface_0 =
+                      (mask & MatrixFreeFunctions::ConstraintKinds::type_y) ==
+                      MatrixFreeFunctions::ConstraintKinds::unconstrained;
+
+                    const Number *weights =
+                      interpolation_matrices[is_subface_0].data();
+
+                    // ... faces
+                    if (is_face_0)
+                      interpolate_3D_face<fe_degree, 1, 0, transpose>(
+                        p0,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 0
+                    else if (is_face_1)
+                      interpolate_3D_face<fe_degree, 1, 0, transpose>(
+                        p1,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 1
+                    if (is_face_4)
+                      interpolate_3D_face<fe_degree, 1, 2, transpose>(
+                        p0,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 4
+                    else if (is_face_5)
+                      interpolate_3D_face<fe_degree, 1, 2, transpose>(
+                        p4,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 5
+
+                    // ... edges
+                    if (is_face_0 || is_face_4 || is_edge_0)
+                      interpolate_3D_edge<fe_degree, 1, transpose>(
+                        p0,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 0
+                    if (is_face_1 || is_face_4 || is_edge_1)
+                      interpolate_3D_edge<fe_degree, 1, transpose>(
+                        p1,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 1
+                    if (is_face_0 || is_face_5 || is_edge_4)
+                      interpolate_3D_edge<fe_degree, 1, transpose>(
+                        p4,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 4
+                    if (is_face_1 || is_face_5 || is_edge_5)
+                      interpolate_3D_edge<fe_degree, 1, transpose>(
+                        p5,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 5
+                  }
+
+                  // direction 2:
+                  {
+                    const bool is_subface_0 =
+                      (mask & MatrixFreeFunctions::ConstraintKinds::type_z) ==
+                      MatrixFreeFunctions::ConstraintKinds::unconstrained;
+
+                    const Number *weights =
+                      interpolation_matrices[is_subface_0].data();
+
+                    // ... faces
+                    if (is_face_0)
+                      interpolate_3D_face<fe_degree, 2, 0, transpose>(
+                        p0,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 0
+                    else if (is_face_1)
+                      interpolate_3D_face<fe_degree, 2, 0, transpose>(
+                        p1,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 1
+                    if (is_face_2)
+                      interpolate_3D_face<fe_degree, 2, 1, transpose>(
+                        p0,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 2
+                    else if (is_face_3)
+                      interpolate_3D_face<fe_degree, 2, 1, transpose>(
+                        p2,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // face 3
+
+                    // ... edges
+                    if (is_face_0 || is_face_2 || is_edge_8)
+                      interpolate_3D_edge<fe_degree, 2, transpose>(
+                        p0,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 8
+                    if (is_face_1 || is_face_2 || is_edge_9)
+                      interpolate_3D_edge<fe_degree, 2, transpose>(
+                        p1,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 9
+                    if (is_face_0 || is_face_3 || is_edge_10)
+                      interpolate_3D_edge<fe_degree, 2, transpose>(
+                        p2,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 10
+                    if (is_face_1 || is_face_3 || is_edge_11)
+                      interpolate_3D_edge<fe_degree, 2, transpose>(
+                        p3,
+                        given_degree,
+                        v,
+                        weights,
+                        values); // edge 11
+                  }
+                }
+              else
+                {
+                  Assert(false, ExcNotImplemented());
+                }
+            }
+
+          values += fe_eval.get_shape_info().dofs_per_component_on_cell;
+        }
+    }
+  };
+
+
+} // end of namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 62fc91bcd0d80b10afad5ced814689f0c04404b5..1c9d8b61efc9ae8e077ffefa7a9ebe46ea32cb2f 100644 (file)
@@ -19,7 +19,9 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/evaluation_kernels.h>
+#include <deal.II/matrix_free/shape_info.h>
 
 DEAL_II_NAMESPACE_OPEN
 
index f242cb5d23fc2089c11edaec8cc2cc93d456af54..b246ee4e445bfcc6a6ba72ac443b3d9c9372d190 100644 (file)
 #include <deal.II/matrix_free/evaluation_kernels.h>
 #include <deal.II/matrix_free/evaluation_selector.h>
 #include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/evaluation_template_factory_internal.h>
 #include <deal.II/matrix_free/fe_evaluation.h>
 
-#ifndef FE_EVAL_FACTORY_DEGREE_MAX
-#  define FE_EVAL_FACTORY_DEGREE_MAX 6
-#endif
-
 DEAL_II_NAMESPACE_OPEN
 
 namespace internal
 {
-  template <int degree, typename EvaluatorType, typename... Args>
-  bool
-  instantiation_helper_run(const unsigned int given_degree,
-                           const unsigned int n_q_points_1d,
-                           Args &...args)
-  {
-    if (given_degree == degree)
-      {
-        if (n_q_points_1d == degree + 1)
-          return EvaluatorType::template run<degree, degree + 1>(args...);
-        else if (n_q_points_1d == degree + 2)
-          return EvaluatorType::template run<degree, degree + 2>(args...);
-        else if (n_q_points_1d == degree)
-          return EvaluatorType::template run<degree, degree>(args...);
-        else if (n_q_points_1d == (3 * degree) / 2 + 1)
-          return EvaluatorType::template run<degree, (3 * degree) / 2 + 1>(
-            args...);
-        else
-          // slow path
-          return EvaluatorType::template run<-1, 0>(args...);
-      }
-    else if (degree < FE_EVAL_FACTORY_DEGREE_MAX)
-      return instantiation_helper_run<
-        (degree < FE_EVAL_FACTORY_DEGREE_MAX ? degree + 1 : degree),
-        EvaluatorType>(given_degree, n_q_points_1d, args...);
-    else
-      // slow path
-      return EvaluatorType::template run<-1, 0>(args...);
-  }
-
   struct FastEvaluationSupported
   {
     template <int fe_degree, int n_q_points_1d>
@@ -420,57 +387,6 @@ namespace internal
       fe_degree, n_q_points_1d, n_components, fe_eval, in_array, out_array);
   }
 
-
-
-  template <int dim, typename Number, typename VectorizedArrayType>
-  void
-  FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::apply(
-    const unsigned int n_components,
-    const unsigned int fe_degree,
-    const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType>
-      &                                            fe_eval,
-    const bool                                     transpose,
-    const std::array<MatrixFreeFunctions::ConstraintKinds,
-                     VectorizedArrayType::size()> &c_mask,
-    VectorizedArrayType *                          values)
-  {
-    instantiation_helper_run<
-      1,
-      FEEvaluationImplHangingNodes<dim, VectorizedArrayType, false>>(
-      fe_degree,
-      fe_degree + 1,
-      n_components,
-      fe_eval,
-      transpose,
-      c_mask,
-      values);
-  }
-
-
-
-  template <int dim, typename Number, typename VectorizedArrayType>
-  void
-  FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::apply(
-    const unsigned int n_components,
-    const unsigned int fe_degree,
-    const FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &fe_eval,
-    const bool                                     transpose,
-    const std::array<MatrixFreeFunctions::ConstraintKinds,
-                     VectorizedArrayType::size()> &c_mask,
-    VectorizedArrayType *                          values)
-  {
-    instantiation_helper_run<
-      1,
-      FEEvaluationImplHangingNodes<dim, VectorizedArrayType, true>>(
-      fe_degree,
-      fe_degree + 1,
-      n_components,
-      fe_eval,
-      transpose,
-      c_mask,
-      values);
-  }
-
 } // end of namespace internal
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h b/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h
new file mode 100644 (file)
index 0000000..cf03c7c
--- /dev/null
@@ -0,0 +1,72 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_evaluation_template_factory_hanging_nodes_templates_h
+#define dealii_matrix_free_evaluation_template_factory_hanging_nodes_templates_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/matrix_free/evaluation_kernels_hanging_nodes.h>
+#include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/evaluation_template_factory_internal.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+  template <int dim, typename Number, typename VectorizedArrayType>
+  void
+  FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::apply(
+    const unsigned int n_components,
+    const unsigned int fe_degree,
+    const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType>
+      &                                            fe_eval,
+    const bool                                     transpose,
+    const std::array<MatrixFreeFunctions::ConstraintKinds,
+                     VectorizedArrayType::size()> &c_mask,
+    VectorizedArrayType *                          values)
+  {
+    instantiation_helper_degree_run<
+      1,
+      FEEvaluationImplHangingNodes<dim, VectorizedArrayType, false>>(
+      fe_degree, n_components, fe_eval, transpose, c_mask, values);
+  }
+
+
+
+  template <int dim, typename Number, typename VectorizedArrayType>
+  void
+  FEEvaluationHangingNodesFactory<dim, Number, VectorizedArrayType>::apply(
+    const unsigned int n_components,
+    const unsigned int fe_degree,
+    const FEEvaluationBaseData<dim, Number, true, VectorizedArrayType> &fe_eval,
+    const bool                                     transpose,
+    const std::array<MatrixFreeFunctions::ConstraintKinds,
+                     VectorizedArrayType::size()> &c_mask,
+    VectorizedArrayType *                          values)
+  {
+    instantiation_helper_degree_run<
+      1,
+      FEEvaluationImplHangingNodes<dim, VectorizedArrayType, true>>(
+      fe_degree, n_components, fe_eval, transpose, c_mask, values);
+  }
+} // end of namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/matrix_free/evaluation_template_factory_internal.h b/include/deal.II/matrix_free/evaluation_template_factory_internal.h
new file mode 100644 (file)
index 0000000..94f66ca
--- /dev/null
@@ -0,0 +1,83 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_evaluation_template_factory_internal_h
+#define dealii_matrix_free_evaluation_template_factory_internal_h
+
+
+#include <deal.II/base/config.h>
+
+#ifndef FE_EVAL_FACTORY_DEGREE_MAX
+#  define FE_EVAL_FACTORY_DEGREE_MAX 6
+#endif
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+  template <int degree, typename EvaluatorType, typename... Args>
+  bool
+  instantiation_helper_run(const unsigned int given_degree,
+                           const unsigned int n_q_points_1d,
+                           Args &...args)
+  {
+    if (given_degree == degree)
+      {
+        if (n_q_points_1d == degree + 1)
+          return EvaluatorType::template run<degree, degree + 1>(args...);
+        else if (n_q_points_1d == degree + 2)
+          return EvaluatorType::template run<degree, degree + 2>(args...);
+        else if (n_q_points_1d == degree)
+          return EvaluatorType::template run<degree, degree>(args...);
+        else if (n_q_points_1d == (3 * degree) / 2 + 1)
+          return EvaluatorType::template run<degree, (3 * degree) / 2 + 1>(
+            args...);
+        else
+          // slow path
+          return EvaluatorType::template run<-1, 0>(args...);
+      }
+    else if (degree < FE_EVAL_FACTORY_DEGREE_MAX)
+      return instantiation_helper_run<
+        (degree < FE_EVAL_FACTORY_DEGREE_MAX ? degree + 1 : degree),
+        EvaluatorType>(given_degree, n_q_points_1d, args...);
+    else
+      // slow path
+      return EvaluatorType::template run<-1, 0>(args...);
+  }
+
+  template <int degree, typename EvaluatorType, typename... Args>
+  bool
+  instantiation_helper_degree_run(const unsigned int given_degree,
+                                  Args &...args)
+  {
+    if (given_degree == degree)
+      {
+        return EvaluatorType::template run<degree, degree + 1>(args...);
+      }
+    else if (degree < FE_EVAL_FACTORY_DEGREE_MAX)
+      return instantiation_helper_degree_run<
+        (degree < FE_EVAL_FACTORY_DEGREE_MAX ? degree + 1 : degree),
+        EvaluatorType>(given_degree, args...);
+    else
+      // slow path
+      return EvaluatorType::template run<-1, 0>(args...);
+  }
+
+} // end of namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index bb85bc75d761ccf4a27d87a61fac34fe1f9a908c..716b9b85117e18a77682c6f9aac7fbac6d7d2e0e 100644 (file)
@@ -23,6 +23,12 @@ SET(_src
   evaluation_template_factory_inst4.cc
   evaluation_template_factory_inst5.cc
   evaluation_template_factory_inst6.cc
+  evaluation_template_factory_hanging_nodes.cc
+  evaluation_template_factory_hanging_nodes_inst2.cc
+  evaluation_template_factory_hanging_nodes_inst3.cc
+  evaluation_template_factory_hanging_nodes_inst4.cc
+  evaluation_template_factory_hanging_nodes_inst5.cc
+  evaluation_template_factory_hanging_nodes_inst6.cc
   fe_point_evaluation.cc
   mapping_info.cc
   mapping_info_inst2.cc
@@ -35,6 +41,7 @@ SET(_src
 
 SET(_inst
   evaluation_template_factory.inst.in
+  evaluation_template_factory_hanging_nodes.inst.in
   fe_point_evaluation.inst.in
   mapping_info.inst.in
   matrix_free.inst.in
index f3e05fbf096b3b4aa4d8e936c65b891adb286bbd..468f6f3724e30c90ebe5a4f7c991b248970fd710 100644 (file)
@@ -32,10 +32,4 @@ for (deal_II_dimension : DIMENSIONS;
       deal_II_dimension,
       deal_II_scalar_vectorized::value_type,
       deal_II_scalar_vectorized>;
-
-    // inverse mass
-    template struct dealii::internal::FEEvaluationHangingNodesFactory<
-      deal_II_dimension,
-      deal_II_scalar_vectorized::value_type,
-      deal_II_scalar_vectorized>;
   }
diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes.cc
new file mode 100644 (file)
index 0000000..0d7a736
--- /dev/null
@@ -0,0 +1,28 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+#define SPLIT_INSTANTIATIONS_COUNT 6
+#ifndef SPLIT_INSTANTIATIONS_INDEX
+#  define SPLIT_INSTANTIATIONS_INDEX 0
+#endif
+
+#include "evaluation_template_factory_hanging_nodes.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes.inst.in b/source/matrix_free/evaluation_template_factory_hanging_nodes.inst.in
new file mode 100644 (file)
index 0000000..0ed6662
--- /dev/null
@@ -0,0 +1,25 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS;
+     deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)
+  {
+    // hanging nodes
+    template struct dealii::internal::FEEvaluationHangingNodesFactory<
+      deal_II_dimension,
+      deal_II_scalar_vectorized::value_type,
+      deal_II_scalar_vectorized>;
+  }
diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst2.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst2.cc
new file mode 100644 (file)
index 0000000..0a69454
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 1
+#include "evaluation_template_factory_hanging_nodes.cc"
diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst3.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst3.cc
new file mode 100644 (file)
index 0000000..4e17425
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 2
+#include "evaluation_template_factory_hanging_nodes.cc"
diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst4.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst4.cc
new file mode 100644 (file)
index 0000000..cab92f1
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 3
+#include "evaluation_template_factory_hanging_nodes.cc"
diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst5.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst5.cc
new file mode 100644 (file)
index 0000000..ea2904f
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 4
+#include "evaluation_template_factory_hanging_nodes.cc"
diff --git a/source/matrix_free/evaluation_template_factory_hanging_nodes_inst6.cc b/source/matrix_free/evaluation_template_factory_hanging_nodes_inst6.cc
new file mode 100644 (file)
index 0000000..e6d89c6
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 5
+#include "evaluation_template_factory_hanging_nodes.cc"
index 4ee4fe50970bfecdd034929c41df77db653a8be2..d93b708336acf370d25c64508a0f532d4fb09ea5 100644 (file)
@@ -24,6 +24,7 @@
 
 #include <deal.II/grid/grid_generator.h>
 
+#include <deal.II/matrix_free/evaluation_kernels_hanging_nodes.h>
 #include <deal.II/matrix_free/fe_evaluation.h>
 #include <deal.II/matrix_free/matrix_free.h>
 

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.