From 91321380462d532da61aa36f2a325c2d6e93bd56 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <martin.kronbichler@uni-a.de>
Date: Wed, 22 Feb 2023 09:26:10 +0100
Subject: [PATCH] Fix bug in MappingQCache(1) with cell similarity

---
 source/fe/mapping_q.cc                   |  4 +-
 tests/mappings/mapping_q_cache_09.cc     | 78 ++++++++++++++++++++++++
 tests/mappings/mapping_q_cache_09.output | 22 +++++++
 3 files changed, 103 insertions(+), 1 deletion(-)
 create mode 100644 tests/mappings/mapping_q_cache_09.cc
 create mode 100644 tests/mappings/mapping_q_cache_09.output

diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc
index c32a86edb1..d945c3a066 100644
--- a/source/fe/mapping_q.cc
+++ b/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)
     {
diff --git a/tests/mappings/mapping_q_cache_09.cc b/tests/mappings/mapping_q_cache_09.cc
new file mode 100644
index 0000000000..937ead3c20
--- /dev/null
+++ b/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>();
+}
diff --git a/tests/mappings/mapping_q_cache_09.output b/tests/mappings/mapping_q_cache_09.output
new file mode 100644
index 0000000000..ba7eaca25a
--- /dev/null
+++ b/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::
-- 
2.39.5