]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add TriaAccessor<>::enclosing_ball() 4320/head
authorvishalkenchan <vishal.boddu@fau.de>
Wed, 26 Apr 2017 12:47:57 +0000 (14:47 +0200)
committervishalkenchan <vishal.boddu@fau.de>
Fri, 28 Apr 2017 12:33:20 +0000 (14:33 +0200)
doc/news/changes/minor/20170426VishalBoddu [new file with mode: 0644]
include/deal.II/grid/tria_accessor.h
include/deal.II/grid/tria_accessor.templates.h
tests/grid/enclosing_sphere_01.cc [new file with mode: 0644]
tests/grid/enclosing_sphere_01.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170426VishalBoddu b/doc/news/changes/minor/20170426VishalBoddu
new file mode 100644 (file)
index 0000000..dc80c93
--- /dev/null
@@ -0,0 +1,5 @@
+New: TriaAccessor::enclosing_ball() computes and return a pair of Point 
+and double corresponding to the center and the radius of a reasonably small 
+enclosing ball of the TriaAccessor object.
+<br>
+(Vishal Boddu, Denis Davydov 2017/04/26)
index 5f2f0a9e63f85f074e82223fd6e8d2854ee48b9d..1c6a7e2d6bb955bf34685f98752fac0166809044 100644 (file)
@@ -1237,6 +1237,34 @@ public:
    */
   double diameter () const;
 
+  /**
+   * Return a pair of Point and double corresponding to the center and
+   * the radius of a reasonably small enclosing ball of the object.
+   *
+   * The function implements Ritter's O(n) algorithm to get a reasonably
+   * small enclosing ball around the vertices of the object.
+   * The initial guess for the enclosing ball is taken to be the ball
+   * which contains the largest diagonal of the object as its diameter.
+   * Starting from such an initial guess, the algorithm tests whether all
+   * the vertices of the object (except the vertices of the largest diagonal)
+   * are geometrically within the ball.
+   * If any vertex (v) is found to be geometrically outside the ball,
+   * a new iterate of the ball is constructed by shifting its center and
+   * increasing the radius so as to geometrically enclose both the previous
+   * ball and the vertex (v). The algorithm terminates when all the vertices
+   * are geometrically inside the ball.
+   *
+   * If a vertex (v) is geometrically inside a particular iterate of the ball,
+   * then it will continue to be so in the subsequent iterates of
+   * the ball (this is true \a by \a construction).
+   *
+   * @note This function assumes d-linear mapping from the reference cell.
+   *
+   * <a href="http://geomalgorithms.com/a08-_containers.html">see this</a> and
+   * [Ritter 1990]
+   */
+  std::pair<Point<spacedim>,double> enclosing_ball () 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
index 37310ca79ab25a51eba9c674ba5051118a82c013..7d1a6a2b17c460eca9eeb3cc0fc520172fd02bee 100644 (file)
@@ -2040,6 +2040,140 @@ TriaAccessor<structdim, dim, spacedim>::diameter () const
 
 
 
