From: heltai Date: Thu, 15 Feb 2007 19:16:08 +0000 (+0000) Subject: Workaround for gmsh output mesh. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=43a3614d7ff679b9351d8653dbdc92a23bd05605;p=dealii-svn.git Workaround for gmsh output mesh. New funcition: GridGenerator::hyper_cube_with_cylindrical_hole. git-svn-id: https://svn.dealii.org/trunk@14477 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index ad85d444af..eb70ac94b2 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -109,12 +109,24 @@ template class SparseMatrix; * modeled by a 2d half shell. * *
  • 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. + * + *
  • 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 + * * * * Some of these functions receive a flag @p colorize. If this is @@ -661,7 +673,39 @@ class GridGenerator 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 + static void hyper_cube_with_cylindrical_hole (Triangulation &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 diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index 803b82ee8c..716c97b841 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -2245,6 +2245,207 @@ void GridGenerator::laplace_transformation (Triangulation &tria, #endif +#if deal_II_dimension == 1 + +template +void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation &, + const double, + const double, + const double, + const unsigned int, + bool){ + Assert(false, ExcNotImplemented()); +} + +#endif + + +#if deal_II_dimension == 2 + +template +void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation &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 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::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + std::vector treated_vertices(triangulation.n_vertices(), false); + for(; cell != endc; ++cell) { + for(unsigned int f=0; f::faces_per_cell; ++f) + if(cell->face(f)->at_boundary()) { + for(unsigned int v=0; v < GeometryInfo::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(outer_radius,outer_radius); + break; + case 3: + cell->face(f)->vertex(v) = center+Point(-outer_radius,outer_radius); + break; + case 5: + cell->face(f)->vertex(v) = center+Point(-outer_radius,-outer_radius); + break; + case 7: + cell->face(f)->vertex(v) = center+Point(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::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 +void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation &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::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + std::vector treated_vertices(triangulation.n_vertices(), false); + for(; cell != endc; ++cell) { + for(unsigned int f=0; f::faces_per_cell; ++f) + if(cell->face(f)->at_boundary()) { + for(unsigned int v=0; v < GeometryInfo::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(outer_radius,outer_radius,d); + break; + case 3: + cell->face(f)->vertex(v) = Point(-outer_radius,outer_radius,d); + break; + case 5: + cell->face(f)->vertex(v) = Point(-outer_radius,-outer_radius,d); + break; + case 7: + cell->face(f)->vertex(v) = Point(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::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::lines_per_face; ++l) + cell->face(f)->line(l)->set_boundary_indicator(6); + } + + } else { + Point 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::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 ( @@ -2324,6 +2525,11 @@ GridGenerator::half_hyper_shell ( const Point&, double, double, unsigned int); +template void +GridGenerator::hyper_cube_with_cylindrical_hole ( + Triangulation &, + const double, const double, const double, const unsigned int, bool); + template void GridGenerator:: diff --git a/deal.II/deal.II/source/grid/grid_out.cc b/deal.II/deal.II/source/grid/grid_out.cc index 91c18d21a1..6dd51f4bf9 100644 --- a/deal.II/deal.II/source/grid/grid_out.cc +++ b/deal.II/deal.II/source/grid/grid_out.cc @@ -743,10 +743,13 @@ void GridOut::write_msh_faces (const Triangulation &tria, Assert (false, ExcNotImplemented()); }; out << static_cast(face->boundary_indicator()) - << " 1 " << GeometryInfo::vertices_per_face; + << " " + << static_cast(face->boundary_indicator()) + << " " << GeometryInfo::vertices_per_face; // note: vertex numbers are 1-base for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) - out << ' ' << face->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1; + out << ' ' + << face->vertex_index(GeometryInfo::ucd_to_deal[vertex])+1; out << '\n'; ++index; diff --git a/deal.II/doc/doxygen/images/cubes_hole.png b/deal.II/doc/doxygen/images/cubes_hole.png new file mode 100644 index 0000000000..e0d8613bc2 Binary files /dev/null and b/deal.II/doc/doxygen/images/cubes_hole.png differ diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index 2e011d7927..c29b934ad1 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -1016,6 +1016,25 @@ inconvenience this causes.

    deal.II

      +
    1. New: Added function GridGenerator::hyper_cube_with_cylindrical_hole that produces + a square with a circular hole in the middle in 2d, and extrudes it along + the z-direction between 0 and L. +
      + (Luca Heltai 2007/02/15) +

      + +
    2. WorkAround: The class GridOut::write_msh produces a mesh which can be visualized in + the Gmsh reader. A bug in Gmsh prevented the boundary indicator to be + properly visualized. The boundary indicator was added to the material flag of the + faces (which is ignored when reading back the mesh) in order for the Gmsh reader + to be able to display the boundary indicator in a proper way. +
      + (Luca Heltai 2007/02/15) +

      + +
    3. Fixed: The function DoFTools::distribute_cell_to_dof_vector produced wrong results if the vector into which the result is to be