</code>
</pre>
<p>
-which does exactly what it says. This function accepts a list of pairs
-of boundary indicators and the according functions and returns a list of
+which does exactly what it says. This function accepts a number
+describing which part of the boundary we want to consider
+and the according function which shall be interpolated on that portion
+of the boundary. It returns a list of
pairs of degrees of freedom numbers and values denoting the respective
Dirichlet boundary values.
</p>
element <code>fe</code>.
</p>
<pre class="example"><code>
+ // this will hold the boundary value for
+ // the degrees of freedom on the boundary
+ // part with number '0'
map<int,double> boundary_values;
-DoFHandler<2>::FunctionMap dirichlet_bc;
-BoundaryFct bfct;
-dirichlet_bc[0]=&bfct;
-VectorTools<2>::interpolate_boundary_values(dof,dirichlet_bc,fe,boundary,boundary_values);
-u.reinit(f.size());
-MatrixTools<2>::apply_boundary_values(boundary_values,A,u,f);
+ // and this is some user-defined function
+ // describing the boundary values
+BoundaryFunction<dim> boundary_function;
+ // interpolate the DoFs on boundary
+ // portion '0':
+VectorTools<dim>::interpolate_boundary_values
+ (dof_handler,
+ 0, boundary_function,
+ boundary_values);
+
+ // set solution vector to the right size
+u.reinit(dof_handler.n_dofs());
+ // we assume that the global matrix and
+ // right hand side have already been
+ // assembled.
+
+ // now insert the boundary values into
+ // the linear system of equations:
+MatrixTools<dim>::apply_boundary_values
+ (boundary_values,
+ A, u, f);
</code></pre>
<p>
-First, we need a few definitions:
-</p>
-<ul>
-<li>
-<code>boundary_values</code> maps DoF indices at the boundary computed by <code>interpolate_boundary_values</code> to their respective values.
-</li>
-<li><code>dirichlet_bc</code> maps boundary indicators to boundary functions, supplied by us. All boundary indicators are zero by default, therefore the
-above statement maps the same function to all the boundaries. The boundary functions compute the boundary values.
-</li>
-<li><code>bfct</code> is a function returning <code>cos(2*PI*x)*sin(2*PI*y)
-</code>, thereby supplying boundary values.
-</ul>
-<p>
-This may seem a bit confusing. What actually happens is the following:
-</p>
-<ol>
-<li><code>interpolate_boundary_values</code> takes the boundary functions
-<code>bfct</code>, its relation to boundaries <code>dirichlet_bc</code> and
-the triangulation <code>dof, fe</code> and returns a
-mapping <code>boundary_values</code> that maps values instead of functions
-to our boundaries. The function looks at all the boundaries the index
-of which is listed in <code>dirichlet_bc</code> (in this example, all
-the boundaries with indicator 0). All we
-ever need to do is specify the initial triangulation.
-</li>
-<li><code>apply_boundary_values</code> subsequently takes that mapping and
-our system of equations <tt>Au=f</tt> and inserts the boundary values into
-the system of equations which can then be solved.
-</li>
-</ol>
+You may have several different parts of the boundary which are
+identified by different boundary indicators. For example, this is
+commonly the case if you have inflow and outflow boundaries, or if
+your domain has inner boundaries, such as obstacles inside your
+domain. You may then want to associate different boundary functions
+with the different parts of the boundary, and you can do so by
+multiply calling <code>interpolate_boundary_values</code> with
+different boundary indicators (the <code>'0'</code> above) and
+respective functions. Since the function does not clear the contents
+of the map, you may use it for each of this successive calls and pass
+only the final object to <code>apply_boundary_values</code>.
<!-- Page Foot -->
<h2>Different kinds of matrices</h2>
<p>
-There are several kinds of matrices we will deal with:
+There are several kinds of matrices which we will deal with:
</p>
<ul>
<li>
-<a href="#full">Full matrices</a>, i.e. matrices where the majority of
-cells have non-zero values.
+<a href="#full">Full matrices</a>, i.e. matrices where we store all
+entries, which is useful if the majority of entries is non-zero.
</li>
<li>
<a href="#sparse">Sparse matrices</a>, i.e. matrices where the majority of
<h3><a name="full">The full matrix</a></h3>
<p>
-The full matrix is assumed to have mostly cells of non-zero value.
-Therefore it consumes a lot of memory, much more than the other matrix types.
-It is initialized using <code>FullMatrix::reinit(const unsigned
-int rows, const unsigned int cols)</code>.
+The full matrix is useful is we can assume that the we have mostly
+entries of non-zero value. Since we then store all entries, it
+consumes a lot of memory, much more than the other matrix types. It
+is initialized using <code>FullMatrix::reinit(const unsigned int rows,
+const unsigned int cols)</code>.
</p>
<p class="Example">
const unsigned int dim=2; // For example
-SparseMatrixStruct<double> sparsity_pattern;
+SparseMatrixStruct sparsity_pattern;
SparseMatrix<double> sparse_matrix;
-DoFHandler<dim> dof;
+DoFHandler<dim> dof_handler;
// Your degrees of freedom must already be distributed
-sparsity_pattern.reinit(dof.n_dofs(),dof.n_dofs(),dof.max_couplings_between_dofs());
+sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
// If your grid is locally refined, <a href="matrix_structure.html#sparsity_pattern">condense the hanging nodes</a> into the
// structure.
sparsity_pattern.compress();
-sparse_matrix.reinit(sparsity_pattern);
+// initialize the actual matrix with the sparsity pattern
+sparse_matrix.reinit (sparsity_pattern);
</code>
</pre>
<h2>Vector Generation</h2>
<p>
-Vector operations are supplied by the class
-<code><a href="http://gaia.iwr.uni-heidelberg.de/~deal/doc/auto/kdoc/lac/Vector.html">Vector</a></code>.
-The first and most important operation on a vector is its initialization
-using <code>void Vector::reinit(const usigned int N, const bool fast=false)</code>.
+Vector operations are supplied by the class <code><a
+href="http://gaia.iwr.uni-heidelberg.de/~deal/doc/auto/kdoc/lac/Vector.html">Vector</a></code>.
+The first and most important operation on a vector is its
+initialization using <code>void Vector::reinit(const usigned int N,
+const bool fast=false)</code>, which resets the vector's size to
+<code>N</code> and sets all elements to zero. If the second parameter
+<code>fast</code> is set true, this means that the entries will not be
+cleared; this is useful if the vector is being assigned to immediately
+after, in which case an initialization with zeroes would be
+redundant. Note that this mode of use is rather uncommon and mostly
+only useful inside linear solvers where it is important to save as
+much computing time as possible.
</p>
<p>
</pre>
+<h2>A word on memory allocation</h2>
+
+Since vectors and matrices are used <i>very</i> often, it is necessary
+to be rather conservative with operations that may use significant
+computing time. One of such operations is memory allocation and
+de-allocation. For this reason, the <code>reinit</code> functions of
+all vector and matrix classes are implemented in a way to only allow
+the allocated memory to grow, since this reduces the number of times
+allocation and de-allocation has to happen, and in addition reduces
+memory fragmentation.
+</p>
+<p>
+One immediate consequence is that if you write
+<pre class="example">
+<code>
+Vector<double> f;
+f.reinit(100);
+f.reinit(50);
+</code>
+</pre>
+then the second <code>reinit</code> will only have the effect of
+reducing the internal size of the vector, but the memory that is no
+more needed will not be released. Likewise, if you later write
+<pre class="example">
+<code>
+f.reinit(100);
+</code>
+</pre>
+this will again only set the size of the vector to a different value,
+but since enough memory is still allocated, no re-allocation has to
+happen. Only if the required number of entries exceeds the memory
+already allocated, will there be a re-allocation.
+</p>
+<p>
+All matrix and vector classes in <acronym>deal.II</acronym> work this
+way, and in fact this is also the way the standard C++ container
+classes like <code>vector</code> or <code>map</code> work.
+However, this is not very useful if you use vectors or matrices
+as member variables, since at times they may not be used at all if
+your program works on other objects or you are in other parts of the
+object. You would then like to let the vectors and matrices sleep with
+no memory presently lying around unused. For this purpose, all
+matrices and vectors support a mode of their <code>reinit</code>
+function which frees all memory that is presently allocated: just call
+it with parameters such that the new size is zero. For example, the
+following lines will release all memory from their objects:
+<pre class="example">
+<code>
+f.reinit(0);
+full_matrix.reinit(0,0);
+sparse_matrix_struct.reinit (0,0);
+sparse_matrix.reinit (sparse_matrix_struct);
+</code>
+</pre>
+Note that in the last line, you have to pass an empty sparsity pattern
+to the sparse matrix to tell it to release all of its memory.
+
+
+<h4>A word of caution</h4>
+
+It is noted that several vector and matrix classes have a
+<code>clear</code> function. These, however, work differently than
+might be expected from the analogous functions in the C++ standard
+library. There, <code>clear</code> removes all data from the objects
+and resets it into a virgin state. In <acronym>deal.II</acronym>,
+<code>clear</code> only resets the values of a metrix or vector to
+zero, but does not change the size and does no re-allocation
+either. This unfortunate name clash is due to historical reasons.
+
<!-- Page Foot -->
<hr>
<table class="navbar">
we shall now discuss how to fill them. You have to:
</p>
<ol>
- <li>Define a quadrature used for approximations. </li>
+ <li>Define a quadrature rule which will be used for integration of
+ the local cell matrices and vectors. </li>
<li><a href="#calcfe">Calculate the trial functions</a> for the
finite elements and their faces.
</li>
</p>
<p class="Example">
-<span class="example">Example:</span>
-The two lines below calculate trial functions for the two-dimensional finite element <code>fe</code> and
-for its faces using Gaussian quadrature. The first line initializes
-an object for the trial function of
-function for the finite element associated with the degree of freedom
-handler <code>dof</code>, telling it to
-update the values of the gradients and of the Jacobi determinant
-multiplied by a
-weight function given by the quadrature <code>qc</code>
-whenever <code>fe_values.reinit(fe)</code> is called. The second line
-does the same for the faces of the finite element, telling it to update the
-<code>JxW</code> values and the quadrature points.
+<span class="example">Example:</span> The two lines below calculate
+trial functions for the two-dimensional finite element <code>fe</code>
+and for its faces using Gaussian quadrature. The first line
+initializes an object for the trial and test functions of the
+finite element associated with the degree of freedom handler
+<code>dof_handler</code>, telling it to update the values of the gradients and
+of the Jacobi determinant multiplied by a weight function given by the
+quadrature <code>quadrature</code> whenever <code>fe_values.reinit(fe)</code>
+is called on a cell. The second line does the same for the faces of the finite
+element, telling it to update the <code>JxW</code> values and the
+quadrature points.
</p>
<pre class="example">
<code>
// Initialize the trial functions on the cell faces.
-FEValues<2> fevalues(fe, qc, UpdateFlags(update_gradients |
- update_JxW_values));
-FEFaceValues<2> ffvalues(fe, qf,
- UpdateFlags(update_JxW_values | update_q_points));
+FEValues<2> fe_values(fe, quadrature,
+ UpdateFlags(update_gradients | update_JxW_values));
+FEFaceValues<2> fe_face_values (fe, face_quadrature,
+ UpdateFlags(update_JxW_values | update_q_points));
</code>
</pre>
<pre class="example">
<code>
- vector<int> indices(fe.total_dofs);
- Vector<double> elvec(fe.total_dofs);
+ vector<int> dof_indices (fe.dofs_per_cell);
+ Vector<double> element_vector (fe.dofs_per_cell);
- FullMatrix<double> elmat(fe.total_dofs);
+ FullMatrix<double> element_matrix (fe.dofs_per_cell, fe.dofs_per_cell);
</code>
</pre>
<p class="example">
Next we traverse all the cells and integrate the Laplace problem using the
-discretized Laplace operator. <tt>qc</tt> is a Gaussian
-<code>Quadrature<dim></code>.
+discretized Laplace operator. <tt>quadrature</tt> is some quadrature
+rule, for example a Gaussian rule.
</p>
<p class="example">
-The outer loop traverses all the points of the quadrature <code>qc</code>.
-The inner two loops traverse the degrees of freedom of the finite element
-<code>fe</code> where <code>du</code> and <code>dv</code> are the gradients
-with respect to the quadrature points. <code>fevalues.JxW(k)</code> gives
-the Jacobi determinant multiplied by the weight of the quadrature point
-<code>k</code>. Taken together the line <br>
+Within the loop over all cells, the outer loop traverses all the
+points of the quadrature <code>quadrature</code>. The inner two loops
+traverse the degrees of freedom of the finite element <code>fe</code>
+where <code>du</code> and <code>dv</code> are the gradients of the
+trial and test functions at the present quadrature points.
+<code>fe_values.JxW(k)</code> gives
+the Jacobi determinant multiplied by the weight of the quadrature
+point <code>k</code>. Taken together the line <br>
<code>
-elmat(i,j) += fevalues.JxW(k) * du * dv;
+element_matrix(i,j) += fevalues.JxW(k) * du * dv;
</code>
<br>
gives the discretized Laplace operator.
<pre class="example">
<code>
// Integrate the problem locally...
-vector<int> indices(fe.total_dofs);
-Vector<double> elvec(fe.total_dofs);
-
-FullMatrix<double> elmat(fe.total_dofs);
-
-for (DoFHandler<2>::active_cell_iterator c = dof.begin_active()
- ; c != dof.end() ; ++c)
-{
- fevalues.reinit(c, stb);
- elmat.clear();
- elvec.clear();
- c->get_dof_indices(indices);
-
- // The loop over the quadrature points starts here. In effect, the
- // integration is performed in this loop.
- for (unsigned k=0;k<qc.n_quadrature_points;++k)
- {
- for (unsigned i=0;i<fe.total_dofs;++i)
- {
- const Point<2> dv = fevalues.shape_grad(i,k);
+vector<int> local_dof_indices (fe.dofs_per_cell);
+Vector<double> element_vector (fe.dofs_per_cell);
+FullMatrix<double> element_matrix (fe.dofs_per_cell);
-
- for (unsigned j=0;j<fe.total_dofs;++j)
+ // traverse all cells
+for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active();
+ cell != dof_handler.end(); ++cell)
+{
+ // compute gradients and some other values
+ // of the shape functions on this cell
+ fe_values.reinit (cell);
+ // set entries of local matrix and vector
+ // to zero
+ element_matrix.clear();
+ element_vector.clear();
+
+ // The loop over the quadrature points
+ // starts here. In effect, the integration
+ // is performed in this loop.
+ for (unsigned q_point=0; q_point<quadrature.n_quadrature_points; ++q_point)
+ for (unsigned i=0; i<fe.dofs_per_cell; ++i)
{
- const Point<2> du = fevalues.shape_grad(j,k);
+ const Point<dim> dv = fe_values.shape_grad (i, q_point);
+
+ for (unsigned j=0; j<fe.dofs_per_cell; ++j)
+ {
+ const Point<2> du = fe_values.shape_grad(j,q_point);
- // Perform the integration using the discretized Laplace operator.
- elmat(i,j) += fevalues.JxW(k)
- * du * dv
- ;
-
+ // Perform the integration using
+ // the weak form of the Laplace operator
+ element_matrix(i,j) += fevalues.JxW(k) * du * dv;
+ }
}
- }
- }
- // ...and insert the local matrix into the global one.
- for (unsigned i=0;i<fe.total_dofs;++i)
- {
- f(indices[i]) += elvec(i);
-
- for (unsigned j=0;j<fe.total_dofs;++j)
- {
- A.add(indices[i], indices[j], elmat(i,j));
- }
- }
-}
+ // get the global indices of the degrees
+ // of freedom on the present cell
+ cell->get_dof_indices (local_dof_indices);
+ // and insert the local matrix
+ // into the global one.
+ for (unsigned i=0; i<fe.dofs_per_cell; ++i)
+ for (unsigned j=0; j<fe.dofs_per_cell; ++j)
+ A.add (local_indices[i], local_indices[j],
+ element_matrix (i,j));
+}
</code>
</pre>
+The right hand side vector is set up analogously, but obviously you
+only need one loop over the shape functions, since the right hand side
+is a linear form instead of a bilinear one.
<!-- Page Foot -->