Reference documentation for deal.II version GIT relicensing-249-g48dc7357c7 2024-03-29 12:30:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The step-49 tutorial program

This tutorial depends on step-1.

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, Wolfgang Bangerth, David Wells, and Sean Ingimarson.

Introduction

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

Note
This tutorial is also available as a Jupyter Python notebook that uses the deal.II python interface. The notebook is available in the same directory as the original C++ program.

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 .vtu 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(), and 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.

One of the issues you will encounter is that deal.II, at least until version 9.2, can only deal with meshes that only consist of quadrilaterals and hexahedra – tetrahedral meshes were not supported and will likely not be supported with all of the features deal.II offers for quadrilateral and hexahedral meshes for several versions following the 9.3 release that introduced support for simplicial and mixed meshes first. As a consequence, let us not show how to work with the tetgen mesh generator that is widely used but can only generate tetrahedral meshes, but instead illustrate how one can work with Gmsh instead.

Gmsh is the smallest and most quickly set up open source tool we are aware of. It can generate unstructured 2d quad meshes. In 3d, it can extrude 2d meshes to get hexahedral meshes; 3D meshing of unstructured geometry into hexahedra is possible, though there are some issues with the quality of these meshes that imply that these meshes only sometimes work in deal.II.

Gmsh has a graphical user interface, but what it does is in essence only to generate a text file that then drives the actual mesh generation. The graphical user interface does allow querying information about the geometry and mesh, however, and so is useful. But it is useful to understand the text format regardless because it allows automating the mesh generation workflow. These text files, with a .geo suffix, can contain computations, loops, variables, etc. This format is quite flexible in allowing the description of complex geometries. The mesh is then generated from a surface representation, which is built from a list of line loops, which is built 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 if you want to write it by hand and see the effects in the graphical user interface of gmsh.

This tutorial contains an example .geo file that describes a box with two circles cut out in the interior and several slits on the sides. For the example.geo file that you can find in the examples/step-49 directory, Gmsh will show the following view (displaying the boundary indicators as well as the mesh are discussed further down below):

We will go through the basics of Gmsh and show how we can obtain a mesh like this.

It's also worth mentioning this video (which can be found here https://www.youtube.com/watch?v=xL2LmDsDLYw&t=921s) which does a good job of explaining how to use gmsh to create basic shapes and meshes. Let us first go through the documented example.geo file that describes the geometry.

We will now go through the main parts of example.geo file that defines the mesh.

Point declaration

In this section we create the points that make up the domain.

The points can either be typed in the .geo file manually or found via Geometry -> Elementary Entities -> Add -> Point where you simply enter the coordinates in the prompt or use the mouse to place them. The latter approach yields a list of the following kind:

Point(1) = {-3.5, -3.36, 0, 1.0};
Point(2) = {3.5, -3.36, 0, 1.0};
Point(3) = {3.5, 3.36, 0, 1.0};
Point(4) = {-3.5, 3.36, 0, 1.0};
Point(5) = {-3.5, 0.2, 0, 1.0};
Point(6) = {-0.98 - b, 0.2, 0, 1.0};
...
Definition point.h:111
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)

It is relevant to note that all points in Gmsh are three-dimensional objects. Since we here want to generate a two-dimensional mesh, the points simply have a zero \(z\) coordinate. The fourth number in the curly braces for each point (equal to 1.0 for all of the points above) indicates the desired mesh size in the vicinity of this point. Gmsh's graphical user interfaces writes this into the .geo file automatically, but it can be omitted and one would probably do that if one were to write this file by hand.

The file contains many more points than just these six. If you look into the file, you will also realize that one does not have to enumerate points consecutively: One can number them in whichever way one wants, which is often useful when constructing complex geometries. In such cases, one might for example want to number all points for one particular part of the geometry starting at zero, and the points for another part at, say, 1000. It does not matter whether all numbers between zero and 1000 are used.

Lines on outer domain

To create lines of the mesh, go to Geometry -> Elementary entities -> Add -> Line. You do not get a prompt to enter in specific coordinates, rather you simply click a starting point and ending point for each line.

The generated code in the geo file then looks like this:

Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 9};
...
constexpr const ReferenceCell Line

