]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a class that provides algorithms working on triangulations. New test case.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 16 Aug 2001 20:53:10 +0000 (20:53 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 16 Aug 2001 20:53:10 +0000 (20:53 +0000)
git-svn-id: https://svn.dealii.org/trunk@4896 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_tools.h [new file with mode: 0644]
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/source/grid/grid_tools.cc [new file with mode: 0644]
deal.II/deal.II/source/grid/tria.cc
deal.II/doc/news/2001/c-3-1.html
tests/deal.II/Makefile.in
tests/deal.II/grid_test.cc
tests/deal.II/grid_tools.cc [new file with mode: 0644]
tests/deal.II/grid_tools.checked [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/grid/grid_tools.h b/deal.II/deal.II/include/grid/grid_tools.h
new file mode 100644 (file)
index 0000000..9d5912d
--- /dev/null
@@ -0,0 +1,58 @@
+//----------------------------  grid_tools.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 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.
+//
+//----------------------------  grid_tools.h  ---------------------------
+#ifndef __deal2__grid_tools_H
+#define __deal2__grid_tools_H
+
+
+template <int dim> class Triangulation;
+
+
+
+/**
+ * This class is a collection of algorithms working on
+ * triangulations. See the descriptions of the individual functions
+ * for more information.
+ *
+ * @author Wolfgang Bangerth, 2001
+ */
+class GridTools
+{
+  public:
+                                    /**
+                                     * Return the diameter of a
+                                     * triangulation. The diameter is
+                                     * computed using only the
+                                     * vertices, i.e. if the diameter
+                                     * should be larger than the
+                                     * maximal distance between
+                                     * boundary vertices due to a
+                                     * higher order mapping, then
+                                     * this function will not catch
+                                     * this.
+                                     */
+    template <int dim>
+    static
+    double diameter (const Triangulation<dim> &tria);
+
+                                    /**
+                                     * Same function, but for 1d.
+                                     */
+    static
+    double diameter (const Triangulation<1> &tria);
+};
+
+
+/*----------------------------   grid_tools.h     ---------------------------*/
+/* end of #ifndef __deal2__grid_tools_H */
+#endif
+/*----------------------------   grid_tools.h     ---------------------------*/
index ddfe2cd1f6b44617282c21354cb645a08e9e3ca8..bd3c02d33152f2a3d532542bbcaf41449d9a0987 100644 (file)
@@ -2841,42 +2841,81 @@ class Triangulation : public TriaDimensionInfo<dim>,
     unsigned int n_levels () const;
 
                                     /**
-                                     * Return the total number of vertices.
-                                     * Some of them may not be used, which
-                                     * usually happens upon coarsening of
-                                     * a triangulation when some vertices are
-                                     * discarded, but we do not want to
-                                     * renumber the remaining one, leading to
-                                     * holes in the numbers of used vertices.
-                                     * You can get the number of used vertices
-                                     * using @p{n_used_vertices} function.
+                                     * Return the total number of
+                                     * vertices.  Some of them may
+                                     * not be used, which usually
+                                     * happens upon coarsening of a
+                                     * triangulation when some
+                                     * vertices are discarded, but we
+                                     * do not want to renumber the
+                                     * remaining one, leading to
+                                     * holes in the numbers of used
+                                     * vertices.  You can get the
+                                     * number of used vertices using
+                                     * @p{n_used_vertices} function.
                                      */
     unsigned int n_vertices () const;
+
+                                    /**
+                                     * Return a constant reference to
+                                     * all the vertices used in this
+                                     * triangulation. Note that not
+                                     * necessarily all vertices in
+                                     * this array are actually used;
+                                     * for example, if you coarsen a
+                                     * mesh, then some vertices are
+                                     * deleted, but their positions
+                                     * in this array are unchanged as
+                                     * the indices of vertices are
+                                     * only allocated once. You can
+                                     * find out about which vertices
+                                     * are actually used by the
+                                     * function
+                                     * @ref{get_used_vertices}.
+                                     */
+    const std::vector<Point<dim> > &
+    get_vertices () const;
     
                                     /**
-                                     * Return the number of vertices that are
-                                     * presently in use, i.e. belong to at least
-                                     * one used element.
+                                     * Return the number of vertices
+                                     * that are presently in use,
+                                     * i.e. belong to at least one
+                                     * used element.
                                      */
     unsigned int n_used_vertices () const;
+
+                                    /**
+                                     * Return a constant reference to
+                                     * the array of @p{bool}s
+                                     * indicating whether an entry in
+                                     * the vertex array is used or
+                                     * not.
+                                     */
+    const std::vector<bool> &
+    get_used_vertices () const;
     
                                     /**
-                                     * Return the maximum number of cells
-                                     * meeting at a common vertex. Since this
-                                     * number is an invariant under refinement,
-                                     * only the cells on
-                                     * the coarsest level are considered. The
-                                     * operation is thus reasonably fast. The
-                                     * invariance is only true for sufficiently
-                                     * many cells in the coarsest triangulation
-                                     * (e.g. for a single cell one would be
-                                     * returned),
-                                     * so a minimum of four is returned in
-                                     * two dimensions, 8 in three dimensions,
-                                     * etc, which is how many cells meet if the
+                                     * Return the maximum number of
+                                     * cells meeting at a common
+                                     * vertex. Since this number is
+                                     * an invariant under refinement,
+                                     * only the cells on the coarsest
+                                     * level are considered. The
+                                     * operation is thus reasonably
+                                     * fast. The invariance is only
+                                     * true for sufficiently many
+                                     * cells in the coarsest
+                                     * triangulation (e.g. for a
+                                     * single cell one would be
+                                     * returned), so a minimum of
+                                     * four is returned in two
+                                     * dimensions, 8 in three
+                                     * dimensions, etc, which is how
+                                     * many cells meet if the
                                      * triangulation is refined.
                                      *
-                                     * In one space dimension, two is returned.
+                                     * In one space dimension, two is
+                                     * returned.
                                      */
     unsigned int max_adjacent_cells () const;
                                     /*@}*/
diff --git a/deal.II/deal.II/source/grid/grid_tools.cc b/deal.II/deal.II/source/grid/grid_tools.cc
new file mode 100644 (file)
index 0000000..c9e04e8
--- /dev/null
@@ -0,0 +1,99 @@
+//----------------------------  grid_tools.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 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.
+//
+//----------------------------  grid_tools.cc  ---------------------------
+
+
+
+#include <grid/grid_tools.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+
+#if deal_II_dimension != 1
+
+template <int dim>
+double
+GridTools::diameter (const Triangulation<dim> &tria)
+{
+                                  // the algorithm used simply
+                                  // traverses all cells and picks
+                                  // out the boundary vertices. it
+                                  // may or may not be faster to
+                                  // simply get all vectors, don't
+                                  // mark boundary vertices, and
+                                  // compute the distances thereof,
+                                  // but at least as the mesh is
+                                  // refined, it seems better to
+                                  // first mark boundary nodes, as
+                                  // marking is O(N) in the number of
+                                  // cells/vertices, while computing
+                                  // the maximal distance is O(N*N)
+  const std::vector<Point<dim> > &vertices = tria.get_vertices ();
+  std::vector<bool> boundary_vertices (vertices.size(), false);
+
+  typename Triangulation<dim>::active_cell_iterator
+    cell = tria.begin_active();
+  const typename Triangulation<dim>::active_cell_iterator
+    endc = tria.end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+      if (cell->face(face)->at_boundary ())
+       for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
+         boundary_vertices[cell->face(face)->vertex_index(i)] = true;
+
+                                  // now traverse the list of
+                                  // boundary vertices and check
+                                  // distances. since distances are
+                                  // symmetric, we only have to check
+                                  // one half
+  double max_distance_sqr = 0;
+  std::vector<bool>::const_iterator pi = boundary_vertices.begin();
+  const unsigned int N = boundary_vertices.size();
+  for (unsigned int i=0; i<N; ++i, ++pi)
+    {
+      std::vector<bool>::const_iterator pj = pi+1;
+      for (unsigned int j=i+1; j<N; ++j, ++pj)
+       if ((*pi==true) && (*pj==true) &&
+           ((vertices[i]-vertices[j]).square() > max_distance_sqr))
+         max_distance_sqr = (vertices[i]-vertices[j]).square();
+    };
+
+  return sqrt(max_distance_sqr);
+};
+
+
+#else
+
+double
+GridTools::diameter (const Triangulation<1> &tria)
+{
+                                  // for 1d, simply check the
+                                  // vertices of the left- and
+                                  // rightmost coarse grid cell
+  Triangulation<1>::cell_iterator leftmost  = tria.begin(0);
+  Triangulation<1>::cell_iterator rightmost = tria.begin(0);
+
+  while (!leftmost->at_boundary(0))  leftmost  = leftmost->neighbor(0);
+  while (!rightmost->at_boundary(1)) rightmost = rightmost->neighbor(1);
+
+  return sqrt((leftmost->vertex(0) - rightmost->vertex(1)).square());
+};
+
+#endif
+
+
+
+#if deal_II_dimension != 1
+template
+double
+GridTools::diameter (const Triangulation<deal_II_dimension> &);
+#endif
index 6a07e50324eb2915672ee50d82a2bf8da15a6cf7..f1d2983dad0f81e45a45fad5b6a74aa0ec958fae 100644 (file)
@@ -3680,6 +3680,7 @@ unsigned int Triangulation<dim>::n_levels () const
 };
 
 
+
 template <int dim>
 unsigned int
 Triangulation<dim>::n_vertices () const 
@@ -3688,6 +3689,16 @@ Triangulation<dim>::n_vertices () const
 };
 
 
