From 114254c6642db892a9006dc236e1b79901480071 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 23 May 2022 12:04:08 +0200 Subject: [PATCH] MatrixFree: Simplify initialization of ShapeInfo --- .../matrix_free/shape_info.templates.h | 127 +++++++----------- 1 file changed, 48 insertions(+), 79 deletions(-) diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 9b44982a6e..2dc95d4e4c 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -81,6 +81,7 @@ namespace internal } + template void get_element_type_specific_information( @@ -131,15 +132,13 @@ namespace internal else Assert(false, ExcNotImplemented()); - // Finally store the renumbering into the member variable of this - // class + // Finally store the renumbering into the respective field if (fe_in.n_components() == 1) lexicographic_numbering = scalar_lexicographic; else { - // have more than one component, get the inverse - // permutation, invert it, sort the components one after one, - // and invert back + // have more than one component, get the inverse permutation, invert + // it, sort the components one by one, and invert back std::vector scalar_inv = Utilities::invert_permutation(scalar_lexicographic); std::vector lexicographic( @@ -171,6 +170,8 @@ namespace internal } } + + template std::unique_ptr> create_fe(const FiniteElement &fe) @@ -739,15 +740,13 @@ namespace internal quadrature_data_on_face[0].resize(quad.size() * 3); quadrature_data_on_face[1].resize(quad.size() * 3); - dealii::FE_DGQArbitraryNodes<1> fe_quad(quad); + const std::vector> poly_coll = + Polynomials::generate_complete_Lagrange_basis(quad.get_points()); for (unsigned int i = 0; i < quad.size(); ++i) { - Point<1> q_point; - q_point[0] = 0; - quadrature_data_on_face[0][i] = fe_quad.shape_value(i, q_point); - q_point[0] = 1; - quadrature_data_on_face[1][i] = fe_quad.shape_value(i, q_point); + quadrature_data_on_face[0][i] = poly_coll[i].value(0.0); + quadrature_data_on_face[1][i] = poly_coll[i].value(1.0); } } @@ -762,81 +761,51 @@ namespace internal auto &subface_interpolation_matrix_scalar_1 = univariate_shape_data.subface_interpolation_matrices_scalar[1]; - const auto fe_1d = create_fe<1>(fe); - const auto fe_2d = create_fe<2>(fe); - - 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_0, - 0); - - fe_2d->get_subface_interpolation_matrix(*fe_2d, - 1, - interpolation_matrix_1, - 0); - - ElementType element_type; - std::vector scalar_lexicographic; - std::vector lexicographic_numbering; - - get_element_type_specific_information(*fe_1d, - *fe_1d, - 0, - element_type, - scalar_lexicographic, - lexicographic_numbering); + const unsigned int nn = fe_degree + 1; + subface_interpolation_matrix_0.resize(nn * nn); + subface_interpolation_matrix_1.resize(nn * nn); + subface_interpolation_matrix_scalar_0.resize(nn * nn); + subface_interpolation_matrix_scalar_1.resize(nn * nn); - 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()); + std::vector> fe_q_points = QGaussLobatto<1>(nn).get_points(); + const std::vector> poly = + Polynomials::generate_complete_Lagrange_basis(fe_q_points); - subface_interpolation_matrix_scalar_0.resize( - fe_1d->n_dofs_per_cell() * fe_1d->n_dofs_per_cell()); - subface_interpolation_matrix_scalar_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) + for (unsigned int i = 0, c = 0; i < nn; ++i) + for (unsigned int j = 0; j < nn; ++j, ++c) { - 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]); - subface_interpolation_matrix_scalar_0[c] = - interpolation_matrix_0(scalar_lexicographic[i], - scalar_lexicographic[j]); + poly[j].value(0.5 * fe_q_points[i][0]); + subface_interpolation_matrix_0[c] = + subface_interpolation_matrix_scalar_0[c]; subface_interpolation_matrix_scalar_1[c] = - interpolation_matrix_1(scalar_lexicographic[i], - scalar_lexicographic[j]); + poly[j].value(0.5 + 0.5 * fe_q_points[i][0]); + subface_interpolation_matrix_1[c] = + subface_interpolation_matrix_scalar_1[c]; } } // get gradient and Hessian transformation matrix for the polynomial // space associated with the quadrature rule (collocation space). We // need to avoid the case with more than a few hundreds of quadrature - // points when the Lagrange polynomials constructed in - // FE_DGQArbitraryNodes underflow. + // points when the Lagrange polynomials might underflow. Note that 200 + // is not an exact value, as different quadrature formulas behave + // slightly differently, but 200 has been observed to be low enough for + // all common quadrature formula types. For QGauss, the actual limit is + // 517 points, for example. if (n_q_points_1d < 200) { shape_gradients_collocation.resize(n_q_points_1d * n_q_points_1d); shape_hessians_collocation.resize(n_q_points_1d * n_q_points_1d); - FE_DGQArbitraryNodes<1> fe_coll(quad.get_points()); + const std::vector> poly_coll = + Polynomials::generate_complete_Lagrange_basis(quad.get_points()); + std::array values; for (unsigned int i = 0; i < n_q_points_1d; ++i) for (unsigned int q = 0; q < n_q_points_1d; ++q) { - shape_gradients_collocation[i * n_q_points_1d + q] = - fe_coll.shape_grad(i, quad.get_points()[q])[0]; - shape_hessians_collocation[i * n_q_points_1d + q] = - fe_coll.shape_grad_grad(i, quad.get_points()[q])[0][0]; + poly_coll[i].value(quad.get_points()[q][0], 2, values.data()); + shape_gradients_collocation[i * n_q_points_1d + q] = values[1]; + shape_hessians_collocation[i * n_q_points_1d + q] = values[2]; } // compute the inverse shape functions in three steps: we first @@ -859,7 +828,7 @@ namespace internal for (unsigned int i = 0; i < n_q_points_1d; ++i) for (unsigned int j = 0; j < n_q_points_1d; ++j) transform_to_gauss(i, j) = - fe_coll.shape_value(j, quad_gauss.point(i)); + poly_coll[j].value(quad_gauss.point(i)[0]); // step 2: computation for the projection (in reference coordinates) // from higher to lower polynomial degree @@ -872,15 +841,17 @@ namespace internal // polynomials where most of the interpolation matrices are unit // matrices when applying the inverse mass matrix, so we do not need // to compute much. - QGauss<1> quad_project(n_dofs_1d); - FE_DGQArbitraryNodes<1> fe_project(quad_project.get_points()); + QGauss<1> quad_project(n_dofs_1d); + const std::vector> poly_project = + Polynomials::generate_complete_Lagrange_basis( + quad_project.get_points()); FullMatrix project_gauss(n_dofs_1d, n_q_points_1d); for (unsigned int i = 0; i < n_dofs_1d; ++i) for (unsigned int q = 0; q < n_q_points_1d; ++q) project_gauss(i, q) = - fe_project.shape_value(i, quad_gauss.get_points()[q]) * + poly_project[i].value(quad_gauss.get_points()[q][0]) * (quad_gauss.weight(q) / quad_project.weight(i)); FullMatrix project_to_dof_space(n_dofs_1d, n_q_points_1d); project_gauss.mmult(project_to_dof_space, transform_to_gauss); @@ -899,11 +870,8 @@ namespace internal { for (unsigned int i = 0; i < n_dofs_1d; ++i) for (unsigned int j = 0; j < n_dofs_1d; ++j) - transform_from_gauss(i, j) = fe_project.shape_value( - j, - Point<1>( - fe.get_unit_support_points()[scalar_lexicographic[i]] - [0])); + transform_from_gauss(i, j) = poly_project[j].value( + fe.get_unit_support_points()[scalar_lexicographic[i]][0]); FullMatrix result(n_dofs_1d, n_q_points_1d); transform_from_gauss.mmult(result, project_to_dof_space); @@ -1178,8 +1146,9 @@ namespace internal shape_hessians_eo = convert_to_eo(shape_hessians, fe_degree + 1, n_q_points_1d); - // FE_DGQArbitraryNodes underflow (see also above where - // shape_gradients_collocation and shape_hessians_collocation is set up). + // Avoid underflow of Lagrange polynomials on typical quadrature + // formulas (see also above where shape_gradients_collocation and + // shape_hessians_collocation is set up). if (n_q_points_1d < 200) { shape_gradients_collocation_eo = -- 2.39.5