Reference documentation for deal.II version Git 4a360c6 2017-05-26 05:30:59 +0200
The step-49 tutorial program

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program
This program was contributed by Timo Heister. Parts of the results section were contributed by Yuhan Zhou and Wolfgang Bangerth.

Introduction

This tutorial is an extension to step-1 and demonstrates several ways to obtain more involved meshes than the ones shown there.

Generating complex geometries is a challenging task, especially in three space dimensions. We will discuss several ways to do this, but this list is not exhaustive. Additionally, there is not one approach that fits all problems.

This example program shows some of ways to create and modify meshes for computations and outputs them as .eps files in much the same way as we do in step-1. No other computations or adaptive refinements are done; the idea is that you can use the techniques used here as building blocks in other, more involved simulators. Please note that the example program does not show all the ways to generate meshes that are discussed in this introduction.

General concerns about meshes

When you use adaptive mesh refinement, you definitely want the initial mesh to be as coarse as possible. The reason is that you can make it as fine as you want using adaptive refinement as long as you have memory and CPU time available. However, this requires that you don't waste mesh cells in parts of the domain where they don't pay off. As a consequence, you don't want to start with a mesh that is too fine to start with, because that takes up a good part of your cell budget already, and because you can't coarsen away cells that are in the initial mesh.

That said, your mesh needs to capture the given geometry adequately.

How to create meshes

There are several ways to create an initial mesh. Meshes can be modified or combined in many ways as discussed later on.

Using GridGenerator

The easiest way to generate meshes is to use the functions in namespace GridGenerator, as already discussed in step-1. There are many different helper functions available, including GridGenerator::hyper_cube(), GridGenerator::hyper_shell(), GridGenerator::hyper_ball(), GridGenerator::hyper_cube_with_cylindrical_hole(), ...

Constructing your own mesh programmatically

If there is no good fit in the GridGenerator namespace for what you want to do, you can always create a Triangulation in your program "by hand". For that, you need a list of vertices with their coordinates and a list of cells referencing those vertices. You can find an example in the function create_coarse_grid in step-14. All the functions in GridGenerator are implemented in this fashion.

We are happy to accept more functions to be added to GridGenerator. So, if you end up writing a function that might be useful for a larger audience, please contribute it.

Importing from external programs

The class GridIn can read many different mesh formats from a file from disk. How this is done is explained in step-5 and can be seen in the function grid_1 in this example, see the code below.

Meshes can be generated from different tools like gmsh, lagrit and cubit. See the documentation of GridIn for more information. The problem is that deal.II needs meshes that only consist of quads and hexas – tetrahedral meshes won't work (this means tools like tetgen can not be used directly).

We will describe a possible workflow using Gmsh. Gmsh is the smallest and most quickly set up open source tool we are aware of. It can generate unstructured 2d quad meshes, but in 3d it can only extrude 2d meshes to get hexahedral meshes. 3D meshing of unstructured geometry into hexahedra is not supported at the time of writing this tutorial (early 2013).

In gmsh, a mesh is described in a text based .geo file, that can contain computations, loops, variables, etc. It is very flexible. The mesh is generated from a surface representation, which is build from a list of line loops, which is build from a list of lines, which are in turn built from points. The .geo script can be written and edited by hand or it can be generated automatically by creating objects graphically inside gmsh. In many cases it is best to combine both approaches. The file can be easily reloaded by pressing "reload" under the "Geometry" tab.

This tutorial contains an example .geo file, that describes a box with two objects cut out in the interior. This is how untitled.geo looks like in gmsh (displaying the boundary indicators as well as the mesh discussed further down below):

You might want to open the untitled.geo file in a text editor (it is located in the same directory as the step-49.cc source file) to see how it is structured. You can see how the boundary of the domain is composed of a number of lines and how later on we combine several lines into "physical lines" (or "physical surfaces") that list the logical lines' numbers. "Physical" object are the ones that carry information about the boundary indicator (see this glossary entry).

Note
It is important that this file contain "physical lines" and "physical surfaces". These give the boundary indicators and material ids for use in deal.II. Without these physical entities, nothing will be imported into deal.II.

deal.II's GridIn class can read the .msh format written by gmsh and that contains a mesh created for the geometry described by the .geo file. You generate the .msh from the .geo by running the commands

gmsh -2 untitled.geo

on the command line, or by clicking "Mesh" and then "2D" inside Gmsh after loading the file. Now this is the mesh read from the .msh file and saved again by deal.II as an image (see the grid_1 function of the current program):

Modifying a Mesh

