From c75fe59d04141a09c5abe877bfc2d5da47370e97 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 27 May 2018 13:55:48 +0800 Subject: [PATCH] Update some of the introduction of step-2. --- examples/step-2/doc/intro.dox | 92 ++++++++++++++++++++++++++--------- 1 file changed, 69 insertions(+), 23 deletions(-) diff --git a/examples/step-2/doc/intro.dox b/examples/step-2/doc/intro.dox index 2d89d08c0d..3d1fc7128a 100644 --- a/examples/step-2/doc/intro.dox +++ b/examples/step-2/doc/intro.dox @@ -14,12 +14,12 @@ with edges, faces, or cells. The term "degree of freedom" is commonly used in the finite element community to indicate two slightly different, but related things. The first is that we'd like to represent the finite element solution as a linear combination of shape -function, in the form $u_h(\mathbf x) = \sum_{j=0}^{N-1} U_j \varphi_j(\mathbf +functions, in the form $u_h(\mathbf x) = \sum_{j=0}^{N-1} U_j \varphi_j(\mathbf x)$. Here, $U_j$ is a vector of expansion coefficients. Because we don't know their values yet (we will compute them as the solution of a linear or nonlinear system), they are called "unknowns" or "degrees of freedom". The second meaning of the term can be explained as follows: A mathematical -description of finite element problem is often to say that we are looking for +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 @@ -32,8 +32,10 @@ element functions that are traditionally defined locally on the cells of the mesh. Describing "degrees of freedom" in this context requires us to simply enumerate 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 elements one also has to enumerate the shape functions that are -associated with edges, faces, or cell interiors of the mesh. The class that +higher order elements, one also has to enumerate the shape functions that are +associated with edges, faces, or cell interiors of the mesh. In other words, +the enumeration of degrees of freedom is an entirely separate thing from the +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 @@ -42,27 +44,37 @@ 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 -to describe the process of enumerating the basis functions as discussed +to describe the process of enumerating the basis functions as discussed above). The DoFHandler is a class that -manages which degrees of freedom live where, i.e., it can answer +knows which degrees of freedom live where, i.e., it can answer questions like "how many degrees of freedom are there globally" and "on this cell, give me the global indices of the shape functions that live here". This is the sort of information you need when determining how big your system matrix should be, and when copying the contributions of a single cell into the global matrix. +

Sparsity

+ The next step would then be to compute a matrix and right hand side corresponding to a particular differential equation using this finite element and mesh. We will keep this step for the step-3 program and rather talk about one practical aspect of a finite element program, namely that finite element -matrices are almost always very sparse, i.e. almost all entries in these -matrices are zero. (To be more precise, we say a discretization leads to a -sparse matrix if the number of nonzero entries per row in the matrix is +matrices are always very sparse: almost all entries in these +matrices are zero. + +To be more precise, we say that a matrix is sparse +if the number of nonzero entries per row in the matrix is bounded by a number that is independent of the overall number of degrees of freedom. For example, the simple 5-point stencil of a finite difference approximation of the Laplace equation leads to a sparse matrix since the number of nonzero entries per row is five, and therefore independent of the -total size of the matrix.) Sparsity is one of the distinguishing feature of +total size of the matrix. For more complicated problems -- say, the Stokes +problem of step-22 -- and in particular in 3d, the number of entries per row +may be several hundred. But the important point is that this number is +independent of the overall size of the problem: If you refine the mesh, the +maximal number of unknowns per row remains the same. + +Sparsity is one of the distinguishing feature of the finite element method compared to, say, approximating the solution of a partial differential equation using a Taylor expansion and matching coefficients, or using a Fourier basis. @@ -77,23 +89,57 @@ number of matrix-vector multiplications to come up with the solution of a linear system with this matrix, then we would have a solver that can find the values of all $N$ unknowns with optimal complexity, i.e., with a total of ${\cal O}(N)$ operations. It is clear that this wouldn't be possible if the -matrix were not sparse, but it also requires very specialized solvers such as +matrix were not sparse (because then the number of entries in the matrix would +have to be ${\cal O}(N^s)$ with some $s>1$, and doing a fixed number of +matrix-vector products would take ${\cal O}(N^s)$ operations), +but it also requires very specialized solvers such as multigrid methods to satisfy the requirement that the solution requires only a fixed number of matrix-vector multiplications. We will frequently look at the question of what solver to use in the remaining programs of this tutorial. The sparsity is generated by the fact that finite element shape -functions are defined locally on individual cells, rather than +functions are defined locally on individual cells, rather than globally, and that the local differential operators in the bilinear -form only couple shape functions that have some overlap. By default, -the DoFHandler class enumerates degrees of freedom on a mesh in a +form only couple shape functions whose support overlaps. (The "support" of +a function is the area where it is nonzero. For the finite element method, +the support of a shape function is generally the cells adjacent to the vertex, +edge, or face it is defined on.) In other words, degrees of freedom $i$ and $j$ +that are not defined on the same cell do not overlap, and consequently +the matrix entry $A_{ij}$ will be zero. (In some cases such +as the Discontinuous Galerkin method, shape functions may also connect +to neighboring cells through face integrals. But finite element +methods do not generally couple shape functions beyond the immediate +neighbors of a cell on which the function is defined.) + + +

How degrees of freedom are enumerated

+ +By default, the DoFHandler class enumerates degrees of freedom on a mesh in a rather random way; consequently, the sparsity pattern is also not -optimized for any particular purpose. However, for -some algorithms, especially for some linear solvers and preconditioners, it is -advantageous to have the degrees of freedom numbered in a certain -order, and we will use the algorithm of Cuthill and McKee to do -so. This can be thought of as choosing a different, permuted basis of the -finite element space. -The results are written to a file and visualized using a simple -visualization program; you get to -see the outcome in the results section below. +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 +elements of a matrix one may build when discretizing a partial differential +equation on a mesh and its DoFHandler. This lack of structure in the sparsity +pattern will be apparent from the pictures we show below. + +For most applications and algorithms, the exact way in which degrees of freedom +are numbered does not matter. For example, the Conjugate Gradient method we +use to solve linear systems does not care. On the other hand, +some algorithms do care: in particular, some preconditioners such as SSOR +will work better if they can walk through degrees of freedom in a particular +order, and it would be nice if we could just sort them in such a way that +SSOR can iterate through them from zero to $N$ in this order. Other examples +include computing incomplete LU or Cholesky factorizations, or if we care +about the block structure of matrices (see step-20 for an example). +deal.II therefore has algorithms that can re-enumerate degrees of freedom +in particular ways in namespace DoFRenumbering. Renumbering can be thought +of as choosing a different, permuted basis of the finite element space. The +sparsity pattern and matrices that result from this renumbering are therefore +also simply a permutation of rows and columns compared to the ones we would +get without explicit renumbering. + +In the program below, we will use the algorithm of Cuthill and McKee to do +so. We will show the sparsity pattern for both the original enumeration of +degrees of freedom and of the renumbered version below, +in the results section. -- 2.39.5