--- /dev/null
+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)
*/
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
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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> ¢er= 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>();
+}
+
+
--- /dev/null
+
+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!
projects / dealii.git / commitdiff