Merge pull request #14822 from kronbichler/fix_mapping_q_cache_1

Fix bug in MappingQCache(1) with cell similarity

doc/news/changes/minor/20230222Kronbichler

@@ -0,0 +1,6 @@
+Fixed: When used with the linear-polynomial MappingQCache(1),
+FEValues::reinit() could previously trigger CellSimilarity in case the
+underlying mesh was affine, but the deformed state described via MappingQCache
+was not. This is now fixed.
+<br>
+(Martin Kronbichler, 2023/02/22)

source/fe/mapping_q.cc

@@ -974,7 +974,9 @@ MappingQ<dim, spacedim>::fill_fe_values(
   // value is computed with just cell vertices and does not take into account
   // cell curvature.
   const CellSimilarity::Similarity computed_cell_similarity =
-    (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
+    (polynomial_degree == 1 && this->preserves_vertex_locations() ?
+       cell_similarity :
+       CellSimilarity::none);
 
   if (dim > 1 && data.tensor_product_quadrature)
     {

tests/mappings/mapping_q_cache_09.cc

@@ -0,0 +1,78 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test that MappingQCache based on MappingQ1 works correctly in case some
+// deformation is applied and CellSimilarity would get activated for the
+// undeformed geometry, but not in the deformed state
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q_cache.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+do_test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, 0, 1);
+  tria.refine_global(4 - dim);
+
+  MappingQ<dim>      mapping(1);
+  MappingQCache<dim> mapping_cache(1);
+  const auto transform = [](const typename Triangulation<dim>::cell_iterator &,
+                            const Point<dim> &in) {
+    Point<dim> out;
+    // move points more densely near zero
+    for (unsigned int d = 0; d < dim; ++d)
+      out[d] = std::pow(in[d], 1.2 + 0.2 * d);
+    return out;
+  };
+  mapping_cache.initialize(mapping, tria, transform, false);
+
+  const Vector<double> aspect_ratios =
+    GridTools::compute_aspect_ratio_of_cells(mapping, tria, QGauss<dim>(2));
+  const Vector<double> aspect_ratios_deformed =
+    GridTools::compute_aspect_ratio_of_cells(mapping_cache,
+                                             tria,
+                                             QGauss<dim>(2));
+
+  deallog.push(std::to_string(dim) + "d");
+  deallog << "Aspect ratios undeformed mesh: " << std::endl;
+  aspect_ratios.print(deallog.get_file_stream());
+  deallog << std::endl;
+  deallog << "Aspect ratios deformed mesh: " << std::endl;
+  aspect_ratios_deformed.print(deallog.get_file_stream());
+  deallog << std::endl;
+  deallog.pop();
+  deallog << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+  MultithreadInfo::set_thread_limit(1);
+  do_test<1>();
+  do_test<2>();
+  do_test<3>();
+}

tests/mappings/mapping_q_cache_09.output

@@ -0,0 +1,22 @@
+
+DEAL:1d::Aspect ratios undeformed mesh: 
+1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 
+DEAL:1d::
+DEAL:1d::Aspect ratios deformed mesh: 
+1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 
+DEAL:1d::
+DEAL::
+DEAL:2d::Aspect ratios undeformed mesh: 
+1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 
+DEAL:2d::
+DEAL:2d::Aspect ratios deformed mesh: 
+1.320e+00 1.712e+00 1.242e+00 1.044e+00 1.900e+00 2.033e+00 1.159e+00 1.240e+00 1.528e+00 1.178e+00 1.750e+00 1.349e+00 1.061e+00 1.008e+00 1.215e+00 1.136e+00 
+DEAL:2d::
+DEAL::
+DEAL:3d::Aspect ratios undeformed mesh: 
+1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 
+DEAL:3d::
+DEAL:3d::Aspect ratios deformed mesh: 
+1.320e+00 1.712e+00 1.883e+00 1.883e+00 1.768e+00 1.768e+00 1.540e+00 1.187e+00 
+DEAL:3d::
+DEAL::