From: danshapero <shapero.daniel@gmail.com>
Date: Tue, 26 Sep 2017 21:49:06 +0000 (-0700)
Subject: Added functions to compute bounding box for triangulation, accessors
X-Git-Tag: v9.0.0-rc1~986^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6219117d71ee761da222df2ed3c1b7e489282ec8;p=dealii.git

Added functions to compute bounding box for triangulation, accessors
---

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index c03c8eb75b..aea8291e4b 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -18,6 +18,7 @@
 
 
 #include <deal.II/base/config.h>
+#include <deal.II/base/bounding_box.h>
 #include <deal.II/base/geometry_info.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/fe/mapping.h>
@@ -170,6 +171,22 @@ namespace GridTools
   template <int dim, typename T>
   double cell_measure (const T &, ...);
 
+  /**
+   * Compute the smallest box containing the entire triangulation.
+   *
+   * If the input triangulation is a `parallel::distributed::Triangulation`,
+   * then each processor will compute a bounding box enclosing all locally
+   * owned, ghost, and artificial cells. In the case of a domain without curved
+   * boundaries, these bounding boxes will all agree between processors because
+   * the union of the areas occupied by artificial and ghost cells equals the
+   * union of the areas occupied by the cells that other processors own.
+   * However, if the domain has curved boundaries, this is no longer the case.
+   * The bounding box returned may be appropriate for the current processor,
+   * but different from the bounding boxes computed on other processors.
+   */
+  template <int dim, int spacedim>
+  BoundingBox<spacedim> compute_bounding_box(const Triangulation<dim, spacedim> &triangulation);
+
   /*@}*/
   /**
    * @name Functions supporting the creation of meshes
diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h
index d45031260f..81828d2a3b 100644
--- a/include/deal.II/grid/tria_accessor.h
+++ b/include/deal.II/grid/tria_accessor.h
@@ -17,6 +17,7 @@
 #define dealii_tria_accessor_h
 
 
+#include <deal.II/base/bounding_box.h>
 #include <deal.II/base/config.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/geometry_info.h>
@@ -1300,6 +1301,11 @@ public:
    */
   std::pair<Point<spacedim>,double> enclosing_ball () const;
 
