@dealiiVideoLecture{9}
-After we have created a grid in the previous tutorial program (in step-1),
-we now show how to define degrees of freedom on this mesh. For this example, we
-will use the lowest order ($Q_1$) finite elements, for which the degrees
-of freedom are associated with the vertices of the mesh. Later
+The finite element method is based on approximating the solution $u$ of a
+differential equation such as $-\Delta u=f$ by a function $u_h$ that is
+"piecewise" polynomial; that is, we subdivide the domain $\Omega$ on which
+the equation is posed into small cells that in the documentation we will
+generally denote by the symbol $K$. On each cell $K$, the approximating
+function $u_h$ we seek is then a polynomial. (Or, strictly speaking, a function
+that is the image of a polynomial from a "reference cell", but let us not make
+things more complicated than necessary for now.)
+
+In the previous tutorial program (in step-1), we showed how we should think of
+the subdivision of the domain into cells as a "mesh" represented by the
+Triangulation class, and how this looks like in code. In the current tutorial program,
+we now show how one represents piecewise polynomial functions through the
+concept of degrees of freedom defined on this mesh. For this example, we
+will use the lowest order ($Q_1$) finite elements, that is the approximating
+function $u_h$ we are looking for will be "bi-linear" on each quadrilateral
+cell $K$ of the mesh. (They would be linear if we would work on triangles.)
+
+In practice, we represent the function as a linear combination of shape
+functions $\varphi_j(\mathbf x)$ with multipliers $U_j$ that we call the
+"degrees of freedom". For the bi-linear functions we consider here, each
+of these shape functions and degrees of freedom is associated with a
+vertex of the mesh. Later
examples will demonstrate higher order elements where degrees of freedom are
not necessarily associated with vertices any more, but can be associated
with edges, faces, or cells.
description of finite element problems is often to say that we are looking for
a finite dimensional function $u_h \in V_h$ that satisfies some set of equations
(e.g. $a(u_h,\varphi_h)=(f,\varphi_h)$ for all test functions $\varphi_h\in
-V_h$). In other words, all we say here that the solution needs to lie in some
+V_h$). In other words, all we say here is that the solution needs to lie in some
space $V_h$. However, to actually solve this problem on a computer we need to
choose a basis of this space; this is the set of shape functions
$\varphi_j(\mathbf x)$ we have used above in the expansion of $u_h(\mathbf x)$
with coefficients $U_j$. There are of course many bases of the space $V_h$,
but we will specifically choose the one that is described by the finite
element functions that are traditionally defined locally on the cells of the
-mesh. Describing "degrees of freedom" in this context requires us to simply
+mesh.
+
+<h3> Enumerating degrees of freedom </h3>
+
+Describing "degrees of freedom" in this context requires us to simply
<i>enumerate</i> the basis functions of the space $V_h$. For $Q_1$ elements
this means simply enumerating the vertices of the mesh in some way, but for
higher order elements, one also has to enumerate the shape functions that are
indices we use for vertices. The class that
provides this enumeration of the basis functions of $V_h$ is called DoFHandler.
-Defining degrees of freedom ("DoF"s in short) on a mesh is a rather
+Defining degrees of freedom ("DoF"s in short) on a mesh is, in practice, a rather
simple task, since the library does all the work for you. Essentially,
all you have to do is create a finite element object (from one of the
many finite element classes deal.II already has, see for example the
@ref fe documentation) and give it to a DoFHandler object through the
-DoFHandler::distribute_dofs function ("distributing DoFs" is the term we use
+DoFHandler::distribute_dofs() function ("distributing DoFs" is the term we use
to describe the process of <i>enumerating</i> the basis functions as discussed
above). The DoFHandler is a class that
knows which degrees of freedom live where, i.e., it can answer
how big your system matrix should be, and when copying the
contributions of a single cell into the global matrix.
+The first task of the current program is therefore to take a mesh
+and a finite element, and enumerate the degrees of freedom. In the
+current context, this means simply giving each vertex of the mesh a
+DoF index. Once that has happened, we will output in a picture which
+vertex ended up with which DoF index. You can find the corresponding
+pictures in the <a href="#Results">results section</a> of this tutorial.
+
+It is probably worth pointing out that where each DoF is geometrically
+located is not a question we typically ask in finite element
+codes. Most often, we only care about the fact that there *is* an
+enumeration of all degrees of freedom, but not which DoF is where.
+(We will also come back to this below
+where we talk about renumbering degrees of freedom.) At the same
+time, it is probably instructive to see this once, and so this program
+shows such a figure.
+
+
<h3> Sparsity </h3>
The next step would then be to compute a matrix and right hand side
<h3> How degrees of freedom are enumerated </h3>
-By default, the DoFHandler class enumerates degrees of freedom on a mesh in a
-rather random way; consequently, the sparsity pattern is also not
+By default, the DoFHandler class enumerates degrees of freedom on a mesh using
+an algorithm that is difficult to describe and leads to results that do look
+right if you know what it is doing but otherwise appears rather random;
+consequently, the sparsity pattern is also not
optimized for any particular purpose. To show this, the code below will
demonstrate a simple way to output the "sparsity pattern" that corresponds to
a DoFHandler, i.e., an object that represents all of the potentially nonzero
<h1>Results</h1>
-The program has, after having been run, produced two sparsity
-patterns. We can visualize them by opening the <code>.svg</code> files in a web browser.
-
-The results then look like this (every point denotes an entry which
-might be nonzero; of course the fact whether the entry actually is
+The program has, after having been run, produced two files of
+DoF locations and sparsity patterns each (once for the original numbering
+and once after renumbering), along with one mesh file.
+
+Let us start with the DoF locations. There is no particularly convenient
+program to visualize this kind of information, but we can resort to
+<a href="http://www.gnuplot.info/">GNUPLOT</a> (one of the simpler visualization
+programs; maybe not the easiest to use since it is command line driven, but
+also universally available on all Linux and other Unix-like systems). The
+command that produces the following pictures reads as follows:
+```
+plot [-0.1:2.1][-1.1:1.1] "mesh.gnuplot" with lines, "dof_locations_1.gnuplot" using 1:2:3 with labels point offset .3,.2 font "4,6"
+```
+This may be cryptic, but what exactly this does is also not particularly important and you shouldn't
+spend too much time understanding what it does. Rather, the important part is
+to look at what we get as output:
+<table style="width:60%" align="center">
+ <tr>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-2.dof-locations-1.png" alt=""></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-2.dof-locations-2.png" alt=""></td>
+ </tr>
+</table>
+What these figures show is (i) a numeric label attached to each vertex -- the
+DoF index, and (ii) that the original enumeration on the left differs from
+the renumbered one on the right. Which of the two is "better" is of course
+a different question (with the answer depending on what we want to do
+with these degrees of freedom); the important point is that for the
+same mesh, one can come up with many different enumerations of the degrees
+of freedom.
+
+As for the sparsity patterns, we can visualize these by opening the
+<code>.svg</code> files in a web browser. The pictures below
+represent the matrix, and every red square denotes an entry which
+might be nonzero. (Whether the entry actually is
zero or not depends on the equation under consideration, but the
indicated positions in the matrix tell us which shape functions can
and which can't couple when discretizing a local, i.e. differential,
-equation):
+equation.)
<table style="width:60%" align="center">
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-2.sparsity-1.svg" alt=""></td>
(that's what the argument "1" to the FE_Q object is). Explore how the
sparsity pattern changes if you use higher order elements, for example
cubic or quintic ones (by using 3 and 5 as the respective arguments).
+You might also want to see where DoFs are now located -- but for that
+you likely want to work with a mesh with fewer cells because DoFs
+are now also located on edges and in the interior of cells.
You could also explore how the sparsity pattern changes by refining
the mesh. You will see that not only the size of the matrix
gnuplot> set style data points
gnuplot> plot "sparsity-pattern-1.gnuplot"
@endcode
-
-Another practice based on
-<a href="http://www.gnuplot.info/">GNUPLOT</a> is trying to
-print out the mesh with locations and numbering of the support
-points. For that, you need to include header files for GridOut and MappingQ1.
-The code for this is:
-@code
- std::ofstream out("gnuplot.gpl");
- out << "plot '-' using 1:2 with lines, "
- << "'-' with labels point pt 2 offset 1,1"
- << std::endl;
- GridOut().write_gnuplot (triangulation, out);
- out << "e" << std::endl;
- const int dim = 2;
- std::map<types::global_dof_index, Point<dim> > support_points;
- DoFTools::map_dofs_to_support_points (MappingQ1<dim>(),
- dof_handler,
- support_points);
- DoFTools::write_gnuplot_dof_support_point_info(out,
- support_points);
- out << "e" << std::endl;
-@endcode
-After we run the code, we get a file called gnuplot.gpl. To view this
-file, we can run the following code in the command line:
-@code
-gnuplot -p gnuplot.gpl
-@endcode.
-With that, you will get a picture similar to
-@image html support_point_dofs1.png
-depending on the mesh you are looking at. For more information, see DoFTools::write_gnuplot_dof_support_point_info.
// require additional comments:
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
// However, the next file is new. We need this include file for the
// association of degrees of freedom ("DoF"s) to vertices, lines, and cells:
// but also 1d and 3d.)
#include <deal.II/fe/fe_q.h>
// In the following file, several tools for manipulating degrees of freedom
-// can be found:
+// can be found, and the one after it is necessary to call one of the
+// functions imported from `dof_tools.h`:
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/mapping_q1.h>
+
// We will use a sparse matrix to visualize the pattern of nonzero entries
// resulting from the distribution of degrees of freedom on the grid. That
// class can be found here:
// scope:
using namespace dealii;
+
// @sect3{Mesh generation}
// This is the function that produced the circular grid in the previous step-1
// example program with fewer refinements steps. The sole difference is that it
// returns the grid it produces via its argument.
+//
+// At the end of the function, we also output this mesh into a
+// file. We will use this as one piece of information when visualizing
+// the location of degrees of freedom. To output a mesh, we use the
+// GridOut class that you have already seen in step-1; the difference
+// is only that we use gnuplot rather than SVG format, because gnuplot
+// is the program we will use to visualize DoF locations.
void make_grid(Triangulation<2> &triangulation)
{
const Point<2> center(1, 0);
triangulation.execute_coarsening_and_refinement();
}
+
+ std::ofstream mesh_file("mesh.gnuplot");
+ GridOut().write_gnuplot(triangulation, mesh_file);
}
+
+// @sect3{Outputting the location of degrees of freedom}
+
+// The next function outputs the locations of degrees of freedom for
+// later visualization. Where each DoF is located is something the
+// DoFHandler object knows, so that is one of the arguments to this
+// function. Since we want to do all of this twice (once for the
+// original enumeration and once for the renumbered set of degrees of
+// freedom), the function also takes as a second argument the name of
+// the file into which we want the output to be written.
+//
+// In order to learn deal.II, it is probably not terribly important to
+// understand exactly what this function does, and you can skip over
+// it. But if you would like to know anyway: We want to call the
+// function DoFTools::map_dofs_to_support_points() that returns a list
+// of locations (in the form of Point objects that store
+// coordinates). It does so in the form of a map through which we can
+// query (in a statement such as `dof_location_map[42]`) where the DoF
+// is located (in the example, where the 42nd DoF is). It puts this
+// information into the `dof_location_map` object.
+//
+// We then use the function
+// DoFTools::write_gnuplot_dof_support_point_info() to write this
+// information into a file in a format that is understandable to the
+// gnuplot program that we will use for visualization in the results
+// section.
+void write_dof_locations(const DoFHandler<2> &dof_handler,
+ const std::string & filename)
+{
+ std::map<types::global_dof_index, Point<2>> dof_location_map;
+ DoFTools::map_dofs_to_support_points(MappingQ1<2>(),
+ dof_handler,
+ dof_location_map);
+
+ std::ofstream dof_location_file(filename);
+ DoFTools::write_gnuplot_dof_support_point_info(dof_location_file,
+ dof_location_map);
+}
+
+
// @sect3{Creation of a DoFHandler}
// Up to now, we only have a grid, i.e. some geometrical (the position of the
const FE_Q<2> finite_element(1);
dof_handler.distribute_dofs(finite_element);
- // Now that we have associated a degree of freedom with a global number to
- // each vertex, we wonder how to visualize this? There is no simple way to
- // directly visualize the DoF number associated with each vertex. However,
- // such information would hardly ever be truly important, since the
- // numbering itself is more or less arbitrary. There are more important
- // factors, of which we will demonstrate one in the following.
+ // Now that we have associated a degree of freedom with a global
+ // number to each vertex, Let us output this information using the
+ // function above:
+ write_dof_locations(dof_handler, "dof_locations_1.gnuplot");
+
+ // In practice, we do not often care about where a degree of freedom
+ // is geometrically located, and so other than seeing it once via
+ // the call above is not practically useful. But where two degrees
+ // of freedom are in relation to each other matters in other ways.
//
// Associated with each vertex of the triangulation is a shape
// function. Assume we want to solve something like Laplace's equation, then
{
DoFRenumbering::Cuthill_McKee(dof_handler);
+ write_dof_locations(dof_handler, "dof_locations_2.gnuplot");
+
+
DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs(),
dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dynamic_sparsity_pattern);
* given map @p support_points. For each support point location, a string
* label containing a list of all DoFs from the map is generated. The map
* can be generated with a call to map_dofs_to_support_points() and is useful
- * to visualize location and global numbering of unknowns.
+ * to visualize location and global numbering of unknowns. This function is
+ * used in step-2.
*
* An example for the format of each line in the output is:
* @code
projects / dealii.git / commitdiff