Add weight_fe_q_dofs_by_entity_shifted() and compute_weights_fe_q_dofs_by_entity_shif...

diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h
index fddbbf3ee78b4902f410564d3e8171f7929b7bb8..9aad0f510085fc738be302224ff91da4f55bf765 100644 (file)
--- a/include/deal.II/matrix_free/tensor_product_kernels.h
+++ b/include/deal.II/matrix_free/tensor_product_kernels.h
@@ -3361,6 +3361,51 @@ namespace internal
   }
 
 
+  template <int dim, int n_points_1d_template, typename Number>
+  inline void
+  weight_fe_q_dofs_by_entity_shifted(const Number *     weights,
+                                     const unsigned int n_components,
+                                     const int n_points_1d_non_template,
+                                     Number *  data)
+  {
+    const int n_points_1d = n_points_1d_template != -1 ?
+                              n_points_1d_template :
+                              n_points_1d_non_template;
+
+    Assert((n_points_1d % 2) == 1,
+           ExcMessage("The function can only with add number of points"));
+    Assert(n_points_1d > 0, ExcNotImplemented());
+    Assert(n_points_1d < 100, ExcNotImplemented());
+
+    const unsigned int n_inside_1d = n_points_1d / 2;
+
+    unsigned int compressed_index[100];
+
+    unsigned int c = 0;
+    for (int i = 0; i < n_inside_1d; ++i)
+      compressed_index[c++] = 0;
+    compressed_index[c++] = 1;
+    for (int i = 0; i < n_inside_1d; ++i)
+      compressed_index[c++] = 2;
+
+    for (unsigned int c = 0; c < n_components; ++c)
+      for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
+        for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
+          {
+            const unsigned int shift =
+              9 * compressed_index[k] + 3 * compressed_index[j];
+
+            unsigned int c = 0;
+            for (int i = 0; i < n_inside_1d; ++i)
+              data[c++] *= weights[shift];
+            data[c++] *= weights[shift + 1];
+            for (int i = 0; i < n_inside_1d; ++i)
+              data[c++] *= weights[shift + 2];
+            data += n_points_1d;
+          }
+  }
+
+
   template <int dim, int n_points_1d_template, typename Number>
   inline bool
   compute_weights_fe_q_dofs_by_entity(const Number *     data,
@@ -3422,6 +3467,78 @@ namespace internal
   }
 
 
+  template <int dim, int n_points_1d_template, typename Number>
+  inline bool
+  compute_weights_fe_q_dofs_by_entity_shifted(
+    const Number *     data,
+    const unsigned int n_components,
+    const int          n_points_1d_non_template,
+    Number *           weights)
+  {
+    const int n_points_1d = n_points_1d_template != -1 ?
+                              n_points_1d_template :
+                              n_points_1d_non_template;
+
+    Assert((n_points_1d % 2) == 1,
+           ExcMessage("The function can only with add number of points"));
+    Assert(n_points_1d > 0, ExcNotImplemented());
+    Assert(n_points_1d < 100, ExcNotImplemented());
+
+    const unsigned int n_inside_1d = n_points_1d / 2;
+
+    unsigned int compressed_index[100];
+
+    unsigned int c = 0;
+    for (int i = 0; i < n_inside_1d; ++i)
+      compressed_index[c++] = 0;
+    compressed_index[c++] = 1;
+    for (int i = 0; i < n_inside_1d; ++i)
+      compressed_index[c++] = 2;
+
+    // Insert the number data into a storage position for weight,
+    // ensuring that the array has either not been touched before
+    // or the previous content is the same. In case the previous
+    // content has a different value, we exit this function and
+    // signal to outer functions that the compression was not possible.
+    const auto check_and_set = [](Number &weight, const Number &data) {
+      if (weight == Number(-1.0) || weight == data)
+        {
+          weight = data;
+          return true; // success for the entry
+        }
+
+      return false; // failure for the entry
+    };
+
+    for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
+      weights[c] = Number(-1.0);
+
+    for (unsigned int comp = 0; comp < n_components; ++comp)
+      for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
+        for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
+             ++j, data += n_points_1d)
+          {
+            const unsigned int shift =
+              9 * compressed_index[k] + 3 * compressed_index[j];
+
+            unsigned int c = 0;
+
+            for (int i = 0; i < n_inside_1d; ++i)
+              if (!check_and_set(weights[shift], data[c++]))
+                return false; // failure
+
+            if (!check_and_set(weights[shift + 1], data[c++]))
+              return false; // failure
+
+            for (int i = 0; i < n_inside_1d; ++i)
+              if (!check_and_set(weights[shift + 2], data[c++]))
+                return false; // failure
+          }
+
+    return true; // success
+  }
+
+
 } // end of namespace internal
 
 

