From 66ccfce9f1678343cb0900cc181554252977550c Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 23 Sep 2021 09:03:43 +0200 Subject: [PATCH] ShapeInfo: store both subface interpolation matrices --- .../matrix_free/cuda_matrix_free.templates.h | 4 +- .../deal.II/matrix_free/evaluation_kernels.h | 167 +++++------------- include/deal.II/matrix_free/shape_info.h | 6 +- .../matrix_free/shape_info.templates.h | 36 ++-- tests/matrix_free/hanging_node_kernels_01.cc | 10 +- 5 files changed, 79 insertions(+), 144 deletions(-) diff --git a/include/deal.II/matrix_free/cuda_matrix_free.templates.h b/include/deal.II/matrix_free/cuda_matrix_free.templates.h index 5039981471..4c64f02188 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h @@ -275,9 +275,9 @@ namespace CUDAWrappers { cudaError_t error_code = cudaMemcpyToSymbol( constraint_weights, - shape_info.data.front().subface_interpolation_matrix.data(), + shape_info.data.front().subface_interpolation_matrices[0].data(), sizeof(double) * - shape_info.data.front().subface_interpolation_matrix.size()); + shape_info.data.front().subface_interpolation_matrices[0].size()); AssertCuda(error_code); local_dof_indices.resize(data->dofs_per_cell); diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 8ca4b63eb9..2fa379d2b9 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -4940,29 +4940,6 @@ namespace internal private: template static void - interpolate_2D(const unsigned int given_degree, - bool is_subface_0, - const unsigned int v, - const Number * weight, - Number * values) - { - if (is_subface_0) - interpolate_2D(given_degree, - v, - weight, - values); - else - interpolate_2D(given_degree, - v, - weight, - values); - } - - template - static void interpolate_2D(const unsigned int given_degree, const unsigned int v, const Number * weight, @@ -4996,10 +4973,8 @@ namespace internal { typename Number::value_type sum = 0.0; for (unsigned int h = 0; h < points; ++h) - sum += weight[((transpose ? 1 : points) * - (is_subface_0 ? k : (points - 1 - k))) + - ((transpose ? points : 1) * - (is_subface_0 ? h : (points - 1 - h)))][v] * + sum += weight[(transpose ? 1 : points) * k + + (transpose ? points : 1) * h][v] * temp[h]; values[i * offset + stride + j][v] = sum; } @@ -5007,30 +4982,9 @@ namespace internal template static void - interpolate_3D_face(const unsigned int dof_offset, - const unsigned int given_degree, - bool is_subface_0, - const unsigned int v, - const Number * weight, - Number * values) - { - if (is_subface_0) - interpolate_3D_face( - dof_offset, given_degree, v, weight, values); - else - interpolate_3D_face( - dof_offset, given_degree, v, weight, values); - } - - template - static void interpolate_3D_face(const unsigned int dof_offset, const unsigned int given_degree, const unsigned int v, @@ -5045,7 +4999,6 @@ namespace internal AssertIndexRange(points, 40); const unsigned int stride = Utilities::pow(points, direction); - const unsigned int d = side / 2; // direction side0 side1 side2 // 0 - p^2 p @@ -5068,10 +5021,8 @@ namespace internal { typename Number::value_type sum = 0.0; for (unsigned int h = 0; h < points; ++h) - sum += weight[((transpose ? 1 : points) * - (is_subface_0 ? k : (points - 1 - k))) + - ((transpose ? points : 1) * - (is_subface_0 ? h : (points - 1 - h)))][v] * + sum += weight[(transpose ? 1 : points) * k + + (transpose ? points : 1) * h][v] * temp[h]; values[dof_offset + k * stride + stride2 * g][v] = sum; } @@ -5080,26 +5031,6 @@ namespace internal template static void - interpolate_3D_edge(const unsigned int p, - const unsigned int given_degree, - bool is_subface_0, - const unsigned int v, - const Number * weight, - Number * values) - { - if (is_subface_0) - interpolate_3D_edge( - p, given_degree, v, weight, values); - else - interpolate_3D_edge( - p, given_degree, v, weight, values); - } - - template - static void interpolate_3D_edge(const unsigned int p, const unsigned int given_degree, const unsigned int v, @@ -5124,10 +5055,8 @@ namespace internal { typename Number::value_type sum = 0.0; for (unsigned int h = 0; h < points; ++h) - sum += weight[((transpose ? 1 : points) * - (is_subface_0 ? k : (points - 1 - k))) + - ((transpose ? points : 1) * - (is_subface_0 ? h : (points - 1 - h)))][v] * + sum += weight[(transpose ? 1 : points) * k + + (transpose ? points : 1) * h][v] * temp[h]; values[p + k * stride][v] = sum; } @@ -5144,14 +5073,13 @@ namespace internal Number::size()> & constraint_mask, Number * values) { - const Number *weights = fe_eval.get_shape_info() - .data.front() - .subface_interpolation_matrix.data(); - 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; }; @@ -5178,20 +5106,22 @@ namespace internal MatrixFreeFunctions::ConstraintKinds::unconstrained) { const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_x) != + (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( given_degree, - is_subface_0, v, weights, values); // face 2 else interpolate_2D( given_degree, - is_subface_0, v, weights, values); // face 3 @@ -5202,20 +5132,22 @@ namespace internal MatrixFreeFunctions::ConstraintKinds::unconstrained) { const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_y) != + (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( given_degree, - is_subface_0, v, weights, values); // face 0 else interpolate_2D( given_degree, - is_subface_0, v, weights, values); // face 1 @@ -5336,39 +5268,38 @@ namespace internal // direction 0: { const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_x) != + (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( + interpolate_3D_face( p0, given_degree, - is_subface_0, v, weights, values); // face 2 else if (is_face_3) - interpolate_3D_face( + interpolate_3D_face( p2, given_degree, - is_subface_0, v, weights, values); // face 3 if (is_face_4) - interpolate_3D_face( + interpolate_3D_face( p0, given_degree, - is_subface_0, v, weights, values); // face 4 else if (is_face_5) - interpolate_3D_face( + interpolate_3D_face( p4, given_degree, - is_subface_0, v, weights, values); // face 5 @@ -5378,7 +5309,6 @@ namespace internal interpolate_3D_edge( p0, given_degree, - is_subface_0, v, weights, values); // edge 2 @@ -5386,7 +5316,6 @@ namespace internal interpolate_3D_edge( p2, given_degree, - is_subface_0, v, weights, values); // edge 3 @@ -5394,7 +5323,6 @@ namespace internal interpolate_3D_edge( p4, given_degree, - is_subface_0, v, weights, values); // edge 6 @@ -5402,7 +5330,6 @@ namespace internal interpolate_3D_edge( p6, given_degree, - is_subface_0, v, weights, values); // edge 7 @@ -5411,39 +5338,38 @@ namespace internal // direction 1: { const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_y) != + (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( p0, given_degree, - is_subface_0, v, weights, values); // face 0 else if (is_face_1) - interpolate_3D_face( + interpolate_3D_face( p1, given_degree, - is_subface_0, v, weights, values); // face 1 if (is_face_4) - interpolate_3D_face( + interpolate_3D_face( p0, given_degree, - is_subface_0, v, weights, values); // face 4 else if (is_face_5) - interpolate_3D_face( + interpolate_3D_face( p4, given_degree, - is_subface_0, v, weights, values); // face 5 @@ -5453,7 +5379,6 @@ namespace internal interpolate_3D_edge( p0, given_degree, - is_subface_0, v, weights, values); // edge 0 @@ -5461,7 +5386,6 @@ namespace internal interpolate_3D_edge( p1, given_degree, - is_subface_0, v, weights, values); // edge 1 @@ -5469,7 +5393,6 @@ namespace internal interpolate_3D_edge( p4, given_degree, - is_subface_0, v, weights, values); // edge 4 @@ -5477,7 +5400,6 @@ namespace internal interpolate_3D_edge( p5, given_degree, - is_subface_0, v, weights, values); // edge 5 @@ -5486,39 +5408,38 @@ namespace internal // direction 2: { const bool is_subface_0 = - (mask & MatrixFreeFunctions::ConstraintKinds::type_z) != + (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( p0, given_degree, - is_subface_0, v, weights, values); // face 0 else if (is_face_1) - interpolate_3D_face( + interpolate_3D_face( p1, given_degree, - is_subface_0, v, weights, values); // face 1 if (is_face_2) - interpolate_3D_face( + interpolate_3D_face( p0, given_degree, - is_subface_0, v, weights, values); // face 2 else if (is_face_3) - interpolate_3D_face( + interpolate_3D_face( p2, given_degree, - is_subface_0, v, weights, values); // face 3 @@ -5528,7 +5449,6 @@ namespace internal interpolate_3D_edge( p0, given_degree, - is_subface_0, v, weights, values); // edge 8 @@ -5536,7 +5456,6 @@ namespace internal interpolate_3D_edge( p1, given_degree, - is_subface_0, v, weights, values); // edge 9 @@ -5544,7 +5463,6 @@ namespace internal interpolate_3D_edge( p2, given_degree, - is_subface_0, v, weights, values); // edge 10 @@ -5552,7 +5470,6 @@ namespace internal interpolate_3D_edge( p3, given_degree, - is_subface_0, v, weights, values); // edge 11 diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index ac641152f8..bab18e73a4 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -262,10 +262,10 @@ namespace internal std::array, 2> hessians_within_subface; /** - * A 1D subface interpolation matrix to the first quadrant. This data - * structure is only set up for FE_Q for dim > 1. + * A 1D subface interpolation matrices to the first and second quadrant. + * This data structure is only set up for FE_Q for dim > 1. */ - AlignedVector subface_interpolation_matrix; + std::array, 2> subface_interpolation_matrices; /** * We store a copy of the one-dimensional quadrature formula diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 3c716a4fe4..e1643d51da 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -627,18 +627,27 @@ namespace internal if (dim > 1 && dynamic_cast *>(&fe)) { - auto &subface_interpolation_matrix = - univariate_shape_data.subface_interpolation_matrix; + auto &subface_interpolation_matrix_0 = + univariate_shape_data.subface_interpolation_matrices[0]; + auto &subface_interpolation_matrix_1 = + univariate_shape_data.subface_interpolation_matrices[1]; const auto fe_1d = create_fe<1>(fe); const auto fe_2d = create_fe<2>(fe); - FullMatrix interpolation_matrix(fe_2d->n_dofs_per_face(0), - fe_2d->n_dofs_per_face(0)); + FullMatrix interpolation_matrix_0(fe_2d->n_dofs_per_face(0), + fe_2d->n_dofs_per_face(0)); + FullMatrix interpolation_matrix_1(fe_2d->n_dofs_per_face(0), + fe_2d->n_dofs_per_face(0)); fe_2d->get_subface_interpolation_matrix(*fe_2d, 0, - interpolation_matrix, + interpolation_matrix_0, + 0); + + fe_2d->get_subface_interpolation_matrix(*fe_2d, + 1, + interpolation_matrix_1, 0); ElementType element_type; @@ -652,14 +661,21 @@ namespace internal scalar_lexicographic, lexicographic_numbering); - subface_interpolation_matrix.resize(fe_1d->n_dofs_per_cell() * - fe_1d->n_dofs_per_cell()); + subface_interpolation_matrix_0.resize(fe_1d->n_dofs_per_cell() * + fe_1d->n_dofs_per_cell()); + subface_interpolation_matrix_1.resize(fe_1d->n_dofs_per_cell() * + fe_1d->n_dofs_per_cell()); for (unsigned int i = 0, c = 0; i < fe_1d->n_dofs_per_cell(); ++i) for (unsigned int j = 0; j < fe_1d->n_dofs_per_cell(); ++j, ++c) - subface_interpolation_matrix[c] = - interpolation_matrix(scalar_lexicographic[i], - scalar_lexicographic[j]); + { + subface_interpolation_matrix_0[c] = + interpolation_matrix_0(scalar_lexicographic[i], + scalar_lexicographic[j]); + subface_interpolation_matrix_1[c] = + interpolation_matrix_1(scalar_lexicographic[i], + scalar_lexicographic[j]); + } } // get gradient and Hessian transformation matrix for the polynomial diff --git a/tests/matrix_free/hanging_node_kernels_01.cc b/tests/matrix_free/hanging_node_kernels_01.cc index ff0999b81b..4ee4fe5097 100644 --- a/tests/matrix_free/hanging_node_kernels_01.cc +++ b/tests/matrix_free/hanging_node_kernels_01.cc @@ -113,8 +113,9 @@ namespace dealii Number::size()> & constraint_mask, Number * values) { - const auto &constraint_weights = - fe_eval.get_shape_info().data.front().subface_interpolation_matrix; + const auto &constraint_weights = fe_eval.get_shape_info() + .data.front() + .subface_interpolation_matrices[0]; const unsigned int fe_degree = fe_degree_ != -1 ? fe_degree_ : @@ -252,8 +253,9 @@ namespace dealii Number::size()> & constraint_mask, Number * values) { - const auto &constraint_weights = - fe_eval.get_shape_info().data.front().subface_interpolation_matrix; + const auto &constraint_weights = fe_eval.get_shape_info() + .data.front() + .subface_interpolation_matrices[0]; const unsigned int fe_degree = fe_degree_ != -1 ? fe_degree_ : -- 2.39.5