What appears on the right are pairs of point indices that define a line. As before, the indices of the lines (here 4 to 8) need not be consecutive.

Creating holes in domain

Here are the points and lines that create the two holes in the domain. We start by defining the relevant points for bottom left hole:

Point(31) = {-2.1, -1.54, 0, 1.0}; // origin
Point(32) = {-2.1, -0.84, 0, 1.0}; // top point
Point(33) = {-2.1, -2.24, 0, 1.0}; // bottom point
Point(103) = {-1.4, -1.54, 0, 1.0}; // right point
Point(104) = {-2.8, -1.54, 0, 1.0}; // left point

Then we use the Circle arc entity, which is found in the same category as Point and Lines. The file reads "Circle" but the difference to note is that this function uses three points to make an arc of a circle, namely the starting point, center, and ending point. The angle of an arc is required to be less than 180 degrees. We therefore split the circle into four arcs, which also allows us to set different boundary IDs. This then looks as follows where point 31 is the common center for all arcs:

Circle(48) = {103, 31, 32};
Circle(49) = {32, 31, 104};
Circle(50) = {104, 31, 33};
Circle(51) = {33, 31, 103};

We follow the same procedure for the second circle.

Plane surfaces

This section describes the "Plane Surfaces", i.e., the 2D surfaces for meshing. This can be found in Elementary entities as Plane Surface. The purpose of this entity is to specify the domain in 2D. Because we have our mesh split into 4 sections, we must have 4 plane surfaces.

To make a plane surface, simply click Plane surface and click on all the relevant borders for the mesh. This includes holes or other objects contained in the mesh. Again, since we have 4 different meshes, we will do this 4 times.

The code below is for the top right mesh. Notice we have 2 instances of Curve Loop, which are constructed from lists of the line segments we have previously built from pairs of two points. The first of the curve loops is the outer boundary of the part of the domain we are describing, whereas the second one is for the hole. This can easily be traced back by following along with the numbers in the braces.

Curve Loop(1) = {40, -43, 10, 57, 59, 60, 58};
Curve Loop(2) = {12, 13, 14, 15};
Plane Surface(1) = {1, 2};

In the description of the first curve loop above, the second line segment (segment 43) is described with a negative sign. This is because we have not paid enough attention when creating the line segment and have specified its two vertices in the wrong order:

Line(14) = {7, 14};
...
Line(40) = {3, 36};
...
Line(43) = {7, 36};

The negative sign then ensures that it is considered in the correct orientation, with an order for vertices of 3 -> 36 -> 7 -> 14.

The other four surfaces are declared in a similar way.

Next we define the physical surface itself. This line is what makes our mesh 2D. The values in the braces on the right hand side are the tags for the Plane Surfaces we declared above. The number in the parentheses of physical surfaces will become the material ID of the mesh as well. For more information about material IDs, check the glossary entry material indicators .

Physical Surface(2) = {1, 2, 3, 4};

Physical IDs

The steps above would be enough to define a two dimensional geometry that can be meshed. In practice, however, one often needs to attach additional information to parts of the geometry. The prototypical example are boundary indicators and material indicators .

To do this, we need to tell Gmsh about all of the "physical entities" and assign them with indices. For example, the physical surface above has been assigned index 2, and so this index "2" will be output as a tag for the cells that make up the surface (i.e., all cells of this mesh) and deal.II will then interpret tags on cells as material ids.

Physical Curve(0) = {38, 39, 52:57, 59, 60, 58, 40, 41, 4:9};

The indices in the list on the right correspond to line segments, with a colon "`:`" denoting a range.

We then assign boundary ID of 1 to the top right circle, where the index 12 corresponds to the circle (a different kind of line segment):

Physical Curve(1) = {12};

Finally, we assign boundary IDs of 2 and 3 to the top and bottom half of the bottom left hole, respectively. We can see here that we needed to pick two boundaries to make both of the physical curves. Recall this circle is comprised of four circle arcs, so Physical Curve(2) is the top left and top right quarter of the circle, and Physical Curve(3) is the bottom left and bottom right quarter of the circle.

Physical Curve(2) = {48, 49};
Physical Curve(3) = {50, 51};
Note
It is important that the .geo file contains some "physical lines" and "physical surfaces", not just for the purposes of material and boundary ids: Without physical lines and surfaces (and volumes in 3d), there is only a collection of geometric objects, but no "domain" that Gmsh can actually define a mesh on.