diff --git a/tests/matrix_free/weight_fe_q_dofs_by_entity.cc b/tests/matrix_free/weight_fe_q_dofs_by_entity.cc
new file mode 100644 (file)
index 0000000..2e5a809
--- /dev/null
+++ b/tests/matrix_free/weight_fe_q_dofs_by_entity.cc
@@ -0,0 +1,104 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test weight_fe_q_dofs_by_entity() and
+// compute_weights_fe_q_dofs_by_entity_shifted()
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+
+#include <iomanip>
+#include <iostream>
+
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  {
+    const unsigned int dim       = 2;
+    const unsigned int fe_degree = 5;
+    using Number                 = double;
+
+    std::vector<Number> weights(Utilities::pow(3, dim));
+    for (unsigned int i = 0; i < weights.size(); ++i)
+      weights[i] = i;
+
+    std::vector<Number> values(Utilities::pow(fe_degree + 1, dim), 1.0);
+
+    internal::weight_fe_q_dofs_by_entity<dim, -1, Number>(weights.data(),
+                                                          1,
+                                                          fe_degree + 1,
+                                                          values.data());
+
+    for (unsigned int i_1 = 0, c = 0; i_1 < fe_degree + 1; ++i_1)
+      {
+        for (unsigned int i_0 = 0; i_0 < fe_degree + 1; ++i_0, ++c)
+          deallog << values[c] << " ";
+        deallog << std::endl;
+      }
+    deallog << std::endl;
+
+    for (auto &i : weights)
+      i = 0.0;
+
+    internal::compute_weights_fe_q_dofs_by_entity<dim, -1, Number>(
+      values.data(), 1, fe_degree + 1, weights.data());
+
+    for (const auto i : weights)
+      deallog << i << " ";
+
+    deallog << std::endl;
+    deallog << std::endl;
+  }
+
+  {
+    const unsigned int dim       = 2;
+    const unsigned int fe_degree = 4;
+    using Number                 = double;
+
+    std::vector<Number> weights(Utilities::pow(3, dim));
+    for (unsigned int i = 0; i < weights.size(); ++i)
+      weights[i] = i;
+
+    std::vector<Number> values(Utilities::pow((2 * fe_degree - 1), dim), 1.0);
+
+    internal::weight_fe_q_dofs_by_entity_shifted<dim, -1, Number>(
+      weights.data(), 1, 2 * fe_degree - 1, values.data());
+
+    for (unsigned int i_1 = 0, c = 0; i_1 < (2 * fe_degree - 1); ++i_1)
+      {
+        for (unsigned int i_0 = 0; i_0 < (2 * fe_degree - 1); ++i_0, ++c)
+          deallog << values[c] << " ";
+        deallog << std::endl;
+      }
+    deallog << std::endl;
+
+    for (auto &i : weights)
+      i = 0.0;
+
+    internal::compute_weights_fe_q_dofs_by_entity_shifted<dim, -1, Number>(
+      values.data(), 1, 2 * fe_degree - 1, weights.data());
+
+    for (const auto i : weights)
+      deallog << i << " ";
+
+    deallog << std::endl;
+    deallog << std::endl;
+  }
+}

diff --git a/tests/matrix_free/weight_fe_q_dofs_by_entity.output b/tests/matrix_free/weight_fe_q_dofs_by_entity.output
new file mode 100644 (file)
index 0000000..e2f83fe
--- /dev/null
+++ b/tests/matrix_free/weight_fe_q_dofs_by_entity.output
@@ -0,0 +1,20 @@
+
+DEAL::0.00000 1.00000 1.00000 1.00000 1.00000 2.00000 
+DEAL::3.00000 4.00000 4.00000 4.00000 4.00000 5.00000 
+DEAL::3.00000 4.00000 4.00000 4.00000 4.00000 5.00000 
+DEAL::3.00000 4.00000 4.00000 4.00000 4.00000 5.00000 
+DEAL::3.00000 4.00000 4.00000 4.00000 4.00000 5.00000 
+DEAL::6.00000 7.00000 7.00000 7.00000 7.00000 8.00000 
+DEAL::
+DEAL::0.00000 1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 
+DEAL::
+DEAL::0.00000 0.00000 0.00000 1.00000 2.00000 2.00000 2.00000 
+DEAL::0.00000 0.00000 0.00000 1.00000 2.00000 2.00000 2.00000 
+DEAL::0.00000 0.00000 0.00000 1.00000 2.00000 2.00000 2.00000 
+DEAL::3.00000 3.00000 3.00000 4.00000 5.00000 5.00000 5.00000 
+DEAL::6.00000 6.00000 6.00000 7.00000 8.00000 8.00000 8.00000 
+DEAL::6.00000 6.00000 6.00000 7.00000 8.00000 8.00000 8.00000 
+DEAL::6.00000 6.00000 6.00000 7.00000 8.00000 8.00000 8.00000 
+DEAL::
+DEAL::0.00000 1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 
+DEAL::