* the triangulation to be refined at the outer boundaries.
*
* <li> Half Hyper ball:
- * You get half of the circle generated by Hyper ball.
+ * You get half of the circle or sphere generated by Hyper ball.
* with center @p p and with radius @p r by calling
* GridGenerator::half_hyper_ball (tria, p, r). The half-circle is
- * triangulated by four cells. The diameter of the center cell is
+ * triangulated by four cells, while the half-cell is triangulated by
+ * six cells. The diameter of the center cell is
* chosen to be the same as for the Hyper ball class.
* To create a half-hyperball in one dimension results in
* an error.
*
* Do not forget to also attach a suitable boundary approximation object
* to the triangulation object you passed to this function if you later want
- * the triangulation to be refined at the outer boundaries. The class
- * HalfHyperBallBoundary will provide a boundary object.
+ * the triangulation to be refined at the outer boundaries. Both classes
+ * HalfHyperBallBoundary and HyperBallBoundary will provide a valid
+ * boundary object.
*
* <li> Hyper shell: A hyper shell is the region between two hyper
* sphere with the same origin. Therefore, it is a ring in two
/**
* This class produces a half hyper-ball,
- * which contains four elements.
+ * which contains four elements in 2d and
+ * 6 in 3d. The orientation is such that
+ * the cutting plane is the x=0.
*
* The triangulation needs to be
* void upon calling this
* function.
*
- * Currently, only a two-dimensional
- * version is implemented. The appropriate
- * boundary class is
- * HalfHyperBallBoundary.
+ * The boundary indicators for the final
+ * triangulation are 0 for the curved boundary and
+ * 1 for the cut plane.
*
+ * The appropriate
+ * boundary class is
+ * HalfHyperBallBoundary, or HyperBallBoundary.
+ *
* @note The triangulation needs to be
* void upon calling this
* function.
}
-//TODO: Colorize edges as circumference and cut plane
-// Implementation for 2D only
template <int dim>
void
GridGenerator::half_hyper_ball (Triangulation<dim> &tria,
std::vector<Point<dim> >(&vertices[0], &vertices[8]),
cells,
SubCellData()); // no boundary information
+
+ typename Triangulation<dim>::cell_iterator cell = tria.begin();
+ typename Triangulation<dim>::cell_iterator end = tria.end();
+
+
+ while (cell != end)
+ {
+ for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
+ {
+ if (cell->face(i)->boundary_indicator() == 255)
+ continue;
+
+ // If x is zero, then this is part of the plane
+ if (cell->face(i)->center()(0) < p(0)+1.e-5)
+ cell->face(i)->set_boundary_indicator(1);
+ }
+ ++cell;
+ }
}
// Implementation for 3D only
template <int dim>
void
-GridGenerator::half_hyper_ball (Triangulation<dim>&,
- const Point<dim>&,
- const double)
+GridGenerator::half_hyper_ball (Triangulation<dim>& tria,
+ const Point<dim>& center,
+ const double radius)
{
- Assert (false, ExcNotImplemented());
-}
+ // These are for the two lower squares
+ const double d = radius/std::sqrt(2.0);
+ const double a = d/(1+std::sqrt(2.0));
+ // These are for the two upper square
+ const double b = a/2.0;
+ const double c = d/2.0;
+ // And so are these
+ const double hb = a*std::sqrt(3.0)/2.0;
+ const double hc = d*std::sqrt(3.0)/2.0;
+
+ Point<dim> vertices[16] = {
+ center+Point<dim>( 0, d, -d),
+ center+Point<dim>( 0, -d, -d),
+ center+Point<dim>( 0, a, -a),
+ center+Point<dim>( 0, -a, -a),
+ center+Point<dim>( 0, a, a),
+ center+Point<dim>( 0, -a, a),
+ center+Point<dim>( 0, d, d),
+ center+Point<dim>( 0, -d, d),
+
+ center+Point<dim>(hc, c, -c),
+ center+Point<dim>(hc, -c, -c),
+ center+Point<dim>(hb, b, -b),
+ center+Point<dim>(hb, -b, -b),
+ center+Point<dim>(hb, b, b),
+ center+Point<dim>(hb, -b, b),
+ center+Point<dim>(hc, c, c),
+ center+Point<dim>(hc, -c, c),
+ };
+ int cell_vertices[6][8] = {
+ {0, 1, 8, 9, 2, 3, 10, 11},
+ {0, 2, 8, 10, 6, 4, 14, 12},
+ {2, 3, 10, 11, 4, 5, 12, 13},
+ {1, 7, 9, 15, 3, 5, 11, 13},
+ {6, 4, 14, 12, 7, 5, 15, 13},
+ {8, 10, 9, 11, 14, 12, 15, 13}
+ };
+
+ std::vector<CellData<dim> > cells (6, CellData<dim>());
+
+ for (unsigned int i=0; i<6; ++i)
+ {
+ for (unsigned int j=0; j<8; ++j)
+ cells[i].vertices[j] = cell_vertices[i][j];
+ cells[i].material_id = 0;
+ };
+
+ tria.create_triangulation (
+ std::vector<Point<dim> >(&vertices[0], &vertices[16]),
+ cells,
+ SubCellData()); // no boundary information
+
+ typename Triangulation<dim>::cell_iterator cell = tria.begin();
+ typename Triangulation<dim>::cell_iterator end = tria.end();
+
+ while (cell != end)
+ {
+ for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
+ {
+ if (cell->face(i)->boundary_indicator() == 255)
+ continue;
+
+ // If the center is on the plane x=0, this is a plane
+ // element
+ if (cell->face(i)->center()(0) < center(0)+1.e-5)
+ cell->face(i)->set_boundary_indicator(1);
+ }
+ ++cell;
+ }
+}
// Implementation for 3D only
2 1 0 3 1 -1 800 0 -1 0.0 0 0 -1 0 0 2
2400 -4800
4945 -3745
-2 1 0 3 1 -1 800 0 -1 0.0 0 0 -1 0 0 2
+2 1 0 3 2 -1 800 0 -1 0.0 0 0 -1 0 0 2
2400 -4800
2400 -2254
2 3 0 1 0 1 900 0 25 0.0 0 0 -1 0 0 5
3454 -145
2400 -145
2400 -2254
-2 1 0 3 1 -1 800 0 -1 0.0 0 0 -1 0 0 2
+2 1 0 3 2 -1 800 0 -1 0.0 0 0 -1 0 0 2
2400 -2254
2400 -145
2 3 0 1 0 1 900 0 25 0.0 0 0 -1 0 0 5
3454 -145
4945 1345
2400 2400
-2 1 0 3 1 -1 800 0 -1 0.0 0 0 -1 0 0 2
+2 1 0 3 2 -1 800 0 -1 0.0 0 0 -1 0 0 2
2400 2400
2400 -145
2 1 0 3 1 -1 800 0 -1 0.0 0 0 -1 0 0 2
member "cells" value "cell data"
end
DEAL:3d::half_hyper_ball
-DEAL:3d::Abort!!!
+object "vertices" class array type float rank 1 shape 3 items 16 data follows
+ 2.00000 1.12132 -2.12132
+ 2.00000 -3.12132 -2.12132
+ 2.00000 -0.121320 -0.878680
+ 2.00000 -1.87868 -0.878680
+ 2.00000 -0.121320 0.878680
+ 2.00000 -1.87868 0.878680
+ 2.00000 1.12132 2.12132
+ 2.00000 -3.12132 2.12132
+ 3.83712 0.0606602 -1.06066
+ 3.83712 -2.06066 -1.06066
+ 2.76096 -0.560660 -0.439340
+ 2.76096 -1.43934 -0.439340
+ 2.76096 -0.560660 0.439340
+ 2.76096 -1.43934 0.439340
+ 3.83712 0.0606602 1.06066
+ 3.83712 -2.06066 1.06066
+object "cells" class array type int rank 1 shape 8 items 6 data follows
+ 0 2 8 10 1 3 9 11
+ 0 6 8 14 2 4 10 12
+ 2 4 10 12 3 5 11 13
+ 1 3 9 11 7 5 15 13
+ 6 7 14 15 4 5 12 13
+ 8 14 9 15 10 12 11 13
+attribute "element type" string "cubes"
+attribute "ref" string "positions"
+
+object "material" class array type int rank 0 items 6 data follows
+ 0 0 0 0 0 0
+attribute "dep" string "connections"
+
+object "level" class array type int rank 0 items 6 data follows
+ 0 0 0 0 0 0
+attribute "dep" string "connections"
+
+object "deal data" class field
+component "positions" value "vertices"
+component "connections" value "cells"
+object "cell data" class field
+component "positions" value "vertices"
+component "connections" value "cells"
+component "material" value "material"
+component "level" value "level"
+
+object "grid data" class group
+member "cells" value "cell data"
+end
DEAL:3d::half_hyper_shell
DEAL:3d::Abort!!!