After acquiring one (or several) meshes in the ways described above, there are many ways to manipulate them before using them in a finite element computation.

Transformations

The GridTools namespace contains a collection of small functions to transform a given mesh in various ways. The usage of the functions GridTools::shift, GridTools::rotate, GridTools::scale is fairly obvious, so we won't discuss those functions here.

The function GridTools::transform allows you to transform the vertices of a given mesh using a smooth function. An example of its use is also given in the results section of step-38 but let us show a simpler example here: In the function grid_5() of the current program, we perturb the y coordinate of a mesh with a sine curve:

regular input mesh
output mesh

Similarly, we can transform a regularly refined unit square to a wall-adapted mesh in y direction using the formula \((x,y) \mapsto (x,\tanh(2*y)/\tanh(2))\). This is done in grid_6() of this tutorial:

regular input mesh
wall-adapted output mesh

Finally, the function GridTools::distort_random allows you to move vertices in the mesh (optionally ignoring boundary nodes) by a random amount. This is demonstrated in grid_7() and the result is as follows:

regular input mesh
perturbed output mesh

This function is primarily intended to negate some of the superconvergence effects one gets when studying convergence on regular meshes, as well as to suppress some optimizations in deal.II that can exploit the fact that cells are similar in shape. In practice, it is of course always better to work with a sequence of unstructured meshes (see possible extensions at the end of the this section).

Merging Meshes

The function GridGenerator::merge_triangulations() allows you to merge two given Triangulation objects into a single one. For this to work, the vertices of the shared edge or face have to match exactly. Lining up the two meshes can be achieved using GridTools::shift and GridTools::scale. In the function grid_2() of this tutorial, we merge a square with a round hole (generated with GridGenerator::hyper_cube_with_cylindrical_hole()) and a rectangle (generated with GridGenerator::subdivided_hyper_rectangle()). The function GridGenerator::subdivided_hyper_rectangle() allows you to specify the number of repetitions and the positions of the corners, so there is no need to shift the triangulation manually here. You should inspect the mesh graphically to make sure that cells line up correctly and no unpaired nodes exist in the merged Triangulation.

These are the input meshes and the output mesh:

input mesh 1
input mesh 2
merged mesh

Moving Vertices

The function grid_3() demonstrates the ability to pick individual vertices and move them around in an existing mesh. Note that this has the potential to produce degenerate or inverted cells and you shouldn't expect anything useful to come of using such meshes. Here, we create a box with a cylindrical hole that is not exactly centered by moving the top vertices upwards:

input mesh

top vertices moved upwards

For the exact way how this is done, see the code below.

Extruding Meshes

If you need a 3d mesh that can be created by extruding a given 2d mesh (that can be created in any of the ways given above), you can use the function GridGenerator::extrude_triangulation(). See the grid_4() function in this tutorial for an example. Note that for this particular case, the given result could also be achieved using the 3d version of GridGenerator::hyper_cube_with_cylindrical_hole(). The main usage is a 2d mesh, generated for example with gmsh, that is read in from a .msh file as described above. This is the output from grid_4():

input mesh

extruded output mesh

After you have a coarse mesh

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 — in fact, there is nothing else to do — 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 results section below.

The commented program

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.

Include files

#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <iostream>
#include <fstream>
#include <map>
using namespace dealii;

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:

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 <int dim>
void print_mesh_info(const Triangulation<dim> &triangulation,
const std::string &filename)
{
std::cout << "Mesh info:" << std::endl
<< " dimension: " << dim << std::endl
<< " no. of cells: " << triangulation.n_active_cells() << std::endl;

Next loop over all faces of all cells and find how often each boundary indicator is used (recall that if you access an element of a std::map object that doesn't exist, it is implicitly created and default initialized – to zero, in the current case – before we then increment it):

{
std::map<unsigned int, unsigned int> boundary_count;
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
{
if (cell->face(face)->at_boundary())
boundary_count[cell->face(face)->boundary_id()]++;
}
}
std::cout << " boundary indicators: ";
for (std::map<unsigned int, unsigned int>::iterator it=boundary_count.begin();
it!=boundary_count.end();
++it)
{
std::cout << it->first << "(" << it->second << " times) ";
}
std::cout << std::endl;
}

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 (triangulation, out);
std::cout << " written to " << filename
<< std::endl
<< std::endl;
}

Main routines

