* This namespace provides a collection of functions for generating
* triangulations for some basic geometries.
*
- * Some of these functions receive a flag @p colorize. If this is set, parts
- * of the boundary receive different
- * @ref GlossBoundaryIndicator "boundary indicators"),
- * allowing them to be distinguished for the purpose of attaching geometry
- * objects and evaluating different boundary conditions.
+ * Some of these functions receive a flag @p colorize. If this is set, parts of
+ * the boundary receive different @ref GlossBoundaryIndicator "boundary
+ * indicators" allowing them to be distinguished for the purpose of evaluating
+ * different boundary conditions.
+ *
+ * If the domain is curved, each of the domain parts that should be
+ * refined by following an appropriate Manifold description will
+ * receive a different @ref GlossManifoldIndicator "manifold
+ * indicator", and the correct Manifold descriptor will be attached to
+ * the Triangulation. Notice that if you later tranform the
+ * triangulation, you have to make sure you attach the correct new Manifold
+ * to the triangulation.
*
* @ingroup grid
*/
* tensor product interval $[left,right]^{\text{dim}}$ in the present number
* of dimensions, where the limits are given as arguments. They default to
* zero and unity, then producing the unit hypercube.
- *
+ *
* If the argument @p colorize is false, all boundary indicators are set to
* zero ("not colorized") for 2d and 3d. If it is true, the boundary is
* colorized as in hyper_rectangle(). In 1d the indicators are always
* GeometryInfo class. Importantly, however, in 3d colorization does not set
* @p boundary_ids of <i>edges</i>, but only of <i>faces</i>, because each
* boundary edge is shared between two faces and it is not clear how the
- * boundary id of an edge should be set in that case. This may later on lead
- * to problems if one wants to assign boundary or manifold objects to parts
- * of the boundary with certain boundary indicators since then the boundary
- * object may not apply to the edges bounding the face it is meant to
- * describe.
+ * boundary id of an edge should be set in that case.
*
* Additionally, if @p colorize is @p true, material ids are assigned to the
* cells according to the octant their center is in: being in the right half
* space for any coordinate direction <i>x<sub>i</sub></i> adds
* 2<sup>i</sup>. For instance, a cell with center point (1,-1,1) yields a
* material id 5, assuming that the center of the hyper rectangle lies at
- * the origin.
+ * the origin. No manifold id is set for the cells.
*
* If @p dim < @p spacedim, this will create a @p dim dimensional object in
* the first @p dim coordinate directions embedded into the @p spacedim
* direction is 1.
*
* If the @p colorize flag is set, the @p boundary_ids of the surfaces are
- * assigned, such that the lower one in @p x-direction is 0, the upper one
- * is 1 (the left and the right vertical face). The indicators for the
- * surfaces in @p y-direction are 2 and 3, the ones for @p z are 4 and 5.
+ * assigned, such that the lower one in @p x-direction is 0, the upper one is
+ * 1 (the left and the right vertical face). The indicators for the surfaces
+ * in @p y-direction are 2 and 3, the ones for @p z are 4 and 5.
* Additionally, material ids are assigned to the cells according to the
* octant their center is in: being in the right half plane for any
* coordinate direction <i>x<sub>i</sub></i> adds 2<sup>i</sup>. For
const std::vector< std::vector<double> > &spacing,
const Point<dim> &p,
const Table<dim,types::material_id> &material_id,
- const bool colorize=false);
+ const bool colorize=false);
/**
* \brief Rectangular domain with rectangular pattern of holes
* This function is declared to exist for triangulations of all space
* dimensions, but throws an error if called in 1d.
*
- * You should attach a SphericalManifold to the cells and faces for correct
- * placement of vertices upon refinement and to be able to use higher order
- * mappings. However, it turns out that creating a mesh for a hyperball is
- * not entirely trivial since the central cell has to be treated
- * differently than the "cap" cells. The "Possibilities for extensions"
- * section of step-6 has an extensive discussion of how one would construct
- * such meshes and what one needs to do for it.
+ * By default, the manifold_id is set to 0 on the boundary faces, 1 on the
+ * the boundary cells, and types::flat_manifold_id on the central cell and on
+ * internal faces.
+ *
+ * A SphericalManifold is attached by default to the boundary faces for
+ * correct placement of boundary vertices upon refinement and to be able to
+ * use higher order mappings. However, it turns out that this strategy may
+ * not be the optimal one to create a good a mesh for a hyperball. The
+ * "Possibilities for extensions" section of step-6 has an extensive
+ * discussion of how one would construct better meshes and what one needs to
+ * do for it. Selecting the argument @p
+ * attach_spherical_manifold_on_boundary_cells to true attaches a
+ * SphericalManifold manifold also to the boundary cells, and not only to the
+ * boundary faces.
*
* @note The triangulation passed as argument needs to be empty when calling this function.
*/
template <int dim>
void hyper_ball (Triangulation<dim> &tria,
const Point<dim> ¢er = Point<dim>(),
- const double radius = 1.);
+ const double radius = 1.,
+ const bool attach_spherical_manifold_on_boundary_cells=false);
/**
* Creates a hyper sphere, i.e., a surface of a ball in @p spacedim
* dimensions. This function only exists for dim+1=spacedim in 2 and 3 space
* dimensions. (To create a mesh of a ball, use GridGenerator::hyper_ball().)
*
- * You should attach a SphericalManifold to the cells and faces for correct
- * placement of vertices upon refinement and to be able to use higher order
- * mappings.
+ * By default, all manifold ids of the triangulation are set to zero, and a
+ * SphericalManifold is attached to the grid.
*
* The following pictures are generated with:
* @code
* static SphericalManifold<2,3> surface_description;
*
* GridGenerator::hyper_sphere(triangulation);
- *
- * triangulation.set_all_manifold_ids(0);
- * triangulation.set_manifold (0, surface_description);
* triangulation.refine_global(3);
* @endcode
*
* relative to @p center, which contains three elements in 2d and four in 3d.
*
* The boundary indicators for the final triangulation are 0 for the curved
- * boundary and 1 for the cut plane.
- *
- * The appropriate boundary class is HyperBallBoundary.
+ * boundary and 1 for the cut plane. The manifold id for the curved boundary is
+ * set to zero, and a SphericalManifold is attached to it.
*
- * @note The triangulation passed as argument needs to be empty when calling this function.
+ * @note The triangulation passed as argument needs to be empty when calling
+ * this function.
*/
template <int dim>
void quarter_hyper_ball (Triangulation<dim> &tria,
* <i>x</i>-axis.
*
* 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.
+ * boundary and 1 for the cut plane. The manifold id for the curved boundary is
+ * set to zero, and a SphericalManifold is attached to it.
*
* @note The triangulation passed as argument needs to be empty when calling this function.
*/
* the GridTools::transform() function using a rotation operator as
* argument.
*
+ * The manifold id for the hull of the cylinder is set to zero, and a
+ * CylindricalManifold is attached to it.
+ *
* @note The triangulation passed as argument needs to be empty when calling this function.
*/
template <int dim>
*
* The boundaries are colored according to the following scheme: 0 for the
* hull of the cone, 1 for the left hand face and 2 for the right hand face.
+ * Both the boundary indicators and the manifold indicators are set.
*
- * In three dimensions, the CylindricalManifold class is an appropriate choice
- * for the description of the hull, with which you probably want to associate
- * boundary indicator 0.
- * In two dimensions the default FlatManifold is sufficient.
+ * In three dimensions, the manifold id of the hull is set to zero, and a
+ * CylindricalManifold is attached to it.
*
* @note The triangulation passed as argument needs to be empty when calling this function.
*
* both the faces and the edges of these boundaries. If the flag is @p
* false, both have indicator zero.
*
- * You should attach a SphericalManifold to the cells and faces for correct
- * placement of vertices upon refinement and to be able to use higher order
- * mappings. Alternatively, it is also possible to attach a
- * HyperShellBoundary to the inner and outer boundary. This will create
- * inferior meshes as described below.
+ * All manifold ids are set to zero, and a SphericalManifold is attached to the
+ * the triangulation.
*
* In 2d, the number <tt>n_cells</tt> of elements for this initial
* triangulation can be chosen arbitrarily. If the number of initial cells
* boundary where $x=0$, get indicator 0, 1, and 2, respectively. Otherwise
* all indicators are set to 0.
*
+ * All manifold ids are set to zero, and a SphericalManifold is attached
+ * to the triangulation.
+ *
* @note The triangulation passed as argument needs to be empty when calling this function.
*/
template <int dim>
/**
- * This class produces a square in the <i>xy</i>-plane with a circular hole
- * in the middle. Square and circle are centered at the origin. In 3d, this
- * geometry is extruded in $z$ direction to the interval $[0,L]$.
+ * This class produces a square in the <i>xy</i>-plane with a cylindrical
+ * hole in the middle. The square and the circle are centered at the
+ * origin. In 3d, this geometry is extruded in $z$ direction to the interval
+ * $[0,L]$.
+ *
+ * The inner boundary has a manifold id of $0$ and a boundary id of
+ * $6$. This function attaches a PolarManifold or CylindricalManifold to the
+ * interior boundary in 2d and 3d respectively. The other faces have
+ * boundary ids of $0, 1, 2, 3, 4$, or $5$ given in the standard order of
+ * faces in 2d or 3d.
*
* @image html cubes_hole.png
*
template <>
void hyper_ball (Triangulation<1> &,
const Point<1> &,
- const double)
+ const double,
+ const bool)
{
Assert (false, ExcNotImplemented());
}
cell->face(0)->set_boundary_id(0);
cell->face(1)->set_boundary_id(4);
cell->face(3)->set_boundary_id(5);
-
}
}
void
hyper_ball (Triangulation<2> &tria,
const Point<2> &p,
- const double radius)
+ const double radius,
+ const bool internal_manifolds)
{
// equilibrate cell sizes at
// transition from the inner part
for (unsigned int j=0; j<4; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
cells[i].material_id = 0;
+ cells[i].manifold_id = i==2 ? 1 : numbers::flat_manifold_id;
};
tria.create_triangulation (
std::vector<Point<2> >(std::begin(vertices), std::end(vertices)),
cells,
SubCellData()); // no boundary information
+ tria.set_all_manifold_ids_on_boundary(0);
+ tria.set_manifold(0, SphericalManifold<2>(p));
+ if (internal_manifolds)
+ tria.set_manifold(1, SphericalManifold<2>(p));
}
if (colorize)
colorize_hyper_shell(tria, center, inner_radius, outer_radius);
+
+ tria.set_all_manifold_ids(0);
+ tria.set_manifold(0, SphericalManifold<2>(center));
}
Triangulation<dim>::cell_iterator cell = tria.begin();
Triangulation<dim>::cell_iterator end = tria.end();
+ tria.set_all_manifold_ids_on_boundary(0);
+
while (cell != end)
{
for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
// component of the center, then this is part of the plane
if (cell->face(i)->center()(0) < p(0)+1.e-5 * radius
|| cell->face(i)->center()(1) < p(1)+1.e-5 * radius)
- cell->face(i)->set_boundary_id(1);
+ {
+ cell->face(i)->set_boundary_id(1);
+ cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
+ }
}
++cell;
}
+ tria.set_manifold(0, SphericalManifold<2>(p));
}
Triangulation<2>::cell_iterator cell = tria.begin();
Triangulation<2>::cell_iterator end = tria.end();
+ tria.set_all_manifold_ids_on_boundary(0);
while (cell != end)
{
// If x is zero, then this is part of the plane
if (cell->face(i)->center()(0) < p(0)+1.e-5 * radius)
- cell->face(i)->set_boundary_id(1);
+ {
+ cell->face(i)->set_boundary_id(1);
+ cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
+ }
}
++cell;
}
+ tria.set_manifold(0, SphericalManifold<2>(p));
}
tria.last()->face(1)->set_boundary_id(2);
}
+ tria.set_all_manifold_ids(0);
+ tria.set_manifold(0, SphericalManifold<2>(center));
}
tria.last()->face(1)->set_boundary_id(2);
}
+
+ tria.set_all_manifold_ids(0);
+ tria.set_manifold(0, SphericalManifold<2>(center));
}
}
triangulation.create_triangulation (vertices, cells, SubCellData ());
+ triangulation.set_all_manifold_ids_on_boundary(0);
for (Triangulation<3>::cell_iterator cell = triangulation.begin ();
cell != triangulation.end (); ++cell)
if (cell->vertex (0) (0) == -half_length)
{
cell->face (4)->set_boundary_id (1);
+ cell->face (4)->set_manifold_id(numbers::flat_manifold_id);
for (unsigned int i = 0; i < 4; ++i)
cell->line (i)->set_boundary_id (0);
if (cell->vertex (4) (0) == half_length)
{
cell->face (5)->set_boundary_id (2);
+ cell->face (5)->set_manifold_id(numbers::flat_manifold_id);
for (unsigned int i = 4; i < 8; ++i)
cell->line (i)->set_boundary_id (0);
for (unsigned int i = 0; i < 4; ++i)
cell->face (i)->set_boundary_id (0);
}
+
+ triangulation.set_manifold(0, CylindricalManifold<3>());
}
void
hyper_ball (Triangulation<3> &tria,
const Point<3> &p,
- const double radius)
+ const double radius,
+ const bool internal_manifold)
{
const double a = 1./(1+std::sqrt(3.0)); // equilibrate cell sizes at transition
// from the inner part to the radial
for (unsigned int j=0; j<GeometryInfo<3>::vertices_per_cell; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
cells[i].material_id = 0;
+ cells[i].manifold_id = i == 0 ? numbers::flat_manifold_id : 1;
};
tria.create_triangulation (
std::vector<Point<3> >(std::begin(vertices), std::end(vertices)),
cells,
SubCellData()); // no boundary information
+ tria.set_all_manifold_ids_on_boundary(0);
+ tria.set_manifold(0, SphericalManifold<3>(p));
+ if (internal_manifold)
+ tria.set_manifold(1, SphericalManifold<3>(p));
}
boundary_ids.insert (0);
GridGenerator::extract_boundary_mesh (volume_mesh, tria,
boundary_ids);
+ tria.set_all_manifold_ids(0);
+ tria.set_manifold(0, SphericalManifold<spacedim-1, spacedim>(p));
}
Triangulation<3>::cell_iterator cell = tria.begin();
Triangulation<3>::cell_iterator end = tria.end();
+ tria.set_all_manifold_ids_on_boundary(0);
+
for (; cell != end; ++cell)
for (unsigned int i=0; i<GeometryInfo<3>::faces_per_cell; ++i)
if (cell->at_boundary(i))
if (cell->face(i)->center()(0) > half_length-1.e-5)
{
cell->face(i)->set_boundary_id(2);
+ cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
for (unsigned int e=0; e<GeometryInfo<3>::lines_per_face; ++e)
if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
(std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
(std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
(std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
- cell->face(i)->line(e)->set_boundary_id(2);
+ {
+ cell->face(i)->line(e)->set_boundary_id(2);
+ cell->face(i)->line(e)->set_manifold_id(numbers::flat_manifold_id);
+ }
}
else if (cell->face(i)->center()(0) < -half_length+1.e-5)
{
cell->face(i)->set_boundary_id(1);
+ cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
for (unsigned int e=0; e<GeometryInfo<3>::lines_per_face; ++e)
if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
(std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
(std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
(std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
- cell->face(i)->line(e)->set_boundary_id(1);
+ {
+ cell->face(i)->line(e)->set_boundary_id(1);
+ cell->face(i)->line(e)->set_manifold_id(numbers::flat_manifold_id);
+ }
}
}
+ tria.set_manifold(0, CylindricalManifold<3>());
}
Triangulation<dim>::cell_iterator cell = tria.begin();
Triangulation<dim>::cell_iterator end = tria.end();
+ tria.set_all_manifold_ids_on_boundary(0);
while (cell != end)
{
for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
|| cell->face(i)->center()(2) < center(2)+1.e-5 * radius)
{
cell->face(i)->set_boundary_id(1);
+ cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
// also set the boundary indicators of the bounding lines,
// unless both vertices are on the perimeter
for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
||
(std::fabs(line_vertices[1].distance(center)-radius) >
1e-5*radius))
- cell->face(i)->line(j)->set_boundary_id(1);
+ {
+ cell->face(i)->line(j)->set_boundary_id(1);
+ cell->face(i)->line(j)->set_manifold_id(numbers::flat_manifold_id);
+ }
}
}
}
++cell;
}
+ tria.set_manifold(0, SphericalManifold<3>(center));
}
Triangulation<3>::cell_iterator cell = tria.begin();
Triangulation<3>::cell_iterator end = tria.end();
+ tria.set_all_manifold_ids_on_boundary(0);
+
// go over all faces. for the ones on the flat face, set boundary
// indicator for face and edges to one; the rest will remain at
// zero but we have to pay attention to those edges that are
if (cell->face(i)->center()(0) < center(0)+1.e-5*radius)
{
cell->face(i)->set_boundary_id(1);
+ cell->face(i)->set_manifold_id(numbers::flat_manifold_id);
for (unsigned int j=0; j<GeometryInfo<3>::lines_per_face; ++j)
{
const Point<3> line_vertices[2]
||
(std::fabs(line_vertices[1].distance(center)-radius) >
1e-5*radius))
- cell->face(i)->line(j)->set_boundary_id(1);
+ {
+ cell->face(i)->line(j)->set_boundary_id(1);
+ cell->face(i)->line(j)->set_manifold_id(numbers::flat_manifold_id);
+ }
}
}
}
++cell;
}
+ tria.set_manifold(0, SphericalManifold<3>(center));
}
if (colorize)
colorize_hyper_shell(tria, p, inner_radius, outer_radius);
+ tria.set_all_manifold_ids(0);
+ tria.set_manifold(0, SphericalManifold<3>(p));
}
}
}
}
+ tria.set_all_manifold_ids(0);
+ tria.set_manifold(0, SphericalManifold<3>(center));
}
if (colorize)
colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
+
+ tria.set_all_manifold_ids(0);
+ tria.set_manifold(0, SphericalManifold<3>(center));
}
tria.create_triangulation (
vertices_3d, cells, SubCellData());
+ tria.set_all_manifold_ids(0);
+ tria.set_manifold(0, CylindricalManifold<3>(2));
}
else if (std::abs(dy - outer_radius) < eps)
cell->face(f)->set_boundary_id(3);
else
- cell->face(f)->set_boundary_id(4);
+ {
+ cell->face(f)->set_boundary_id(4);
+ cell->face(f)->set_manifold_id(0);
+ }
}
else
{
double d = (cell->face(f)->center() - center).norm();
if (d-inner_radius < 0)
- cell->face(f)->set_boundary_id(1);
+ {
+ cell->face(f)->set_boundary_id(1);
+ cell->face(f)->set_manifold_id(0);
+ }
else
cell->face(f)->set_boundary_id(0);
}
}
}
+ triangulation.set_manifold(0, PolarManifold<2>(center));
}
else
{
- cell->face(f)->set_boundary_id(6);
- for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
- cell->face(f)->line(l)->set_boundary_id(6);
+ cell->face(f)->set_all_boundary_ids(6);
+ cell->face(f)->set_all_manifold_ids(0);
}
}
double d = c.norm();
if (d-inner_radius < 0)
{
- cell->face(f)->set_boundary_id(1);
- for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
- cell->face(f)->line(l)->set_boundary_id(1);
+ cell->face(f)->set_all_boundary_ids(1);
+ cell->face(f)->set_all_manifold_ids(0);
}
else
cell->face(f)->set_boundary_id(0);
}
}
}
+ triangulation.set_manifold(0, CylindricalManifold<3>(2));
}
template <int dim, int spacedim1, int spacedim2>
projects / dealii.git / commitdiff