From c26a9ec145f2c3d2c8207051537677db6b400d7e Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 15 Feb 2011 18:31:55 +0000 Subject: [PATCH] Fix a problem with colorize_hyper_shell which wasn't written in a particularly defensive way and produced garbage with the new 96-cell hyper shell. git-svn-id: https://svn.dealii.org/trunk@23370 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/grid/grid_generator.cc | 64 +++++++-- tests/deal.II/grid_hyper_shell_05.cc | 82 +++++++++++ tests/deal.II/grid_hyper_shell_05/cmp/generic | 136 ++++++++++++++++++ 3 files changed, 270 insertions(+), 12 deletions(-) create mode 100644 tests/deal.II/grid_hyper_shell_05.cc create mode 100644 tests/deal.II/grid_hyper_shell_05/cmp/generic diff --git a/deal.II/source/grid/grid_generator.cc b/deal.II/source/grid/grid_generator.cc index 9666862117..dbafe0257d 100644 --- a/deal.II/source/grid/grid_generator.cc +++ b/deal.II/source/grid/grid_generator.cc @@ -1442,6 +1442,7 @@ GridGenerator::colorize_hyper_shell ( for (Triangulation<2>::cell_iterator cell = tria.begin(); cell != tria.end(); ++cell) { + Assert (cell->face(2)->at_boundary(), ExcInternalError()); cell->face(2)->set_boundary_indicator(1); } } @@ -2236,29 +2237,68 @@ colorize_hyper_shell (Triangulation<3>& tria, const double, const double) { - // Inspite of receiving geometrical - // data, we do this only based on - // topology. - - // For the mesh based on cube, - // this is highly irregular + // the following uses a good amount + // of knowledge about the + // orientation of cells. this is + // probably not good style... if (tria.n_cells() == 6) { Triangulation<3>::cell_iterator cell = tria.begin(); + cell->face(4)->set_boundary_indicator(1); + Assert (cell->face(4)->at_boundary(), ExcInternalError()); + (++cell)->face(2)->set_boundary_indicator(1); + Assert (cell->face(2)->at_boundary(), ExcInternalError()); + (++cell)->face(2)->set_boundary_indicator(1); + Assert (cell->face(2)->at_boundary(), ExcInternalError()); + (++cell)->face(0)->set_boundary_indicator(1); + Assert (cell->face(0)->at_boundary(), ExcInternalError()); + (++cell)->face(2)->set_boundary_indicator(1); + Assert (cell->face(2)->at_boundary(), ExcInternalError()); + (++cell)->face(0)->set_boundary_indicator(1); + Assert (cell->face(0)->at_boundary(), ExcInternalError()); } - else - // For higher polyhedra, this is regular. + else if (tria.n_cells() == 12) { + // again use some internal + // knowledge for (Triangulation<3>::cell_iterator cell = tria.begin(); cell != tria.end(); ++cell) - cell->face(5)->set_boundary_indicator(1); + { + Assert (cell->face(5)->at_boundary(), ExcInternalError()); + cell->face(5)->set_boundary_indicator(1); + } } + else if (tria.n_cells() == 96) + { + // the 96-cell hypershell is + // based on a once refined + // 12-cell mesh. consequently, + // since the outer faces all + // are face_no==5 above, so + // they are here (unless they + // are in the interior). Use + // this to assign boundary + // indicators, but also make + // sure that we encounter + // exactly 48 such faces + unsigned int count = 0; + for (Triangulation<3>::cell_iterator cell = tria.begin(); + cell != tria.end(); ++cell) + if (cell->face(5)->at_boundary()) + { + cell->face(5)->set_boundary_indicator(1); + ++count; + } + Assert (count == 48, ExcInternalError()); + } + else + Assert (false, ExcNotImplemented()); } @@ -2422,11 +2462,11 @@ GridGenerator::hyper_shell (Triangulation<3>& tria, const double r = inner_radius / outer_radius; const double pi = numbers::PI; const double rho = 2*r/(1+r); - + // then this is the distance of the // interior nodes from the center: const double middle_radius = rho * outer_radius; - + // mark vertices we've already moved or // that we want to ignore: we don't // want to move vertices at the inner @@ -2457,7 +2497,7 @@ GridGenerator::hyper_shell (Triangulation<3>& tria, const Point<3> old_distance = cell->vertex(v) - p; const double old_radius = cell->vertex(v).distance(p); cell->vertex(v) = p + old_distance * (middle_radius / old_radius); - + vertex_already_treated[cell->vertex_index(v)] = true; } diff --git a/tests/deal.II/grid_hyper_shell_05.cc b/tests/deal.II/grid_hyper_shell_05.cc new file mode 100644 index 0000000000..a287aeefb8 --- /dev/null +++ b/tests/deal.II/grid_hyper_shell_05.cc @@ -0,0 +1,82 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2007, 2008, 2009, 2010, 2011 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. +// +//---------------------------------------------------------------------- + +// generate and refine a hyper shell in 3d with 96 cells. +// +// this mesh ran into an interesting problem (somewhat hard to find +// :-( ) in that it violates the condition +// cell->face(f)->at_boundary() == cell->at_boundary(f) +// when looping over all cells. The underlying cause was that the +// colorize option for the mesh was causing trouble. Verify that this +// problem has been fixed + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +std::ofstream logfile("grid_hyper_shell_05/output"); + + +template +void check (const unsigned int n) +{ + deallog << "n=" << n << std::endl; + + Point center; + Triangulation tria (Triangulation::none); + GridGenerator::hyper_shell (tria, center, 0.5, 1, n, true); + + // this is the test that failed + // before + for (typename Triangulation::active_cell_iterator + cell = tria.begin_active(); + cell != tria.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + Assert (cell->face(f)->at_boundary() == cell->at_boundary(f), + ExcInternalError()); + + // also output something slightly + // more useful + for (typename Triangulation::active_cell_iterator + cell = tria.begin_active(); + cell != tria.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->at_boundary(f)) + deallog << cell->face(f) << ' ' << (int)cell->face(f)->boundary_indicator() + << ' ' << cell->face(f)->center().norm() << std::endl; +} + + +int main() +{ + deallog << std::setprecision(3); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + check<3> (6); + check<3> (12); + check<3> (96); +} diff --git a/tests/deal.II/grid_hyper_shell_05/cmp/generic b/tests/deal.II/grid_hyper_shell_05/cmp/generic new file mode 100644 index 0000000000..c7e57732bb --- /dev/null +++ b/tests/deal.II/grid_hyper_shell_05/cmp/generic @@ -0,0 +1,136 @@ + +DEAL::n=6 +DEAL::10 1 0.577 +DEAL::1 0 0.289 +DEAL::13 1 0.577 +DEAL::3 0 0.289 +DEAL::20 1 0.577 +DEAL::5 0 0.289 +DEAL::11 1 0.577 +DEAL::2 0 0.289 +DEAL::9 1 0.577 +DEAL::0 0 0.289 +DEAL::17 1 0.577 +DEAL::4 0 0.289 +DEAL::n=12 +DEAL::24 0 0.381 +DEAL::44 1 0.762 +DEAL::11 0 0.381 +DEAL::40 1 0.762 +DEAL::25 0 0.381 +DEAL::45 1 0.762 +DEAL::3 0 0.381 +DEAL::37 1 0.762 +DEAL::20 0 0.381 +DEAL::42 1 0.762 +DEAL::16 0 0.381 +DEAL::41 1 0.762 +DEAL::21 0 0.381 +DEAL::43 1 0.762 +DEAL::1 0 0.381 +DEAL::36 1 0.762 +DEAL::33 0 0.381 +DEAL::47 1 0.762 +DEAL::8 0 0.381 +DEAL::39 1 0.762 +DEAL::31 0 0.381 +DEAL::46 1 0.762 +DEAL::4 0 0.381 +DEAL::38 1 0.762 +DEAL::n=96 +DEAL::24 0 0.467 +DEAL::125 0 0.467 +DEAL::120 0 0.467 +DEAL::252 0 0.467 +DEAL::44 1 0.935 +DEAL::174 1 0.935 +DEAL::172 1 0.935 +DEAL::296 1 0.935 +DEAL::11 0 0.467 +DEAL::86 0 0.467 +DEAL::82 0 0.467 +DEAL::215 0 0.467 +DEAL::40 1 0.935 +DEAL::165 1 0.935 +DEAL::164 1 0.935 +DEAL::292 1 0.935 +DEAL::25 0 0.467 +DEAL::129 0 0.467 +DEAL::123 0 0.467 +DEAL::255 0 0.467 +DEAL::45 1 0.935 +DEAL::175 1 0.935 +DEAL::173 1 0.935 +DEAL::297 1 0.935 +DEAL::3 0 0.467 +DEAL::54 0 0.467 +DEAL::55 0 0.467 +DEAL::190 0 0.467 +DEAL::37 1 0.935 +DEAL::157 1 0.935 +DEAL::158 1 0.935 +DEAL::289 1 0.935 +DEAL::20 0 0.467 +DEAL::112 0 0.467 +DEAL::107 0 0.467 +DEAL::238 0 0.467 +DEAL::42 1 0.935 +DEAL::171 1 0.935 +DEAL::169 1 0.935 +DEAL::294 1 0.935 +DEAL::16 0 0.467 +DEAL::95 0 0.467 +DEAL::96 0 0.467 +DEAL::229 0 0.467 +DEAL::41 1 0.935 +DEAL::166 1 0.935 +DEAL::167 1 0.935 +DEAL::293 1 0.935 +DEAL::21 0 0.467 +DEAL::106 0 0.467 +DEAL::110 0 0.467 +DEAL::242 0 0.467 +DEAL::43 1 0.935 +DEAL::168 1 0.935 +DEAL::170 1 0.935 +DEAL::295 1 0.935 +DEAL::1 0 0.467 +DEAL::62 0 0.467 +DEAL::50 0 0.467 +DEAL::181 0 0.467 +DEAL::36 1 0.935 +DEAL::161 1 0.935 +DEAL::156 1 0.935 +DEAL::288 1 0.935 +DEAL::33 0 0.467 +DEAL::145 0 0.467 +DEAL::146 0 0.467 +DEAL::282 0 0.467 +DEAL::47 1 0.935 +DEAL::178 1 0.935 +DEAL::179 1 0.935 +DEAL::299 1 0.935 +DEAL::8 0 0.467 +DEAL::73 0 0.467 +DEAL::74 0 0.467 +DEAL::209 0 0.467 +DEAL::39 1 0.935 +DEAL::162 1 0.935 +DEAL::163 1 0.935 +DEAL::291 1 0.935 +DEAL::31 0 0.467 +DEAL::140 0 0.467 +DEAL::141 0 0.467 +DEAL::274 0 0.467 +DEAL::46 1 0.935 +DEAL::176 1 0.935 +DEAL::177 1 0.935 +DEAL::298 1 0.935 +DEAL::4 0 0.467 +DEAL::57 0 0.467 +DEAL::60 0 0.467 +DEAL::193 0 0.467 +DEAL::38 1 0.935 +DEAL::159 1 0.935 +DEAL::160 1 0.935 +DEAL::290 1 0.935 -- 2.39.5