Meshing Parameters

Finally, the .geo file contains meshing parameters that can all be adjusted in the GMSH GUI. We include these parameters at the bottom of our .geo file to keep the changes when creating our mesh (with suffix .msh) file.

Mesh.Algorithm = 3;
Mesh.RecombineAll = 1;
Mesh.CharacteristicLengthFactor = .6;
Mesh.SubdivisionAlgorithm = 1;
Mesh.Smoothing = 20;
Show "*";

To view what each of these do, press Tools -> Options -> Mesh.

You may change the algorithm used for meshing in 2D algorithm.

CharacteristicLengthFactor can be thought of the distance between nodes on the boundaries of the mesh. This sets the initial size of the nodes of the mesh. If we want a finer mesh we want a smaller characteristic length. For this examples purposes we would like to do more refinement in deal.II so we pick a relatively coarse mesh. This can also be adjusted in the Options menu.

SubdivisionAlgorithm is set to 1 in our file because we want to use "All Quads" for subdivision. In the options menu the default is "None", but we use "All Quads" for the reasons mentioned earlier.

Smoothing is a post-processing step that iteratively improves mesh quality (for example by moving inner vertices). The number of steps is not crucial in this example. For our example we pick 20, however this value can be adjusted accordingly as it's not the most important step in mesh building.

Note
That ends the code for example.geo. It's important to note that the .geo file did not come out as organized as it appears here. It was edited after creating the mesh and geometry to be more readable.

Before loading this into deal.II, we need to create the .msh file. If everything was set up correctly, we can press Mesh -> 2D. If all this goes well, then we should see something similar to the following output:

It behooves us at this point to also make sure our boundaries were declared correctly. So in the Options menu go to Mesh -> Visibility. In the dropdown menu labeled "Label type" we may choose a label for our boundaries. We wish to check physical IDs, so we would like to have "Physical tag" picked. Now we check "1D element labels" and we should see the numbers assigned to the specified boundaries.

Finally, we can choose Mesh -> Save to save the .msh file, which deal.II's GridIn class can read. You can also generate the .msh from the .geo file by running the command

gmsh -2 example.geo
Note
Gmsh has a number of other interfaces by which one can describe geometries to it. In particular, it has the ability to interface with scripting languages like Python and Julia, but it can also be scripted from C++. These interfaces are useful if one doesn't just want to generate a mesh for a single geometry (in which case the graphical interface or, in simple cases, a hand-written .geo file is probably the simplest approach), but instead wants to do parametric studies over the geometry for which it is necessary to generate many meshes for geometries that differ in certain parameters. Another case where this is useful is if there is already a CAD geometry for which one only needs a mesh; indeed, this can be done from within deal.II using the Gmsh::create_triangulation_from_boundary_curve() function.

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. (Superconvergence refers to the fact that if a mesh has certain symmetries – for example, if the edges running into a vertex are symmetric to this vertex, and if this is so for all vertices of a cell – that the solution is then often convergent with a higher order than one would have expected from the usual error analysis. In the end, this is a result of the fact that if one were to make a Taylor expansion of the error, the symmetry leads to the fact that the expected next term of the expansion happens to be zero, and the error order is determined by the second next* term. A distorted mesh does not have these symmetries and consequently the error reflects what one will see when solving the equation on any kind of mesh, rather than showing something that is only reflective of a particular situation.)

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:

  • 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 VTU format that can easily be visualized in Paraview or VisIt.

  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;
 
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

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<types::boundary_id, unsigned int> boundary_count;
  for (const auto &face : triangulation.active_face_iterators())
  if (face->at_boundary())
  boundary_count[face->boundary_id()]++;
 
  std::cout << " boundary indicators: ";
  for (const std::pair<const types::boundary_id, unsigned int> &pair :
  boundary_count)
  {
  std::cout << pair.first << '(' << pair.second << " times) ";
  }
  std::cout << std::endl;
  }
 
Point< 2 > second
Definition grid_out.cc:4614
Point< 2 > first
Definition grid_out.cc:4613
STL namespace.

