]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix a problem with colorize_hyper_shell which wasn't written in a particularly defens...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 15 Feb 2011 18:31:55 +0000 (18:31 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 15 Feb 2011 18:31:55 +0000 (18:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@23370 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/grid/grid_generator.cc
tests/deal.II/grid_hyper_shell_05.cc [new file with mode: 0644]
tests/deal.II/grid_hyper_shell_05/cmp/generic [new file with mode: 0644]

index 9666862117e7dce8e5bbde07db3b68b6c772192b..dbafe0257d2f5097421ebb839709e1aa3a01429f 100644 (file)
@@ -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 (file)
index 0000000..a287aee
--- /dev/null
@@ -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 <grid/tria_boundary_lib.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_out.h>
+#include <grid/grid_tools.h>
+#include <base/logstream.h>
+#include <cmath>
+#include <cstdlib>
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+
+std::ofstream logfile("grid_hyper_shell_05/output");
+
+
+template<int dim>
+void check (const unsigned int n)
+{
+  deallog << "n=" << n << std::endl;
+
+  Point<dim> center;
+  Triangulation<dim> tria (Triangulation<dim>::none);
+  GridGenerator::hyper_shell (tria, center, 0.5, 1, n, true);
+
+                                  // this is the test that failed
+                                  // before
+  for (typename Triangulation<dim>::active_cell_iterator
+        cell = tria.begin_active();
+       cell != tria.end(); ++cell)
+    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+      Assert (cell->face(f)->at_boundary() == cell->at_boundary(f),
+             ExcInternalError());
+
+                                  // also output something slightly
+                                  // more useful
+  for (typename Triangulation<dim>::active_cell_iterator
+        cell = tria.begin_active();
+       cell != tria.end(); ++cell)
+    for (unsigned int f=0; f<GeometryInfo<dim>::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 (file)
index 0000000..c7e5773
--- /dev/null
@@ -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

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.