* modeled by a 2d half shell.
*
* <li> Slit domain: The slit domain is a variant of the hyper cube
- * domain. In two spatial dimensions, it is a square into which a slit
- * is sawed; if the initial square is though to be composed of four
- * smaller squares, then two of them are not connected even though
- * they are neighboring each other. Analogously, into the cube in
- * three spatial dimensions, a half-plane is sawed, disconnecting four
- * of the eight child-cubes from one of their neighbors.
+ * domain. In two spatial dimensions, it is a square into which a
+ * slit is sawed; if the initial square is though to be composed
+ * of four smaller squares, then two of them are not connected
+ * even though they are neighboring each other. Analogously, into
+ * the cube in three spatial dimensions, a half-plane is sawed,
+ * disconnecting four of the eight child-cubes from one of their
+ * neighbors.
+ *
+ * <li> Hyper cube with cylindrical hole: This domain is a square on
+ * the xy plane times the interval [0,L] with a cylindrical hole
+ * in the middle. The parameters that can be set are the internal
+ * radius, the external radius (inteded as the radius of the
+ * biggest enclosed cylinder), the depth of the structure (used
+ * only in three dimensions), and the number of repetitions along
+ * the z direction.
+ *
+ * @image html cubes_hole.png
+ *
* </ul>
*
* Some of these functions receive a flag @p colorize. If this is
const double inner_radius,
const double outer_radius,
const unsigned int n_cells = 0);
-
+
+ /**
+ * This class produces a square
+ * on the @p xy-plane with a
+ * circular hole in the middle,
+ * times the interval @p [0.L]
+ * (only in 3d). It is
+ * implemented in 2d and 3d, and
+ * takes the following arguments:
+ *
+ * - @p inner_radius: size of the
+ * internal hole
+ * - @p outer_radius: size of the
+ * biggest enclosed cylinder
+ * - @p L: extension on the @p z-direction
+ * - @p repetitions: number of subdivisions
+ * along the @z-direction
+ * - @p colorize: wether to assign different
+ * boundary indicators to different faces.
+ * The colors are given in lexycografical
+ * ordering for the flat faces (0 to 3 in 2d,
+ * 0 to 5 in 3d) plus the curved hole
+ * (4 in 2d, and 6 in 3d).
+ * If @colorize is set to false, then flat faces
+ * get the number 0 and the hole gets number 1.
+ */
+ template<int dim>
+ static void hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
+ const double inner_radius = .25,
+ const double outer_radius = .5,
+ const double L = .5,
+ const unsigned int repetition = 1,
+ const bool colorize = false);
/**
* Produce a ring of cells in 3D that is
#endif
+#if deal_II_dimension == 1
+
+template<int dim>
+void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &,
+ const double,
+ const double,
+ const double,
+ const unsigned int,
+ bool){
+ Assert(false, ExcNotImplemented());
+}
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template<int dim>
+void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<dim> &triangulation,
+ const double inner_radius,
+ const double outer_radius,
+ const double, // width,
+ const unsigned int, // width_repetition,
+ bool colorize)
+{
+ Assert(inner_radius < outer_radius,
+ ExcMessage("outer_radius has to be bigger than inner_radius."));
+
+ Point<dim> center;
+ // We create an hyper_shell in two dimensions, and then we modify it.
+ GridGenerator::hyper_shell (triangulation,
+ center, inner_radius, outer_radius,
+ 8);
+ typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active(),
+ endc = triangulation.end();
+ std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
+ for(; cell != endc; ++cell) {
+ for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if(cell->face(f)->at_boundary()) {
+ for(unsigned int v=0; v < GeometryInfo<dim>::vertices_per_face; ++v) {
+ unsigned int vv = cell->face(f)->vertex_index(v);
+ if(treated_vertices[vv] == false) {
+ treated_vertices[vv] = true;
+ switch(vv) {
+ case 1:
+ cell->face(f)->vertex(v) = center+Point<dim>(outer_radius,outer_radius);
+ break;
+ case 3:
+ cell->face(f)->vertex(v) = center+Point<dim>(-outer_radius,outer_radius);
+ break;
+ case 5:
+ cell->face(f)->vertex(v) = center+Point<dim>(-outer_radius,-outer_radius);
+ break;
+ case 7:
+ cell->face(f)->vertex(v) = center+Point<dim>(outer_radius,-outer_radius);
+ default:
+ break;
+ }
+ }
+ }
+ }
+ }
+ double eps = 1e-3 * outer_radius;
+ cell = triangulation.begin_active();
+ for(; cell != endc; ++cell) {
+ for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if(cell->face(f)->at_boundary()) {
+ double dx = cell->face(f)->center()(0) - center(0);
+ double dy = cell->face(f)->center()(1) - center(1);
+ if(colorize) {
+ if(std::abs(dx + outer_radius) < eps)
+ cell->face(f)->set_boundary_indicator(0);
+ else if(std::abs(dx - outer_radius) < eps)
+ cell->face(f)->set_boundary_indicator(1);
+ else if(std::abs(dy + outer_radius) < eps)
+ cell->face(f)->set_boundary_indicator(2);
+ else if(std::abs(dy - outer_radius) < eps)
+ cell->face(f)->set_boundary_indicator(3);
+ else
+ cell->face(f)->set_boundary_indicator(4);
+ } else {
+ double d = (cell->face(f)->center() - center).norm();
+ if(d-inner_radius < 0)
+ cell->face(f)->set_boundary_indicator(1);
+ else
+ cell->face(f)->set_boundary_indicator(0);
+ }
+ }
+ }
+}
+
+#endif
+
+#if deal_II_dimension == 3
+
+template<int dim>
+void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<dim> &triangulation,
+ const double inner_radius,
+ const double outer_radius,
+ const double L,
+ const unsigned int Nz,
+ bool colorize)
+{
+ Assert(inner_radius < outer_radius,
+ ExcMessage("outer_radius has to be bigger than inner_radius."));
+ Assert(L > 0,
+ ExcMessage("Must give positive extension L"));
+ Assert(Nz >= 1, ExcLowerRange(1, Nz));
+
+ GridGenerator::cylinder_shell (triangulation,
+ L, inner_radius, outer_radius,
+ 8,
+ Nz);
+
+ typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active(),
+ endc = triangulation.end();
+ std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
+ for(; cell != endc; ++cell) {
+ for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if(cell->face(f)->at_boundary()) {
+ for(unsigned int v=0; v < GeometryInfo<dim>::vertices_per_face; ++v) {
+ unsigned int vv = cell->face(f)->vertex_index(v);
+ if(treated_vertices[vv] == false) {
+ treated_vertices[vv] = true;
+ for(unsigned int i=0; i<=Nz; ++i) {
+ double d = ((double) i)*L/((double) Nz);
+ switch(vv-i*16) {
+ case 1:
+ cell->face(f)->vertex(v) = Point<dim>(outer_radius,outer_radius,d);
+ break;
+ case 3:
+ cell->face(f)->vertex(v) = Point<dim>(-outer_radius,outer_radius,d);
+ break;
+ case 5:
+ cell->face(f)->vertex(v) = Point<dim>(-outer_radius,-outer_radius,d);
+ break;
+ case 7:
+ cell->face(f)->vertex(v) = Point<dim>(outer_radius,-outer_radius,d);
+ break;
+ default:
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ double eps = 1e-3 * outer_radius;
+ cell = triangulation.begin_active();
+ for(; cell != endc; ++cell) {
+ for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if(cell->face(f)->at_boundary()) {
+ double dx = cell->face(f)->center()(0);
+ double dy = cell->face(f)->center()(1);
+ double dz = cell->face(f)->center()(2);
+
+ if(colorize) {
+ if(std::abs(dx + outer_radius) < eps)
+ cell->face(f)->set_boundary_indicator(0);
+
+ else if(std::abs(dx - outer_radius) < eps)
+ cell->face(f)->set_boundary_indicator(1);
+
+ else if(std::abs(dy + outer_radius) < eps)
+ cell->face(f)->set_boundary_indicator(2);
+
+ else if(std::abs(dy - outer_radius) < eps)
+ cell->face(f)->set_boundary_indicator(3);
+
+ else if(std::abs(dz) < eps)
+ cell->face(f)->set_boundary_indicator(4);
+
+ else if(std::abs(dz - L) < eps)
+ cell->face(f)->set_boundary_indicator(5);
+
+ else {
+ cell->face(f)->set_boundary_indicator(6);
+ for(unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+ cell->face(f)->line(l)->set_boundary_indicator(6);
+ }
+
+ } else {
+ Point<dim> c = cell->face(f)->center();
+ c(2) = 0;
+ double d = c.norm();
+ if(d-inner_radius < 0) {
+ cell->face(f)->set_boundary_indicator(1);
+ for(unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+ cell->face(f)->line(l)->set_boundary_indicator(1);
+ } else
+ cell->face(f)->set_boundary_indicator(0);
+ }
+ }
+ }
+}
+
+#endif
+
+
// explicit instantiations
template void
GridGenerator::hyper_cube<deal_II_dimension> (
const Point<deal_II_dimension>&, double, double, unsigned int);
+template void
+GridGenerator::hyper_cube_with_cylindrical_hole (
+ Triangulation<deal_II_dimension> &,
+ const double, const double, const double, const unsigned int, bool);
+
template void
GridGenerator::