+  /**
+   * Return the smallest bounding box that encloses the object.
+   */
+  BoundingBox<spacedim> bounding_box () const;
+
   /**
    * Length of an object in the direction of the given axis, specified in the
    * local coordinate system. See the documentation of GeometryInfo for the
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index e492b2fef4..6fd3160354 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -364,6 +364,22 @@ namespace GridTools
   }
 
 
+
+  template <int dim, int spacedim>
+  BoundingBox<spacedim>
+  compute_bounding_box(const Triangulation<dim, spacedim> &tria)
+  {
+    using iterator = typename Triangulation<dim, spacedim>::active_cell_iterator;
+    const auto predicate = [](const iterator &)
+    {
+      return true;
+    };
+
+    return compute_bounding_box(tria, std::function<bool(const iterator &)>(predicate));
+  }
+
+
+
   template <int dim, int spacedim>
   void
   delete_unused_vertices (std::vector<Point<spacedim> >    &vertices,
diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index 2b3f29b8e1..df32b3a31c 100644
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -159,6 +159,11 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
     (const Triangulation<deal_II_dimension, deal_II_space_dimension> &,
      const Mapping<deal_II_dimension, deal_II_space_dimension> &);
 
+    template
+    BoundingBox<deal_II_space_dimension>
+    compute_bounding_box
+    (const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
+
     template
     void delete_unused_vertices (std::vector<Point<deal_II_space_dimension> > &,
                                  std::vector<CellData<deal_II_dimension> > &,
diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc
index 106b9615a9..74cf9b5fa9 100644
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -1110,6 +1110,29 @@ measure () const
 
 
 
+template <int structdim, int dim, int spacedim>
+BoundingBox<spacedim>
+TriaAccessor<structdim, dim, spacedim>::
+bounding_box () const
+{
+  std::pair<Point<spacedim>, Point<spacedim> > boundary_points =
+    std::make_pair(this->vertex(0), this->vertex(0));
+
+  for (unsigned int v=1; v<GeometryInfo<structdim>::vertices_per_cell; ++v)
+    {
+      const Point<spacedim> &x = this->vertex(v);
+      for (unsigned int k=0; k<spacedim; ++k)
+        {
+          boundary_points.first[k] = std::min(boundary_points.first[k], x[k]);
+          boundary_points.second[k] = std::max(boundary_points.second[k], x[k]);
+        }
+    }
+
+  return BoundingBox<spacedim>(boundary_points);
+}
+
+
+
 template <>
 double TriaAccessor<1,1,1>::extent_in_direction(const unsigned int axis) const
 {
diff --git a/tests/grid/grid_tools_bounding_box_03.cc b/tests/grid/grid_tools_bounding_box_03.cc
new file mode 100644
index 0000000000..a5939dce2d
--- /dev/null
+++ b/tests/grid/grid_tools_bounding_box_03.cc
@@ -0,0 +1,56 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check computing bounding boxes of entire triangulations and individual cells
+
+#include "../tests.h"
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+template <int dim, int spacedim>
+void test_tria_bounding_box()
+{
+  dealii::Point<dim> p1, p2;
+  for (unsigned int k = 0; k < dim; ++k)
+    {
+      p1[k] = k;
+      p2[k] = 2 * k + 1;
+    }
+
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::hyper_rectangle(tria, p1, p2);
+  tria.refine_global(1);
+  const BoundingBox<spacedim> bounding_box = GridTools::compute_bounding_box(tria);
+  const std::pair<Point<spacedim>, Point<spacedim> > &boundary_points =
+    bounding_box.get_boundary_points();
+
+  deallog << boundary_points.first << ", " << boundary_points.second << std::endl;
+}
+
+int main(void)
+{
+  initlog();
+
+  test_tria_bounding_box<1, 1>();
+  test_tria_bounding_box<1, 2>();
+  test_tria_bounding_box<2, 2>();
+  test_tria_bounding_box<1, 3>();
+  test_tria_bounding_box<2, 3>();
+  test_tria_bounding_box<3, 3>();
+
+  return 0;
+}
+
diff --git a/tests/grid/grid_tools_bounding_box_03.output b/tests/grid/grid_tools_bounding_box_03.output
new file mode 100644
index 0000000000..e7c415a7d1
--- /dev/null
+++ b/tests/grid/grid_tools_bounding_box_03.output
@@ -0,0 +1,7 @@
+
+DEAL::0.00000, 1.00000
+DEAL::0.00000 0.00000, 1.00000 0.00000
+DEAL::0.00000 1.00000, 1.00000 3.00000
+DEAL::0.00000 0.00000 0.00000, 1.00000 0.00000 0.00000
+DEAL::0.00000 1.00000 0.00000, 1.00000 3.00000 0.00000
+DEAL::0.00000 1.00000 2.00000, 1.00000 3.00000 5.00000
diff --git a/tests/grid/measure_et_al_02.cc b/tests/grid/measure_et_al_02.cc
index 464871e186..31ce5f90bf 100644
--- a/tests/grid/measure_et_al_02.cc
+++ b/tests/grid/measure_et_al_02.cc
@@ -100,6 +100,11 @@ void test()
               << tria.begin_active()->extent_in_direction(0) << std::endl;
       deallog << "dim" << dim << ":case" << case_no << ":minimum_vertex_distance="
               << tria.begin_active()->minimum_vertex_distance() << std::endl;
+
+      const BoundingBox<dim> box = tria.begin_active()->bounding_box();
+      const std::pair<Point<dim>, Point<dim> > &pts = box.get_boundary_points();
+      deallog << "dim" << dim << ":case" << case_no << ":bounding_box="
+              << pts.first << ", " << pts.second << std::endl;
       tria.clear();
     }
 }
diff --git a/tests/grid/measure_et_al_02.output b/tests/grid/measure_et_al_02.output
index e1db82bb4a..d7ec0015ed 100644
--- a/tests/grid/measure_et_al_02.output
+++ b/tests/grid/measure_et_al_02.output
@@ -2,18 +2,24 @@
 DEAL::dim2:case0:diameter=2.8284
 DEAL::dim2:case0:extent_in_direction=2.0000
 DEAL::dim2:case0:minimum_vertex_distance=2.0000
+DEAL::dim2:case0:bounding_box=1.0000 1.0000, 3.0000 3.0000
 DEAL::dim2:case1:diameter=3.6056
 DEAL::dim2:case1:extent_in_direction=3.0000
 DEAL::dim2:case1:minimum_vertex_distance=2.0000
+DEAL::dim2:case1:bounding_box=0.0000 1.0000, 3.0000 3.0000
 DEAL::dim2:case2:diameter=4.4721
 DEAL::dim2:case2:extent_in_direction=3.0000
 DEAL::dim2:case2:minimum_vertex_distance=2.2361
+DEAL::dim2:case2:bounding_box=0.0000 1.0000, 4.0000 3.0000
 DEAL::dim3:case0:diameter=3.4641
 DEAL::dim3:case0:extent_in_direction=2.0000
 DEAL::dim3:case0:minimum_vertex_distance=2.0000
+DEAL::dim3:case0:bounding_box=1.0000 1.0000 1.0000, 3.0000 3.0000 3.0000
 DEAL::dim3:case1:diameter=4.1231
 DEAL::dim3:case1:extent_in_direction=3.0000
 DEAL::dim3:case1:minimum_vertex_distance=2.0000
+DEAL::dim3:case1:bounding_box=0.0000 1.0000 1.0000, 3.0000 3.0000 3.0000
 DEAL::dim3:case2:diameter=4.1231
 DEAL::dim3:case2:extent_in_direction=3.0000
 DEAL::dim3:case2:minimum_vertex_distance=2.0000
+DEAL::dim3:case2:bounding_box=0.0000 1.0000 1.0000, 3.0000 3.0000 3.0000