grid_1: Loading a mesh 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;
GridIn<2> gridin;
gridin.attach_triangulation(triangulation);
std::ifstream f("untitled.msh");
gridin.read_msh(f);
print_mesh_info (triangulation, "grid-1.eps");
}

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 ()
{
std::vector< unsigned int > repetitions(2);
repetitions[0]=3;
repetitions[1]=2;
Point<2>(1.0,-1.0),
Point<2>(4.0,1.0));
Triangulation<2> triangulation;
GridGenerator::merge_triangulations (tria1, tria2, triangulation);
print_mesh_info (triangulation, "grid-2.eps");
}

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;
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
{
Point<2> &v = cell->vertex(i);
if (std::abs(v(1)-1.0)<1e-5)
v(1) += 0.5;
}
}

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 spherical manifold (i.e., a hyper ball) with appropriate center and assign it to the triangulation. Notice that the function that generates the triangulation sets the boundary indicators of the inner mesh, but leaves unchanged the manifold indicator. We copy the boundary indicator to the manifold indicators in order for the object to be refined accordingly. We can then refine twice:

const SphericalManifold<2> boundary_description(Point<2>(0,0));
triangulation.set_manifold (1, boundary_description);
triangulation.refine_global(2);

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).

print_mesh_info (triangulation, "grid-3.eps");
triangulation.set_manifold (1);
}

There is one snag to doing things as shown above: If one moves the nodes on the boundary as shown here, one often ends up with cells in the interior that are badly distorted since the interior nodes were not moved around. This is not that much of a problem in the current case since the mesh did not contain any internal nodes when the nodes were moved – it was the coarse mesh and it so happened that all vertices are at the boundary. It's also the case that the movement we had here was, compared to the average cell size not overly dramatic. Nevertheless, sometimes one does want to move vertices by a significant distance, and in that case one needs to move internal nodes as well. One way to do that automatically is to call the function GridTools::laplace_transform that takes a set of transformed vertex coordinates and moves all of the other vertices in such a way that the resulting mesh has, in some sense, a small distortion.

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;
GridGenerator::extrude_triangulation (triangulation, 3, 2.0, out);
print_mesh_info (out, "grid-4.eps");
}

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<Point<2>(const Point<2>)> object one can get via std::bind in more complex cases.

Point<2> grid_5_transform (const Point<2> &in)
{
return Point<2>(in(0),
in(1) + std::sin(in(0)/5.0*3.14159));
}
void grid_5()
{
Triangulation<2> triangulation;
std::vector<unsigned int> repetitions(2);
repetitions[0] = 14;
repetitions[1] = 2;
repetitions,
Point<2>(0.0,0.0),
Point<2>(10.0,1.0));
GridTools::transform (&grid_5_transform, triangulation);
print_mesh_info (triangulation, "grid-5.eps");
}

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(const double y) const
{
return std::tanh(2*y)/tanh(2);
}
Point<2> operator() (const Point<2> &in) const
{
return Point<2> (in(0),
trans(in(1)));
}
};
void grid_6()
{
Triangulation<2> triangulation;
std::vector< unsigned int > repetitions(2);
repetitions[0] = repetitions[1] = 40;
repetitions,
Point<2>(0.0,0.0),
Point<2>(1.0,1.0));
GridTools::transform(Grid6Func(), triangulation);
print_mesh_info (triangulation, "grid-6.eps");
}

grid_7: Demonstrating 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> triangulation;
std::vector<unsigned int> repetitions(2);
repetitions[0] = repetitions[1] = 16;
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions,
Point<2>(0.0,0.0),
Point<2>(1.0,1.0));
GridTools::distort_random (0.3, triangulation, true);
print_mesh_info (triangulation, "grid-7.eps");
}

The main function

Finally, the main function. There isn't much to do here, only to call the subfunctions.

int main ()
{
try
{
grid_1 ();
grid_2 ();
grid_3 ();
grid_4 ();
grid_5 ();
grid_6 ();
grid_7 ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
}

Results

The program produces a series of .eps files of the triangulations. The methods are discussed above.

Next steps: Curved boundaries

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:

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:

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:

To illustrate this process in more detail, let us consider an example created by Yuhan Zhou as part of a 2013 semester project at Texas A&M University. The goal was to generate (and use) a geometry that describes a microstructured electric device. In a CAD program, the geometry looks like this:

In the following, we will walk you through the entire process of creating a mesh for this geometry, including a number of common pitfalls by showing the things that can go wrong.

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:

// 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<2>::vertices_per_cell;
++j)
cells[i].vertices[j] = cell_vertices[i][j];
}
coarse_grid.create_triangulation (vertices,
cells,
}
// 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<2> 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<3> triangulation_3d_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<2> 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<3> triangulation_3d_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
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<3> triangulation_3d_tmp;
triangulation_3d_2,
triangulation_3d_tmp);
triangulation_3d_1);
GridGenerator::merge_triangulations (triangulation_3d_tmp,
triangulation_3d_1,
triangulation);
}