Finally, produce a graphical representation of the mesh to an output file :

  std::ofstream out(filename);
  GridOut grid_out;
  grid_out.write_vtu(triangulation, out);
  std::cout << " written to " << filename << std::endl << std::endl;
  }
 
void write_vtu(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:3544

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

It's worth noting that it is possible to save manifold ids when using the gmsh api. If we specify

GMSH_INCLUDE_DIR
GMSH_LIBRARY

when building deal.II, then DEAL_II_GMSH_WITH_API gets defined and and we can use GridIn::read_msh(). More details on the function can be found in its deal.II documentation.

We will be utilizing the SphericalManifold class for the holes. We need to assign manifold IDs for this purpose. As physical IDs from Gmsh are assigned to boundary IDs in deal.II, we will assign manifold IDs based on the boundary IDs loaded from the file.

  void grid_1()
  {
  const Point<2> Top_right_hole_origin(0.42, 2.0);
  const Point<2> Bottom_left_hole_origin(-2.1, -1.54);
 
  const SphericalManifold<2> Top_right_manifold(Top_right_hole_origin);
  const SphericalManifold<2> Bottom_left_manifold(Bottom_left_hole_origin);
 
 
  GridIn<2> gridin;
  std::ifstream f("example.msh");
  gridin.read_msh(f);
 
void attach_triangulation(Triangulation< dim, spacedim > &tria)
Definition grid_in.cc:154

Here is where we get the boundary IDs made in gmsh, which are in the first coordinate position, and assign them to manifold ids. With our example, we have boundary ID 1 on the top right hole and 2 and 3 for the bottom left hole. We assign both of these boundary IDs 2 because together they make a circle to match the manifold we assign it later.

  triangulation.set_all_manifold_ids_on_boundary(1, 1); // top right hole
  triangulation.set_all_manifold_ids_on_boundary(2,
  2); // top of bottom left hole
  triangulation.set_all_manifold_ids_on_boundary(
  3, 2); // bottom of bottom left hole
 
  triangulation.set_manifold(1, Top_right_manifold);
  triangulation.set_manifold(2, Bottom_left_manifold);
 
  triangulation.refine_global(2);
 
  print_mesh_info(triangulation, "grid-1.vtu");
  }
 
 

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;
  repetitions,
  Point<2>(1.0, -1.0),
  Point<2>(4.0, 1.0));
 
 
  print_mesh_info(triangulation, "grid-2.vtu");
  }
 
 
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
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 repetitions=1, const bool colorize=false)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)

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 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()
  {
 
  for (const auto &cell : triangulation.active_cell_iterators())
  {
  for (const auto i : cell->vertex_indices())
  {
  Point<2> &v = cell->vertex(i);
  if (std::abs(v[1] - 1.0) < 1e-5)
  v[1] += 0.5;
  }
  }
 
unsigned int vertex_indices[2]
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)

In the second step we will refine the mesh twice. To do this correctly, we should place new points on the interior boundary along the surface of a circle centered at the origin. Fortunately, GridGenerator::hyper_cube_with_cylindrical_hole already attaches a Manifold object to the interior boundary, so we do not need to do anything but refine the mesh (see the results section for a fully worked example where we do attach a Manifold object).

  triangulation.refine_global(2);
  print_mesh_info(triangulation, "grid-3.vtu");
  }
 

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()
  {
 
  print_mesh_info(out, "grid-4.vtu");
  }
 
 
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})

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 an argument that can be called like a function taking a Point and returning a Point. There are different ways of providing such an argument: It could be a pointer to a function; it could be an object of a class that has an operator(); it could be a lambda function; or it could be anything that is described via a std::function<Point<2>(const Point<2>)> object.

Decidedly the more modern way is to use a lambda function that takes a Point and returns a Point, and that is what we do in the following:

  void grid_5()
  {
  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(
  [](const Point<2> &in) {
  return Point<2>(in[0], in[1] + std::sin(numbers::PI * in[0] / 5.0));
  },
  print_mesh_info(triangulation, "grid-5.vtu");
  }
 
 
 
static constexpr double PI
Definition numbers.h:259
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)

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 {in[0], trans(in[1])};
  }
  };
 
 
  void grid_6()
  {
  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.vtu");
  }
 
 

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 (because results are generally better on meshes with "nicely shaped" cells than on the deformed cells produced by GridTools::distort_random()), 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 superconvergence properties.

  void grid_7()
  {
  std::vector<unsigned int> repetitions(2);
  repetitions[0] = repetitions[1] = 16;
  repetitions,
  Point<2>(0.0, 0.0),
  Point<2>(1.0, 1.0));
 
  print_mesh_info(triangulation, "grid-7.vtu");
  }
 
 
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)

