From b45d3e6118aa9a11bbab550c3067537bcf5c9b97 Mon Sep 17 00:00:00 2001 From: vishalkenchan Date: Wed, 26 Apr 2017 14:47:57 +0200 Subject: [PATCH] add TriaAccessor<>::enclosing_ball() --- doc/news/changes/minor/20170426VishalBoddu | 5 + include/deal.II/grid/tria_accessor.h | 28 +++ .../deal.II/grid/tria_accessor.templates.h | 134 +++++++++++ tests/grid/enclosing_sphere_01.cc | 218 ++++++++++++++++++ tests/grid/enclosing_sphere_01.output | 11 + 5 files changed, 396 insertions(+) create mode 100644 doc/news/changes/minor/20170426VishalBoddu create mode 100644 tests/grid/enclosing_sphere_01.cc create mode 100644 tests/grid/enclosing_sphere_01.output diff --git a/doc/news/changes/minor/20170426VishalBoddu b/doc/news/changes/minor/20170426VishalBoddu new file mode 100644 index 0000000000..dc80c93763 --- /dev/null +++ b/doc/news/changes/minor/20170426VishalBoddu @@ -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. +
+(Vishal Boddu, Denis Davydov 2017/04/26) diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index 5f2f0a9e63..1c6a7e2d6b 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -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. + * + * see this and + * [Ritter 1990] + */ + std::pair,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 diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index 37310ca79a..7d1a6a2b17 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -2040,6 +2040,140 @@ TriaAccessor::diameter () const +template +std::pair,double> +TriaAccessor::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::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 center; + double radius = 0; + + switch (structdim) + { + case 2: + { + const Point p30( this->vertex(3)-this->vertex(0)); + const Point 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 p70( this->vertex(7)-this->vertex(0)); + const Point p61( this->vertex(6)-this->vertex(1)); + const Point p25( this->vertex(2)-this->vertex(5)); + const Point p34( this->vertex(3)-this->vertex(4)); + const std::vector diagonals= { p70.norm(), + p61.norm(), + p25.norm(), + p34.norm() + }; + const std::vector::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,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::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 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::vertices_per_cell; ++v) + if (center.distance(this->vertex(v)) > radius + 100. *std::numeric_limits::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 double TriaAccessor::minimum_vertex_distance () const diff --git a/tests/grid/enclosing_sphere_01.cc b/tests/grid/enclosing_sphere_01.cc new file mode 100644 index 0000000000..36537c561c --- /dev/null +++ b/tests/grid/enclosing_sphere_01.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + +#define PRECISION 5 + + +template +void create_triangulation(const unsigned int, + Triangulation &) +{ + 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 +void test() +{ + Triangulation tria; + for (unsigned int case_no=0; case_no<4; ++case_no) + { + create_triangulation(case_no, tria); + const std::pair, double> + smallest_sphere = tria.begin_active()->enclosing_ball(); + const double &radius = smallest_sphere.second; + const Point ¢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::epsilon()) + for ( int v = 0; v < GeometryInfo::vertices_per_cell; ++v) + AssertThrow( std::fabs(center.distance(tria.begin_active()->vertex(v))) + < radius + 100. *std::numeric_limits::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 index 0000000000..f223dbdf09 --- /dev/null +++ b/tests/grid/enclosing_sphere_01.output @@ -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! -- 2.39.5