From 2b837b18705a75e09c0e3af00dfe7b95278e9c97 Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 4 Apr 2013 03:22:04 +0000
Subject: [PATCH] First part of discussing the use of curved boundaries.

git-svn-id: https://svn.dealii.org/trunk@29175 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/doc/doxygen/headers/glossary.h        |   5 +-
 deal.II/examples/step-49/doc/intro.dox        |  13 +-
 deal.II/examples/step-49/doc/results.dox      | 295 +++++++++++++++++-
 .../include/deal.II/grid/tria_boundary_lib.h  |   3 +-
 4 files changed, 311 insertions(+), 5 deletions(-)

diff --git a/deal.II/doc/doxygen/headers/glossary.h b/deal.II/doc/doxygen/headers/glossary.h
index 2b1c70e65f..688339494a 100644
--- a/deal.II/doc/doxygen/headers/glossary.h
+++ b/deal.II/doc/doxygen/headers/glossary.h
@@ -3,7 +3,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -339,7 +339,8 @@
  * object with a particular boundary indicator. This allows the Triangulation
  * object to use a different method of finding new points on faces and edges
  * to be refined; the default is to use a StraightBoundary object for all
- * faces and edges.
+ * faces and edges. The results section of step-49 has a worked example that
+ * shows all of this in action.
  *
  * The second use of boundary indicators is to describe not only which geometry
  * object to use on a particular boundary but to select a part of the boundary
diff --git a/deal.II/examples/step-49/doc/intro.dox b/deal.II/examples/step-49/doc/intro.dox
index 1a04b14a49..977ae59b9d 100644
--- a/deal.II/examples/step-49/doc/intro.dox
+++ b/deal.II/examples/step-49/doc/intro.dox
@@ -275,9 +275,20 @@ mesh, generated for example with gmsh, that is read in from a
 </table>
 
 
+<h3> After you have a coarse mesh </h3>
+
+Creating a coarse mesh using the methods discussed above is only the first
+step. When you have it, it will typically serve as the basis for further mesh
+refinement. This is not difficult &mdash; in fact, there is nothing else to do
+&mdash; if your geometry consists of only straight faces. However, this is
+often not the case if you have a more complex geometry and more steps than
+just creating the mesh are necessary. We will go over some of these steps in
+the <a href="#Results">results section</a> below.
+
+
 <!--
 
-<h4>Possible Extensions</h4>
+Possible Extensions
 
 - Modify a mesh:
   - change boundary indicators
diff --git a/deal.II/examples/step-49/doc/results.dox b/deal.II/examples/step-49/doc/results.dox
index cc5601ce56..2539e25399 100644
--- a/deal.II/examples/step-49/doc/results.dox
+++ b/deal.II/examples/step-49/doc/results.dox
@@ -1,3 +1,296 @@
 <h1>Results</h1>
 
