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);
}
}
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());
}
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
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;
}
--- /dev/null
+//----------------------------------------------------------------------
+// $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);
+}