--- /dev/null
+//---------------------------- 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 ---------------------------*/
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;
/*@}*/
--- /dev/null
+//---------------------------- 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
};
+
template <int dim>
unsigned int
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
};
+
+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;
};
<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
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
############################################################
// 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)
};
-switch (test_case)
+ switch (test_case)
{
case 1:
{
};
-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();
};
--- /dev/null
+//---------------------------- 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;
+};
--- /dev/null
+
+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
projects / dealii-svn.git / commitdiff