From 9982d375ffbdf8abe4b6612517de707ae2631e5c Mon Sep 17 00:00:00 2001
From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Mon, 4 Jul 2022 11:57:14 -0400
Subject: [PATCH] Use correct tolerance in MappingCartesian check.

---
 source/fe/mapping_cartesian.cc             | 17 +++-
 tests/mappings/mapping_cartesian_03.cc     | 94 ++++++++++++++++++++++
 tests/mappings/mapping_cartesian_03.output |  5 ++
 3 files changed, 112 insertions(+), 4 deletions(-)
 create mode 100644 tests/mappings/mapping_cartesian_03.cc
 create mode 100644 tests/mappings/mapping_cartesian_03.output

diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc
index 63020f1516..47349a5dd1 100644
--- a/source/fe/mapping_cartesian.cc
+++ b/source/fe/mapping_cartesian.cc
@@ -55,12 +55,21 @@ is_cartesian(const CellType &cell)
   if (!cell->reference_cell().is_hyper_cube())
     return false;
 
-  const double tolerance    = 1e-14 * cell->diameter();
-  const auto   bounding_box = cell->bounding_box();
+  const double tolerance         = 1e-14;
+  const auto   bounding_box      = cell->bounding_box();
+  const auto & bounding_vertices = bounding_box.get_boundary_points();
+  const auto   cell_measure =
+    bounding_vertices.first.distance_square(bounding_vertices.second);
 
   for (const unsigned int v : cell->vertex_indices())
-    if (cell->vertex(v).distance(bounding_box.vertex(v)) > tolerance)
-      return false;
+    {
+      const double vertex_tolerance =
+        tolerance * std::max(cell->vertex(v).norm_square(), cell_measure);
+
+      if (cell->vertex(v).distance_square(bounding_box.vertex(v)) >
+          vertex_tolerance)
+        return false;
+    }
 
   return true;
 }
diff --git a/tests/mappings/mapping_cartesian_03.cc b/tests/mappings/mapping_cartesian_03.cc
new file mode 100644
index 0000000000..78c79b1721
--- /dev/null
+++ b/tests/mappings/mapping_cartesian_03.cc
@@ -0,0 +1,94 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2003 - 2018 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// there used to be a bug in MappingCartesian, where an assert was triggered
+// for valid cells if the cells where too small compared to the dimension
+// of the triangulation. The bug would not occur for coarse grids,
+// but would reliably be triggered with higher refinement level.
+// In this test we move the origin of the box to trigger the bug
+// easier, but it would also occur with a zero origin (just requiring higher
+// refinement levels).
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_cartesian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+check(const Triangulation<dim> &tria)
+{
+  MappingCartesian<dim> mapping;
+  FE_Q<dim>             fe(1);
+  DoFHandler<dim>       dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  QGauss<dim> quadrature(1);
+
+  UpdateFlags update_flags = update_quadrature_points | update_values;
+
+  FEValues<dim> fe_values(mapping, fe, quadrature, update_flags);
+
+  for (auto cell : dof_handler.active_cell_iterators())
+    fe_values.reinit(cell);
+
+  deallog << std::endl;
+
+  deallog << "OK" << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  {
+    Triangulation<2>          coarse_grid;
+    std::vector<unsigned int> rep_vec({9, 2});
+    Point<2>                  box_origin(3000000., 660000.);
+    Point<2>                  extents(3000000., 660000.);
+    GridGenerator::subdivided_hyper_rectangle(
+      coarse_grid, rep_vec, box_origin, box_origin + extents, true);
+    coarse_grid.refine_global(2);
+    check(coarse_grid);
+  }
+
+  {
+    Triangulation<3>          coarse_grid;
+    std::vector<unsigned int> rep_vec({9, 9, 2});
+    Point<3>                  box_origin(3000000., 3000000., 660000.);
+    Point<3>                  extents(3000000., 3000000., 660000.);
+    GridGenerator::subdivided_hyper_rectangle(
+      coarse_grid, rep_vec, box_origin, box_origin + extents, true);
+    coarse_grid.refine_global(2);
+    check(coarse_grid);
+  }
+}
diff --git a/tests/mappings/mapping_cartesian_03.output b/tests/mappings/mapping_cartesian_03.output
new file mode 100644
index 0000000000..5a6edd9d09
--- /dev/null
+++ b/tests/mappings/mapping_cartesian_03.output
@@ -0,0 +1,5 @@
+
+DEAL::
+DEAL::OK
+DEAL::
+DEAL::OK
-- 
2.39.5