+
+template <int dim>
+const std::vector<Point<dim> > &
+Triangulation<dim>::get_vertices () const 
+{
+  return vertices;
+};
+
+
+
 template <int dim>
 unsigned int
 Triangulation<dim>::n_used_vertices () const 
@@ -3697,10 +3708,20 @@ Triangulation<dim>::n_used_vertices () const
 };
 
 
+
+template <int dim>
+const std::vector<bool> &
+Triangulation<dim>::get_used_vertices () const 
+{
+  return vertices_used;
+};
+
+
 #if deal_II_dimension == 1
 
 template <>
-unsigned int Triangulation<1>::max_adjacent_cells () const {
+unsigned int Triangulation<1>::max_adjacent_cells () const
+{
   return 2;
 };
 
index 2b4b10178f3121c9f1c2428ac71a24371343ab48..eee8321b95547fde4809082c056467519980b9de 100644 (file)
@@ -545,6 +545,15 @@ documentation, etc</a>.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       New: There is now a class <code class="class">GridTools</code>
+       which provides algorithms working on triangulations. At
+       present, it offers a function computing the diameter of a
+       triangulation.
+       <br>
+       (WB 2001/08/16)
+       </p>
+
   <li> <p>
        Changed: The <code class="class">MatrixCreator</code> and <code
        class="class">MatrixTools</code> class have lost their template
index 8e0fb1a07952d0952d65ae8478e491bffa0c0386..81b8d2ff7496e405fd0af24e5a2b6d86d33ca2fb 100644 (file)
@@ -43,12 +43,13 @@ support_point_map.exe        : support_point_map.go        $(lib-1d) $(lib-2d) $
 filtered_matrix.exe          : filtered_matrix.go          $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 boundaries.exe               : boundaries.go               $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 sparsity_pattern.exe         : sparsity_pattern.go         $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
+grid_tools.exe               : grid_tools.go               $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 
 tests = grid_test grid_transform dof_test data_out derivatives gradients constraints mg \
         mglocal block_matrices second_derivatives derivative_approximation \
         matrices error_estimator intergrid_constraints intergrid_map \
         wave-test-3 dof_renumbering support_point_map filtered_matrix \
-       boundaries sparsity_pattern
+       boundaries sparsity_pattern grid_tools
 
 ############################################################
 
index b0fceeb12b7f6bd627ed608d8d4fd05c60a96c39..5934586395bf17c16f6c0543fa89b670721baeba 100644 (file)
@@ -93,23 +93,23 @@ CurvedLine<dim>::get_new_point_on_line (const typename Triangulation<dim>::line_
                                   // vertices of the line is like that
   if (dim>=3)
     if (((middle(2) == 0) || (middle(2) == 1))
-                                      // find out, if the line is in the
-                                      // interior of the top or bottom face
-                                      // of the domain, or at the edge.
-                                      // lines at the edge need to undergo
-                                      // the usual treatment, while for
-                                      // interior lines taking the midpoint
-                                      // is sufficient
-                                      //
-                                      // note: the trick with the boundary
-                                      // id was invented after the above was
-                                      // written, so we are not very strict
-                                      // here with using these flags
+                                        // find out, if the line is in the
+                                        // interior of the top or bottom face
+                                        // of the domain, or at the edge.
+                                        // lines at the edge need to undergo
+                                        // the usual treatment, while for
+                                        // interior lines taking the midpoint
+                                        // is sufficient
+                                        //
+                                        // note: the trick with the boundary
+                                        // id was invented after the above was
+                                        // written, so we are not very strict
+                                        // here with using these flags
        && (line->boundary_indicator() == 1))
       return middle;
 
 
-double x=middle(0),
+  double x=middle(0),
         y=middle(1);
   
   if (y<x)
@@ -183,7 +183,7 @@ void test (const int test_case)
     };
 
 
-switch (test_case) 
+  switch (test_case) 
     {
       case 1: 
       {
@@ -256,10 +256,10 @@ switch (test_case)
     };
 
 
-GridOut().write_ucd (tria, logfile);
+  GridOut().write_ucd (tria, logfile);
     
   deallog << "     Total number of cells        = " << tria.n_cells() << std::endl
-       << "     Total number of active cells = " << tria.n_active_cells() << std::endl;
+         << "     Total number of active cells = " << tria.n_active_cells() << std::endl;
 
   deallog.pop();
 };
diff --git a/tests/deal.II/grid_tools.cc b/tests/deal.II/grid_tools.cc
new file mode 100644 (file)
index 0000000..21f7b86
--- /dev/null
@@ -0,0 +1,79 @@
+//----------------------------  grid_test.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001 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.
+//
+//----------------------------  grid_test.cc  ---------------------------
+
+
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_tools.h>
+
+#include <fstream>
+
+
+std::ofstream logfile("grid_tools.output");
+
+
+
+template <int dim>
+void test ()
+{
+                                  // test 1: hypercube
+  if (true)
+    {
+      Triangulation<dim> tria;
+      GridGenerator::hyper_cube(tria);
+
+      for (unsigned int i=0; i<2; ++i)
+       {
+         tria.refine_global(2);
+         deallog << dim << "d, "
+                 << "hypercube diameter, "
+                 << i*2
+                 << " refinements: "
+                 << GridTools::diameter (tria)
+                 << std::endl;
+       };
+    };
+
+                                  // test 2: hyperball
+  if (dim == 2)
+    {
+      Triangulation<dim> tria;
+      GridGenerator::hyper_ball(tria, Point<dim>(), 1);
+
+      for (unsigned int i=0; i<2; ++i)
+       {
+         tria.refine_global(2);
+         deallog << dim << "d, "
+                 << "hyperball diameter, "
+                 << i*2
+                 << " refinements: "
+                 << GridTools::diameter (tria)
+                 << std::endl;
+       };
+    };
+};
+
+
+int main ()
+{
+  logfile.precision(4);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<1> ();
+  test<2> ();
+  test<3> ();
+  
+  return 0;
+};
diff --git a/tests/deal.II/grid_tools.checked b/tests/deal.II/grid_tools.checked
new file mode 100644 (file)
index 0000000..a96beeb
--- /dev/null
@@ -0,0 +1,9 @@
+
+DEAL::1d, hypercube diameter, 0 refinements: 1.000
+DEAL::1d, hypercube diameter, 2 refinements: 1.000
+DEAL::2d, hypercube diameter, 0 refinements: 1.414
+DEAL::2d, hypercube diameter, 2 refinements: 1.414
+DEAL::2d, hyperball diameter, 0 refinements: 2.000
+DEAL::2d, hyperball diameter, 2 refinements: 2.000
+DEAL::3d, hypercube diameter, 0 refinements: 1.732
+DEAL::3d, hypercube diameter, 2 refinements: 1.732

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.