From 347beb386ea8cb74f8532ab8ce716553a60feb1a Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 2 Oct 2020 12:15:09 +0200 Subject: [PATCH] Multiple points in Mapping::transform_real_to_unit_cell --- include/deal.II/fe/mapping.h | 17 ++++++ include/deal.II/fe/mapping_q_generic.h | 7 +++ source/fe/mapping.cc | 24 ++++++++ source/fe/mapping_q_generic.cc | 84 ++++++++++++++++++++------ 4 files changed, 114 insertions(+), 18 deletions(-) diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index 8b302d8ba7..ac4f6652b0 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -452,6 +452,23 @@ public: const typename Triangulation::cell_iterator &cell, const Point & p) const = 0; + /** + * Map multiple points from the real point locations to points in reference + * locations. The functionality is essentially the same as looping over all + * points and calling the Mapping::transform_real_to_unit_cell() function + * for each point individually, but it can be much faster for certain + * mappings that implement a more specialized version such as + * MappingQGeneric. The only difference in behavior is that this function + * will never throw an ExcTransformationFailed() exception. If the + * transformation fails for `real_points[i]`, the returned `unit_points[i]` + * contains std::numeric_limits::infinity() as the first entry. + */ + virtual void + transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> & real_points, + ArrayView> &unit_points) const; + /** * Transform the point @p p on the real @p cell to the corresponding point * on the reference cell, and then project this point to a (dim-1)-dimensional diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index 6dec094643..501ef16289 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -181,6 +181,13 @@ public: const typename Triangulation::cell_iterator &cell, const Point &p) const override; + // for documentation, see the Mapping base class + virtual void + transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> & real_points, + ArrayView> &unit_points) const override; + /** * @} */ diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc index 5433bf7471..a70503a51b 100644 --- a/source/fe/mapping.cc +++ b/source/fe/mapping.cc @@ -76,6 +76,30 @@ Mapping::get_bounding_box( +template +void +Mapping::transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> & real_points, + ArrayView> & unit_points) const +{ + AssertDimension(real_points.size(), unit_points.size()); + for (unsigned int i = 0; i < real_points.size(); ++i) + { + try + { + unit_points[i] = transform_real_to_unit_cell(cell, real_points[i]); + } + catch (typename Mapping::ExcTransformationFailed &) + { + unit_points[i] = Point(); + unit_points[i][0] = std::numeric_limits::infinity(); + } + } +} + + + template Point Mapping::project_real_point_to_unit_point_on_face( diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 9e4eab6b04..ab1058c97e 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -893,17 +894,15 @@ namespace internal /** * Implementation of transform_real_to_unit_cell for dim==spacedim */ - template + template Point do_transform_real_to_unit_cell_internal( - const Point & p, + const Point & p, const Point & initial_p_unit, - const std::vector> & points, + const std::vector> & points, const std::vector> &polynomials_1d, const std::vector & renumber) { - const unsigned int spacedim = dim; - AssertDimension(points.size(), Utilities::pow(polynomials_1d.size(), dim)); @@ -917,8 +916,8 @@ namespace internal // The shape values and derivatives of the mapping at this point are // previously computed. - Point p_unit = initial_p_unit; - std::pair, Tensor<2, dim>> p_real = + Point p_unit = initial_p_unit; + std::pair, Tensor<2, spacedim>> p_real = compute_mapped_location_of_point(points, polynomials_1d, renumber, @@ -967,6 +966,9 @@ namespace internal const double eps = 1.e-11; const unsigned int newton_iteration_limit = 20; + Point invalid_point; + invalid_point[0] = std::numeric_limits::infinity(); + unsigned int newton_iteration = 0; double last_f_weighted_norm_square; do @@ -979,9 +981,9 @@ namespace internal const Tensor<2, spacedim> &df = p_real.second; // Solve [f'(x)]d=f(x) - AssertThrow( - determinant(df) > 0, - (typename Mapping::ExcTransformationFailed())); + if (determinant(df) <= 0) + return invalid_point; + const Tensor<2, spacedim> df_inverse = invert(df); const Tensor<1, spacedim> delta = df_inverse * f; @@ -1035,18 +1037,13 @@ namespace internal else if (step_length > 0.05) step_length /= 2; else - AssertThrow( - false, - (typename Mapping::ExcTransformationFailed())); + return invalid_point; } while (true); ++newton_iteration; if (newton_iteration > newton_iteration_limit) - AssertThrow( - false, - (typename Mapping::ExcTransformationFailed())); + return invalid_point; last_f_weighted_norm_square = (df_inverse * f).norm_square(); } while (last_f_weighted_norm_square > eps * eps); @@ -2295,7 +2292,58 @@ MappingQGeneric::transform_real_to_unit_cell( // perform the Newton iteration and return the result. note that this // statement may throw an exception, which we simply pass up to the caller - return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit); + const Point p_unit = + this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit); + if (p_unit[0] == std::numeric_limits::infinity()) + AssertThrow(false, + (typename Mapping::ExcTransformationFailed())); + return p_unit; +} + + + +template +void +MappingQGeneric::transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> & real_points, + ArrayView> & unit_points) const +{ + AssertDimension(real_points.size(), unit_points.size()); + const std::vector> support_points = + this->compute_mapping_support_points(cell); + + // From the chosen (high-order) support points, now only pick the first + // 2^dim points and construct an affine approximation from those. + const std::pair, Tensor<1, spacedim>> + affine_factors = GridTools::affine_cell_approximation( + ArrayView>(support_points.data(), + GeometryInfo::vertices_per_cell)); + const DerivativeForm<1, spacedim, dim> A_inv = + affine_factors.first.covariant_form().transpose(); + + for (unsigned int i = 0; i < real_points.size(); ++i) + { + try + { + // Compute an initial guess by inverting the affine approximation + // A * x_unit + b = x_real + Point initial_guess( + apply_transformation(A_inv, + real_points[i] - affine_factors.second)); + unit_points[i] = internal::MappingQGenericImplementation:: + do_transform_real_to_unit_cell_internal( + real_points[i], + GeometryInfo::project_to_unit_cell(initial_guess), + support_points, + polynomials_1d, + renumber_lexicographic_to_hierarchic); + } + catch (typename Mapping::ExcTransformationFailed &) + { + unit_points[i][0] = std::numeric_limits::infinity(); + } + } } -- 2.39.5