The step-3 tutorial program

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

Introduction

This is the first example where we actually use finite elements to compute something. We will solve a simple version of Laplace's equation with zero boundary values, but a nonzero right hand side. This example is still quite simple, but it already shows the basic structure of most finite element programs, which are along the following lines:

The commented program

Many new include files

These include files are already known to you. They declare the classes which handle triangulations and enumeration of degrees of freedom:

 #include <grid/tria.h>
 #include <dofs/dof_handler.h>

And this is the file in which the functions are declared that create grids:

 #include <grid/grid_generator.h>

The next three files contain classes which are needed for loops over all cells and to get the information from the cell objects. The first two have been used before to get geometric information from cells; the last one is new and provides information about the degrees of freedom local to a cell:

 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>

In this file contains the description of the Lagrange interpolation finite element:

 #include <fe/fe_q.h>

And this file is needed for the creation of sparsity patterns of sparse matrices, as shown in previous examples:

 #include <dofs/dof_tools.h>

The next two file are needed for assembling the matrix using quadrature on each cell. The classes declared in them will be explained below:

 #include <fe/fe_values.h>
 #include <base/quadrature_lib.h>

The following three include files we need for the treatment of boundary values:

 #include <base/function.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>

We're now almost to the end. The second to last group of include files is for the linear algebra which we employ to solve the system of equations arising from the finite element discretization of the Laplace equation. We will use vectors and full matrices for assembling the system of equations locally on each cell, and transfer the results into a sparse matrix. We will then use a Conjugate Gradient solver to solve the problem, for which we need a preconditioner (in this program, we use the identity preconditioner which does nothing, but we need to include the file anyway):

 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>

Finally, this is for output to a file and to the console:

 #include <numerics/data_out.h>
 #include <fstream>
 #include <iostream>

...and this is to import the deal.II namespace into the global scope:

 using namespace dealii;

The LaplaceProblem class

Instead of the procedural programming of previous examples, we encapsulate everything into a class for this program. The class consists of functions which each perform certain aspects of a finite element program, a `main' function which controls what is done first and what is done next, and a list of member variables.

The public part of the class is rather short: it has a constructor and a function `run' that is called from the outside and acts as something like the `main' function: it coordinates which operations of this class shall be run in which order. Everything else in the class, i.e. all the functions that actually do anything, are in the private section of the class:

 class LaplaceProblem 
 {
   public:
     LaplaceProblem ();
 
     void run ();

Then there are the member functions that mostly do what their names suggest. Since they do not need to be called from outside, they are made private to this class.

   private:
     void make_grid_and_dofs ();
     void assemble_system ();
     void solve ();
     void output_results () const;

And finally we have some member variables. There are variables describing the triangulation and the global numbering of the degrees of freedom (we will specify the exact polynomial degree of the finite element in the constructor of this class)...

     Triangulation<2>     triangulation;
     FE_Q<2>              fe;
     DoFHandler<2>        dof_handler;

...variables for the sparsity pattern and values of the system matrix resulting from the discretization of the Laplace equation...

     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;

...and variables which will hold the right hand side and solution vectors.

     Vector<double>       solution;
     Vector<double>       system_rhs;
 };

LaplaceProblem::LaplaceProblem

Here comes the constructor. It does not much more than first to specify that we want bi-linear elements (denoted by the parameter to the finite element object, which indicates the polynomial degree), and to associate the dof_handler variable to the triangulation we use. (Note that the triangulation isn't set up with a mesh at all at the present time, but the DoFHandler doesn't care: it only wants to know which triangulation it will be associated with, and it only starts to care about an actual mesh once you try to distribute degree of freedom on the mesh using the distribute_dofs() function.) All the other member variables of the LaplaceProblem class have a default constructor which does all we want.

 LaplaceProblem::LaplaceProblem () :
                 fe (1),
                 dof_handler (triangulation)
 {}

LaplaceProblem::make_grid_and_dofs

