From: Matthias Maier <tamiko@43-1.org> Date: Mon, 16 May 2022 17:36:16 +0000 (-0500) Subject: TransfiniteInterpolationManifold: fix codimension one case X-Git-Tag: v9.4.0-rc1~208^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bddd1966f8fc2ddac7dce358c5c1bfd947d1b605;p=dealii.git TransfiniteInterpolationManifold: fix codimension one case The codimension one case currently fails with an assertion: ``` An error occurred in line <290> of file <source/grid/grid_tools.cc> in function std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>> dealii::GridTools::affine_cell_approximation(const ArrayView<const Point<spacedim>> &) [dim = 2, spacedim = 3] The violated condition was: ::dealii::deal_II_exceptions::internals::compare_for_equality(vertices.size(), GeometryInfo<dim>::vertices_per_cell) Additional information: Dimension 9 not equal to 4. ``` This is triggered because we take a compatibility branch in case of `dim < spacedim` (in `include/deal.II/fe/mapping_q_internal.h`): ``` 908 if (real_support_points.size() == 909 GeometryInfo<dim>::vertices_per_cell || 910 dim < spacedim) 911 { ``` but still initialize the vector `unit_points` by a subdivided quadrature (which is needed for the default quadratic approximation; in `source/grid/manifold_lib.cc`): ``` 1686 std::vector<Point<dim>> unit_points = 1687 QIterated<dim>(QTrapez<1>(), 2).get_points(); 1688 std::vector<Point<spacedim>> real_points(unit_points.size()); ``` Fix this by simply initializing `unit_points` with the right number of interpolation points for the `dim < spacedim` variant. --- diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index ade514c553..2174e1480b 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -901,13 +901,17 @@ namespace internal 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 (real_support_points.size() == - GeometryInfo<dim>::vertices_per_cell || - dim < spacedim) + // polynomial due to a lack of points (interpolation matrix would + // get singular). Similarly, it is not entirely clear how to gather + // enough information for the case dim < spacedim. + // + // In both cases we require the vector real_support_points to + // contain the vertex positions and fall back to an affine + // approximation: + Assert(dim == spacedim || real_support_points.size() == + GeometryInfo<dim>::vertices_per_cell, + ExcInternalError()); + if (real_support_points.size() == GeometryInfo<dim>::vertices_per_cell) { const auto affine = GridTools::affine_cell_approximation<dim>( make_array_view(real_support_points)); diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 1b9ad3730c..085f8e8248 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -1673,7 +1673,7 @@ TransfiniteInterpolationManifold<dim, spacedim>::initialize( const Triangulation<dim, spacedim> &triangulation) { this->triangulation = &triangulation; - // in case the triangulation is cleared, remove the pointers by a signal + // In case the triangulation is cleared, remove the pointers by a signal: clear_signal.disconnect(); clear_signal = triangulation.signals.clear.connect([&]() -> void { this->triangulation = nullptr; @@ -1683,8 +1683,16 @@ TransfiniteInterpolationManifold<dim, spacedim>::initialize( coarse_cell_is_flat.resize(triangulation.n_cells(level_coarse), false); quadratic_approximation.clear(); + // In case of dim == spacedim we perform a quadratic approximation in + // InverseQuadraticApproximation(), thus initialize the unit_points + // vector with one subdivision to get 3^dim unit_points. + // + // In the co-dimension one case (meaning dim < spacedim) we have to fall + // back to a simple GridTools::affine_cell_approximation<dim>() which + // requires 2^dim points, instead. Thus, initialize the QIteraded + // quadrature with no subdivisions. std::vector<Point<dim>> unit_points = - QIterated<dim>(QTrapez<1>(), 2).get_points(); + QIterated<dim>(QTrapez<1>(), (dim == spacedim ? 2 : 1)).get_points(); std::vector<Point<spacedim>> real_points(unit_points.size()); for (const auto &cell : triangulation.active_cell_iterators())