With this code, you get a mesh that looks like this:

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. A first attempt looks like this:

Triangulation<3> triangulation;
create_3d_grid (triangulation);
// Create the objects that describe the boundaries and attach them
// to the triangulation as the ones to use on boundaries marked
// with boundary indicators 8 and 9
const double inner_radius = 1.5;
const double outer_radius = 2.5;
static const CylinderBoundary<3> inner_cylinder(inner_radius, 2);
static const CylinderBoundary<3> outer_cylinder(outer_radius, 2);
triangulation.set_boundary (8, inner_cylinder);
triangulation.set_boundary (9, outer_cylinder);
// Then loop over all faces of the domain and, if for the position
// of the center of a face the following holds then set boundary
// indicators:
// - if y>3 and z<=2.5 or z>=5 then use boundary indicator 8
// - if y>3 and 2.5<=z<=5 then use boundary indicator 9
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
for (unsigned int f=0;
f < GeometryInfo<3>::faces_per_cell;
++f)
{
const Point<3> face_center = cell->face(f)->center();
if (cell->face(f)->at_boundary())
{
if ((face_center[2] <= 2.5 || face_center[2] >= 5) &&
face_center[1] >= 3)
cell->face(f)->set_boundary_id(8);
if (face_center[2] >= 2.5 &&
face_center[2] <= 5
&& face_center[1] >= 3)
cell->face(f)->set_boundary_id(9);
}
}
// Then refine the mesh once
triangulation.refine_global (1);

With this code, we get a mesh that looks like this:

This is clearly not correct: The new vertices that have been entered at mid-edge and mid-face points are not where they should have been. Upon some reflection, it turns out that while the radii of the cylinders are correct, the axes of the two cylinder objects should not have been along coordinate axes but shifted. This can be corrected by creating them as follows, the two points given as arguments indicating the direction and a point on the axis:

static const CylinderBoundary<3> inner_cylinder (inner_radius,
Point<3>(0,0,1),
Point<3>(0,3,0));
static const CylinderBoundary<3> outer_cylinder (outer_radius,
Point<3>(0,0,1),
Point<3>(0,3,0));
triangulation.set_boundary (9, outer_cylinder);

This yields an improvement, though it is still not quite correct:

Looking closely at this mesh, we realize that the new points on mid-face vertices are where they should be, though the new vertices inserted at mid-edge points are in the wrong place (you see this by comparing the picture with the one of the coarse mesh). What is happening is that we are only telling the triangulation to use these geometry objects for the faces but not for the adjacent edges as well. This is easily fixed by using the function TriaAccessor::set_all_boundary_ids() instead of TriaAccessor::set_boundary_id() used above. With this change, the grid now looks like this:

This is already better. However, something is still going wrong on the front left face. On second look, we can also see that the faces where the geometry widens have been refined at the bottom, that there is one transition face that looks wrong because it has a triangle rather than a quadrilateral, and that finally the transition faces in the cylindrical region appear not to have been refined at all in radial direction.

This due to the fact that we have (erroneously) marked all boundary faces between \(0\le z\le 2.5\) with the boundary indicator for the small cylinder and similarly for the other regions. This condition includes the faces parallel to the x-y plane. To fix it, we need to exclude faces whose center points have \(z\) values equal to (or at least close to, since we should not compare for equality in floating point arithmetic) 0, 2.5, 5 or 7.5. This replacement code does the trick:

// Then loop over all faces of the domain and, if for the position
// of the center of a face the following holds then set boundary
// indicators:
// - if y>3 and z<2.5 or z>5 then use boundary indicator 8
// - if y>3 and 2.5<z<5 then use boundary indicator 9
// In this process, exclude faces whose z-coordinates are
// within a small distance of z=0, z=2.5, z=5 or z=7.5.
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
for (unsigned int f=0;
f < GeometryInfo<3>::faces_per_cell;
++f)
{
const Point<3> face_center = cell->face(f)->center();
if (cell->face(f)->at_boundary())
if ((face_center[2]>1e-6) &&
(face_center[2]<7.5-1e-6) &&
(std::fabs(face_center[2]-2.5)>1e-6) &&
(std::fabs(face_center[2]-5.0)>1e-6))
{
if ((face_center[2] < 2.5 || face_center[2] > 5)
&& face_center[1] >= 3)
cell->face(f)->set_all_boundary_ids(8);
if (face_center[2] > 2.5 && face_center[2] < 5
&& face_center[1] >= 3)
cell->face(f)->set_all_boundary_ids(9);
}
}
}

