From: wolf Date: Thu, 30 Oct 2003 14:23:51 +0000 (+0000) Subject: GridTools::find_active_cell_around_point X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8c9b77b201edf0ca4dce37fae25c7dccf4159c7f;p=dealii-svn.git GridTools::find_active_cell_around_point git-svn-id: https://svn.dealii.org/trunk@8174 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/grid_tools.h b/deal.II/deal.II/include/grid/grid_tools.h index be9414dfe0..6623d7ae3f 100644 --- a/deal.II/deal.II/include/grid/grid_tools.h +++ b/deal.II/deal.II/include/grid/grid_tools.h @@ -23,11 +23,12 @@ /** - * This class is a collection of algorithms working on - * triangulations. See the descriptions of the individual functions - * for more information. + * This class is a collection of algorithms working on triangulations, + * such as shifting or rotating triangulations, but also finding a + * cell that contains a given point. See the descriptions of the + * individual functions for more information. * - * @author Wolfgang Bangerth, 2001 + * @author Wolfgang Bangerth, 2001, 2003 */ class GridTools { @@ -142,18 +143,125 @@ class GridTools static void scale (const double scaling_factor, Triangulation &triangulation); + + /** + * Find and return an iterator to + * the active cell that surrounds + * a given point @p{ref}. The + * type of the first parameter + * may be either + * @ref{Triangulation}, + * @ref{DoFHandler}, or + * @ref{MGDoFHandler}, i.e. we + * can find the cell around a + * point for iterators into each + * of these classes. + * + * The algorithm used in this + * function proceeds by first + * looking for the surrounding + * cell on the coarse grid, and + * then recursively checking its + * sibling cells. The complexity + * is thus @p{O(M+log N)} where + * @p{M} is the number of coarse + * grid cells, and @p{N} the + * total number of cells. + * + * There are cases where this + * function will not found a + * given point in space + * dimensions higher than one, + * even though it is inside the + * domain being discretized, or + * will find a point that is + * actually outside the + * domain. The reason for this is + * that we use piecewise d-linear + * mappings of the unit cell to + * real cells. Thus, if a point + * is close to a convex boundary + * or on it, it may not be inside + * any of the cells since they + * have straight boundaries that + * lie entirely inside the + * domain. + * + * Another case for this is that + * a point may not be found even + * though it is actually in one + * of the cells. This may happen, + * if the point is not in one of + * the coarse grid cells, even + * though it is in one of the + * cells on finer levels of the + * triangulation. Note that this + * of course implies that mother + * and child cells do not exactly + * overlap, a case that is + * frequent along curved + * boundaries. In this latter + * case, a different algorithm + * may be used instead that uses + * a linear search over all + * active cells, rather than + * first searchin for a coarse + * grid cell. Note, however, that + * such an algorithm has a + * significantly higher numerical + * cost than the logarithmic + * algorithm used here. + * + * Lastly, if a point lies on the + * boundary of two or more cells, + * then the algorithm may return + * with any of these cells. While + * this is in general not really + * problem, if may be a nuisance + * if the point lies at the + * boundary of cells with + * different refinement levels + * and one would rather like to + * evaluate a solution on the + * cell with more refinement. For + * this, more sophisticated + * algorithms would be necessary, + * though. + */ + template + static + typename Container::active_cell_iterator + find_active_cell_around_point (const Container &container, + const Point &p); /** * Exception */ DeclException0 (ExcTriangulationHasBeenRefined); - /** * Exception */ DeclException1 (ExcScalingFactorNotPositive, double, << "The scaling factor must be positive, but is " << arg1); + /** + * Exception + */ + template + DeclException1 (ExcPointNotFoundInCoarseGrid, + Point, + << "The point <" << arg1 + << "> could not be found inside any of the " + << "coarse grid cells."); + /** + * Exception + */ + template + DeclException1 (ExcPointNotFound, + Point, + << "The point <" << arg1 + << "> could not be found inside any of the " + << "subcells of a coarse grid cell."); }; diff --git a/deal.II/deal.II/source/grid/grid_tools.cc b/deal.II/deal.II/source/grid/grid_tools.cc index 1fa69bf0e6..08032d0252 100644 --- a/deal.II/deal.II/source/grid/grid_tools.cc +++ b/deal.II/deal.II/source/grid/grid_tools.cc @@ -17,6 +17,10 @@ #include #include #include +#include +#include +#include +#include #include @@ -186,6 +190,55 @@ GridTools::scale (const double scaling_factor, +template +typename Container::active_cell_iterator +GridTools::find_active_cell_around_point (const Container &container, + const Point &p) +{ + // first find the coarse grid cell + // that contains the point. we can + // only do this by a linear search + typename Container::cell_iterator cell = container.begin(0); + for (; cell!=container.end(0); ++cell) + if (cell->point_inside (p)) + break; + + // make sure that we found a cell + // in the coarse grid that contains + // this point. for cases where this + // might happen unexpectedly, see + // the documentation of this + // function + AssertThrow (cell != container.end(0), + ExcPointNotFoundInCoarseGrid (p)); + + // now do the logarithmic part of + // the algorithm: go from child to + // grandchild + while (cell->has_children()) + { + unsigned int c=0; + for (; c::children_per_cell; ++c) + if (cell->point_inside (p)) + break; + + // make sure we found a child + // cell + AssertThrow (c != GeometryInfo::children_per_cell, + ExcPointNotFound (p)); + + // then reset cell to the child + cell = cell->child(c); + } + + // now that we have a terminal + // cell, return it + return cell; +} + + + + #if deal_II_dimension != 1 template double @@ -199,3 +252,18 @@ void GridTools::shift (const Point &, template void GridTools::scale (const double, Triangulation &); + +template +Triangulation::active_cell_iterator +GridTools::find_active_cell_around_point (const Triangulation &, + const Point &p); + +template +DoFHandler::active_cell_iterator +GridTools::find_active_cell_around_point (const DoFHandler &, + const Point &p); + +template +MGDoFHandler::active_cell_iterator +GridTools::find_active_cell_around_point (const MGDoFHandler &, + const Point &p); diff --git a/deal.II/doc/news/c-4-0.html b/deal.II/doc/news/c-4-0.html index 7910f97102..3008cb1e20 100644 --- a/deal.II/doc/news/c-4-0.html +++ b/deal.II/doc/news/c-4-0.html @@ -295,6 +295,15 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

