From 898e8785ce45efbd9d87641e7a682870ef38f88e Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 28 Feb 2025 09:27:07 +0100 Subject: [PATCH] FEFaceEvaluation: enable FE_DGP --- .../deal.II/matrix_free/evaluation_kernels.h | 70 ++++-------- .../matrix_free/evaluation_kernels_common.h | 101 ++++++++++++++++++ .../matrix_free/evaluation_kernels_face.h | 80 +++++++++----- .../deal.II/matrix_free/fe_evaluation_data.h | 2 +- tests/matrix_free/matrix_vector_faces_36.cc | 42 ++++++++ .../matrix_free/matrix_vector_faces_36.output | 17 +++ 6 files changed, 236 insertions(+), 76 deletions(-) create mode 100644 include/deal.II/matrix_free/evaluation_kernels_common.h create mode 100644 tests/matrix_free/matrix_vector_faces_36.cc create mode 100644 tests/matrix_free/matrix_vector_faces_36.output diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 236c157b1c..4817976afe 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -242,36 +243,19 @@ namespace internal (type == MatrixFreeFunctions::truncated_tensor) ? Utilities::pow(shape_data.front().fe_degree + 1, dim) : fe_eval.get_shape_info().dofs_per_component_on_cell; - const Number *values_dofs = values_dofs_actual; + const Number *values_dofs = + (type == MatrixFreeFunctions::truncated_tensor) ? + temp1 + 2 * (std::max( + fe_eval.get_shape_info().dofs_per_component_on_cell, + n_q_points)) : + values_dofs_actual; + if (type == MatrixFreeFunctions::truncated_tensor) - { - const std::size_t n_dofs_per_comp = - fe_eval.get_shape_info().dofs_per_component_on_cell; - Number *values_dofs_tmp = - temp1 + 2 * (std::max(n_dofs_per_comp, n_q_points)); - const int degree = - fe_degree != -1 ? fe_degree : shape_data.front().fe_degree; - for (unsigned int c = 0; c < n_components; ++c) - for (int i = 0, count_p = 0, count_q = 0; - i < (dim > 2 ? degree + 1 : 1); - ++i) - { - for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j) - { - for (int k = 0; k < degree + 1 - j - i; - ++k, ++count_p, ++count_q) - values_dofs_tmp[c * dofs_per_comp + count_q] = - values_dofs_actual[c * n_dofs_per_comp + count_p]; - for (int k = degree + 1 - j - i; k < degree + 1; - ++k, ++count_q) - values_dofs_tmp[c * dofs_per_comp + count_q] = Number(); - } - for (int j = degree + 1 - i; j < degree + 1; ++j) - for (int k = 0; k < degree + 1; ++k, ++count_q) - values_dofs_tmp[c * dofs_per_comp + count_q] = Number(); - } - values_dofs = values_dofs_tmp; - } + embed_truncated_into_full_tensor_product( + n_components, + const_cast(values_dofs), + values_dofs_actual, + fe_eval); Number *values_quad = fe_eval.begin_values(); Number *gradients_quad = fe_eval.begin_gradients(); @@ -436,6 +420,7 @@ namespace internal (type == MatrixFreeFunctions::truncated_tensor) ? Utilities::fixed_power(shape_data.front().fe_degree + 1) : fe_eval.get_shape_info().dofs_per_component_on_cell; + // expand dof_values to tensor product for truncated tensor products Number *values_dofs = (type == MatrixFreeFunctions::truncated_tensor) ? @@ -582,28 +567,11 @@ namespace internal } if (type == MatrixFreeFunctions::truncated_tensor) - { - const std::size_t n_dofs_per_comp = - fe_eval.get_shape_info().dofs_per_component_on_cell; - values_dofs -= dofs_per_comp * n_components; - const int degree = - fe_degree != -1 ? fe_degree : shape_data.front().fe_degree; - for (unsigned int c = 0; c < n_components; ++c) - for (int i = 0, count_p = 0, count_q = 0; - i < (dim > 2 ? degree + 1 : 1); - ++i) - { - for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j) - { - for (int k = 0; k < degree + 1 - j - i; - ++k, ++count_p, ++count_q) - values_dofs_actual[c * n_dofs_per_comp + count_p] = - values_dofs[c * dofs_per_comp + count_q]; - count_q += j + i; - } - count_q += i * (degree + 1); - } - } + truncate_tensor_product_to_complete_degrees( + n_components, + values_dofs_actual, + values_dofs - dofs_per_comp * n_components, + fe_eval); } diff --git a/include/deal.II/matrix_free/evaluation_kernels_common.h b/include/deal.II/matrix_free/evaluation_kernels_common.h new file mode 100644 index 0000000000..07a5d15537 --- /dev/null +++ b/include/deal.II/matrix_free/evaluation_kernels_common.h @@ -0,0 +1,101 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2025 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +#ifndef dealii_matrix_free_evaluation_kernels_common_h +#define dealii_matrix_free_evaluation_kernels_common_h + +#include + +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace internal +{ + template + inline void + embed_truncated_into_full_tensor_product( + const unsigned int n_components, + Number *values_dofs, + const Number *values_dofs_actual, + FEEvaluationData &fe_eval) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + + const std::size_t dofs_per_comp = + Utilities::pow(shape_data.fe_degree + 1, dim); + const std::size_t n_dofs_per_comp = shape_info.dofs_per_component_on_cell; + const int degree = fe_degree != -1 ? fe_degree : shape_data.fe_degree; + + for (unsigned int c = 0; c < n_components; ++c) + for (int i = 0, count_p = 0, count_q = 0; i < (dim > 2 ? degree + 1 : 1); + ++i) + { + for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j) + { + for (int k = 0; k < degree + 1 - j - i; ++k, ++count_p, ++count_q) + values_dofs[c * dofs_per_comp + count_q] = + values_dofs_actual[c * n_dofs_per_comp + count_p]; + for (int k = degree + 1 - j - i; k < degree + 1; ++k, ++count_q) + values_dofs[c * dofs_per_comp + count_q] = Number(); + } + for (int j = degree + 1 - i; j < degree + 1; ++j) + for (int k = 0; k < degree + 1; ++k, ++count_q) + values_dofs[c * dofs_per_comp + count_q] = Number(); + } + } + + + + template + inline void + truncate_tensor_product_to_complete_degrees( + const unsigned int n_components, + Number *values_dofs_actual, + const Number *values_dofs, + FEEvaluationData &fe_eval) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + + const unsigned int dofs_per_comp = + Utilities::fixed_power(shape_data.fe_degree + 1); + const std::size_t n_dofs_per_comp = shape_info.dofs_per_component_on_cell; + const int degree = fe_degree != -1 ? fe_degree : shape_data.fe_degree; + for (unsigned int c = 0; c < n_components; ++c) + for (int i = 0, count_p = 0, count_q = 0; i < (dim > 2 ? degree + 1 : 1); + ++i) + { + for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j) + { + for (int k = 0; k < degree + 1 - j - i; ++k, ++count_p, ++count_q) + values_dofs_actual[c * n_dofs_per_comp + count_p] = + values_dofs[c * dofs_per_comp + count_q]; + count_q += j + i; + } + count_q += i * (degree + 1); + } + } +} // end of namespace internal + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/matrix_free/evaluation_kernels_face.h b/include/deal.II/matrix_free/evaluation_kernels_face.h index e3f2e118bd..6214c6585c 100644 --- a/include/deal.II/matrix_free/evaluation_kernels_face.h +++ b/include/deal.II/matrix_free/evaluation_kernels_face.h @@ -25,6 +25,7 @@ #include #include +#include #include #include #include @@ -874,16 +875,20 @@ namespace internal interpolate_raviart_thomas( n_components, input, output, flags, face_no, shape_info); else - interpolate_generic( - n_components, - input, - output, - flags, - face_no, - shape_info.data.front().fe_degree + 1, - shape_info.data.front().shape_data_on_face, - shape_info.dofs_per_component_on_cell, - 3 * shape_info.dofs_per_component_on_face); + { + const unsigned int fe_degree_ = shape_info.data.front().fe_degree; + + interpolate_generic( + n_components, + input, + output, + flags, + face_no, + fe_degree_ + 1, + shape_info.data.front().shape_data_on_face, + Utilities::pow(fe_degree_ + 1, dim), + 3 * Utilities::pow(fe_degree_ + 1, dim - 1)); + } } /** @@ -1751,7 +1756,7 @@ namespace internal static bool evaluate_tensor(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number *values_dofs, + const Number *values_dofs_actual, FEEvaluationData &fe_eval) { const auto &shape_info = fe_eval.get_shape_info(); @@ -1764,8 +1769,22 @@ namespace internal // Note: we always keep storage of values, 1st and 2nd derivatives in an // array, so reserve space for all three here - Number *temp = fe_eval.get_scratch_data().begin(); - Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face; + Number *temp1 = fe_eval.get_scratch_data().begin(); + Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face; + + const Number *values_dofs = + (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ? + temp2 + + 2 * (std::max(fe_eval.get_shape_info().dofs_per_component_on_cell, + shape_info.n_q_points)) : + values_dofs_actual; + + if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) + embed_truncated_into_full_tensor_product( + n_components, + const_cast(values_dofs), + values_dofs_actual, + fe_eval); bool use_vectorization = true; if (fe_eval.get_dof_access_index() == @@ -1781,15 +1800,15 @@ namespace internal values_dofs, fe_eval, use_vectorization, - temp, - scratch_data); + temp1, + temp2); evaluate_in_face( - n_components, evaluation_flag, fe_eval, temp, scratch_data); + n_components, evaluation_flag, fe_eval, temp1, temp2); if (dim == 3) adjust_quadrature_for_face_orientation( - n_components, evaluation_flag, fe_eval, use_vectorization, temp); + n_components, evaluation_flag, fe_eval, use_vectorization, temp1); return false; } @@ -2290,7 +2309,7 @@ namespace internal static bool integrate_tensor(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, - Number *values_dofs, + Number *values_dofs_actual, FEEvaluationData &fe_eval, const bool sum_into_values) { @@ -2302,8 +2321,16 @@ namespace internal Utilities::pow(fe_degree + 1, dim - 1) : Utilities::fixed_power(shape_data.fe_degree + 1); - Number *temp = fe_eval.get_scratch_data().begin(); - Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face; + Number *temp1 = fe_eval.get_scratch_data().begin(); + Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face; + + // expand dof_values to tensor product for truncated tensor products + Number *values_dofs = + (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ? + temp2 + 2 * (std::max( + fe_eval.get_shape_info().dofs_per_component_on_cell, + fe_eval.get_shape_info().n_q_points)) : + values_dofs_actual; bool use_vectorization = true; @@ -2321,20 +2348,25 @@ namespace internal if (dim == 3) adjust_quadrature_for_face_orientation( - n_components, integration_flag, fe_eval, use_vectorization, temp); + n_components, integration_flag, fe_eval, use_vectorization, temp1); integrate_in_face( - n_components, integration_flag, fe_eval, temp, scratch_data); + n_components, integration_flag, fe_eval, temp1, temp2); collect_from_face(n_components, integration_flag, values_dofs, fe_eval, use_vectorization, - temp, - scratch_data, + temp1, + temp2, sum_into_values); + + if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) + truncate_tensor_product_to_complete_degrees( + n_components, values_dofs_actual, values_dofs, fe_eval); + return false; } diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h index 7a54ba8c56..517e6cd9a7 100644 --- a/include/deal.II/matrix_free/fe_evaluation_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_data.h @@ -1201,7 +1201,7 @@ FEEvaluationData::set_data_pointers( const unsigned int size_scratch_data = std::max(tensor_dofs_per_component + 1, dofs_per_component) * n_components * - 3 + + 4 + 2 * n_quadrature_points; const unsigned int size_data_arrays = n_components * dofs_per_component + diff --git a/tests/matrix_free/matrix_vector_faces_36.cc b/tests/matrix_free/matrix_vector_faces_36.cc new file mode 100644 index 0000000000..d562b65532 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_36.cc @@ -0,0 +1,42 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2025 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// Similar to matrix_vector_faces_36 but with FE_DGP. + +#include + +#include + +#include "../tests.h" + +#include "matrix_vector_faces_common.h" + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(5 - dim); + + FE_DGP fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + do_test(dof, constraints, true); +} diff --git a/tests/matrix_free/matrix_vector_faces_36.output b/tests/matrix_free/matrix_vector_faces_36.output new file mode 100644 index 0000000000..c4fe177ff1 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_36.output @@ -0,0 +1,17 @@ + +DEAL:2d::Testing FE_DGP<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d::Norm of difference parallel: 0 +DEAL:2d:: +DEAL:2d::Testing FE_DGP<2>(2) +DEAL:2d::Norm of difference: 0 +DEAL:2d::Norm of difference parallel: 0 +DEAL:2d:: +DEAL:3d::Testing FE_DGP<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d::Norm of difference parallel: 0 +DEAL:3d:: +DEAL:3d::Testing FE_DGP<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d::Norm of difference parallel: 0 +DEAL:3d:: -- 2.39.5