<h1>Assembling the problem</h1>
<h2>What's to be done</h2>
+<p>
In order to assemble the matrices we basically need to:
+</p>
<ol>
-<li><a href=#matrix">Generate the matrices</a>, i.e. call the <em>DEAL</em> functions that reserve
-storage space for us.
+<li><a href="#matrix">Generate the matrices</a>, i.e. call the <em>DEAL</em> functions that reserve storage space for us.
</li>
<li>
<a href="#calcfe">Calculate the finite element trial functions</a>
<pre>
<code>
void
-Laplace::assemble_primal(const Function<2>& exact, const Function<2>&)
+Laplace::assemble_primal(const Function<2>&exact, const Function<2>&)
{
</code>
</pre>
-<h3><a name="matrix">Generating the matrix structures<a></h3>
+<h3><a name="matrix">Generating the matrix structures</a></h3>
<p>
First we generate an n times n square matrix where n is the number
</code>
</pre>
-<h3><a name+"calcfe">Calculatinginite element trial functions</a></h3>
+<h3><a name="calcfe">Calculatinginite element trial functions</a></h3>
<p>
The two lines below calculate trial functions for the finite elements and
for their faces using Gaussian quadrature.
</p>
<pre>
<code>
- FEValues<2> fevalues(fe_primal, qc_primal, UpdateFlags(update_gradients |
+ FEValues<2> fevalues(fe_primal, qc_primal, UpdateFlags(update_gradients |
update_JxW_values));
- FEFaceValues<2> ffvalues(fe_primal, qf_primal,
+ FEFaceValues<2> ffvalues(fe_primal, qf_primal,
UpdateFlags(update_JxW_values | update_q_points));
</code>
</pre>
<h3><a name="integrate">Integrating the problem</a></h3>
<p>
Integration is done locally. Therefore we need appropriate definitions for
+</p>
<ul>
<li>
an index vector that will allow us to reassemble the global matrix later on
<li>a vector of doubles with the dimension of the total number of degrees of freedom</li>
<li>and a square matrix of doubles with the same dimension</li>
</ul>
-</p>
+
<pre>
<code>
- vector<int> indices(fe_primal.total_dofs);
+ vector<int> indices(fe_primal.total_dofs);
dVector elvec(fe_primal.total_dofs);
dFMatrix elmat(fe_primal.total_dofs);
</p>
<pre>
<code>
- for (DoFHandler<2>::active_cell_iterator c = dof_primal.begin_active()
+ for (DoFHandler<2>::active_cell_iterator c = dof_primal.begin_active()
; c != dof_primal.end() ; ++c)
{
fevalues.reinit(c, stb);
elvec.clear();
c->get_dof_indices(indices);
- for (unsigned k=0;k<qc_primal.n_quadrature_points;++k)
+ for (unsigned k=0;k<qc_primal.n_quadrature_points;++k)
{
- for (unsigned i=0;i<fe_primal.total_dofs;++i)
+ for (unsigned i=0;i<fe_primal.total_dofs;++i)
{
- const Point<2> dv = fevalues.shape_grad(i,k);
+ const Point<2> dv = fevalues.shape_grad(i,k);
- for (unsigned j=0;j<fe_primal.total_dofs;++j)
+ for (unsigned j=0;j<fe_primal.total_dofs;++j)
{
- const Point<2> du = fevalues.shape_grad(j,k);
+ const Point<2> du = fevalues.shape_grad(j,k);
elmat(i,j) += fevalues.JxW(k)
* du * dv
<p>
There are two <em>DEAL</em> functions relevant for us at the moment:
+</p>
<pre>
<code>
static_void interpolate_boundary_values(...)
</code>
</pre>
+<p>
which does exactly what it says. This function returns a list of pairs
of boundary indicators and the according functions denoting the respective
Dirichlet boundary values.
</p>
<p>
This output is used by
+</p>
<pre>
<code>
static void apply_boundary_values(...)
</code>
</pre>
+<p>
that inserts the proper boundary conditions into the equation system:
</p>
<P>
<pre><code>
- map<int,double> boundary_values;
- DoFHandler<2>::FunctionMap dirichlet_bc;
+ map<int,double> boundary_values;
+ DoFHandler<2>::FunctionMap dirichlet_bc;
BoundaryFct bfkt;
- dirichlet_bc[0]=&bfkt;
- VectorTools<2>::interpolate_boundary_values(dof_primal,dirichlet_bc,fe_primal,boundary,boundary_values);
+ dirichlet_bc[0]=&bfkt;
+ VectorTools<2>::interpolate_boundary_values(dof_primal,dirichlet_bc,fe_primal,boundary,boundary_values);
u.reinit(f);
- MatrixTools<2>::apply_boundary_values(boundary_values,A,u,f);
+ MatrixTools<2>::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 boundary values computed by <code>interpolate_boundary_values</code> to boundary indicators,i.e. to boundaries.
</li>
<li><code>dirichlet_bc</code> maps boundary functions, supplied by us, to boundary indicators. The boundary functions compute the boundary values.
</li>
-<li><code>bfkt</bfkt> is a function returning <code>sin(x)*sin(y)
+<li><code>bfkt</code> is a function returning <code>sin(x)*sin(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>bfkt</code>, its relation to boundaries <code>dirichlet_bc</code> and
ever need to do is specify the initial triangulation.
</li>
<li><code>apply_boundary_values</code> subsequently takes that mapping and
-our equation system <code>Au=f</code> and inserts the boundary values into
+our equation system <tt>Au=f</tt> and inserts the boundary values into
the equation system which can then be solved.
</li>
</ol>
-</p>
+
<hr>
<p>
<body>
<h1>Generating a triangulation</h1>
+<p>
The initial triangulation is generated by the constructor of the class <em>Laplace</em>.
It is a hypercube with the coordinates along each axis ranging from -1 to 1. This first
triangulation consists of only one finite element.
+</p>
<pre>
<code>
Laplace::Laplace(Function<2>& solution)
tr.create_hypercube(-1.,1.);
</code>
</pre>
-<p>
-Right afterwards an iterator for the cells is defined (in effect a pointer to the first cell
-in this case) and all the cell's faces are numbered.
-</p>
-<pre>
-<code>
- Triangulation<2>::cell_iterator c = tr.begin();
- for (unsigned int i=0;i<GeometryInfo<2>::faces_per_cell;++i)
- {
- Triangulation<2>::face_iterator f = c->face(i);
- f->set_boundary_indicator(i+1);
- }
-</code>
-</pre>
+
+
<p>
Afterwards all the degrees of freedom, i.e. all the cell vertices, are renumbered.
This step is necessary whenever the triangulation has changed, e.g. after each refinement.
</p>
<hr>
-<p>
+
<address><a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
-</p>
<p>
Last modified: Fri Feb 12, 1999