-The program produces a series of .eps files of the Triangulations. The methods are discussed above.
+The program produces a series of <code>.eps</code> files of the
+Triangulations. The methods are discussed above.
+
+
+<h3>Next steps</h3>
+
+As mentioned in the introduction,
+creating a coarse mesh using the methods discussed here is only the first
+step. In order to refine a mesh, the Triangulation needs to know where to put
+new vertices on the mid-points of edges and faces. By default, these new
+points will be placed at the centers of the old edge but this isn't what you
+want if you need curved boundaries that aren't already adequately resolved by
+the coarse mesh. Several of the meshes shown in the introduction section fall
+into this category. For example, for this mesh the central hole is supposed to
+be round:
+
+<img
+ src="http://www.dealii.org/images/steps/developer/step-49.grid-2a.png"
+ alt="" height="200px">
+
+On the other hand, if you simply refine it, the Triangulation class can not
+know whether you wanted the hole to be round or to be an octagon. The default
+is to place new points along existing edges. After two mesh refinement steps,
+this would yield the following mesh, which is not what we wanted:
+
+<img
+ src="http://www.dealii.org/images/steps/developer/step-49.grid-2d-refined.png"
+ alt="" height="200px">
+
+What needs to happen is that you tell the triangulation that you in fact want
+to use a curved boundary. The way to do this requires three steps:
+- Create an object that describes the boundary in terms that allow the
+  Triangulation::execute_coarsening_and_refinement function to ask the
+  boundary description where a new point should be located upon mesh
+  refinement.
+- Tell the triangulation object that you want this object to be used for all
+  boundaries with boundary indicates equal to a particular value (for more
+  information on boundary indicators, see the
+  @ref GlossBoundaryIndicator "glossary entry on this topic".)
+- Mark those parts of the boundary of the domain for which you want the
+  boundary to be so treated with the value of the boundary indicator used in
+  the previous step. (The order of this step and the previous one does not
+  matter.)
+
+To illustrate this process in more detail, let us consider an example created
+by Yuhan Zhou as part of a semester project at Texas A&amp;M University in
+2013. The goal was to generate (and use) a geometry that describes a
+microstructured electric device. In a CAD program, the geometry looks like
+this:
+
+<img
+ src="http://www.dealii.org/images/steps/developer/step-49.yuhan.1.png"
+ alt="" height="200px">
+
+The first step in getting there was to create a coarse mesh, which was done
+by creating a 2d coarse mesh for each of the two cross section, extruding them
+into the third direction, and gluing them together. The following code does
+this, using the techniques previously described:
+
+@code
+// Given a list of points and how vertices connect to cells,
+// create a mesh. This is in the same way as we do in step-14.
+void create_2d_grid (const Point<2> vertices_1[],
+	  	     const unsigned int n_vertices,
+		     const int cell_vertices[][4],
+		     const unsigned int n_cells,
+                     Triangulation<2> &coarse_grid)
+{
+  const std::vector<Point<2> > vertices (&vertices_1[0],
+			    	         &vertices_1[n_vertices]);
+
+  std::vector<CellData<2> > cells (n_cells, CellData<2>());
+  for (unsigned int i=0; i<n_cells; ++i)
+    {
+      for (unsigned int j=0;
+	   j<GeometryInfo<dimm>::vertices_per_cell;
+	   ++j)
+	cells[i].vertices[j] = cell_vertices[i][j];
+    }
+
+  coarse_grid.create_triangulation (vertices,
+				    cells,
+				    SubCellData());
+}
+
+
+// Create a triangulation that covers the entire volume
+void create_3d_grid (Triangulation<3> &triangulation)
+{
+  // Generate first cross section
+  const Point<2> vertices_1[]
+    = {  Point<2> (-1.5,   0.),
+	 Point<2> (-0.5,   0.),
+	 Point<2> ( 0.5,   0.),
+	 Point<2> ( 1.5,   0.),
+
+	 Point<2> (-1.5,  1.5),
+	 Point<2> (-0.5,  1.5),
+	 Point<2> ( 0.5,  1.5),
+	 Point<2> ( 1.5,  1.5),
+
+	 Point<2> (-1.5,   3.),
+	 Point<2> (-0.5,   3.),
+	 Point<2> ( 0.5,   3.),
+	 Point<2> ( 1.5,   3.),
+
+	 Point<2> (-0.5,   3+0.5*sqrt(3)),
+	 Point<2> ( 0.5,   3+0.5*sqrt(3)),
+
+	 Point<2> (-0.75,  3+0.75*sqrt(3)),
+	 Point<2> ( 0.75,  3+0.75*sqrt(3))
+  };
+  const int cell_vertices_1[][GeometryInfo<2>::vertices_per_cell]
+    = {{0, 1, 4, 5},
+       {1, 2, 5, 6},
+       {3, 7, 2, 6},
+       {4, 5, 8, 9},
+       {5, 6, 9, 10},
+       {7,11, 6,10},
+       {8, 9, 14,12},
+       {9, 10,12,13},
+       {11,15,10,13},
+       {14,12,15,13}
+  };
+
+  // Copy vertices into a 2d triangulation
+  Triangulation<dim-1> triangulation_2d_1;
+  create_2d_grid (vertices_1,
+                  sizeof(vertices_1)/sizeof(vertices_1[0]),
+		  cell_vertices_1,
+                  sizeof(cell_vertices_1)/sizeof(cell_vertices_1[0]),
+                  triangulation_2d_1);
+
+  // Then extrude it into a 3d piece
+  Triangulation<dim> triangulation_3d_1;
+  GridGenerator::extrude_triangulation (triangulation_2d_1,
+                                        5, 2.5,
+                                        triangulation_3d_1);
+
+  // Now do the same with the second volume
+  const Point<2> vertices_2[]
+    = {  Point<2> (-2.5,   0.),
+	 Point<2> (-1.5,   0.),
+	 Point<2> (-0.5,   0.),
+	 Point<2> ( 0.5,   0.),
+	 Point<2> ( 1.5,   0.),
+	 Point<2> ( 2.5,   0.),
+
+	 Point<2> (-2.5,  1.5),
+	 Point<2> (-1.5,  1.5),
+	 Point<2> (-0.5,  1.5),
+	 Point<2> ( 0.5,  1.5),
+	 Point<2> ( 1.5,  1.5),
+	 Point<2> ( 2.5,  1.5),
+
+	 Point<2> (-2.5,  3.),
+	 Point<2> (-1.5,  3.),
+	 Point<2> (-0.5,  3.),
+	 Point<2> ( 0.5,  3.),
+	 Point<2> ( 1.5,  3.),
+	 Point<2> ( 2.5,  3.),
+
+	 Point<2> (-0.5,   3.+0.5*sqrt(3)),
+	 Point<2> ( 0.5,   3.+0.5*sqrt(3)),
+
+	 Point<2> (-0.75,  3.+0.75*sqrt(3)),
+	 Point<2> ( 0.75,  3.+0.75*sqrt(3)),
+
+	 Point<2> (-1.25,  3.+1.25*sqrt(3)),
+	 Point<2> ( 1.25,  3.+1.25*sqrt(3))
+  };
+  const int cell_vertices_2[][GeometryInfo<2>::vertices_per_cell]
+    = {{0, 1, 6, 7},
+       {1, 2, 7, 8},
+       {2, 3, 8, 9},
+       {4, 10, 3, 9},
+       {5, 11, 4, 10},
+       {6, 7, 12, 13},
+       {7, 8, 13, 14},
+       {8, 9, 14, 15},
+       {10, 16, 9, 15},
+       {11, 17, 10, 16},
+       {12, 13, 22, 20},
+       {13, 14, 20, 18},
+       {14, 15, 18, 19},
+       {16, 21, 15, 19},
+       {17, 23, 16, 21},
+       {20, 18, 21, 19},
+       {22, 20, 23, 21}
+  };
+
+  Triangulation<dim-1> triangulation_2d_2;
+  create_2d_grid (vertices_2,
+                  sizeof(vertices_2)/sizeof(vertices_2[0]),
+		  cell_vertices_2,
+                  sizeof(cell_vertices_2)/sizeof(cell_vertices_2[0]),
+                  triangulation_2d_2);
+
+  Triangulation<dim> triangulation_3d_2;
+  GridGenerator::extrude_triangulation (triangulation_2d_2,
+                                        5, 2.5,
+                                        triangulation_3d_2);
+
+  // Also shift this triangulation in the z-direction so
+  // that it matches the end face of the first part
+  GridTools::shift (Point<3>(0,0,2.5),
+                    triangulation_3d_2);
+
+
+  // Now first merge these two pieces, then shift the
+  // first piece in z-direction beyond the second, and
+  // merge the shifted piece with the two previously
+  // merged one into the final one:
+  Triangulation<dim> triangulation_3d_tmp;
+  GridGenerator::merge_triangulations (triangulation_3d_1,
+                                       triangulation_3d_2,
+                                       triangulation_3d_tmp);
+
+  GridTools::shift (Point<3>(0,0,5),
+                    triangulation_3d_1);
+
+  GridGenerator::merge_triangulations (triangulation_3d_tmp,
+                                       triangulation_3d_1,
+                                       triangulation);
+@endcode
+
+With this code, you get a mesh that looks like this:
+
+<img
+ src="http://www.dealii.org/images/steps/developer/step-49.yuhan.2.png"
+ alt="" height="200px">
+
+The next step is to teach each of the top surfaces that they should be
+curved. We can do this by creating CylinderBoundary objects that
+describe this, where after some playing with the arguments one finds the
+correct values for the axes, offsets, and radii:
+
+@code
+  const double inner_radius = 1.5;
+  const double outer_radius = 2.5;
+
+  typename Triangulation<dim>::active_cell_iterator
+    cell = triangulation.begin_active(),
+    endc = triangulation.end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int f=0;
+	 f < GeometryInfo<dim>::faces_per_cell;
+	 ++f)
+      {
+	const Point<dim> face_center = cell->face(f)->center();
+
+	if (cell->face(f)->at_boundary())
+	  {
+	    const double dist = sqrt(pow(face_center[1]-3,2)+pow(face_center[0],2));
+
+	    if((face_center[2] <= 2.5 || face_center[2] >= 5) && face_center[1] >= 3 && dist <= (inner_radius+outer_radius)/2)
+	      cell->face(f)->set_boundary_indicator(8);
+
+	    if(face_center[2] >= 2.5 && face_center[2] <= 5 && face_center[1] >= 3 && dist >= (inner_radius+outer_radius)/2)
+	      cell->face(f)->set_boundary_indicator(9);
+	  }
+      }
+  static const CylinderBoundary<dim> inner_cylinder(inner_radius,2);
+  static const CylinderBoundary<dim> outer_cylinder(outer_radius,2);
+  triangulation.set_boundary (8, inner_cylinder);
+  triangulation.set_boundary (9, outer_cylinder);
+
+
+  for (typename Triangulation<dim>::active_face_iterator
+	 face=triangulation.begin_active_face();
+       face!=triangulation.end_face(); ++face)
+    if (face->at_boundary())
+      if ((face->boundary_indicator() == 8)
+	  ||
+	  (face->boundary_indicator() == 9))
+	for (unsigned int edge = 0; edge<GeometryInfo<dim>::lines_per_face;
+	     ++edge)
+	  face->line(edge)
+	    ->set_boundary_indicator (face->boundary_indicator());
+
+
+
+  //std::string filename = "grid-";
+  //filename += ('0');
+  //filename += ".eps";
+  //std::ofstream output (filename.c_str());
+
+  //GridOut grid_out;
+  //grid_out.write_eps (triangulation_2d_1, output);
+
+
+}
+
+@endcode
diff --git a/deal.II/include/deal.II/grid/tria_boundary_lib.h b/deal.II/include/deal.II/grid/tria_boundary_lib.h
index c89a796bc4..b482962bf1 100644
--- a/deal.II/include/deal.II/grid/tria_boundary_lib.h
+++ b/deal.II/include/deal.II/grid/tria_boundary_lib.h
@@ -31,7 +31,8 @@ DEAL_II_NAMESPACE_OPEN
  *
  * This class was developed to be used in conjunction with the
  * @p cylinder function of GridGenerator. It should be used for
- * the hull of the cylinder only (boundary indicator 0).
+ * the hull of the cylinder only (boundary indicator 0). Its use is
+ * discussed in detail in the results section of step-49.
  *
  * This class is derived from StraightBoundary rather than from
  * Boundary, which would seem natural, since this way we can use the
-- 
2.39.5