/* $Id$ */
/* Version: $Name$ */
/* */
-/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors */
+/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2006 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. */
+
+ // @sect3{Many new include files}
+
// These include files are already
// known to you. They declare the
// classes which handle
- // triangulations and enumerate the
- // degrees of freedom.
+ // 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 which
- // create grids.
+ // 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 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 are the finite
- // element descriptions.
+ // In this file contains the description of
+ // the Lagrange interpolation finite element:
#include <fe/fe_q.h>
// And this file is needed for the
// assembling the matrix using
// quadrature on each cell. The
// classes declared in them will be
- // explained below.
+ // explained below:
#include <fe/fe_values.h>
#include <base/quadrature_lib.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
- // These include files are 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), and a class which
- // provides the solver with some
- // memory for temporary vectors.
+ // 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 <iostream>
- // Instead of the procedural
- // programming of previous examples,
- // we encapsulate everything into a
- // class for this program. The class
- // consists of functions which do
- // 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.
+ // @sect3{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:
- // This is the constructor:
LaplaceProblem ();
- // And the top-level function,
- // which is called from the
- // outside to start the whole
- // program (see the `main'
- // function at the bottom of this
- // file):
void run ();
- // Then there are some member
- // functions that mostly do what
- // their names suggest. Since
- // they do not need to be called
- // from outside, they are made
+ // 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 solve ();
void output_results () const;
- // And then we have the member
+ // And finally we have some member
// variables. There are variables
// describing the triangulation
- // and the numbering of the
+ // and the global numbering of the
// degrees of freedom (we will
// specify the exact polynomial
// degree of the finite element
Vector<double> system_rhs;
};
-
- // 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 specifies the polynomial
- // degree), and to associate the
- // dof_handler variable to the
- // triangulation we use. All the
- // other member variables of the
- // LaplaceProblem class have a
- // default constructor which does all
- // we want.
+ // @sect4{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)
{}
- // Now, the first thing we've got to
+ // @sect4{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
// 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 terminal cells. By terminal
- // we mean the cells on the finest
- // grid.
+ // 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;
- // We stress the adjective
- // `terminal' or `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:
+ // 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
- // previous examples. Since we use
- // the FEQ1 class, i.e. bilinear
- // elements, this associates one
- // degree of freedom with each
- // vertex.
+ // 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);
-
- // Now that we have the degrees of
- // freedom, we can take a look at
- // how many there are:
std::cout << "Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
// 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 sparse
- // matrix for the system matrix and
- // tag those entries that might be
- // nonzero. Since that has already
- // been done, we won't discuss the
- // next few lines:
+ // 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());
// `compress' has been called, you
// can't add nonzero entries
// anymore; the sparsity pattern is
- // `sealed', so to say), and we can
+ // `sealed', so to say), we can
// initialize the matrix itself
// with it. Note that the
// SparsityPattern object does
system_rhs.reinit (dof_handler.n_dofs());
}
+ // @sect4{LaplaceProblem::assemble_system}
+
// Now comes the difficult part:
// assembling matrices and
// 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
- // idea now is that since we only
- // need the finite element shape
- // functions on the quadrature points
- // of each cell, we don't need the
- // shape functions of the finite
- // element themselves any
- // more. Therefore, we won't deal
- // with the finite element object
- // `fe' (which was of type FEQ1), but
- // with another object which only
- // provides us with the values,
- // gradients, etc of the shape
- // functions at the quadrature
- // points. The objects which do this
- // are of type FEValues.
+ // 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). Since the integrands in
- // the matrix entries are quadratic
- // (in 1D), this is sufficient. The
- // same holds for 2D.
+ // 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 the finite element is that
- // we want to use, the quadrature
- // points and their
- // weights. 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, their gradients, and
- // also the weights of the
- // quadrature points and the
+ // 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. The
- // values of the shape functions
- // computed by specifying
- // update_values; the gradients are
- // done alike, using
- // update_gradients. The
- // determinants of the Jacobians
- // and the weights are always used
- // together, so only the products
- // (Jacobians times weights, or
- // short JxW) are computed; since
- // we also need them, we have to
- // list them as well. The advantage
- // of this proceeding is that we
- // calculate only what we
- // need. This optimatizes the
- // process of solving:
+ // 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,
- UpdateFlags(update_values |
- update_gradients |
- update_JxW_values));
-
- // For use further down below, we
- // define two short cuts 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). We also
- // define an abbreviation for the
- // number of quadrature points
- // (here that should be nine). 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 the matrix
- // assemblage.
+ 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 standard names
- // for these purposes.
+ // 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.n_quadrature_points;
- // 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
- // ell 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,
- // although access times are not so
- // problematic for them.
+ // 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);
// we transfer the result into the
// global matrix, we have to know
// the global numbers of the
- // degrees of freedom. When we get
+ // 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. You have seen before how
- // this works, so this should be
- // familiar to you:
- DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ // 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 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:
+ // 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);
- // Reset the values of the
- // contributions of this cell
- // to global matrix and global
- // right hand side to zero,
- // before we fill them.
+ // 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;
- // 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
+ // 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 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)
// 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 take into account
- // boundary values.
+ }
+
+
+ // 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 them onto
- // the boundary. There is a
- // function in the library which
- // does exactly this:
- // interpolate_boundary_values. Its
- // parameters are (omitting
- // parameters for which default
- // values exist which are
- // sufficient here): 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.
+ // 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,
- // are 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 got different
- // boundaries, you have to number
- // them differently and have to be
- // attentive when using the method
- // above.
+ // 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
- // described a function which is
- // zero everywhere. We create such
- // an object in-place and pass it
- // to the
- // interpolate_boundary_values
+ // 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
// here for all entries). This
// mapping of DoF numbers to
// boundary values is done by the
- // `map' class.
+ // ``std::map'' class.
std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
}
- // The following function simply
+ // @sect4{LaplaceProblem::solve}
+
+ // The following function simply
// solves the discretized
// equation. As the system is quite a
// large one for direct solvers such
// 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,
+ // 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 ()
{
- // We need to tell the algorithm
- // where to stop. This is done by
- // using a SolverControl object,
- // and as stopping criterion we
- // say: maximally 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.
+ // 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
- // here are the matrix type and the
- // type of the vectors, but the
- // empty angle brackets indicate
- // that we simply take the default
- // arguments.
+ // 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, but we don't
- // want to use one, so we tell it
- // to use the identity operation as
- // preconditioner.
+ // 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
}
+ // @sect4{LaplaceProblem::output_results}
+
// The last part of a typical finite
// element program is to output the
// results and maybe do some
// To write the output to a file,
// we need an object which knows
// about output formats and the
- // like. This is the DataOut class,
+ // 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 we
- // add the solution vector (and the
- // name by which it shall be
- // written to disk) to the list of
- // data that is to be written. 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:
+ // 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
// handle. The reason is that we
// have separated the frontend
// (which knows about how to treat
- // DoFHandler objects and data
+ // ``DoFHandler'' objects and data
// vectors) from the back end (which
- // knows several output formats)
+ // knows many different output formats)
// and use an intermediate data
// format to transfer data from the
// front- to the backend. The data
}
- // The following function 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:
+ // @sect4{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 ();
output_results ();
}
-
+
+ // @sect3{The ``main'' function}
// This is the main function of the
// program. Since the concept of a
// programming, it often does not
// much more than creating an object
// of the top-level class and calling
- // it principle function. This is
- // what is done here as well.
+ // its principle function. This is
+ // what is done here as well:
int main ()
{
LaplaceProblem laplace_problem;
laplace_problem.run ();
+
return 0;
}