From e1dd4e9a25420bbf3c545611c883a422973e5ead Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 20 Dec 2021 10:20:13 +0100 Subject: [PATCH] Bug fix of MF face eval for hanging nodes and non-standard orientation --- .../deal.II/matrix_free/evaluation_kernels.h | 25 ++++--- .../deal.II/matrix_free/face_setup_internal.h | 15 +++++ include/deal.II/matrix_free/shape_info.h | 7 ++ .../matrix_free/shape_info.templates.h | 67 +++++++++---------- 4 files changed, 66 insertions(+), 48 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index d4101092fa..d26e9ce2b4 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3783,19 +3783,18 @@ namespace internal } } else if (fe_eval.get_face_orientation() != 0) - for (unsigned int c = 0; c < n_components; ++c) - adjust_for_face_orientation( - dim, - n_components, - evaluation_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(), 0), - false, - Utilities::pow(n_q_points_1d, dim - 1), - temp, - fe_eval.begin_values(), - fe_eval.begin_gradients(), - fe_eval.begin_hessians()); + adjust_for_face_orientation( + dim, + n_components, + evaluation_flag, + &fe_eval.get_shape_info().face_orientations_quad( + fe_eval.get_face_orientation(), 0), + false, + Utilities::pow(n_q_points_1d, dim - 1), + temp, + fe_eval.begin_values(), + fe_eval.begin_gradients(), + fe_eval.begin_hessians()); } return false; diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h index bcfd049961..bdcc6516e6 100644 --- a/include/deal.II/matrix_free/face_setup_internal.h +++ b/include/deal.II/matrix_free/face_setup_internal.h @@ -28,6 +28,7 @@ #include #include +#include #include #include @@ -1041,6 +1042,20 @@ namespace internal } else info.face_orientation = exterior_face_orientation; + + // make sure to select correct subface index in case of non-standard + // orientation of the coarser neighbor face + if (cell->level() > neighbor->level() && exterior_face_orientation > 0) + { + const Table<2, unsigned int> orientation = + ShapeInfo::compute_orientation_table(2); + const std::array inverted_orientations{ + {0, 1, 2, 3, 6, 5, 4, 7}}; + info.subface_index = + orientation[inverted_orientations[exterior_face_orientation]] + [info.subface_index]; + } + return info; } diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index d7269738da..1187d0708c 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -375,6 +375,13 @@ namespace internal static bool is_supported(const FiniteElement &fe); + /** + * Compute a table with numbers of re-orientation for all versions of + * face flips, orientation, and rotation (relating only to 3D elements). + */ + static Table<2, unsigned int> + compute_orientation_table(const unsigned int n_points_per_dim); + /** * Return data of univariate shape functions which defines the * dimension @p dimension of tensor product shape functions diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index fd13f5fb0f..9f1a084a81 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -854,41 +854,8 @@ namespace internal // (similar to MappingInfoStorage::QuadratureDescriptor::initialize) if (dim == 3) { - const auto compute_orientations = - [](const unsigned int n, - Table<2, unsigned int> &face_orientations) { - face_orientations.reinit(8, n * n); - for (unsigned int j = 0, i = 0; j < n; ++j) - for (unsigned int k = 0; k < n; ++k, ++i) - { - // face_orientation=true, face_flip=false, - // face_rotation=false - face_orientations[0][i] = i; - // face_orientation=false, face_flip=false, - // face_rotation=false - face_orientations[1][i] = j + k * n; - // face_orientation=true, face_flip=true, - // face_rotation=false - face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n; - // face_orientation=false, face_flip=true, - // face_rotation=false - face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n; - // face_orientation=true, face_flip=false, - // face_rotation=true - face_orientations[4][i] = j + (n - 1 - k) * n; - // face_orientation=false, face_flip=false, - // face_rotation=true - face_orientations[5][i] = k + (n - 1 - j) * n; - // face_orientation=true, face_flip=true, - // face_rotation=true - face_orientations[6][i] = (n - 1 - j) + k * n; - // face_orientation=false, face_flip=true, - // face_rotation=true - face_orientations[7][i] = (n - 1 - k) + j * n; - } - }; - compute_orientations(fe_degree + 1, face_orientations_dofs); - compute_orientations(n_q_points_1d, face_orientations_quad); + face_orientations_dofs = compute_orientation_table(fe_degree + 1); + face_orientations_quad = compute_orientation_table(n_q_points_1d); } else { @@ -1163,6 +1130,36 @@ namespace internal + template + Table<2, unsigned int> + ShapeInfo::compute_orientation_table(const unsigned int n) + { + Table<2, unsigned int> face_orientations(8, n * n); + for (unsigned int j = 0, i = 0; j < n; ++j) + for (unsigned int k = 0; k < n; ++k, ++i) + { + // face_orientation=true, face_flip=false, face_rotation=false + face_orientations[0][i] = i; + // face_orientation=false, face_flip=false, face_rotation=false + face_orientations[1][i] = j + k * n; + // face_orientation=true, face_flip=true, face_rotation=false + face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n; + // face_orientation=false, face_flip=true, face_rotation=false + face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n; + // face_orientation=true, face_flip=false, face_rotation=true + face_orientations[4][i] = j + (n - 1 - k) * n; + // face_orientation=false, face_flip=false, face_rotation=true + face_orientations[5][i] = k + (n - 1 - j) * n; + // face_orientation=true, face_flip=true, face_rotation=true + face_orientations[6][i] = (n - 1 - j) + k * n; + // face_orientation=false, face_flip=true, face_rotation=true + face_orientations[7][i] = (n - 1 - k) + j * n; + } + return face_orientations; + } + + + template std::size_t ShapeInfo::memory_consumption() const -- 2.39.5