With this, we finally get a mesh that looks good:

We can then refine the mesh two more times to see in more detail what happens to the curved part of the boundary:

So, yes!, this is finally what we were looking for!

The plain program

/* ---------------------------------------------------------------------
*
* Copyright (C) 2013 - 2017 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*
* Author: Timo Heister, Texas A&M University, 2013
*/
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <iostream>
#include <fstream>
#include <map>
using namespace dealii;
template <int dim>
void print_mesh_info(const Triangulation<dim> &triangulation,
const std::string &filename)
{
std::cout << "Mesh info:" << std::endl
<< " dimension: " << dim << std::endl
<< " no. of cells: " << triangulation.n_active_cells() << std::endl;
{
std::map<unsigned int, unsigned int> boundary_count;
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
{
if (cell->face(face)->at_boundary())
boundary_count[cell->face(face)->boundary_id()]++;
}
}
std::cout << " boundary indicators: ";
for (std::map<unsigned int, unsigned int>::iterator it=boundary_count.begin();
it!=boundary_count.end();
++it)
{
std::cout << it->first << "(" << it->second << " times) ";
}
std::cout << std::endl;
}
std::ofstream out (filename.c_str());
GridOut grid_out;
grid_out.write_eps (triangulation, out);
std::cout << " written to " << filename
<< std::endl
<< std::endl;
}
void grid_1 ()
{
Triangulation<2> triangulation;
GridIn<2> gridin;
gridin.attach_triangulation(triangulation);
std::ifstream f("untitled.msh");
gridin.read_msh(f);
print_mesh_info (triangulation, "grid-1.eps");
}
void grid_2 ()
{
std::vector< unsigned int > repetitions(2);
repetitions[0]=3;
repetitions[1]=2;
Point<2>(1.0,-1.0),
Point<2>(4.0,1.0));
Triangulation<2> triangulation;
GridGenerator::merge_triangulations (tria1, tria2, triangulation);
print_mesh_info (triangulation, "grid-2.eps");
}
void grid_3 ()
{
Triangulation<2> triangulation;
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
{
Point<2> &v = cell->vertex(i);
if (std::abs(v(1)-1.0)<1e-5)
v(1) += 0.5;
}
}
const SphericalManifold<2> boundary_description(Point<2>(0,0));
triangulation.set_manifold (1, boundary_description);
triangulation.refine_global(2);
print_mesh_info (triangulation, "grid-3.eps");
triangulation.set_manifold (1);
}
void grid_4()
{
Triangulation<2> triangulation;
GridGenerator::extrude_triangulation (triangulation, 3, 2.0, out);
print_mesh_info (out, "grid-4.eps");
}
Point<2> grid_5_transform (const Point<2> &in)
{
return Point<2>(in(0),
in(1) + std::sin(in(0)/5.0*3.14159));
}
void grid_5()
{
Triangulation<2> triangulation;
std::vector<unsigned int> repetitions(2);
repetitions[0] = 14;
repetitions[1] = 2;
repetitions,
Point<2>(0.0,0.0),
Point<2>(10.0,1.0));
GridTools::transform (&grid_5_transform, triangulation);
print_mesh_info (triangulation, "grid-5.eps");
}
struct Grid6Func
{
double trans(const double y) const
{
return std::tanh(2*y)/tanh(2);
}
Point<2> operator() (const Point<2> &in) const
{
return Point<2> (in(0),
trans(in(1)));
}
};
void grid_6()
{
Triangulation<2> triangulation;
std::vector< unsigned int > repetitions(2);
repetitions[0] = repetitions[1] = 40;
repetitions,
Point<2>(0.0,0.0),
Point<2>(1.0,1.0));
GridTools::transform(Grid6Func(), triangulation);
print_mesh_info (triangulation, "grid-6.eps");
}
void grid_7()
{
Triangulation<2> triangulation;
std::vector<unsigned int> repetitions(2);
repetitions[0] = repetitions[1] = 16;
GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions,
Point<2>(0.0,0.0),
Point<2>(1.0,1.0));
GridTools::distort_random (0.3, triangulation, true);
print_mesh_info (triangulation, "grid-7.eps");
}
int main ()
{
try
{
grid_1 ();
grid_2 ();
grid_3 ();
grid_4 ();
grid_5 ();
grid_6 ();
grid_7 ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
}