Added unit_to_real and real_to_unit to BoundingBox

diff --git a/doc/news/changes/minor/20200629LucaHeltai b/doc/news/changes/minor/20200629LucaHeltai
new file mode 100644 (file)
index 0000000..e187db4
--- /dev/null
+++ b/doc/news/changes/minor/20200629LucaHeltai
@@ -0,0 +1,5 @@
+New: BoundingBox::real_to_unit() and BoundingBox::unit_to_real() allow one to
+apply both the direct and the inverse transformation that are needed to map the
+unit bounding box to the current box, and viceversa.
+<br>
+(Luca Heltai, 2020/06/29)

diff --git a/include/deal.II/base/bounding_box.h b/include/deal.II/base/bounding_box.h
index 62bdbbc9d7ef7380e54ad296e6366498c9df7af2..36d515d40e07a0c67929ab50c34d69c32059581d 100644 (file)
--- a/include/deal.II/base/bounding_box.h
+++ b/include/deal.II/base/bounding_box.h
@@ -263,6 +263,28 @@ public:
   BoundingBox<spacedim - 1, Number>
   cross_section(const unsigned int direction) const;
 
+  /**
+   * Apply the affine transformation that transforms this BoundingBox to a unit
+   * BoundingBox object.
+   *
+   * If $B$ is this bounding box, and $\hat{B}$ is the unit bounding box,
+   * compute the affine mapping that satisfies $G(B) = \hat{B}$ and apply it to
+   * @p point.
+   */
+  Point<spacedim, Number>
+  real_to_unit(const Point<spacedim, Number> &point) const;
+
+  /**
+   * Apply the affine transformation that transforms the unit BoundingBox object
+   * to this object.
+   *
+   * If $B$ is this bounding box, and $\hat{B}$ is the unit bounding box,
+   * compute the affine mapping that satisfies $F(\hat{B}) = B$ and apply it to
+   * @p point.
+   */
+  Point<spacedim, Number>
+  unit_to_real(const Point<spacedim, Number> &point) const;
+
   /**
    * Boost serialization function
    */

diff --git a/source/base/bounding_box.cc b/source/base/bounding_box.cc
index c071e52b32b485398ef99a2de6b57f4c79b1377a..73ff6d088aea7f4814c0d340dbde6eff48c49e78 100644 (file)
--- a/source/base/bounding_box.cc
+++ b/source/base/bounding_box.cc
@@ -315,6 +315,35 @@ BoundingBox<spacedim, Number>::cross_section(const unsigned int direction) const
 
 
 
+template <int spacedim, typename Number>
+Point<spacedim, Number>
+BoundingBox<spacedim, Number>::real_to_unit(
+  const Point<spacedim, Number> &point) const
+{
+  auto       unit = point;
+  const auto diag = boundary_points.second - boundary_points.first;
+  unit -= boundary_points.first;
+  for (unsigned int d = 0; d < spacedim; ++d)
+    unit[d] /= diag[d];
+  return unit;
+}
+
+
+
+template <int spacedim, typename Number>
+Point<spacedim, Number>
+BoundingBox<spacedim, Number>::unit_to_real(
+  const Point<spacedim, Number> &point) const
+{
+  auto       real = boundary_points.first;
+  const auto diag = boundary_points.second - boundary_points.first;
+  for (unsigned int d = 0; d < spacedim; ++d)
+    real[d] += diag[d] * point[d];
+  return real;
+}
+
+
+
 template <int dim, typename Number>
 BoundingBox<dim, Number>
 create_unit_bounding_box()

diff --git a/tests/base/bounding_box_7.cc b/tests/base/bounding_box_7.cc
new file mode 100644 (file)
index 0000000..b0a2fe2
--- /dev/null
+++ b/tests/base/bounding_box_7.cc
@@ -0,0 +1,85 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2020 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 the following functions of the BoundingBox class
+// unit_to_real
+// real_to_unit
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/std_cxx20/iota_view.h>
+
+#include <deal.II/boost_adaptors/bounding_box.h>
+#include <deal.II/boost_adaptors/point.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+using iota = std_cxx20::ranges::iota_view<unsigned int, unsigned int>;
+
+template <int dim>
+void
+test()
+{
+  const unsigned int n_boxes  = 3;
+  const unsigned int n_points = 3;
+
+  const auto unit_box = create_unit_bounding_box<dim>();
+
+  for (auto i : iota(0, n_boxes))
+    {
+      const auto box = random_box<dim>();
+      // Check that all vertices get transformed correctly
+      for (auto v : iota(0, GeometryInfo<dim>::vertices_per_cell))
+        {
+          const auto vertex      = unit_box.vertex(v);
+          const auto real_vertex = box.unit_to_real(vertex);
+          if (real_vertex.distance(box.vertex(v)) > 1e-10)
+            deallog << "Error: " << vertex << " -> " << real_vertex
+                    << " != " << box.vertex(v) << std::endl;
+
+          const auto unit_vertex = box.real_to_unit(real_vertex);
+          if (unit_vertex.distance(vertex) > 1e-10)
+            deallog << "Error: " << real_vertex << " -> " << unit_vertex
+                    << " != " << vertex << std::endl;
+        }
+
+      // Now check random points
+      for (auto j : iota(0, n_points))
+        {
+          const auto unit_point   = random_point<dim>();
+          const auto real_point   = box.unit_to_real(unit_point);
+          const auto real_to_unit = box.real_to_unit(real_point);
+
+          if (real_to_unit.distance(unit_point) > 1e-10)
+            deallog << "Error: " << unit_point << " -> " << real_point << " -> "
+                    << real_to_unit << " != " << unit_point << std::endl;
+        }
+    }
+  deallog << "Spacedim " << dim << " OK" << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+}

diff --git a/tests/base/bounding_box_7.output b/tests/base/bounding_box_7.output
new file mode 100644 (file)
index 0000000..a44332f
--- /dev/null
+++ b/tests/base/bounding_box_7.output
@@ -0,0 +1,4 @@
+
+DEAL::Spacedim 1 OK
+DEAL::Spacedim 2 OK
+DEAL::Spacedim 3 OK