Now, the first thing we've got to do is to generate the triangulation on which we would like to do our computation and number each vertex with a degree of freedom. We have seen this in the previous examples before. Then we have to set up space for the system matrix and right hand side of the discretized problem. This is what this function does:

 void LaplaceProblem::make_grid_and_dofs ()
 {

First create the grid and refine all cells five times. Since the initial grid (which is the square [-1,1]x[-1,1]) consists of only one cell, the final grid has 32 times 32 cells, for a total of 1024.

   GridGenerator::hyper_cube (triangulation, -1, 1);
   triangulation.refine_global (5);

Unsure that 1024 is the correct number? Let's see: n_active_cells return the number of active cells:

   std::cout << "Number of active cells: "
             << triangulation.n_active_cells()
             << std::endl;

Here, by active we mean the cells on the finest level, i.e. cells that aren't refined any further. We stress the adjective `active', since there are more cells, namely the parent cells of the finest cells, their parents, etc, up to the one cell which made up the initial grid. Of course, on the next coarser level, the number of cells is one quarter that of the cells on the finest level, i.e. 256, then 64, 16, 4, and 1. We can get the total number of cells like this:

   std::cout << "Total number of cells: "
             << triangulation.n_cells()
             << std::endl;

Note the distinction between n_active_cells() and n_cells().

Next we enumerate all the degrees of freedom. This is done by using the distribute_dofs function, as we have seen in the step-2 example. Since we use the FE_Q class with a polynomial degree of 1, i.e. bilinear elements, this associates one degree of freedom with each vertex. While we're at generating output, let us also take a look at how many degrees of freedom are generated:

   dof_handler.distribute_dofs (fe);
   std::cout << "Number of degrees of freedom: "
             << dof_handler.n_dofs()
             << std::endl;

There should be one DoF for each vertex. Since we have a 32 times 32 grid, the number of DoFs should be 33 times 33, or 1089.

As we have seen in the previous example, we set up a sparsity pattern for the system matrix and tag those entries that might be nonzero. Compared to what we did in step-2, the only change is that instead of giving a magically obtained maximal number of nonzero entries per row, we now use a function in the DoFHandler class that can compute this number for us:

   sparsity_pattern.reinit (dof_handler.n_dofs(),
                            dof_handler.n_dofs(),
                            dof_handler.max_couplings_between_dofs());
   DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
   sparsity_pattern.compress();

Now the sparsity pattern is built and fixed (after `compress' has been called, you can't add nonzero entries anymore; the sparsity pattern is `sealed', so to say), we can initialize the matrix itself with it. Note that the SparsityPattern object does not hold the values of the matrix, it only stores the places where entries are. The entries are themselves stored in objects of type SparseMatrix, of which our variable system_matrix is one.

The distinction between sparsity pattern and matrix was made to allow several matrices to use the same sparsity pattern. This may not seem relevant, but when you consider the size which matrices can have, and that it may take some time to build the sparsity pattern, this becomes important in large-scale problems.

   system_matrix.reinit (sparsity_pattern);

The last thing to do in this function is to set the sizes of the right hand side vector and the solution vector to the right values:

   solution.reinit (dof_handler.n_dofs());
   system_rhs.reinit (dof_handler.n_dofs());
 }

LaplaceProblem::assemble_system

Now comes the difficult part: assembling matrices and vectors. In fact, this is not overly difficult, but it is something that the library can't do for you as for most of the other things in the functions above and below.

The general way to assemble matrices and vectors is to loop over all cells, and on each cell compute the contribution of that cell to the global matrix and right hand side by quadrature. The point to realize now is that we need the values of the shape functions at the locations of quadrature points on the real cell. However, both the finite element shape functions as well as the quadrature points are only defined on the unit cell. They are therefore of little help to us, and we will in fact hardly ever query information about finite element shape functions or quadrature points from these objects directly.

Rather, what is required is a way to map this data from the unit cell to the real cell. Classes that can do that are derived from the Mapping class, though one again often does not have to deal with them directly: many functions in the library can take a mapping object as argument, but when it is omitted they simply resort to the standard bilinear Q1 mapping. We will go this route, and not bother with it for the moment (we come back to this in step-10, step-11, and step-12).

So what we now have is a collection of three classes to deal with: finite element, quadrature, and mapping objects. That's too much, so there is one type of class that orchestrates information exchange between these three: the FEValues class. If given one instance of each three of these objects, it will be able to provide you with information about values and gradients of shape functions at quadrature points on a real cell.

Using all this, we will assemble the linear system for this problem in the following function:

 void LaplaceProblem::assemble_system () 
 {

Ok, let's start: we need a quadrature formula for the evaluation of the integrals on each cell. Let's take a Gauss formula with two quadrature points in each direction, i.e. a total of four points since we are in 2D. This quadrature formula integrates polynomials of degrees up to three exactly (in 1D). It is easy to check that this is sufficient for the present problem:

   QGauss<2>  quadrature_formula(2);

And we initialize the object which we have briefly talked about above. It needs to be told which finite element we want to use, and the quadrature points and their weights (jointly described by a Quadrature object). As mentioned, we use the implied Q1 mapping, rather than specifying one ourselves explicitly. Finally, we have to tell it what we want it to compute on each cell: we need the values of the shape functions at the quadrature points (for the right hand side (f,phi)), their gradients (for the matrix entries (grad phi_i, grad phi_j)), and also the weights of the quadrature points and the determinants of the Jacobian transformations from the unit cell to the real cells.

This list of what kind of information we actually need is given as a bitwise connection of flags as the third argument to the constructor of FEValues. Since these values have to be recomputed, or updated, every time we go to a new cell, all of these flags start with the prefix update_ and then indicate what it actually is that we want updated. The flag to give if we want the values of the shape functions computed is update_values; for the gradients it is update_gradients. The determinants of the Jacobians and the quadrature weights are always used together, so only the products (Jacobians times weights, or short JxW) are computed; since we need them, we have to list update_JxW_values as well:

   FEValues<2> fe_values (fe, quadrature_formula, 
                          update_values | update_gradients | update_JxW_values);

The advantage of this proceeding is that we can specify what kind of information we actually need on each cell. It is easily understandable that this approach can significant speed up finite element computations, compared to approaches where everything, including second derivatives, normal vectors to cells, etc are computed on each cell, regardless whether they are needed or not.

For use further down below, we define two short cuts for values that will be used very frequently. First, an abbreviation for the number of degrees of freedom on each cell (since we are in 2D and degrees of freedom are associated with vertices only, this number is four, but we rather want to write the definition of this variable in a way that does not preclude us from later choosing a different finite element that has a different number of degrees of freedom per cell, or work in a different space dimension).

Secondly, we also define an abbreviation for the number of quadrature points (here that should be four). In general, it is a good idea to use their symbolic names instead of hard-coding these number even if you know them, since you may want to change the quadrature formula and/or finite element at some time; the program will just work with these changes, without the need to change anything in this function.

The shortcuts, finally, are only defined to make the following loops a bit more readable. You will see them in many places in larger programs, and `dofs_per_cell' and `n_q_points' are more or less by convention the standard names for these purposes:

   const unsigned int   dofs_per_cell = fe.dofs_per_cell;
   const unsigned int   n_q_points    = quadrature_formula.size();

Now, we said that we wanted to assemble the global matrix and vector cell-by-cell. We could write the results directly into the global matrix, but this is not very efficient since access to the elements of a sparse matrix is slow. Rather, we first compute the contribution of each cell in a small matrix with the degrees of freedom on the present cell, and only transfer them to the global matrix when the computations are finished for this cell. We do the same for the right hand side vector. So let's first allocate these objects (these being local objects, all degrees of freedom are coupling with all others, and we should use a full matrix object rather than a sparse one for the local operations; everything will be transferred to a global sparse matrix later on):

   FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
   Vector<double>       cell_rhs (dofs_per_cell);

When assembling the contributions of each cell, we do this with the local numbering of the degrees of freedom (i.e. the number running from zero through dofs_per_cell-1). However, when we transfer the result into the global matrix, we have to know the global numbers of the degrees of freedom. When we query them, we need a scratch (temporary) array for these numbers:

   std::vector<unsigned int> local_dof_indices (dofs_per_cell);

Now for the loop over all cells. We have seen before how this works, so this should be familiar including the conventional names for these variables:

   DoFHandler<2>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (; cell!=endc; ++cell)
     {

We are now sitting on one cell, and we would like the values and gradients of the shape functions be computed, as well as the determinants of the Jacobian matrices of the mapping between unit cell and true cell, at the quadrature points. Since all these values depend on the geometry of the cell, we have to have the FEValues object re-compute them on each cell:

       fe_values.reinit (cell);

Next, reset the local cell's contributions contributions to global matrix and global right hand side to zero, before we fill them:

       cell_matrix = 0;
       cell_rhs = 0;

Then finally assemble the matrix: For the Laplace problem, the matrix on each cell is the integral over the gradients of shape function i and j. Since we do not integrate, but rather use quadrature, this is the sum over all quadrature points of the integrands times the determinant of the Jacobian matrix at the quadrature point times the weight of this quadrature point. You can get the gradient of shape function i at quadrature point q_point by using fe_values.shape_grad(i,q_point); this gradient is a 2-dimensional vector (in fact it is of type Tensor<1,dim>, with here dim=2) and the product of two such vectors is the scalar product, i.e. the product of the two shape_grad function calls is the dot product. This is in turn multiplied by the Jacobian determinant and the quadrature point weight (that one gets together by the call to fe_values.JxW). Finally, this is repeated for all shape functions phi_i and phi_j:

       for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int j=0; j<dofs_per_cell; ++j)
           for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
             cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
                                  fe_values.shape_grad (j, q_point) *
                                  fe_values.JxW (q_point));

We then do the same thing for the right hand side. Here, the integral is over the shape function i times the right hand side function, which we choose to be the function with constant value one (more interesting examples will be considered in the following programs). Again, we compute the integral by quadrature, which transforms the integral to a sum over all quadrature points of the value of the shape function at that point times the right hand side function (i.e. 1) times the Jacobian determinant times the weight of that quadrature point:

       for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
           cell_rhs(i) += (fe_values.shape_value (i, q_point) *
                           1 *
                           fe_values.JxW (q_point));

Now that we have the contribution of this cell, we have to transfer it to the global matrix and right hand side. To this end, we first have to find out which global numbers the degrees of freedom on this cell have. Let's simply ask the cell for that information:

       cell->get_dof_indices (local_dof_indices);

Then again loop over all shape functions i and j and transfer the local elements to the global matrix. The global numbers can be obtained using local_dof_indices[i]:

       for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int j=0; j<dofs_per_cell; ++j)
           system_matrix.add (local_dof_indices[i],
                              local_dof_indices[j],
                              cell_matrix(i,j));