The main function

Finally, the main function. There isn't much to do here, only to call all the various functions we wrote above.

  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 .vtu files of the triangulations. The methods are discussed above.

Next steps: Curved Cells

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, faces, and cells. By default, these new points will be placed at the arithmetic mean of the surrounding points, but this isn't what you want if you need curved boundaries that aren't already adequately resolved by the coarse mesh. For example, for this mesh the central hole is supposed to be round:

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 straight lines. 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 geometry. The way to do this requires three steps:

All functions in the GridGenerator namespace which create a mesh where some cells should be curved also attach the correct Manifold object to the provided Triangulation: i.e., for those functions we get the correct behavior by default. For a hand-generated mesh, however, the situation is much more interesting.

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 cross sections, 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 @ref step_14 "step-14".
void create_2d_grid(
const std::vector<Point<2>> &vertices,
const std::vector<
std::array<unsigned int, GeometryInfo<2>::vertices_per_cell>>
Triangulation<2> &coarse_grid)
{
std::vector<CellData<2>> cells(vertex_indices.size());
for (unsigned int i = 0; i < cells.size(); ++i)
{
for (unsigned int j = 0; j < vertex_indices[i].size(); ++j)
cells[i].vertices[j] = vertex_indices[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 std::vector<Point<2>> vertices_1{{-1.5, 0.},
{-0.5, 0.},
{0.5, 0.},
{1.5, 0.},
{-1.5, 1.5},
{-0.5, 1.5},
{0.5, 1.5},
{1.5, 1.5},
{-1.5, 3.},
{-0.5, 3.},
{0.5, 3.},
{1.5, 3.},
{-0.5, 3 + 0.5 * sqrt(3)},
{0.5, 3 + 0.5 * sqrt(3)},
{-0.75, 3 + 0.75 * sqrt(3)},
{0.75, 3 + 0.75 * sqrt(3)}};
const std::vector<std::array<unsigned int, GeometryInfo<2>::vertices_per_cell>>
cell_vertices_1 = {{{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, cell_vertices_1, 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 std::vector<Point<2>> vertices_2{{-2.5, 0.},
{-1.5, 0.},
{-0.5, 0.},
{0.5, 0.},
{1.5, 0.},
{2.5, 0.},
{-2.5, 1.5},
{-1.5, 1.5},
{-0.5, 1.5},
{0.5, 1.5},
{1.5, 1.5},
{2.5, 1.5},
{-2.5, 3.},
{-1.5, 3.},
{-0.5, 3.},
{0.5, 3.},
{1.5, 3.},
{2.5, 3.},
{-0.5, 3. + 0.5 * sqrt(3)},
{0.5, 3. + 0.5 * sqrt(3)},
{-0.75, 3. + 0.75 * sqrt(3)},
{0.75, 3. + 0.75 * sqrt(3)},
{-1.25, 3. + 1.25 * sqrt(3)},
{1.25, 3. + 1.25 * sqrt(3)}};
const std::vector<std::array<unsigned int, GeometryInfo<2>::vertices_per_cell>>
cell_vertices_2 = {{{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, cell_vertices_2, 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
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<3> triangulation_3d_tmp;
triangulation_3d_2,
triangulation_3d_tmp);
GridTools::shift(Point<3>(0, 0, 5), triangulation_3d_1);
triangulation_3d_1,
}
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
Point< 3 > vertices[4]
void consistently_order_cells(std::vector< CellData< dim > > &cells)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)

This creates the following mesh:

This mesh has the right general shape, but the top cells are now polygonal: their edges are no longer along circles and we do not have a very accurate representation of the original geometry. The next step is to teach the top part of the domain that it should be curved. Put another way, all calculations done on the top boundary cells should be done in cylindrical coordinates rather than Cartesian coordinates. We can do this by creating a CylindricalManifold object and associating it with the cells above \(y = 3\). This way, when we refine the cells on top, we will place new points along concentric circles instead of straight lines.

In deal.II we describe all geometries with classes that inherit from Manifold. The default geometry is Cartesian and is implemented in the FlatManifold class. As the name suggests, Manifold and its inheriting classes provide a way to describe curves and curved cells in a general way with ideas and terminology from differential geometry: for example, CylindricalManifold inherits from ChartManifold, which describes a geometry through pull backs and push forwards. In general, one should think that the Triangulation class describes the topology of a domain (in addition, of course, to storing the locations of the vertices) while the Manifold classes describe the geometry of a domain (e.g., whether or not a pair of vertices lie along a circular arc or a straight line). A Triangulation will refine cells by doing computations with the Manifold associated with that cell regardless of whether or not the cell is on the boundary. Put another way: the Manifold classes do not need any information about where the boundary of the Triangulation actually is: it is up to the Triangulation to query the right Manifold for calculations on a cell. Most Manifold functions (e.g., Manifold::get_intermediate_point) know nothing about the domain itself and just assume that the points given to it lie along a geodesic. In this case, with the CylindricalManifold constructed below, the geodesics are arcs along circles orthogonal to the \(z\)-axis centered along the line \((0, 3, z)\).

Since all three top parts of the domain use the same geodesics, we will mark all cells with centers above the \(y = 3\) line as being cylindrical in nature:

const Tensor<1, 3> axis({0.0, 0.0, 1.0});
const Point<3> axial_point(0, 3.0, 0.0);
const CylindricalManifold<3> cylinder(axis, axial_point);
const types::manifold_id cylinder_id = 8;
create_3d_grid(triangulation);
triangulation.set_manifold(cylinder_id, cylinder);
for (auto &cell : triangulation.active_cell_iterators())
if (cell->center()[1] >= 3.0)
cell->set_all_manifold_ids(cylinder_id);
triangulation.refine_global(1);
Point< 3 > center

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

This change fixes the boundary but creates a new problem: the cells adjacent to the cylinder's axis are badly distorted. We should use Cartesian coordinates for calculations on these central cells to avoid this issue. The cells along the center line all have a face that touches the line \((0, 3, z)\) so, to implement this, we go back and overwrite the manifold_ids on these cells to be zero (which is the default):

const Tensor<1, 3> axis({0.0, 0.0, 1.0});
const Point<3> axial_point(0, 3.0, 0.0);
const CylindricalManifold<3> cylinder(axis, axial_point);
const types::manifold_id cylinder_id = 8;
create_3d_grid(triangulation);
triangulation.set_manifold(cylinder_id, cylinder);
for (auto &cell : triangulation.active_cell_iterators())
if (cell->center()[1] >= 3.0)
cell->set_all_manifold_ids(cylinder_id);
for (const auto &cell : triangulation.active_cell_iterators())
for (const auto &face : cell->face_iterators())
{
const Point<3> face_center = face->center();
if (std::abs(face_center[0]) < 1.0e-5 &&
std::abs(face_center[1] - 3.0) < 1.0e-5)
cell->set_all_manifold_ids(numbers::flat_manifold_id);
}
triangulation.refine_global(1);
const types::manifold_id flat_manifold_id
Definition types.h:325

This gives us the following grid:

This gives us a good mesh, where cells at the center of each circle are still Cartesian and cells around the boundary lie along a circle. We can really see the nice detail of the boundary fitted mesh if we refine two more times:

Possibilities for extensions

Assigning different boundary ids

It is often useful to assign different boundary ids to a mesh that is generated in one form or another as described in this tutorial to apply different boundary conditions.

For example, you might want to apply a different boundary condition for the right boundary of the first grid in this program. To do this, iterate over the cells and their faces and identify the correct faces (for example using cell->center() to query the coordinates of the center of a cell as we do in step-1, or using cell->face(f)->get_boundary_id() to query the current boundary indicator of the \(f\)th face of the cell). You can then use cell->face(f)->set_boundary_id() to set the boundary id to something different. You can take a look back at step-1 how iteration over the meshes is done there.

Extracting a boundary mesh

Computations on manifolds, like they are done in step-38, require a surface mesh embedded into a higher dimensional space. While some can be constructed using the GridGenerator namespace or loaded from a file, it is sometimes useful to extract a surface mesh from a volume mesh.

Use the function GridGenerator::extract_boundary_mesh() to extract the surface elements of a mesh. Using the function on a 3d mesh (a Triangulation<3,3>, for example from grid_4()), this will return a Triangulation<2,3> that you can use in step-38. Also try extracting the boundary mesh of a Triangulation<2,2>.

The plain program

/* ------------------------------------------------------------------------
*
* SPDX-License-Identifier: LGPL-2.1-or-later
* Copyright (C) 2013 - 2024 by the deal.II authors
*
* This file is part of the deal.II library.
*
* Part of the source code is dual licensed under Apache-2.0 WITH
* LLVM-exception OR LGPL-2.1-or-later. Detailed license information
* governing the source code and code contributions can be found in
* LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
*
* ------------------------------------------------------------------------
*
* Author: Timo Heister, Texas A&M University, 2013
*/
#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<types::boundary_id, unsigned int> boundary_count;
for (const auto &face : triangulation.active_face_iterators())
if (face->at_boundary())
boundary_count[face->boundary_id()]++;
std::cout << " boundary indicators: ";
for (const std::pair<const types::boundary_id, unsigned int> &pair :
boundary_count)
{
std::cout << pair.first << '(' << pair.second << " times) ";
}
std::cout << std::endl;
}
std::ofstream out(filename);
GridOut grid_out;
grid_out.write_vtu(triangulation, out);
std::cout << " written to " << filename << std::endl << std::endl;
}
void grid_1()
{
const Point<2> Top_right_hole_origin(0.42, 2.0);
const Point<2> Bottom_left_hole_origin(-2.1, -1.54);
const SphericalManifold<2> Top_right_manifold(Top_right_hole_origin);
const SphericalManifold<2> Bottom_left_manifold(Bottom_left_hole_origin);
GridIn<2> gridin;
std::ifstream f("example.msh");
gridin.read_msh(f);
triangulation.set_all_manifold_ids_on_boundary(1, 1); // top right hole
triangulation.set_all_manifold_ids_on_boundary(2,
2); // top of bottom left hole
triangulation.set_all_manifold_ids_on_boundary(
3, 2); // bottom of bottom left hole
triangulation.set_manifold(1, Top_right_manifold);
triangulation.set_manifold(2, Bottom_left_manifold);
triangulation.refine_global(2);
print_mesh_info(triangulation, "grid-1.vtu");
}
void grid_2()
{
std::vector<unsigned int> repetitions(2);
repetitions[0] = 3;
repetitions[1] = 2;
repetitions,
Point<2>(1.0, -1.0),
Point<2>(4.0, 1.0));
print_mesh_info(triangulation, "grid-2.vtu");
}
void grid_3()
{
for (const auto &cell : triangulation.active_cell_iterators())
{
for (const auto i : cell->vertex_indices())
{
Point<2> &v = cell->vertex(i);
if (std::abs(v[1] - 1.0) < 1e-5)
v[1] += 0.5;
}
}
triangulation.refine_global(2);
print_mesh_info(triangulation, "grid-3.vtu");
}
void grid_4()
{
print_mesh_info(out, "grid-4.vtu");
}
void grid_5()
{
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(
[](const Point<2> &in) {
return Point<2>(in[0], in[1] + std::sin(numbers::PI * in[0] / 5.0));
},
print_mesh_info(triangulation, "grid-5.vtu");
}
struct Grid6Func
{
double trans(const double y) const
{
return std::tanh(2 * y) / tanh(2);
}
Point<2> operator()(const Point<2> &in) const
{
return {in[0], trans(in[1])};
}
};
void grid_6()
{
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.vtu");
}
void grid_7()
{
std::vector<unsigned int> repetitions(2);
repetitions[0] = repetitions[1] = 16;
repetitions,
Point<2>(0.0, 0.0),
Point<2>(1.0, 1.0));
print_mesh_info(triangulation, "grid-7.vtu");
}
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;
}
}
void read_msh(std::istream &in)
Definition grid_in.cc:2006
Expression tanh(const Expression &x)
const InputIterator OutputIterator out
Definition parallel.h:167