From: Martin Kronbichler Date: Wed, 28 Oct 2020 08:47:10 +0000 (+0100) Subject: Construct approximation by two sets of points. X-Git-Tag: v9.3.0-rc1~959^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11085%2Fhead;p=dealii.git Construct approximation by two sets of points. Also better document some variables and mark constant variables as const. --- diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index 8e21cecfdb..1512df3227 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -606,26 +606,38 @@ protected: /* * The default line support points. These are used when computing the * location in real space of the support points on lines and quads, which - * are asked to the Manifold class. + * are needed by the Manifold class. * * The number of points depends on the degree of this class, and it matches * the number of degrees of freedom of an FE_Q<1>(this->degree). */ - std::vector> line_support_points; + const std::vector> line_support_points; /* * The one-dimensional polynomials defined as Lagrange polynomials from the * line support points. These are used for point evaluations and match the * polynomial space of an FE_Q<1>(this->degree). */ - std::vector> polynomials_1d; + const std::vector> polynomials_1d; /* * The numbering from the lexicographic to the hierarchical ordering used * when expanding the tensor product with the mapping support points (which * come in hierarchical numbers). */ - std::vector renumber_lexicographic_to_hierarchic; + const std::vector renumber_lexicographic_to_hierarchic; + + /* + * The support points in reference coordinates. These are used for + * constructing approximations of the output of + * compute_mapping_support_points() when evaluating the mapping on the fly, + * rather than going through the FEValues interface provided by + * InternalData. + * + * The number of points depends on the degree of this class, and it matches + * the number of degrees of freedom of an FE_Q(this->degree). + */ + const std::vector> unit_cell_support_points; /** * A vector of tables of weights by which we multiply the locations of the @@ -646,7 +658,8 @@ protected: * For the definition of this table see equation (8) of the `mapping' * report. */ - std::vector> support_point_weights_perimeter_to_interior; + const std::vector> + support_point_weights_perimeter_to_interior; /** * A table of weights by which we multiply the locations of the vertex @@ -660,7 +673,7 @@ protected: * in 2D, 8 in 3D), and as many rows as there are additional support points * in the mapping, i.e., (degree+1)^dim - 2^dim. */ - Table<2, double> support_point_weights_cell; + const Table<2, double> support_point_weights_cell; /** * Return the locations of support points for the mapping. For example, for diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index 88f2efbb52..9ea8c9c231 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -206,6 +206,34 @@ namespace internal */ namespace MappingQGenericImplementation { + /** + * This function generates the reference cell support points from the 1d + * support points by expanding the tensor product. + */ + template + std::vector> + unit_support_points(const std::vector> & line_support_points, + const std::vector &renumbering) + { + AssertDimension(Utilities::pow(line_support_points.size(), dim), + renumbering.size()); + std::vector> points(renumbering.size()); + const unsigned int n1 = line_support_points.size(); + for (unsigned int q2 = 0, q = 0; q2 < (dim > 2 ? n1 : 1); ++q2) + for (unsigned int q1 = 0; q1 < (dim > 1 ? n1 : 1); ++q1) + for (unsigned int q0 = 0; q0 < n1; ++q0, ++q) + { + points[renumbering[q]][0] = line_support_points[q0][0]; + if (dim > 1) + points[renumbering[q]][1] = line_support_points[q1][0]; + if (dim > 2) + points[renumbering[q]][2] = line_support_points[q2][0]; + } + return points; + } + + + /** * This function is needed by the constructor of * MappingQ for dim= 2 and 3. @@ -844,30 +872,36 @@ namespace internal /** * Constructor. + * + * @param real_support_points The position of the mapping support points + * in real space, queried by + * MappingQGeneric::compute_mapping_support_points(). + * + * @param unit_support_points The location of the support points in + * reference coordinates $[0, 1]^d$ that map to the mapping support + * points in real space by a polynomial map. */ InverseQuadraticApproximation( - const std::vector> &mapping_support_points, - const std::vector> & line_support_points, - const std::vector & renumber) - : p_shift(mapping_support_points[0]) - , scale(1. / - mapping_support_points[0].distance(mapping_support_points[1])) + const std::vector> &real_support_points, + const std::vector> & unit_support_points) + : normalization_shift(real_support_points[0]) + , normalization_length( + 1. / real_support_points[0].distance(real_support_points[1])) + , is_affine(true) { - AssertDimension(mapping_support_points.size(), renumber.size()); - AssertDimension(mapping_support_points.size(), - Utilities::pow(line_support_points.size(), dim)); - - const unsigned int n1 = line_support_points.size(); + AssertDimension(real_support_points.size(), unit_support_points.size()); // For the bi-/trilinear approximation, we cannot build a quadratic // polynomial due to a lack of points (interpolation matrix would get // singular), so pick the affine approximation. Similarly, it is not // entirely clear how to gather enough information for the case dim < // spacedim - if (n1 == 2 || dim < spacedim) + if (real_support_points.size() == + GeometryInfo::vertices_per_cell || + dim < spacedim) { const auto affine = GridTools::affine_cell_approximation( - make_array_view(mapping_support_points)); + make_array_view(real_support_points)); DerivativeForm<1, spacedim, dim> A_inv = affine.first.covariant_form().transpose(); coefficients[0] = apply_transformation(A_inv, affine.second); @@ -880,40 +914,32 @@ namespace internal SymmetricTensor<2, n_functions> matrix; std::array shape_values; - for (unsigned int q2 = 0, q = 0; q2 < (dim > 2 ? n1 : 1); ++q2) - for (unsigned int q1 = 0; q1 < (dim > 1 ? n1 : 1); ++q1) - for (unsigned int q0 = 0; q0 < n1; ++q0, ++q) - { - // Evaluate quadratic shape functions in point, shifted to the - // first support point and scaled by the length between the - // first two support points to avoid roundoff issues with - // scaling far away from 1. - shape_values[0] = 1.; - const Tensor<1, spacedim> p_scaled = - (mapping_support_points[renumber[q]] - p_shift) * scale; - for (unsigned int d = 0; d < spacedim; ++d) - shape_values[1 + d] = p_scaled[d]; - for (unsigned int d = 0, c = 0; d < spacedim; ++d) - for (unsigned int e = 0; e <= d; ++e, ++c) - shape_values[1 + spacedim + c] = p_scaled[d] * p_scaled[e]; - - // Build lower diagonal of least squares matrix and rhs, the - // essential part being that we construct the matrix with the - // real points and the right hand side by comparing to the - // reference point positions which sets up an inverse - // interpolation. - for (unsigned int i = 0; i < n_functions; ++i) - for (unsigned int j = 0; j <= i; ++j) - matrix[i][j] += shape_values[i] * shape_values[j]; - Point reference_point; - reference_point[0] = line_support_points[q0][0]; - if (dim > 1) - reference_point[1] = line_support_points[q1][0]; - if (dim > 2) - reference_point[2] = line_support_points[q2][0]; - for (unsigned int i = 0; i < n_functions; ++i) - coefficients[i] += shape_values[i] * reference_point; - } + for (unsigned int q = 0; q < unit_support_points.size(); ++q) + { + // Evaluate quadratic shape functions in point, with the + // normalization applied in order to avoid roundoff issues with + // scaling far away from 1. + shape_values[0] = 1.; + const Tensor<1, spacedim> p_scaled = + normalization_length * + (real_support_points[q] - normalization_shift); + for (unsigned int d = 0; d < spacedim; ++d) + shape_values[1 + d] = p_scaled[d]; + for (unsigned int d = 0, c = 0; d < spacedim; ++d) + for (unsigned int e = 0; e <= d; ++e, ++c) + shape_values[1 + spacedim + c] = p_scaled[d] * p_scaled[e]; + + // Build lower diagonal of least squares matrix and rhs, the + // essential part being that we construct the matrix with the + // real points and the right hand side by comparing to the + // reference point positions which sets up an inverse + // interpolation. + for (unsigned int i = 0; i < n_functions; ++i) + for (unsigned int j = 0; j <= i; ++j) + matrix[i][j] += shape_values[i] * shape_values[j]; + for (unsigned int i = 0; i < n_functions; ++i) + coefficients[i] += shape_values[i] * unit_support_points[q]; + } // Factorize the matrix A = L * L^T in-place with the // Cholesky-Banachiewicz algorithm. The implementation is similar to @@ -981,11 +1007,13 @@ namespace internal for (unsigned int d = 0; d < dim; ++d) result[d] = coefficients[0][d]; - // Shift point to avoid roundoff problems. Spell out the loop - // because Number might be a vectorized array. + // Apply the normalization to ensure a good conditioning. Since Number + // might be a vectorized array whereas the normalization is a point of + // doubles, we cannot use the overload of operator- and must instead + // loop over the components of the point. Point p_scaled; for (unsigned int d = 0; d < spacedim; ++d) - p_scaled[d] = (p[d] - p_shift[d]) * scale; + p_scaled[d] = (p[d] - normalization_shift[d]) * normalization_length; for (unsigned int d = 0; d < spacedim; ++d) result += coefficients[1 + d] * p_scaled[d]; @@ -1001,10 +1029,31 @@ namespace internal } private: - const Point p_shift; - const double scale; + /** + * In order to guarantee a good conditioning, we need to apply a + * transformation to the points in real space that is computed by a + * shift vector normalization_shift (first point of the mapping support + * points in real space) and an inverse length scale called + * `length_normalization` as the distance between the first two points. + */ + const Point normalization_shift; + + /** + * See the documentation of `normalization_shift` above. + */ + const double normalization_length; + + /** + * The vector of coefficients in the quadratic approximation. + */ std::array, n_functions> coefficients; - bool is_affine; + + /** + * In case the quadratic approximation is not possible due to an + * insufficient number of support points, we switch to an affine + * approximation that always works but is less accurate. + */ + bool is_affine; }; diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 4a2d3d3292..f9c772999a 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -368,6 +368,10 @@ MappingQGeneric::MappingQGeneric(const unsigned int p) Polynomials::generate_complete_Lagrange_basis(line_support_points)) , renumber_lexicographic_to_hierarchic( FETools::lexicographic_to_hierarchic_numbering(p)) + , unit_cell_support_points( + internal::MappingQGenericImplementation::unit_support_points( + line_support_points, + renumber_lexicographic_to_hierarchic)) , support_point_weights_perimeter_to_interior( internal::MappingQGenericImplementation:: compute_support_point_weights_perimeter_to_interior( @@ -741,9 +745,7 @@ MappingQGeneric::transform_points_real_to_unit_cell( // 2^dim points and construct an affine approximation from those. internal::MappingQGenericImplementation:: InverseQuadraticApproximation - inverse_approximation(support_points, - line_support_points, - renumber_lexicographic_to_hierarchic); + inverse_approximation(support_points, unit_cell_support_points); const unsigned int n_points = real_points.size(); const unsigned int n_lanes = VectorizedArray::size(); diff --git a/tests/mappings/mapping_q_inverse_quadratic_approximation.cc b/tests/mappings/mapping_q_inverse_quadratic_approximation.cc index 8972746fe7..411db4c8f7 100644 --- a/tests/mappings/mapping_q_inverse_quadratic_approximation.cc +++ b/tests/mappings/mapping_q_inverse_quadratic_approximation.cc @@ -59,6 +59,9 @@ print_result(const unsigned int mapping_degree, QGaussLobatto<1>(mapping_degree + 1).get_points()); std::vector renumber = FETools::lexicographic_to_hierarchic_numbering(mapping_degree); + std::vector> mapping_unit_support_points = + internal::MappingQGenericImplementation::unit_support_points( + mapping_points, renumber); for (const auto &cell : tria.active_cell_iterators()) { @@ -76,8 +79,7 @@ print_result(const unsigned int mapping_degree, internal::MappingQGenericImplementation:: InverseQuadraticApproximation approx(fe_values.get_quadrature_points(), - mapping_points, - renumber); + mapping_unit_support_points); deallog << "Inverse quadratic approximation: " << approx.compute(p) << std::endl << std::endl;