/**
- * 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
{
static
void scale (const double scaling_factor,
Triangulation<dim> &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 <int dim, typename Container>
+ static
+ typename Container::active_cell_iterator
+ find_active_cell_around_point (const Container &container,
+ const Point<dim> &p);
/**
* Exception
*/
DeclException0 (ExcTriangulationHasBeenRefined);
-
/**
* Exception
*/
DeclException1 (ExcScalingFactorNotPositive,
double,
<< "The scaling factor must be positive, but is " << arg1);
+ /**
+ * Exception
+ */
+ template <int N>
+ DeclException1 (ExcPointNotFoundInCoarseGrid,
+ Point<N>,
+ << "The point <" << arg1
+ << "> could not be found inside any of the "
+ << "coarse grid cells.");
+ /**
+ * Exception
+ */
+ template <int N>
+ DeclException1 (ExcPointNotFound,
+ Point<N>,
+ << "The point <" << arg1
+ << "> could not be found inside any of the "
+ << "subcells of a coarse grid cell.");
};
#include <grid/tria.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
#include <cmath>
+template <int dim, typename Container>
+typename Container::active_cell_iterator
+GridTools::find_active_cell_around_point (const Container &container,
+ const Point<dim> &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<dim> (p));
+
+ // now do the logarithmic part of
+ // the algorithm: go from child to
+ // grandchild
+ while (cell->has_children())
+ {
+ unsigned int c=0;
+ for (; c<GeometryInfo<dim>::children_per_cell; ++c)
+ if (cell->point_inside (p))
+ break;
+
+ // make sure we found a child
+ // cell
+ AssertThrow (c != GeometryInfo<dim>::children_per_cell,
+ ExcPointNotFound<dim> (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
template
void GridTools::scale<deal_II_dimension> (const double,
Triangulation<deal_II_dimension> &);
+
+template
+Triangulation<deal_II_dimension>::active_cell_iterator
+GridTools::find_active_cell_around_point (const Triangulation<deal_II_dimension> &,
+ const Point<deal_II_dimension> &p);
+
+template
+DoFHandler<deal_II_dimension>::active_cell_iterator
+GridTools::find_active_cell_around_point (const DoFHandler<deal_II_dimension> &,
+ const Point<deal_II_dimension> &p);
+
+template
+MGDoFHandler<deal_II_dimension>::active_cell_iterator
+GridTools::find_active_cell_around_point (const MGDoFHandler<deal_II_dimension> &,
+ const Point<deal_II_dimension> &p);
<h3>deal.II</h3>
<ol>
+ <li> <p>
+ New: Long requested but never implemented before in the library: there
+ is now a function <code
+ <class="member">GridToold::find_active_cell_around_point</code> that,
+ given a point, finds the active cell in which this point lies.
+ <br>
+ (WB 2003/10/30)
+ </p>
+
<li> <p>
Changed: The <code
<class="member">SolutionTransfer::interpolate</code>
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 \
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
############################################################
--- /dev/null
+//---------------------------- 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 <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_tools.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <fstream>
+
+
+
+
+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<GeometryInfo<2>::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);
+ }
+}
+
+
+
--- /dev/null
+//---------------------------- 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 <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_tools.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <fstream>
+
+
+
+
+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<GeometryInfo<3>::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);
+ }
+}
+
+
+
--- /dev/null
+//---------------------------- 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 <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_tools.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <fstream>
+
+
+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<GeometryInfo<3>::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);
+ }
+
+}
+
+
+