]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Workaround for gmsh output mesh.
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 15 Feb 2007 19:16:08 +0000 (19:16 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 15 Feb 2007 19:16:08 +0000 (19:16 +0000)
New funcition:
GridGenerator::hyper_cube_with_cylindrical_hole.

git-svn-id: https://svn.dealii.org/trunk@14477 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_generator.h
deal.II/deal.II/source/grid/grid_generator.cc
deal.II/deal.II/source/grid/grid_out.cc
deal.II/doc/doxygen/images/cubes_hole.png [new file with mode: 0644]
deal.II/doc/news/changes.html

index ad85d444afb9211424b83ea41a7d21708a14e2a4..eb70ac94b286f61b862d23b5a7198c128b2e2ce3 100644 (file)
@@ -109,12 +109,24 @@ template <typename number> class SparseMatrix;
  *      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
@@ -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<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
index 803b82ee8c27a4d5f7e6543179e99c034534b922..716c97b841c3b4aefb17ac7f6af5bcdbc74bfc14 100644 (file)
@@ -2245,6 +2245,207 @@ void GridGenerator::laplace_transformation (Triangulation<dim> &tria,
 #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> (
@@ -2324,6 +2525,11 @@ GridGenerator::half_hyper_shell (
   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::
index 91c18d21a17e2fff62962f54f4b6e061ccdd128e..6dd51f4bf993309e9e21ff8124902827e3f16d6f 100644 (file)
@@ -743,10 +743,13 @@ void GridOut::write_msh_faces (const Triangulation<dim> &tria,
                  Assert (false, ExcNotImplemented());
          };
        out << static_cast<unsigned int>(face->boundary_indicator()) 
-           << " 1 " << GeometryInfo<dim>::vertices_per_face;
+           << " " 
+           << static_cast<unsigned int>(face->boundary_indicator()) 
+           << " " << GeometryInfo<dim>::vertices_per_face;
                                         // note: vertex numbers are 1-base
        for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
-         out << ' ' << face->vertex_index(GeometryInfo<dim-1>::ucd_to_deal[vertex])+1;
+         out << ' ' 
+             << face->vertex_index(GeometryInfo<dim-1>::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 (file)
index 0000000..e0d8613
Binary files /dev/null and b/deal.II/doc/doxygen/images/cubes_hole.png differ
index 2e011d7927a1c6005c63cf099a4ef57bfc3e2279..c29b934ad1b5dc1c24cd7de45b0319a27e64494a 100644 (file)
@@ -1016,6 +1016,25 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+   <li> <p>New: Added function <code class="class">GridGenerator</code>::<code 
+       class="member">hyper_cube_with_cylindrical_hole</code> that produces 
+          a square with a circular hole in the middle in 2d, and extrudes it along
+          the z-direction between 0 and L. 
+       <br>
+       (Luca Heltai 2007/02/15)
+       </p>
+
+   <li> <p>WorkAround: The class <code class="class">GridOut</code>::<code 
+       class="member">write_msh</code> 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.
+       <br>
+       (Luca Heltai 2007/02/15)
+       </p>
+
+
    <li> <p>Fixed: The function <code class="class">DoFTools</code>::<code 
        class="member">distribute_cell_to_dof_vector</code> produced
        wrong results if the vector into which the result is to be

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.