From 36338ab0ac8df3c547c8fa3501b6416c2b042594 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 21 Aug 2023 23:23:54 +0200 Subject: [PATCH] MatrixFree ShapeInfo: Move some code to UnivariateShapeData --- include/deal.II/matrix_free/shape_info.h | 48 ++- .../matrix_free/shape_info.templates.h | 351 ++++++++---------- 2 files changed, 194 insertions(+), 205 deletions(-) diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index f280d555ef..258d8f23ce 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -144,6 +144,35 @@ namespace internal std::size_t memory_consumption() const; + /** + * Evaluate the auxiliary polynomial space associated with the Lagrange + * polynomials in points of the given quadrature formula, filling the + * fields shape_[gradients,hessians]_collocation and related + * information. + */ + template + void + evaluate_collocation_space( + const FiniteElement &fe, + const Quadrature<1> & quad, + const std::vector & lexicographic); + + /** + * Check whether we have symmetries in the shape values. In that case, + * also fill the shape_???_eo fields. + */ + bool + check_and_set_shapes_symmetric(); + + /** + * Check whether symmetric 1d basis functions are such that the shape + * values form a diagonal matrix, i.e., the nodal points are collocated + * with the quadrature points. This allows for specialized algorithms + * that save some operations in the evaluation. + */ + bool + check_shapes_collocation() const; + /** * Encodes the type of element detected at construction. FEEvaluation * will select the most efficient algorithm based on the given element @@ -578,25 +607,6 @@ namespace internal * quadrature points to represent the correct order. */ dealii::Table<2, unsigned int> face_orientations_quad; - - private: - /** - * Check whether we have symmetries in the shape values. In that case, - * also fill the shape_???_eo fields. - */ - bool - check_1d_shapes_symmetric( - UnivariateShapeData &univariate_shape_data); - - /** - * Check whether symmetric 1d basis functions are such that the shape - * values form a diagonal matrix, i.e., the nodal points are collocated - * with the quadrature points. This allows for specialized algorithms - * that save some operations in the evaluation. - */ - bool - check_1d_shapes_collocation( - const UnivariateShapeData &univariate_shape_data) const; }; diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 8f3c72cffa..c743f6831e 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -286,7 +286,7 @@ namespace internal lex_normal.push_back(i); lex_normal.push_back(dofs_per_face_normal); - // 'direction' distingusishes between normal and tangential direction + // 'direction' distinguishes between normal and tangential direction for (unsigned int direction = 0; direction < 2; ++direction) { UnivariateShapeData &univariate_shape_data = @@ -610,17 +610,10 @@ namespace internal return; // grant write access to common univariate shape data - auto &shape_values = univariate_shape_data.shape_values; - auto &shape_gradients = univariate_shape_data.shape_gradients; - auto &shape_hessians = univariate_shape_data.shape_hessians; - auto &shape_gradients_collocation = - univariate_shape_data.shape_gradients_collocation; - auto &shape_hessians_collocation = - univariate_shape_data.shape_hessians_collocation; - auto &inverse_shape_values = univariate_shape_data.inverse_shape_values; - auto &shape_data_on_face = univariate_shape_data.shape_data_on_face; - auto &quadrature_data_on_face = - univariate_shape_data.quadrature_data_on_face; + auto &shape_values = univariate_shape_data.shape_values; + auto &shape_gradients = univariate_shape_data.shape_gradients; + auto &shape_hessians = univariate_shape_data.shape_hessians; + auto &shape_data_on_face = univariate_shape_data.shape_data_on_face; auto &values_within_subface = univariate_shape_data.values_within_subface; auto &gradients_within_subface = univariate_shape_data.gradients_within_subface; @@ -737,26 +730,6 @@ namespace internal fe.shape_grad_grad(my_i, q_point)[0][0]; } - if (n_q_points_1d < 200) - { - quadrature_data_on_face[0].resize(quad.size() * 3); - quadrature_data_on_face[1].resize(quad.size() * 3); - - const std::vector> poly_coll = - Polynomials::generate_complete_Lagrange_basis(quad.get_points()); - - for (unsigned int i = 0; i < quad.size(); ++i) - { - std::array values; - poly_coll[i].value(0.0, 2, values.data()); - for (unsigned int d = 0; d < 3; ++d) - quadrature_data_on_face[0][i + d * quad.size()] = values[d]; - poly_coll[i].value(1.0, 2, values.data()); - for (unsigned int d = 0; d < 3; ++d) - quadrature_data_on_face[1][i + d * quad.size()] = values[d]; - } - } - if (dim > 1 && (dynamic_cast *>(&fe) || dynamic_cast *>(&fe))) { @@ -807,136 +780,17 @@ namespace internal } } - // 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 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); - 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) - { - 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 - // change from the given quadrature formula and the associated - // Lagrange polynomials to the Lagrange polynomials at quadrature - // points. in this basis, we can then perform the second step, which - // is the computation of a projection matrix from the potentially - // higher polynomial degree associated to the quadrature points to a - // polynomial space of degree equal to the degree of the given - // elements. in the third step, we change from the Lagrange - // polynomials in the Gauss quadrature points to the polynomial - // space of the given element - - // step 1: change basis from the Lagrange polynomials at the given - // quadrature points to the Lagrange basis at Gauss quadrature - // points. this is often the identity operation as we often compute - // with Gaussian quadrature, but not necessarily so - QGauss<1> quad_gauss(n_q_points_1d); - FullMatrix transform_to_gauss(n_q_points_1d, n_q_points_1d); - 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) = - poly_coll[j].value(quad_gauss.point(i)[0]); - - // step 2: computation for the projection (in reference coordinates) - // from higher to lower polynomial degree - // - // loop over quadrature points, multiply by q-weight on high degree - // integrate loop going from high degree to low degree loop over new - // points, multiply by inverse q-weight on low degree - // - // This projection step is for the special case of Lagrange - // 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); - 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) = - 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); - - // step 3: change the basis back to the given finite element - // space. we can use a shortcut for elements that define support - // points, in which case we can evaluate the Lagrange polynomials of - // the Gauss quadrature in those points. this will give more - // accurate results than the inversion of a matrix. for more general - // polynomial spaces, we must invert a matrix of a Vandermonde type, - // which we do by a Householder transformation to keep roundoff - // errors low. - inverse_shape_values.resize_fast(array_size); - FullMatrix transform_from_gauss(n_dofs_1d, n_dofs_1d); - if (fe.has_support_points()) - { - 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) = 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); - - // set very small entries to zero - we are in reference space - // with normalized numbers, so this is straight-forward to check - // here - for (unsigned int i = 0; i < n_dofs_1d; ++i) - for (unsigned int q = 0; q < n_q_points_1d; ++q) - inverse_shape_values[i * n_q_points_1d + q] = - std::abs(result(i, q)) < 1e-15 ? 0 : result(i, q); - } - else - { - for (unsigned int i = 0; i < n_dofs_1d; ++i) - for (unsigned int j = 0; j < n_dofs_1d; ++j) - { - Point q_point = unit_point; - q_point[0] = quad_project.point(i)[0]; - - transform_from_gauss(i, j) = - fe.shape_value(scalar_lexicographic[j], q_point); - } - Householder H(transform_from_gauss); - Vector in(n_dofs_1d), out(n_dofs_1d); - for (unsigned int q = 0; q < n_q_points_1d; ++q) - { - for (unsigned int i = 0; i < n_dofs_1d; ++i) - in(i) = project_to_dof_space(i, q); - H.least_squares(out, in); - for (unsigned int i = 0; i < n_dofs_1d; ++i) - inverse_shape_values[i * n_q_points_1d + q] = - std::abs(out(i)) < 1e-15 ? 0. : out(i); - } - } - } + univariate_shape_data.evaluate_collocation_space(fe, + quad, + scalar_lexicographic); if (element_type == tensor_general && - check_1d_shapes_symmetric(univariate_shape_data)) + univariate_shape_data.check_and_set_shapes_symmetric()) { if (dynamic_cast *>(&fe) && fe.tensor_degree() > 1) element_type = tensor_symmetric_no_collocation; - else if (check_1d_shapes_collocation(univariate_shape_data)) + else if (univariate_shape_data.check_shapes_collocation()) element_type = tensor_symmetric_collocation; else element_type = tensor_symmetric; @@ -956,7 +810,7 @@ namespace internal } } else if (element_type == tensor_symmetric_plus_dg0) - check_1d_shapes_symmetric(univariate_shape_data); + univariate_shape_data.check_and_set_shapes_symmetric(); nodal_at_cell_boundaries = true; for (unsigned int i = 1; i < n_dofs_1d; ++i) @@ -1057,34 +911,163 @@ namespace internal + template + template + void + UnivariateShapeData::evaluate_collocation_space( + const FiniteElement &fe, + const Quadrature<1> & quad, + const std::vector & lexicographic) + { + const unsigned int n_q_points_1d = quad.size(); + const unsigned int n_dofs_1d = + std::min(fe.n_dofs_per_cell(), fe_degree + 1); + + // 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 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) + return; + + shape_gradients_collocation.resize(n_q_points_1d * n_q_points_1d); + shape_hessians_collocation.resize(n_q_points_1d * n_q_points_1d); + 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) + { + 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 + // change from the given quadrature formula and the associated + // Lagrange polynomials to the Lagrange polynomials at quadrature + // points. in this basis, we can then perform the second step, which + // is the computation of a projection matrix from the potentially + // higher polynomial degree associated to the quadrature points to a + // polynomial space of degree equal to the degree of the given + // elements. in the third step, we change from the Lagrange + // polynomials in the Gauss quadrature points to the polynomial + // space of the given element + + // step 1: change basis from the Lagrange polynomials at the given + // quadrature points to the Lagrange basis at Gauss quadrature + // points. this is often the identity operation as we often compute + // with Gaussian quadrature, but not necessarily so + QGauss<1> quad_gauss(n_q_points_1d); + FullMatrix transform_to_gauss(n_q_points_1d, n_q_points_1d); + 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) = poly_coll[j].value(quad_gauss.point(i)[0]); + + // step 2: computation for the projection (in reference coordinates) + // from higher to lower polynomial degree + // + // loop over quadrature points, multiply by q-weight on high degree + // integrate loop going from high degree to low degree loop over new + // points, multiply by inverse q-weight on low degree + // + // This projection step is for the special case of Lagrange + // 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); + 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) = + 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); + + // step 3: change the basis back to the given finite element + // space. we can use a shortcut for elements that define support + // points, in which case we can evaluate the Lagrange polynomials of + // the Gauss quadrature in those points. this will give more + // accurate results than the inversion of a matrix. for more general + // polynomial spaces, we must invert a matrix of a Vandermonde type, + // which we do by a Householder transformation to keep roundoff + // errors low. + inverse_shape_values.resize_fast(shape_values.size()); + FullMatrix transform_from_gauss(n_dofs_1d, n_dofs_1d); + if (fe.has_support_points()) + { + 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) = poly_project[j].value( + fe.get_unit_support_points()[lexicographic[i]][0]); + FullMatrix result(n_dofs_1d, n_q_points_1d); + transform_from_gauss.mmult(result, project_to_dof_space); + + // set very small entries to zero - we are in reference space + // with normalized numbers, so this is straight-forward to check + // here + for (unsigned int i = 0; i < n_dofs_1d; ++i) + for (unsigned int q = 0; q < n_q_points_1d; ++q) + inverse_shape_values[i * n_q_points_1d + q] = + std::abs(result(i, q)) < 1e-15 ? 0 : result(i, q); + } + else + { + for (unsigned int i = 0; i < n_dofs_1d; ++i) + for (unsigned int j = 0; j < n_dofs_1d; ++j) + { + Point q_point; + q_point[0] = quad_project.point(i)[0]; + + transform_from_gauss(i, j) = + fe.shape_value(lexicographic[j], q_point); + } + Householder H(transform_from_gauss); + Vector in(n_dofs_1d), out(n_dofs_1d); + for (unsigned int q = 0; q < n_q_points_1d; ++q) + { + for (unsigned int i = 0; i < n_dofs_1d; ++i) + in(i) = project_to_dof_space(i, q); + H.least_squares(out, in); + for (unsigned int i = 0; i < n_dofs_1d; ++i) + inverse_shape_values[i * n_q_points_1d + q] = + std::abs(out(i)) < 1e-15 ? 0. : out(i); + } + } + quadrature_data_on_face[0].resize(quad.size() * 3); + quadrature_data_on_face[1].resize(quad.size() * 3); + + for (unsigned int i = 0; i < quad.size(); ++i) + { + std::array values; + poly_coll[i].value(0.0, 2, values.data()); + for (unsigned int d = 0; d < 3; ++d) + quadrature_data_on_face[0][i + d * quad.size()] = values[d]; + poly_coll[i].value(1.0, 2, values.data()); + for (unsigned int d = 0; d < 3; ++d) + quadrature_data_on_face[1][i + d * quad.size()] = values[d]; + } + } + + + template bool - ShapeInfo::check_1d_shapes_symmetric( - UnivariateShapeData &univariate_shape_data) + UnivariateShapeData::check_and_set_shapes_symmetric() { - if (dofs_per_component_on_cell == 0) + if (fe_degree == 0) return false; - const auto n_q_points_1d = univariate_shape_data.n_q_points_1d; - const auto fe_degree = univariate_shape_data.fe_degree; - auto & shape_values = univariate_shape_data.shape_values; - auto & shape_gradients = univariate_shape_data.shape_gradients; - auto & shape_hessians = univariate_shape_data.shape_hessians; - auto & shape_gradients_collocation = - univariate_shape_data.shape_gradients_collocation; - auto &shape_hessians_collocation = - univariate_shape_data.shape_hessians_collocation; - auto &shape_values_eo = univariate_shape_data.shape_values_eo; - auto &shape_gradients_eo = univariate_shape_data.shape_gradients_eo; - auto &shape_hessians_eo = univariate_shape_data.shape_hessians_eo; - auto &shape_gradients_collocation_eo = - univariate_shape_data.shape_gradients_collocation_eo; - auto &shape_hessians_collocation_eo = - univariate_shape_data.shape_hessians_collocation_eo; - auto &inverse_shape_values = univariate_shape_data.inverse_shape_values; - auto &inverse_shape_values_eo = - univariate_shape_data.inverse_shape_values_eo; - const double zero_tol = std::is_same_v == true ? 1e-12 : 1e-7; // symmetry for values @@ -1196,15 +1179,11 @@ namespace internal template bool - ShapeInfo::check_1d_shapes_collocation( - const UnivariateShapeData &univariate_shape_data) const + UnivariateShapeData::check_shapes_collocation() const { - if (dofs_per_component_on_cell != n_q_points) + if (fe_degree + 1 != n_q_points_1d) return false; - const auto fe_degree = univariate_shape_data.fe_degree; - auto & shape_values = univariate_shape_data.shape_values; - const double zero_tol = std::is_same_v == true ? 1e-12 : 1e-7; // check: identity operation for shape values -- 2.39.5