+template <int structdim, int dim, int spacedim>
+std::pair<Point<spacedim>,double>
+TriaAccessor<structdim, dim, spacedim>::enclosing_ball () const
+{
+  // If the object is one dimensional,
+  // the enclosing ball is the initial iterate
+  // i.e., the ball's center and diameter are
+  // the center and the diameter of the object.
+  if (structdim==1)
+    return std::make_pair( (this->vertex(1)+this->vertex(0))*0.5,
+                           (this->vertex(1)-this->vertex(0)).norm()*0.5);
+
+  // The list is_initial_guess_vertex contains bool values and has
+  // the same size as the number of vertices per object.
+  // The entries of is_initial_guess_vertex are set true only for those
+  // two vertices corresponding to the largest diagonal which is being used
+  // to construct the initial ball.
+  // We employ this mask to skip these two vertices while enlarging the ball.
+  std::array<bool, GeometryInfo<structdim>::vertices_per_cell> is_initial_guess_vertex;
+
+  //First let all the vertices be outside
+  std::fill(is_initial_guess_vertex.begin(),
+            is_initial_guess_vertex.end(),
+            false);
+
+  // Get an initial guess by looking at the largest diagonal
+  Point<spacedim> center;
+  double radius = 0;
+
+  switch (structdim)
+    {
+    case 2:
+    {
+      const Point<spacedim> p30( this->vertex(3)-this->vertex(0));
+      const Point<spacedim> p21( this->vertex(2)-this->vertex(1));
+      if (p30.norm() > p21.norm())
+        {
+          center = this->vertex(0) + 0.5* p30;
+          radius = p30.norm()/2.;
+          is_initial_guess_vertex[3] = true;
+          is_initial_guess_vertex[0] = true;
+        }
+      else
+        {
+          center = this->vertex(1) + 0.5* p21;
+          radius = p21.norm()/2.;
+          is_initial_guess_vertex[2] = true;
+          is_initial_guess_vertex[1] = true;
+        }
+      break;
+    }
+    case 3:
+    {
+      const Point<spacedim> p70( this->vertex(7)-this->vertex(0));
+      const Point<spacedim> p61( this->vertex(6)-this->vertex(1));
+      const Point<spacedim> p25( this->vertex(2)-this->vertex(5));
+      const Point<spacedim> p34( this->vertex(3)-this->vertex(4));
+      const std::vector<double> diagonals= { p70.norm(),
+                                             p61.norm(),
+                                             p25.norm(),
+                                             p34.norm()
+                                           };
+      const std::vector<double>::const_iterator
+      it = std::max_element( diagonals.begin(), diagonals.end());
+      if (it == diagonals.begin())
+        {
+          center = this->vertex(0) + 0.5* p70;
+          is_initial_guess_vertex[7] = true;
+          is_initial_guess_vertex[0] = true;
+        }
+      else if (it == diagonals.begin()+1)
+        {
+          center = this->vertex(1) + 0.5* p61;
+          is_initial_guess_vertex[6] = true;
+          is_initial_guess_vertex[1] = true;
+        }
+      else if (it == diagonals.begin()+2)
+        {
+          center = this->vertex(5) + 0.5* p25;
+          is_initial_guess_vertex[2] = true;
+          is_initial_guess_vertex[5] = true;
+        }
+      else
+        {
+          center = this->vertex(4) + 0.5* p34;
+          is_initial_guess_vertex[3] = true;
+          is_initial_guess_vertex[4] = true;
+        }
+      radius = *it * 0.5;
+      break;
+    }
+    default:
+      Assert (false, ExcNotImplemented());
+      return std::pair<Point<spacedim>,double>();
+    }
+
+  // For each vertex that is found to be geometrically outside the ball
+  // enlarge the ball  so that the new ball contains both the previous ball
+  // and the given vertex.
+  for (unsigned int v = 0; v < GeometryInfo<structdim>::vertices_per_cell; ++v)
+    if (!is_initial_guess_vertex[v])
+      {
+        const double distance = center.distance(this->vertex(v));
+        if (distance > radius)
+          {
+            // we found a vertex which is outside of the ball
+            // extend it (move center and change radius)
+            const Point<spacedim> pCV (center - this->vertex(v));
+            radius = (distance + radius) * 0.5;
+            center = this->vertex(v) + pCV * (radius / distance);
+
+            // Now the new ball constructed in this block
+            // encloses the vertex (v) that was found to be geometrically
+            // outside the old ball.
+          }
+      }
+#ifdef DEBUG
+  bool all_vertices_within_ball = true;
+
+  // Set all_vertices_within_ball false if any of the vertices of the object
+  // are geometrically outside the ball
+  for (unsigned int v = 0; v < GeometryInfo<structdim>::vertices_per_cell; ++v)
+    if (center.distance(this->vertex(v)) > radius + 100. *std::numeric_limits<double>::epsilon())
+      {
+        all_vertices_within_ball = false;
+        break;
+      }
+  // If all the vertices are not within the ball throw error
+  Assert (all_vertices_within_ball, ExcInternalError());
+#endif
+  return std::make_pair(center, radius);
+}
+
+
 template <int structdim, int dim, int spacedim>
 double
 TriaAccessor<structdim, dim, spacedim>::minimum_vertex_distance () const
