From: Martin Kronbichler Date: Sun, 1 Nov 2020 09:09:04 +0000 (+0100) Subject: Transfinite interpolation: Use inverse quadratic approximation for pull_back X-Git-Tag: v9.3.0-rc1~948^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5022ca913fbb5703a5c44e1cd45b9ca18154510f;p=dealii.git Transfinite interpolation: Use inverse quadratic approximation for pull_back --- diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index bb3c28ca59..d4d86f9080 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -996,12 +996,18 @@ namespace internal } } + /** + * Copy constructor. + */ + InverseQuadraticApproximation(const InverseQuadraticApproximation &) = + default; + /** * Evaluate the quadratic approximation. */ template Point - compute(const Point &p) + compute(const Point &p) const { Point result; for (unsigned int d = 0; d < dim; ++d) diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index e8b4d86339..ba1fe83e7b 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -26,6 +26,17 @@ DEAL_II_NAMESPACE_OPEN +// forward declaration +namespace internal +{ + namespace MappingQGenericImplementation + { + template + class InverseQuadraticApproximation; + } +} // namespace internal + + /** * Manifold description for a polar coordinate system. * @@ -1133,6 +1144,14 @@ private: */ FlatManifold chart_manifold; + /** + * A vector of quadratic approximations to the inverse map from real points + * to chart points for each of the coarse mesh cells. + */ + std::vector> + quadratic_approximation; + /** * The connection to Triangulation::signals::clear that must be reset once * this class goes out of scope. diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index a2110f226b..91f03bed62 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -17,6 +17,7 @@ #include #include +#include #include #include @@ -1682,10 +1683,13 @@ TransfiniteInterpolationManifold::initialize( }); level_coarse = triangulation.last()->level(); coarse_cell_is_flat.resize(triangulation.n_cells(level_coarse), false); - typename Triangulation::active_cell_iterator - cell = triangulation.begin(level_coarse), - endc = triangulation.end(level_coarse); - for (; cell != endc; ++cell) + quadratic_approximation.clear(); + + std::vector> unit_points = + QIterated(QTrapez<1>(), 2).get_points(); + std::vector> real_points(unit_points.size()); + + for (const auto &cell : triangulation.active_cell_iterators()) { bool cell_is_flat = true; for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) @@ -1700,6 +1704,11 @@ TransfiniteInterpolationManifold::initialize( AssertIndexRange(static_cast(cell->index()), coarse_cell_is_flat.size()); coarse_cell_is_flat[cell->index()] = cell_is_flat; + + // build quadratic interpolation + for (unsigned int i = 0; i < unit_points.size(); ++i) + real_points[i] = push_forward(cell, unit_points[i]); + quadratic_approximation.emplace_back(real_points, unit_points); } } @@ -2281,7 +2290,7 @@ TransfiniteInterpolationManifold:: for (unsigned int i = 0; i < points.size(); ++i) { Point point = - cell->real_to_unit_cell_affine_approximation(points[i]); + quadratic_approximation[cell->index()].compute(points[i]); current_distance += GeometryInfo::distance_to_unit_cell(point); } distances_and_cells.push_back( @@ -2434,11 +2443,11 @@ TransfiniteInterpolationManifold::compute_chart_points( [&](const typename Triangulation::cell_iterator &cell, const unsigned int point_index) { Point guess; - // an optimization: keep track of whether or not we used the affine + // an optimization: keep track of whether or not we used the quadratic // approximation so that we don't call pull_back with the same // initial guess twice (i.e., if pull_back fails the first time, // don't try again with the same function arguments). - bool used_affine_approximation = false; + bool used_quadratic_approximation = false; // if we have already computed three points, we can guess the fourth // to be the missing corner point of a rectangle if (point_index == 3 && surrounding_points.size() >= 8) @@ -2468,9 +2477,9 @@ TransfiniteInterpolationManifold::compute_chart_points( } else { - guess = cell->real_to_unit_cell_affine_approximation( + guess = quadratic_approximation[cell->index()].compute( surrounding_points[point_index]); - used_affine_approximation = true; + used_quadratic_approximation = true; } chart_points[point_index] = pull_back(cell, surrounding_points[point_index], guess); @@ -2480,9 +2489,9 @@ TransfiniteInterpolationManifold::compute_chart_points( // than the cheap methods used above) if (chart_points[point_index][0] == internal::invalid_pull_back_coordinate && - !used_affine_approximation) + !used_quadratic_approximation) { - guess = cell->real_to_unit_cell_affine_approximation( + guess = quadratic_approximation[cell->index()].compute( surrounding_points[point_index]); chart_points[point_index] = pull_back(cell, surrounding_points[point_index], guess);