]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
GridTools::find_active_cell_around_point
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 30 Oct 2003 14:23:51 +0000 (14:23 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 30 Oct 2003 14:23:51 +0000 (14:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@8174 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_tools.h
deal.II/deal.II/source/grid/grid_tools.cc
deal.II/doc/news/c-4-0.html
tests/bits/Makefile
tests/bits/find_cell_1.cc [new file with mode: 0644]
tests/bits/find_cell_2.cc [new file with mode: 0644]
tests/bits/find_cell_3.cc [new file with mode: 0644]

index be9414dfe0acc372c435ea25c08be69e2085c7be..6623d7ae3fb038b9a9cd1d01840225335a01c5b2 100644 (file)
 
 
 /**
- * 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<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.");
 };
 
 
index 1fa69bf0e6bb3ad5ce79c0aa513a6c13608f67b4..08032d0252df4bbb04d861035e62be17074a163b 100644 (file)
 #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>
 
@@ -186,6 +190,55 @@ GridTools::scale (const double        scaling_factor,
 
 
 
+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
@@ -199,3 +252,18 @@ void GridTools::shift<deal_II_dimension> (const Point<deal_II_dimension> &,
 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);
index 7910f9710283e7ea63b2ab422aa966982a96daa0..3008cb1e20677da8a02cfc847829d9e05f3268cc 100644 (file)
@@ -295,6 +295,15 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <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>      
index 7dcd1780d39e6e0982b036145d2f980c72643a25..736568b1d4e7abd333cfc57ec44beb5bf4584d7c 100644 (file)
@@ -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 (file)
index 0000000..a838942
--- /dev/null
@@ -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 <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);
+  }
+}
+
+  
+  
diff --git a/tests/bits/find_cell_2.cc b/tests/bits/find_cell_2.cc
new file mode 100644 (file)
index 0000000..5efdc88
--- /dev/null
@@ -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 <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);
+  }
+}
+
+  
+  
diff --git a/tests/bits/find_cell_3.cc b/tests/bits/find_cell_3.cc
new file mode 100644 (file)
index 0000000..e88fb72
--- /dev/null
@@ -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 <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);
+  }
+  
+}
+
+  
+  

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.