diff --git a/tests/grid/enclosing_sphere_01.cc b/tests/grid/enclosing_sphere_01.cc
new file mode 100644 (file)
index 0000000..36537c5
--- /dev/null
@@ -0,0 +1,218 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Computes a reasonably small enclosing ball on a variety of cells.
+// The design of this test and part of the blessed file are taken
+// from measure_et_al_02 test.
+
+
+#include "../tests.h"
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <fstream>
+#include <iomanip>
+#include <limits>
+
+#define PRECISION 5
+
+
+template<int dim>
+void create_triangulation(const unsigned int,
+                          Triangulation<dim> &)
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template<>
+void create_triangulation(const unsigned int case_no,
+                          Triangulation<2> &tria)
+{
+  switch (case_no)
+    {
+    case 0:
+      GridGenerator::hyper_cube(tria, 1., 3.);
+      break;
+    case 1:
+    {
+      GridGenerator::hyper_cube(tria, 1., 3.);
+      Point<2> &v0=tria.begin_active()->vertex(0);
+      v0(0) = 0.;
+      break;
+    }
+    case 2:
+    {
+      GridGenerator::hyper_cube(tria, 1., 3.);
+      Point<2> &v0=tria.begin_active()->vertex(0);
+      v0(0) = 0.;
+      Point<2> &v3=tria.begin_active()->vertex(3);
+      v3(0) = 4.;
+      break;
+    }
+    case 3:
+    {
+      GridGenerator::hyper_cube(tria, 1., 3.);
+      Point<2> &v0=tria.begin_active()->vertex(0);
+      v0 = Point<2>(1.9,1.9);
+
+      Point<2> &v3=tria.begin_active()->vertex(3);
+      v3 = Point<2>(3.1,3.);
+      break;
+
+      //
+      //  y^  2-------3
+      //   |  |       |
+      //   |  |       |
+      //   |  |       |
+      //   |  0-------1
+      //   *------------>x
+      //
+      //        |
+      //        v
+      //
+      // vertices 0 and 3 are moved
+      // such that the initial guess for the ball
+      // (with its diameter as the largest diagonal (between vertices 1-2))
+      // doesnot enclose vertex 3
+      //
+      //  y^  2--------+3
+      //   |  |        |
+      //   |  |        |
+      //   |  | 0      |
+      //   |  +--------1
+      //   *------------>x
+    }
+    default:
+      Assert(false, ExcNotImplemented());
+    };
+}
+
+
+template<>
+void create_triangulation(const unsigned int case_no,
+                          Triangulation<3> &tria)
+{
+  switch (case_no)
+    {
+    case 0:
+      GridGenerator::hyper_cube(tria, 1., 3.);
+      break;
+    case 1:
+    case 2:// like case 1
+    {
+      GridGenerator::hyper_cube(tria, 1., 3.);
+      Point<3> &v0=tria.begin_active()->vertex(0);
+      v0(0) = 0.;
+      break;
+    }
+    case 3:
+    {
+
+      //       6-------7        6-------7
+      //      /|       |       /       /|
+      //     / |       |      /       / |
+      //    /  |       |     /       /  |
+      //   4   |       |    4-------5   |
+      //   |   2-------3    |       |   3
+      //   |  /       /     |       |  /
+      //   | /       /      |       | /
+      //   |/       /       |       |/
+      //   0-------1        0-------1
+      //
+      //
+      //                |
+      //                v
+      //
+      // moving vertices 1 and 6 such that diagonal 16 is smaller than all other
+      // diagonals but vertex 6 lies outside all the balls constructed
+      // using largest diagonals
+      //
+      //
+      //      6                6
+      //       +-------7        +-------7
+      //      /|       |       /       /|
+      //     / |       |      /       / |
+      //    /  |       |     /       /  |
+      //   4   |       |    4-------5   |
+      //   |   +-------3    |       |   3
+      //   |  /       /     |       |  /
+      //   | /       /      |       | /
+      //   |/     1 /       |       |/
+      //   0-------+        0-------+
+      //
+
+      GridGenerator::hyper_cube(tria, 1., 3.);
+      Point<3> &v1=tria.begin_active()->vertex(1);
+      Point<3> &v6=tria.begin_active()->vertex(6);
+      v1 += Point<3>(-0.9,0.9,0.9); // v1 was (3.,1.,1.)
+      v6 += Point<3>(-0.2,0.2,0.2); // v7 was (1.,3.,3.)
+      const Point<3> &v0 =tria.begin_active()->vertex(0);
+      const Point<3> &v7 =tria.begin_active()->vertex(7);
+      AssertThrow( v1.distance(v6) < v0.distance(v7),
+                   ExcMessage("Vertices not moved as drawn above"));
+      break;
+    }
+    default:
+      Assert(false, ExcNotImplemented());
+    };
+}
+
+
+template<int dim>
+void test()
+{
+  Triangulation<dim> tria;
+  for (unsigned int case_no=0; case_no<4; ++case_no)
+    {
+      create_triangulation(case_no, tria);
+      const std::pair<Point<dim>, double>
+      smallest_sphere = tria.begin_active()->enclosing_ball();
+      const double &radius  = smallest_sphere.second;
+      const Point<dim> &center= smallest_sphere.first;
+
+      deallog << "dim" << dim << ":case" << case_no << ":diameter="
+              << radius << ":center="
+              << center   << std::endl;
+
+      // Check that all the vertices are within the sphere
+      // (sphere with thickness 100. *std::numeric_limits<double>::epsilon())
+      for ( int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
+        AssertThrow( std::fabs(center.distance(tria.begin_active()->vertex(v)))
+                     < radius + 100. *std::numeric_limits<double>::epsilon(),
+                     ExcInternalError());
+
+      tria.clear();
+    }
+  deallog << "PASSED!" << std::endl;
+}
+
+
+int main()
+{
+  std::ofstream logfile ("output");
+  deallog << std::setprecision (PRECISION);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+}
+
+
diff --git a/tests/grid/enclosing_sphere_01.output b/tests/grid/enclosing_sphere_01.output
new file mode 100644 (file)
index 0000000..f223dbd
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::dim2:case0:diameter=1.4142:center=2.0000 2.0000
+DEAL::dim2:case1:diameter=1.8028:center=1.5000 2.0000
+DEAL::dim2:case2:diameter=2.2361:center=2.0000 2.0000
+DEAL::dim2:case3:diameter=1.4504:center=2.0268 2.0243
+DEAL::PASSED!
+DEAL::dim3:case0:diameter=1.7321:center=2.0000 2.0000 2.0000
+DEAL::dim3:case1:diameter=2.0616:center=1.5000 2.0000 2.0000
+DEAL::dim3:case2:diameter=2.0616:center=1.5000 2.0000 2.0000
+DEAL::dim3:case3:diameter=1.9053:center=1.9000 2.1000 2.1000
+DEAL::PASSED!

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.