From: wolf Date: Thu, 16 Aug 2001 20:53:10 +0000 (+0000) Subject: Add a class that provides algorithms working on triangulations. New test case. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c9f08a153dbab9e07b6874dfadac1baefeb98f38;p=dealii-svn.git Add a class that provides algorithms working on triangulations. New test case. git-svn-id: https://svn.dealii.org/trunk@4896 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 new file mode 100644 index 0000000000..9d5912d9e3 --- /dev/null +++ b/deal.II/deal.II/include/grid/grid_tools.h @@ -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 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 + static + double diameter (const Triangulation &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 ---------------------------*/ diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index ddfe2cd1f6..bd3c02d331 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -2841,42 +2841,81 @@ class Triangulation : public TriaDimensionInfo, 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 > & + 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 & + 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 index 0000000000..c9e04e827c --- /dev/null +++ b/deal.II/deal.II/source/grid/grid_tools.cc @@ -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 +#include +#include +#include + +#if deal_II_dimension != 1 + +template +double +GridTools::diameter (const Triangulation &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 > &vertices = tria.get_vertices (); + std::vector boundary_vertices (vertices.size(), false); + + typename Triangulation::active_cell_iterator + cell = tria.begin_active(); + const typename Triangulation::active_cell_iterator + endc = tria.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary ()) + for (unsigned int i=0; i::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::const_iterator pi = boundary_vertices.begin(); + const unsigned int N = boundary_vertices.size(); + for (unsigned int i=0; i::const_iterator pj = pi+1; + for (unsigned int j=i+1; j 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 &); +#endif diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 6a07e50324..f1d2983dad 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -3680,6 +3680,7 @@ unsigned int Triangulation::n_levels () const }; + template unsigned int Triangulation::n_vertices () const @@ -3688,6 +3689,16 @@ Triangulation::n_vertices () const }; + +template +const std::vector > & +Triangulation::get_vertices () const +{ + return vertices; +}; + + + template unsigned int Triangulation::n_used_vertices () const @@ -3697,10 +3708,20 @@ Triangulation::n_used_vertices () const }; + +template +const std::vector & +Triangulation::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; }; diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 2b4b10178f..eee8321b95 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -545,6 +545,15 @@ documentation, etc.

deal.II

    +

  1. + New: There is now a class GridTools + which provides algorithms working on triangulations. At + present, it offers a function computing the diameter of a + triangulation. +
    + (WB 2001/08/16) +

    +

  2. Changed: The MatrixCreator and MatrixTools class have lost their template diff --git a/tests/deal.II/Makefile.in b/tests/deal.II/Makefile.in index 8e0fb1a079..81b8d2ff74 100644 --- a/tests/deal.II/Makefile.in +++ b/tests/deal.II/Makefile.in @@ -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 ############################################################ diff --git a/tests/deal.II/grid_test.cc b/tests/deal.II/grid_test.cc index b0fceeb12b..5934586395 100644 --- a/tests/deal.II/grid_test.cc +++ b/tests/deal.II/grid_test.cc @@ -93,23 +93,23 @@ CurvedLine::get_new_point_on_line (const typename Triangulation::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 +void test () +{ + // test 1: hypercube + if (true) + { + Triangulation 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 tria; + GridGenerator::hyper_ball(tria, Point(), 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 index 0000000000..a96beeb8cf --- /dev/null +++ b/tests/deal.II/grid_tools.checked @@ -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