From: Luca Heltai Date: Tue, 17 Apr 2018 08:22:05 +0000 (+0200) Subject: Have all GridGenerator functions attach curved manifolds. X-Git-Tag: v9.0.0-rc1~154^2~7 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c6ff1dfb4459df7ccf82444c72f0725a9e1ca80b;p=dealii.git Have all GridGenerator functions attach curved manifolds. Co-authored-by: David Wells --- diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 064ff48464..ac34516103 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -33,11 +33,18 @@ DEAL_II_NAMESPACE_OPEN * 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 */ @@ -54,7 +61,7 @@ namespace GridGenerator * 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 @@ -138,18 +145,14 @@ namespace GridGenerator * GeometryInfo class. Importantly, however, in 3d colorization does not set * @p boundary_ids of edges, but only of faces, 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 xi adds * 2i. 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 @@ -176,9 +179,9 @@ namespace GridGenerator * 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 xi adds 2i. For @@ -260,7 +263,7 @@ namespace GridGenerator const std::vector< std::vector > &spacing, const Point &p, const Table &material_id, - const bool colorize=false); + const bool colorize=false); /** * \brief Rectangular domain with rectangular pattern of holes @@ -452,29 +455,36 @@ namespace GridGenerator * 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 void hyper_ball (Triangulation &tria, const Point ¢er = Point(), - 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 @@ -483,9 +493,6 @@ namespace GridGenerator * 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 * @@ -509,11 +516,11 @@ namespace GridGenerator * 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 void quarter_hyper_ball (Triangulation &tria, @@ -526,10 +533,8 @@ namespace GridGenerator * x-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. */ @@ -556,6 +561,9 @@ namespace GridGenerator * 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 @@ -578,11 +586,10 @@ namespace GridGenerator * * 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. * @@ -702,11 +709,8 @@ namespace GridGenerator * 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 n_cells of elements for this initial * triangulation can be chosen arbitrarily. If the number of initial cells @@ -781,6 +785,9 @@ namespace GridGenerator * 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 @@ -876,9 +883,16 @@ namespace GridGenerator /** - * This class produces a square in the xy-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 xy-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 * diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 9f461d62b3..b4e09ed732 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -1985,7 +1985,8 @@ namespace GridGenerator template <> void hyper_ball (Triangulation<1> &, const Point<1> &, - const double) + const double, + const bool) { Assert (false, ExcNotImplemented()); } @@ -2265,7 +2266,6 @@ namespace GridGenerator cell->face(0)->set_boundary_id(0); cell->face(1)->set_boundary_id(4); cell->face(3)->set_boundary_id(5); - } } @@ -2277,7 +2277,8 @@ namespace GridGenerator 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 @@ -2307,12 +2308,17 @@ namespace GridGenerator 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 >(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)); } @@ -2379,6 +2385,9 @@ namespace GridGenerator if (colorize) colorize_hyper_shell(tria, center, inner_radius, outer_radius); + + tria.set_all_manifold_ids(0); + tria.set_manifold(0, SphericalManifold<2>(center)); } @@ -2473,6 +2482,8 @@ namespace GridGenerator Triangulation::cell_iterator cell = tria.begin(); Triangulation::cell_iterator end = tria.end(); + tria.set_all_manifold_ids_on_boundary(0); + while (cell != end) { for (unsigned int i=0; i::faces_per_cell; ++i) @@ -2484,10 +2495,14 @@ namespace GridGenerator // 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)); } @@ -2534,6 +2549,7 @@ namespace GridGenerator 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) { @@ -2544,10 +2560,14 @@ namespace GridGenerator // 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)); } @@ -2632,6 +2652,8 @@ namespace GridGenerator tria.last()->face(1)->set_boundary_id(2); } + tria.set_all_manifold_ids(0); + tria.set_manifold(0, SphericalManifold<2>(center)); } @@ -2712,6 +2734,9 @@ namespace GridGenerator tria.last()->face(1)->set_boundary_id(2); } + + tria.set_all_manifold_ids(0); + tria.set_manifold(0, SphericalManifold<2>(center)); } @@ -2896,6 +2921,7 @@ namespace GridGenerator } 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) @@ -2903,6 +2929,7 @@ namespace GridGenerator 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); @@ -2911,6 +2938,7 @@ namespace GridGenerator 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); @@ -2919,6 +2947,8 @@ namespace GridGenerator for (unsigned int i = 0; i < 4; ++i) cell->face (i)->set_boundary_id (0); } + + triangulation.set_manifold(0, CylindricalManifold<3>()); } @@ -3002,7 +3032,8 @@ namespace GridGenerator 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 @@ -3052,12 +3083,17 @@ namespace GridGenerator for (unsigned int j=0; j::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 >(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)); } @@ -3074,6 +3110,8 @@ namespace GridGenerator boundary_ids.insert (0); GridGenerator::extract_boundary_mesh (volume_mesh, tria, boundary_ids); + tria.set_all_manifold_ids(0); + tria.set_manifold(0, SphericalManifold(p)); } @@ -3163,6 +3201,8 @@ namespace GridGenerator 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::faces_per_cell; ++i) if (cell->at_boundary(i)) @@ -3170,26 +3210,35 @@ namespace GridGenerator 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::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::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>()); } @@ -3245,6 +3294,7 @@ namespace GridGenerator Triangulation::cell_iterator cell = tria.begin(); Triangulation::cell_iterator end = tria.end(); + tria.set_all_manifold_ids_on_boundary(0); while (cell != end) { for (unsigned int i=0; i::faces_per_cell; ++i) @@ -3258,6 +3308,7 @@ namespace GridGenerator || 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::lines_per_face; ++j) @@ -3271,13 +3322,17 @@ namespace GridGenerator || (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)); } @@ -3346,6 +3401,8 @@ namespace GridGenerator 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 @@ -3364,6 +3421,7 @@ namespace GridGenerator 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::lines_per_face; ++j) { const Point<3> line_vertices[2] @@ -3375,12 +3433,16 @@ namespace GridGenerator || (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)); } @@ -3608,6 +3670,8 @@ namespace GridGenerator if (colorize) colorize_hyper_shell(tria, p, inner_radius, outer_radius); + tria.set_all_manifold_ids(0); + tria.set_manifold(0, SphericalManifold<3>(p)); } @@ -3719,6 +3783,8 @@ namespace GridGenerator } } } + tria.set_all_manifold_ids(0); + tria.set_manifold(0, SphericalManifold<3>(center)); } @@ -3783,6 +3849,9 @@ namespace GridGenerator 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)); } @@ -3865,6 +3934,8 @@ namespace GridGenerator tria.create_triangulation ( vertices_3d, cells, SubCellData()); + tria.set_all_manifold_ids(0); + tria.set_manifold(0, CylindricalManifold<3>(2)); } @@ -4263,18 +4334,25 @@ namespace GridGenerator 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)); } @@ -4373,9 +4451,8 @@ namespace GridGenerator else { - cell->face(f)->set_boundary_id(6); - for (unsigned int l=0; l::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); } } @@ -4386,15 +4463,15 @@ namespace GridGenerator double d = c.norm(); if (d-inner_radius < 0) { - cell->face(f)->set_boundary_id(1); - for (unsigned int l=0; l::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