From 5fcd5adf99d056c428e2b0bf4808bcd51dc2f036 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 12 Dec 2022 13:48:21 +0100 Subject: [PATCH] Add weight_fe_q_dofs_by_entity_shifted() and compute_weights_fe_q_dofs_by_entity_shifted() --- .../matrix_free/tensor_product_kernels.h | 117 ++++++++++++++++++ .../matrix_free/weight_fe_q_dofs_by_entity.cc | 104 ++++++++++++++++ .../weight_fe_q_dofs_by_entity.output | 20 +++ 3 files changed, 241 insertions(+) create mode 100644 tests/matrix_free/weight_fe_q_dofs_by_entity.cc create mode 100644 tests/matrix_free/weight_fe_q_dofs_by_entity.output diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index fddbbf3ee7..9aad0f5100 100644 --- 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 + 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 inline bool compute_weights_fe_q_dofs_by_entity(const Number * data, @@ -3422,6 +3467,78 @@ namespace internal } + template + 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(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 index 0000000000..2e5a809762 --- /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 + +#include +#include + +#include "../tests.h" + +int +main() +{ + initlog(); + + { + const unsigned int dim = 2; + const unsigned int fe_degree = 5; + using Number = double; + + std::vector weights(Utilities::pow(3, dim)); + for (unsigned int i = 0; i < weights.size(); ++i) + weights[i] = i; + + std::vector values(Utilities::pow(fe_degree + 1, dim), 1.0); + + internal::weight_fe_q_dofs_by_entity(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( + 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 weights(Utilities::pow(3, dim)); + for (unsigned int i = 0; i < weights.size(); ++i) + weights[i] = i; + + std::vector values(Utilities::pow((2 * fe_degree - 1), dim), 1.0); + + internal::weight_fe_q_dofs_by_entity_shifted( + 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( + 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 index 0000000000..e2f83fe932 --- /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:: -- 2.39.5