Merge pull request #16713 from nmuch/colorize-hyper-rectangle-with-si…

…mplices

Colorize hyper rectangle with simplices

source/grid/grid_generator.cc

@@ -2634,7 +2634,7 @@ namespace GridGenerator
         // use a large epsilon to
         // compare numbers to avoid
         // roundoff problems.
-        double epsilon = 10;
+        double epsilon = std::numeric_limits<double>::max();
         for (unsigned int i = 0; i < dim; ++i)
           epsilon = std::min(epsilon, 0.01 * delta[i][i]);
         Assert(epsilon > 0,
@@ -8645,8 +8645,6 @@ namespace GridGenerator
   {
     AssertDimension(dim, spacedim);
 
-    AssertThrow(colorize == false, ExcNotImplemented());
-
     std::vector<Point<spacedim>> vertices;
     std::vector<CellData<dim>>   cells;
 
@@ -8782,11 +8780,34 @@ namespace GridGenerator
       }
     else
       {
-        AssertThrow(colorize == false, ExcNotImplemented());
+        AssertThrow(false, ExcNotImplemented());
       }
 
     // actually create triangulation
     tria.create_triangulation(vertices, cells, SubCellData());
+
+    if (colorize)
+      {
+        // to colorize, run through all
+        // faces of all cells and set
+        // boundary indicator to the
+        // correct value if it was 0.
+
+        // use a large epsilon to
+        // compare numbers to avoid
+        // roundoff problems.
+        double epsilon = std::numeric_limits<double>::max();
+        for (unsigned int i = 0; i < dim; ++i)
+          epsilon = std::min(epsilon,
+                             0.01 * (std::abs(p2[i] - p1[i]) / repetitions[i]));
+        Assert(epsilon > 0,
+               ExcMessage(
+                 "The distance between corner points must be positive."));
+
+        // actual code is external since
+        // 1-D is different from 2/3d.
+        colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
+      }
   }
 
 

tests/grid/grid_generator_11.cc

@@ -0,0 +1,78 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2016 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// Test GridGenerator::subdivided_hyper_rectangle_with_simplices
+
+#include <deal.II/base/point.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test(std::ostream &out)
+{
+  Point<dim> p1;
+  p1[0] = 2.;
+  if (dim > 1)
+    p1[1] = -1.;
+  if (dim > 2)
+    p1[2] = 0.;
+  Point<dim> p2;
+  p2[0] = 3.;
+  if (dim > 1)
+    p2[1] = 2.;
+  if (dim > 2)
+    p2[2] = 4.;
+
+  GridOut go;
+  go.set_flags(GridOutFlags::Msh(true));
+
+  if (dim > 1)
+    {
+      deallog << "subdivided_hyper_rectangle_with_simplices" << std::endl;
+      Triangulation<dim>        tr;
+      std::vector<unsigned int> sub(dim);
+      sub[0] = 2;
+      if (dim > 1)
+        sub[1] = 3;
+      if (dim > 2)
+        sub[2] = 4;
+      GridGenerator::subdivided_hyper_rectangle_with_simplices(
+        tr, sub, p1, p2, true);
+      if (tr.n_cells() > 0)
+        go.write_msh(tr, out);
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("2d");
+  test<2>(deallog.get_file_stream());
+  deallog.pop();
+  deallog.push("3d");
+  test<3>(deallog.get_file_stream());
+  deallog.pop();
+}