And again, we do the same thing for the right hand side vector.

       for (unsigned int i=0; i<dofs_per_cell; ++i)
         system_rhs(local_dof_indices[i]) += cell_rhs(i);
     }

Now almost everything is set up for the solution of the discrete system. However, we have not yet taken care of boundary values (in fact, Laplace's equation without Dirichlet boundary values is not even uniquely solvable, since you can add an arbitrary constant to the discrete solution). We therefore have to do something about the situation.

For this, we first obtain a list of the degrees of freedom on the boundary and the value the shape function shall have there. For simplicity, we only interpolate the boundary value function, rather than projecting it onto the boundary. There is a function in the library which does exactly this: VectorTools::interpolate_boundary_values. Its parameters are (omitting parameters for which default values exist and that we don't care about): the DoFHandler object to get the global numbers of the degrees of freedom on the boundary; the component of the boundary where the boundary values shall be interpolated; the boundary value function itself; and the output object.

The component of the boundary is meant as follows: in many cases, you may want to impose certain boundary values only on parts of the boundary. For example, you may have inflow and outflow boundaries in fluid dynamics, or clamped and free parts of bodies in deformation computations of bodies. Then you will want to denote these different parts of the boundary by different numbers and tell the interpolate_boundary_values function to only compute the boundary values on a certain part of the boundary (e.g. the clamped part, or the inflow boundary). By default, all boundaries have the number `0', and since we have not changed that, this is still so; therefore, if we give `0' as the desired portion of the boundary, this means we get the whole boundary. If you have boundaries with kinds of boundaries, you have to number them differently. The function call below will then only determine boundary values for parts of the boundary.

The function describing the boundary values is an object of type Function or of a derived class. One of the derived classes is ZeroFunction, which describes (not unexpectedly) a function which is zero everywhere. We create such an object in-place and pass it to the interpolate_boundary_values function.

Finally, the output object is a list of pairs of global degree of freedom numbers (i.e. the number of the degrees of freedom on the boundary) and their boundary values (which are zero here for all entries). This mapping of DoF numbers to boundary values is done by the std::map class.

   std::map<unsigned int,double> boundary_values;
   VectorTools::interpolate_boundary_values (dof_handler,
                                             0,
                                             ZeroFunction<2>(),
                                             boundary_values);

Now that we got the list of boundary DoFs and their respective boundary values, let's use them to modify the system of equations accordingly. This is done by the following function call:

   MatrixTools::apply_boundary_values (boundary_values,
                                       system_matrix,
                                       solution,
                                       system_rhs);
 }

LaplaceProblem::solve

The following function simply solves the discretized equation. As the system is quite a large one for direct solvers such as Gauss elimination or LU decomposition, we use a Conjugate Gradient algorithm. You should remember that the number of variables here (only 1089) is a very small number for finite element computations, where 100.000 is a more usual number. For this number of variables, direct methods are no longer usable and you are forced to use methods like CG.

 void LaplaceProblem::solve () 
 {

First, we need to have an object that knows how to tell the CG algorithm when to stop. This is done by using a SolverControl object, and as stopping criterion we say: stop after a maximum of 1000 iterations (which is far more than is needed for 1089 variables; see the results section to find out how many were really used), and stop if the norm of the residual is below 1e-12. In practice, the latter criterion will be the one which stops the iteration:

   SolverControl           solver_control (1000, 1e-12);

Then we need the solver itself. The template parameters to the SolverCG class are the matrix type and the type of the vectors, but the empty angle brackets indicate that we simply take the default arguments (which are SparseMatrix<double> and Vector<double>):

   SolverCG<>              cg (solver_control);

Now solve the system of equations. The CG solver takes a preconditioner as its fourth argument. We don't feel ready to delve into this yet, so we tell it to use the identity operation as preconditioner:

   cg.solve (system_matrix, solution, system_rhs,
             PreconditionIdentity());

Now that the solver has done its job, the solution variable contains the nodal values of the solution function.

 }

LaplaceProblem::output_results

The last part of a typical finite element program is to output the results and maybe do some postprocessing (for example compute the maximal stress values at the boundary, or the average flux across the outflow, etc). We have no such postprocessing here, but we would like to write the solution to a file.

 void LaplaceProblem::output_results () const
 {

To write the output to a file, we need an object which knows about output formats and the like. This is the DataOut class, and we need an object of that type:

   DataOut<2> data_out;

Now we have to tell it where to take the values from which it shall write. We tell it which DoFHandler object to use, and the solution vector (and the name by which the solution variable shall appear in the output file). If we had more than one vector which we would like to look at in the output (for example right hand sides, errors per cell, etc) we would add them as well:

   data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (solution, "solution");

After the DataOut object knows which data it is to work on, we have to tell it to process them into something the back ends can handle. The reason is that we have separated the frontend (which knows about how to treat DoFHandler objects and data vectors) from the back end (which knows many different output formats) and use an intermediate data format to transfer data from the front- to the backend. The data is transformed into this intermediate format by the following function:

   data_out.build_patches ();

Now we have everything in place for the actual output. Just open a file and write the data into it, using GNUPLOT format (there are other functions which write their data in postscript, AVS, GMV, or some other format):

   std::ofstream output ("solution.gpl");
   data_out.write_gnuplot (output);
 }

LaplaceProblem::run

Finally, the last function of this class is the main function which calls all the other functions of the LaplaceProblem class. The order in which this is done resembles the order in which most finite element programs work. Since the names are mostly self-explanatory, there is not much to comment about:

 void LaplaceProblem::run () 
 {
   make_grid_and_dofs ();
   assemble_system ();
   solve ();
   output_results ();
 }

The main function

This is the main function of the program. Since the concept of a main function is mostly a remnant from the pre-object era in C/C++ programming, it often does not much more than creating an object of the top-level class and calling its principle function. This is what is done here as well:

 int main () 
 {
   LaplaceProblem laplace_problem;
   laplace_problem.run ();
 
   return 0;
 }

Results

The output of the program looks as follows:

Number of active cells: 1024
Total number of cells: 1365
Number of degrees of freedom: 1089
DEAL:cg::Starting value 0.121094
DEAL:cg::Convergence step 48 value 5.33692e-13

The first three lines is what we wrote to cout. The last two lines were generated without our intervention by the CG solver. The first two lines state the residual at the start of the iteration, while the last line tells us that the solver needed 47 iterations to bring the norm of the residual to 5.3e-13, i.e. below the threshold 1e-12 which we have set in the `solve' function. We will show in the next program how to suppress this output, which is sometimes useful for debugging purposes, but often clutters up the screen display.

Apart from the output shown above, the program generated the file solution.gpl, which is in GNUPLOT format. It can be viewed as follows: invoke GNUPLOT and enter the following sequence of commands at its prompt:

examples/step-3> gnuplot

        G N U P L O T
        Version 3.7 patchlevel 3
        last modified Thu Dec 12 13:00:00 GMT 2002
        System: Linux 2.6.11.4-21.10-default

        Copyright(C) 1986 - 1993, 1998 - 2002
        Thomas Williams, Colin Kelley and many others

        Type `help` to access the on-line reference manual
        The gnuplot FAQ is available from
        http://www.gnuplot.info/gnuplot-faq.html

        Send comments and requests for help to <info-gnuplot@dartmouth.edu>
        Send bugs, suggestions and mods to <bug-gnuplot@dartmouth.edu>


Terminal type set to 'x11'
gnuplot> set style data lines
gnuplot> splot "solution.gpl"
This produces the picture of the solution below left. Alternatively, you can order GNUPLOT to do some hidden line removal by the command
gnuplot> set hidden3d
to get the result at the right:

step-3.solution-1.png

step-3.solution-2.png

Possibilities for extensions

If you want to play around a little bit with this program, here are a few suggestions:

The plain program

(If you are looking at a locally installed deal.II version, then the program can be found at /build /buildd /deal.ii-6.2.1 /examples /step-3 /step-3.cc . Otherwise, this is only the path on some remote server.)

 / * Subversion Id:  step-3.cc 16263 2008-05-28 16:41:01Z kanschat  * /
 / * Author: Wolfgang Bangerth, University of Heidelberg, 1999 * /
 
 / *    Subversion Id:  step-3.cc 16263 2008-05-28 16:41:01Z kanschat        * /
 / *                                                                * /
 / *    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008 by the deal.II authors * /
 / *                                                                * /
 / *    This file is subject to QPL and may not be  distributed     * /
 / *    without copyright and license information. Please refer     * /
 / *    to the file deal.II/doc/license.html for the  text  and     * /
 / *    further information on this license.                        * /
 #include <grid/tria.h>
 #include <dofs/dof_handler.h>
 #include <grid/grid_generator.h>
 
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
 
 #include <fe/fe_q.h>
 
 #include <dofs/dof_tools.h>
 
 #include <fe/fe_values.h>
 #include <base/quadrature_lib.h>
 
 #include <base/function.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
 
 #include <numerics/data_out.h>
 #include <fstream>
 #include <iostream>
 
 using namespace dealii;
 class LaplaceProblem 
 {
   public:
     LaplaceProblem ();
 
     void run ();
     
   private:
     void make_grid_and_dofs ();
     void assemble_system ();
     void solve ();
     void output_results () const;
 
     Triangulation<2>     triangulation;
     FE_Q<2>              fe;
     DoFHandler<2>        dof_handler;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
 
     Vector<double>       solution;
     Vector<double>       system_rhs;
 };
 LaplaceProblem::LaplaceProblem () :
                 fe (1),
                dof_handler (triangulation)
 {}
 void LaplaceProblem::make_grid_and_dofs ()
 {
   GridGenerator::hyper_cube (triangulation, -1, 1);
   triangulation.refine_global (5);
   std::cout << "Number of active cells: "
            << triangulation.n_active_cells()
            << std::endl;
   std::cout << "Total number of cells: "
            << triangulation.n_cells()
            << std::endl;
   
   dof_handler.distribute_dofs (fe);
   std::cout << "Number of degrees of freedom: "
            << dof_handler.n_dofs()
            << std::endl;
 
   sparsity_pattern.reinit (dof_handler.n_dofs(),
                           dof_handler.n_dofs(),
                           dof_handler.max_couplings_between_dofs());
   DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
   sparsity_pattern.compress();
 
   system_matrix.reinit (sparsity_pattern);
 
   solution.reinit (dof_handler.n_dofs());
   system_rhs.reinit (dof_handler.n_dofs());
 }
 void LaplaceProblem::assemble_system () 
 {
   QGauss<2>  quadrature_formula(2);
   FEValues<2> fe_values (fe, quadrature_formula, 
                         update_values | update_gradients | update_JxW_values);
   
   const unsigned int   dofs_per_cell = fe.dofs_per_cell;
   const unsigned int   n_q_points    = quadrature_formula.size();
 
   FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
   Vector<double>       cell_rhs (dofs_per_cell);
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
   DoFHandler<2>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (; cell!=endc; ++cell)
     {
       fe_values.reinit (cell);
 
       cell_matrix = 0;
       cell_rhs = 0;
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
            cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
                                 fe_values.shape_grad (j, q_point) *
                                 fe_values.JxW (q_point));
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
          cell_rhs(i) += (fe_values.shape_value (i, q_point) *
                          1 *
                          fe_values.JxW (q_point));
 
       cell->get_dof_indices (local_dof_indices);
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
          system_matrix.add (local_dof_indices[i],
                             local_dof_indices[j],
                             cell_matrix(i,j));
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        system_rhs(local_dof_indices[i]) += cell_rhs(i);
     }
 
 
   std::map<unsigned int,double> boundary_values;
   VectorTools::interpolate_boundary_values (dof_handler,
                                            0,
                                            ZeroFunction<2>(),
                                            boundary_values);
   MatrixTools::apply_boundary_values (boundary_values,
                                      system_matrix,
                                      solution,
                                      system_rhs);
 }
 void LaplaceProblem::solve () 
 {
   SolverControl           solver_control (1000, 1e-12);
   SolverCG<>              cg (solver_control);
 
   cg.solve (system_matrix, solution, system_rhs,
            PreconditionIdentity());
 }
 void LaplaceProblem::output_results () const
 {
   DataOut<2> data_out;
   data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (solution, "solution");
   data_out.build_patches ();
 
   std::ofstream output ("solution.gpl");
   data_out.write_gnuplot (output);
 }
 void LaplaceProblem::run () 
 {
   make_grid_and_dofs ();
   assemble_system ();
   solve ();
   output_results ();
 }
 int main () 
 {
   LaplaceProblem laplace_problem;
   laplace_problem.run ();
 
   return 0;
 }

deal.II documentation generated on Sat Aug 15 16:51:06 2009 by doxygen 1.5.9