deal.II

    +
  1. + New: Long requested but never implemented before in the library: there + is now a function GridToold::find_active_cell_around_point that, + given a point, finds the active cell in which this point lies. +
    + (WB 2003/10/30) +

    +
  2. Changed: The SolutionTransfer::interpolate diff --git a/tests/bits/Makefile b/tests/bits/Makefile index 7dcd1780d3..736568b1d4 100644 --- a/tests/bits/Makefile +++ b/tests/bits/Makefile @@ -134,6 +134,10 @@ volume_1.exe : volume_1.g.$(OBJEXT) $(libraries) volume_2.exe : volume_2.g.$(OBJEXT) $(libraries) volume_3.exe : volume_3.g.$(OBJEXT) $(libraries) volume_4.exe : volume_4.g.$(OBJEXT) $(libraries) +find_cell_1.exe : find_cell_1.g.$(OBJEXT) $(libraries) +find_cell_2.exe : find_cell_2.g.$(OBJEXT) $(libraries) +find_cell_3.exe : find_cell_3.g.$(OBJEXT) $(libraries) + tests = anna_1 anna_2 anna_3 anna_4 anna_5 anna_6 \ geometry_info_1 point_inside_1 point_inside_2 \ @@ -171,7 +175,7 @@ tests = anna_1 anna_2 anna_3 anna_4 anna_5 anna_6 \ volume_1 volume_2 volume_3 volume_4 \ mapping_cartesian_1 \ mapping_q4_3d \ - q_points + q_points find_cell_1 find_cell_2 find_cell_3 ############################################################ diff --git a/tests/bits/find_cell_1.cc b/tests/bits/find_cell_1.cc new file mode 100644 index 0000000000..a838942352 --- /dev/null +++ b/tests/bits/find_cell_1.cc @@ -0,0 +1,68 @@ +//---------------------------- find_cell_1.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- find_cell_1.cc --------------------------- + + +// take a 2d mesh and check that we can find an arbitrary point's cell +// in it + +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +void check (Triangulation<2> &tria) +{ + Triangulation<2>::active_cell_iterator cell + = GridTools::find_active_cell_around_point (tria, + Point<2>(1./3., 1./2.)); + + deallog << cell << std::endl; + for (unsigned int v=0; v::vertices_per_cell; ++v) + deallog << "<" << cell->vertex(v) << "> "; + deallog << std::endl; +} + + +int main () +{ + std::ofstream logfile("find_cell_1.output"); + deallog.attach(logfile); + deallog.depth_console(0); + + { + Triangulation<2> coarse_grid; + GridGenerator::hyper_cube (coarse_grid); + coarse_grid.refine_global (2); + check (coarse_grid); + } + + { + Triangulation<2> coarse_grid; + GridGenerator::hyper_ball (coarse_grid); + static const HyperBallBoundary<2> boundary; + coarse_grid.set_boundary (0, boundary); + coarse_grid.refine_global (2); + check (coarse_grid); + } +} + + + diff --git a/tests/bits/find_cell_2.cc b/tests/bits/find_cell_2.cc new file mode 100644 index 0000000000..5efdc8857d --- /dev/null +++ b/tests/bits/find_cell_2.cc @@ -0,0 +1,69 @@ +//---------------------------- find_cell_2.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- find_cell_2.cc --------------------------- + + +// same as find_cell_2_1, but in 3d + +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +void check (Triangulation<3> &tria) +{ + Triangulation<3>::active_cell_iterator cell + = GridTools::find_active_cell_around_point (tria, + Point<3>(1./3., + 1./2., + 2./3.)); + + deallog << cell << std::endl; + for (unsigned int v=0; v::vertices_per_cell; ++v) + deallog << "<" << cell->vertex(v) << "> "; + deallog << std::endl; +} + + +int main () +{ + std::ofstream logfile("find_cell_2.output"); + deallog.attach(logfile); + deallog.depth_console(0); + + { + Triangulation<3> coarse_grid; + GridGenerator::hyper_cube (coarse_grid); + coarse_grid.refine_global (2); + check (coarse_grid); + } + + { + Triangulation<3> coarse_grid; + GridGenerator::hyper_ball (coarse_grid); + static const HyperBallBoundary<3> boundary; + coarse_grid.set_boundary (0, boundary); + coarse_grid.refine_global (2); + check (coarse_grid); + } +} + + + diff --git a/tests/bits/find_cell_3.cc b/tests/bits/find_cell_3.cc new file mode 100644 index 0000000000..e88fb72985 --- /dev/null +++ b/tests/bits/find_cell_3.cc @@ -0,0 +1,75 @@ +//---------------------------- find_cell_3.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- find_cell_3.cc --------------------------- + + +// like find_cell_2, but with the strange meshes from the mesh_3d_* tests + +#include "mesh_3d.h" + +#include +#include +#include +#include +#include +#include +#include + +#include + + +void check (Triangulation<3> &tria) +{ + Triangulation<3>::active_cell_iterator cell + = GridTools::find_active_cell_around_point (tria, + Point<3>(1./3., + 1./2., + 2./3.)); + + deallog << cell << std::endl; + for (unsigned int v=0; v::vertices_per_cell; ++v) + deallog << "<" << cell->vertex(v) << "> "; + deallog << std::endl; +} + + +int main () +{ + std::ofstream logfile("find_cell_3.output"); + deallog.attach(logfile); + deallog.depth_console(0); + + { + Triangulation<3> coarse_grid; + create_two_cubes (coarse_grid); + coarse_grid.refine_global (1); + check (coarse_grid); + } + + { + Triangulation<3> coarse_grid; + create_L_shape (coarse_grid); + coarse_grid.refine_global (1); + check (coarse_grid); + } + + { + Triangulation<3> coarse_grid; + GridGenerator::hyper_ball (coarse_grid); + coarse_grid.refine_global (1); + check (coarse_grid); + } + +} + + +