From 8484eaf5fc2822a50e95816d77cd6bbade35dc77 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 25 Feb 2013 04:21:05 +0000 Subject: [PATCH] Add more documentation. git-svn-id: https://svn.dealii.org/trunk@28543 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-49/step-49.cc | 198 ++++++++++++++++++++-------- 1 file changed, 145 insertions(+), 53 deletions(-) diff --git a/deal.II/examples/step-49/step-49.cc b/deal.II/examples/step-49/step-49.cc index e8b9937687..499aff5b37 100644 --- a/deal.II/examples/step-49/step-49.cc +++ b/deal.II/examples/step-49/step-49.cc @@ -1,4 +1,4 @@ -/* $Id: $ +/* $Id$ * * Copyright (C) 2013 by the deal.II authors * @@ -8,6 +8,15 @@ * further information on this license. */ + +// This tutorial program is odd in the sense that, unlike for most other +// steps, the introduction already provides most of the information on how to +// use the various strategies to generate meshes. Consequently, there is +// little that remains to be commented on here, and we intersperse the code +// with relatively little text. In essence, the code here simply provides a +// reference implementation of what has already been described in the +// introduction. + // @sect3{Include files} #include @@ -25,17 +34,29 @@ using namespace dealii; -// @sect3{Generate some output for a given mesh} - +// @sect3{Generating output for a given mesh} + +// The following function generates some output for any of the meshes we will +// be generating in the remainder of this program. In particular, it generates +// the following information: +// +// - Some general information about the number of space dimensions in which +// this mesh lives and its number of cells. +// - The number of boundary faces that use each boundary indicator, so that +// it can be compared with what we expect. +// +// Finally, the function outputs the mesh in encapsulated postscript (EPS) +// format that can easily be visualized in the same way as was done in step-1. template -void mesh_info(Triangulation &tria, const char *filename) +void mesh_info(const Triangulation &tria, + const std::string &filename) { - // some general info: std::cout << "Mesh info:" << std::endl << " dimension: " << dim << std::endl << " no. of cells: " << tria.n_active_cells() << std::endl; - // loop over all the cells and find how often each boundary indicator is used: + // Next loop over all faces of all cells and find how often each boundary + // indicator is used: { std::map boundary_count; typename Triangulation::active_cell_iterator @@ -60,9 +81,9 @@ void mesh_info(Triangulation &tria, const char *filename) std::cout << std::endl; } - // Now we want to write a graphical representation of the mesh to an output - // file. - std::ofstream out (filename); + // Finally, produce a graphical representation of the mesh to an output + // file: + std::ofstream out (filename.c_str()); GridOut grid_out; grid_out.write_eps (tria, out); std::cout << " written to " << filename << std::endl; @@ -70,8 +91,12 @@ void mesh_info(Triangulation &tria, const char *filename) // @sect3{Main routines} +// @sect4{grid_1: Loading a mesh generated by gmsh} -// load from a msh generated by gmsh +// In this first example, we show how to load the mesh for which we have +// discussed in the introduction how to generate it. This follows the same +// pattern as used in step-5 to load a mesh, although there it was written in +// a different file format (UCD instead of MSH). void grid_1 () { Triangulation<2> triangulation; @@ -86,28 +111,52 @@ void grid_1 () } -// merge triangulations -// note: the vertices have to be exactly on top of each other! +// @sect4{grid_2: Merging triangulations} + +// Here, we first create two triangulations and then merge them into one. As +// discussed in the introduction, it is important to ensure that the vertices +// at the common interface are located at the same coordinates. void grid_2 () { - Triangulation<2> triangulation; - - Triangulation<2> tria1, tria2; - - + Triangulation<2> tria1; GridGenerator::hyper_cube_with_cylindrical_hole (tria1, 0.25, 1.0); + + Triangulation<2> tria2; std::vector< unsigned int > repetitions(2); repetitions[0]=3; repetitions[1]=2; GridGenerator::subdivided_hyper_rectangle (tria2, repetitions, - Point<2>(1.0,-1.0), Point<2>(4.0,1.0)); + Point<2>(1.0,-1.0), + Point<2>(4.0,1.0)); + Triangulation<2> triangulation; GridGenerator::merge_triangulations (tria1, tria2, triangulation); mesh_info(triangulation, "grid-2.eps"); } -// move vertices + +// @sect4{grid_3: Moving vertices} + +// In this function, we move vertices of a mesh. This is simpler than one +// usually expects: if you ask a cell using cell-@>vertex(i) for +// the coordinates of its ith vertex, it doesn't just provide the +// location of this vertex but in fact a reference to the location where these +// coordinates are stored. We can then modify the value stored there. +// +// So this is what we do in the first part of this function: We create a +// square of geometry $[-1,1]^2$ with a circular hole with radius 0.25 located +// at the origin. We then loop over all cells and all vertices and if a vertex +// has a $y$ coordinate equal to one, we move it upward by 0.5. +// +// Note that this sort of procedure does not usually work this way because one +// will typically encounter the same vertices multiple times and may move them +// more than once. It works here because we select the vertices we want to use +// based on their geometric location, and a vertex moved once will fail this +// test in the future. A more general approach to this problem would have been +// to keep a std::set of of those vertex indices that we have already moved +// (which we can obtain using cell-@>vertex_index(i) and only +// move those vertices whose index isn't in the set yet. void grid_3 () { Triangulation<2> triangulation; @@ -122,74 +171,110 @@ void grid_3 () { Point<2> &v = cell->vertex(i); if (std::abs(v(1)-1.0)<1e-5) - v(1)+=0.5; + v(1) += 0.5; } - } - // here we are going to set a boundary descriptor for the round hole - // and refine the mesh twice. - const HyperBallBoundary<2> boundary_description(Point<2>(0,0),0.25); + // In the second step we will refine the mesh twice. To do this correctly, + // we have to associate a geometry object with the boundary of the hole; + // since the boundary of the hole has boundary indicator 1 (see the + // documentation of the function that generates the mesh), we need to create + // an object that describes a circle (i.e., a hyper ball) with appropriate + // center and radius and assign it to the triangulation. We can then refine + // twice: + const HyperBallBoundary<2> boundary_description(Point<2>(0,0), 0.25); triangulation.set_boundary (1, boundary_description); triangulation.refine_global(2); - mesh_info(triangulation, "grid-3.eps"); - - // remove boundary object from Triangulation again: + // The mesh so generated is then passed to the function that generates + // output. In a final step we remove the boundary object again so that it is + // no longer in use by the triangulation when it is destroyed (the boundary + // object is destroyed first in this function since it was declared after + // the triangulation). + mesh_info (triangulation, "grid-3.eps"); triangulation.set_boundary (1); } -// demonstrate extrude_triangulation + +// @sect4{grid_4: Demonstrating extrude_triangulation} + +// This example takes the initial grid from the previous function and simply extrudes it into the third space dimension: void grid_4() { Triangulation<2> triangulation; Triangulation<3> out; GridGenerator::hyper_cube_with_cylindrical_hole (triangulation, 0.25, 1.0); - GridGenerator::extrude_triangulation(triangulation, 3, 2.0, out); + GridGenerator::extrude_triangulation (triangulation, 3, 2.0, out); mesh_info(out, "grid-4.eps"); - + } + +// @sect4{grid_5: Demonstrating GridTools::transform, part 1} + +// This and the next example first create a mesh and then transform it by +// moving every node of the mesh according to a function that takes a point +// and returns a mapped point. In this case, we transform $(x,y) \mapsto +// (x,y+\sin(\pi x/5))$. +// +// GridTools::transform takes a triangulation and any kind of object that can +// be called like a function as arguments. This function-like argument can be +// simply the address of a function as in the current case, or an object that +// has an operator() as in the next example, or for example a +// std::function@(const Point@<2@>) object one can get +// via std::bind in more complex cases. struct Grid5Func { Point<2> operator() (const Point<2> & in) const { - return Point<2>(in(0), in(1)+sin(in(0)/5.0*3.14159)); - } + return Point<2>(in(0), + in(1) + std::sin(in(0)/5.0*3.14159)); + } }; -// demonstrate GridTools::transform void grid_5() { Triangulation<2> tria; - std::vector< unsigned int > repetitions(2); - repetitions[0]=14; - repetitions[1]=2; + std::vector repetitions(2); + repetitions[0] = 14; + repetitions[1] = 2; GridGenerator::subdivided_hyper_rectangle (tria, repetitions, - Point<2>(0.0,0.0), Point<2>(10.0,1.0)); + Point<2>(0.0,0.0), + Point<2>(10.0,1.0)); - GridTools::transform(Grid5Func(), tria); - mesh_info(tria, "grid-5.eps"); + mesh_info(tria, "grid-5.eps"); } + + +// @sect4{grid_6: Demonstrating GridTools::transform, part 2} + +// In this second example of transforming points from an original to a new +// mesh, we will use the mapping $(x,y) \mapsto (x,\tanh(2y)/\tanh(2))$. To +// make things more interesting, rather than doing so in a single function as +// in the previous example, we here create an object with an +// operator() that will be called by GridTools::transform. Of +// course, this object may in reality be much more complex: the object may +// have member variables that play a role in computing the new locations of +// vertices. struct Grid6Func { - double trans(double x) const + double trans(const double x) const { - return tanh(2*x)/tanh(2); + return std::tanh(2*x)/tanh(2); } - + Point<2> operator() (const Point<2> & in) const { - return Point<2>((in(0)), trans(in(1))); - } + return Point<2> (in(0), + trans(in(1))); + } }; - -// demonstrate GridTools::transform + void grid_6() { Triangulation<2> tria; @@ -200,25 +285,32 @@ void grid_6() Point<2>(0.0,0.0), Point<2>(1.0,1.0)); GridTools::transform(Grid6Func(), tria); - mesh_info(tria, "grid-6.eps"); + mesh_info(tria, "grid-6.eps"); } +// @sect4{grid_7: Demonstrating distort_random} -//demonstrate distort_random +// In this last example, we create a mesh and then distort its (interior) +// vertices by a random perturbation. This is not something you want to do for +// production computations, but it is a useful tool for testing +// discretizations and codes to make sure they don't work just by accident +// because the mesh happens to be uniformly structured and supporting +// super-convergence properties. void grid_7() { Triangulation<2> tria; - std::vector< unsigned int > repetitions(2); - repetitions[0]=16; - repetitions[1]=16; + std::vector repetitions(2); + repetitions[0] = repetitions[1] = 16; GridGenerator::subdivided_hyper_rectangle (tria, repetitions, - Point<2>(0.0,0.0), Point<2>(1.0,1.0)); + Point<2>(0.0,0.0), + Point<2>(1.0,1.0)); tria.distort_random(0.3, true); - mesh_info(tria, "grid-7.eps"); + mesh_info(tria, "grid-7.eps"); } + // @sect3{The main function} // Finally, the main function. There isn't much to do here, only to call the -- 2.39.5