<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
"http://www.w3.org/TR/REC-html40/strict.dtd">
<html>
-<head>
- <title>DEAL tutorial: the Laplace problem</title>
- <meta name="keyword" content="DEAL,DEAL tutorial">
- <link href="../dealtut.css" rel="StyleSheet" media="screen" type="text/css" title="DEAL tutorial">
- <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de>">
- <meta name="keywords" content="DEAL,DEAL tutorial">
-</head>
+ <head>
+ <title>deal.II tutorial: the Laplace problem</title>
+ <meta name="keyword" content="deal.II,deal.II tutorial">
+ <link rel="stylesheet" href="../dealtut.css" type="text/css">
+
+ <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de>">
+ <meta name="keywords" content="deal.II,deal.II tutorial,deal II">
+ </head>
+
<body>
<h1>Assembling the problem</h1>
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.II</em>
+ functions that reserve storage space for us.
</li>
<li>
<a href="#calcfe">Calculate the finite element trial functions</a>
</li>
<li>
Traverse all existing cells and <a href="#integrate">integrate the problem</a>
-using the discretized laplace operator
+using the discretized Laplace operator
</li>
<li>
-Traverse all the cell faces and <a href="#boundary">set the appropriate boundary conditions</a>
+<a href="#insert">Insert</a> the local matrices we have used into the global matrix using the
+appropriate <acronym>deal.II</acronym> functions
</li>
<li>
-Insert the local matrices we have used into the global matrix using the
-appropriate <em>DEAL</em> functions
+Traverse all the cell faces and <a href="#boundary">set the appropriate boundary conditions</a>
</li>
</ol>
<h2>...and how to do it</h2>
-<h3>Function parameters</h3>
-<pre>
-<code>
-void
-Laplace::assemble_primal(const Function<2>&exact, const Function<2>&)
-{
-</code>
-</pre>
-
<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
-of the degrees of freedom, i.e. the number of points of our discretization.
-The parameter <em>max_couplings_between_dofs()</em> returns the maximum
-number of couplings between degrees of freedom and allows <em>DEAL</em>
-to generate the matrix structure more efficiently, for most of its
-elements are zero.
+First we generate a structure for the storage of a sparse n times n matrix
+where n is the number of the degrees of freedom.
+The parameter <code>max_couplings_between_dofs()</code> returns the maximum
+number of couplings between degrees of freedom, i.e. the number of elements in the
+matrix and allows <acronym>deal.II</acronym> to generate the matrix structure more
+efficiently, which in effect is done in the next line of code.
</p>
<p>
-Afterwards the <em>hanging nodes</em> are copied into the matrix, i.e.
-the matrix is generated.
+Afterwards a constraint matrix for the <em>hanging nodes</em> is created and they
+are copied into the matrix structure.
</p>
<pre>
<code>
- matrix_structure.reinit(dof_primal.n_dofs(),dof_primal.n_dofs(),
- dof_primal.max_couplings_between_dofs());
- dof_primal.make_sparsity_pattern(matrix_structure);
+ matrix_structure.reinit(dof.n_dofs(),dof.n_dofs(),
+ dof.max_couplings_between_dofs());
+ dof.make_sparsity_pattern(matrix_structure);
hanging_nodes.clear();
- dof_primal.make_constraint_matrix(hanging_nodes);
+ dof.make_constraint_matrix(hanging_nodes);
hanging_nodes.condense(matrix_structure);
</code>
</pre>
<p>
-The problem is of the form <tt>Au=f</tt>:
+The problem is of the form <tt>Au=f</tt>, we generate the matrix <tt>A</tt> with the
+structure given by <code>matrix_structure</code>:
</p>
<pre>
<code>
A.reinit(matrix_structure);
- f.reinit(dof_primal.n_dofs());
+ f.reinit(dof.n_dofs());
</code>
</pre>
-<h3><a name="calcfe">Calculatinginite element trial functions</a></h3>
+<h3><a name="calcfe">Calculating finite element trial functions</a></h3>
<p>
The two lines below calculate trial functions for the finite elements and
-for their faces using Gaussian quadrature.
+for their faces using Gaussian quadrature. The first line calculates the trial
+function for the finite element associated with the degree of freedom <code>dof</code>
+updating the values of the gradients and of the Jacobi determinant multiplied by a
+weight function given by the quadrature <code>qc</code>.
</p>
<pre>
<code>
- FEValues<2> fevalues(fe_primal, qc_primal, UpdateFlags(update_gradients |
+ FEValues<dim> fevalues(dof->get_fe(), qc, UpdateFlags(update_gradients |
update_JxW_values));
- FEFaceValues<2> ffvalues(fe_primal, qf_primal,
+ FEFaceValues<dim> ffvalues(dof->get_fe(), qf,
UpdateFlags(update_JxW_values | update_q_points));
</code>
</pre>
<li>
an index vector that will allow us to reassemble the global matrix later on
</li>
-<li>a vector of doubles with the dimension of the total number of degrees of freedom</li>
+<li>a vector of doubles with the dimension of the number of degrees of freedom per cell
+</li>
<li>and a square matrix of doubles with the same dimension</li>
</ul>
<pre>
<code>
- vector<int> indices(fe_primal.total_dofs);
- dVector elvec(fe_primal.total_dofs);
+ vector<int> indices(fe.total_dofs);
+ dVector elvec(fe.total_dofs);
- dFMatrix elmat(fe_primal.total_dofs);
+ dFMatrix elmat(fe.total_dofs);
</code>
</pre>
<p>
Next we traverse all the cells and integrate the Laplace problem using the
-discretized Laplace operator. <tt>qc_primal</tt> is a Gaussian quadrature.
+discretized Laplace operator. <tt>qc</tt> is a
+<code>Quadrature<dim&rt;</code>.
</p>
+<p>
+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>
+<code>
+elmat(i,j) += fevalues.JxW(k) * du * dv;
+</code>
+<br>
+gives the discretized Laplace operator.
+</p>
+
<pre>
<code>
- for (DoFHandler<2>::active_cell_iterator c = dof_primal.begin_active()
- ; c != dof_primal.end() ; ++c)
+ 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);
+ c->get_dofindices(indices);
- for (unsigned k=0;k<qc_primal.n_quadrature_points;++k)
+ for (unsigned k=0;k<qc.n_quadrature_points;++k)
{
- for (unsigned i=0;i<fe_primal.total_dofs;++i)
+ for (unsigned i=0;i<fe.total_dofs;++i)
{
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.total_dofs;++j)
{
const Point<2> du = fevalues.shape_grad(j,k);
}
</code>
</pre>
+<p>
+<a name="insert">The insertion of the local matrix into the global one</a>
+happens in the following piece of code (<code>f</code> is the right hand
+vector of the problem, <code>A</code> the problem matrix):
+</p>
+<pre>
+<code>
+ 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));
+ }
+ }
+</code>
+</pre>
<h3><a name="boundary">Setting boundary conditions</a></h3>
<p>
-There are two <em>DEAL</em> functions relevant for us at the moment:
+There are two <acronym>deal.II</acronym> functions relevant for us at the moment:
</p>
<pre>
<code>
-static_void interpolate_boundary_values(...)
+void VectorTools::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
+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
+pairs of DoF numbers and values denoting the respective
Dirichlet boundary values.
</p>
<p>
</p>
<pre>
<code>
-static void apply_boundary_values(...)
+void MatrixTools::apply_boundary_values(...)
</code>
</pre>
<p>
-that inserts the proper boundary conditions into the equation system:
+that inserts the proper boundary conditions into the system of equations:
</p>
<P>
<pre><code>
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);
+ BoundaryFct bfct;
+ dirichlet_bc[0]=&bfct;
+ VectorTools<2>::interpolate_boundary_values(dof,dirichlet_bc,fe,boundary,boundary_values);
u.reinit(f);
MatrixTools<2>::apply_boundary_values(boundary_values,A,u,f);
</code></pre>
</p>
<ul>
<li>
-<code>boundary_values</code> maps boundary values computed by <code>interpolate_boundary_values</code> to boundary indicators,i.e. to boundaries.
+<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 functions, supplied by us, to boundary indicators. The boundary functions compute the boundary values.
+<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>bfkt</code> is a function returning <code>sin(x)*sin(y)
+<li><code>bfct</code> is a function returning <code>cos(2*PI*x)*sin(2*PI*y)
</code>, thereby supplying boundary values.
</ul>
<p>
</p>
<ol>
<li><code>interpolate_boundary_values</code> takes the boundary functions
-<code>bfkt</code>, its relation to boundaries <code>dirichlet_bc</code> and
-the triangulation <code>dof_primal, fe_primal</code> and returns a
+<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 <em>all</em> the boundaries. All we
ever need to do is specify the initial triangulation.
</li>
<li><code>apply_boundary_values</code> subsequently takes that mapping and
-our equation system <tt>Au=f</tt> and inserts the boundary values into
-the equation system which can then be solved.
+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>
<hr>
<address><a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
<p>
-Last modified: Fri Feb 12, 1999
+Last modified: Tue Mar 9, 1999
</p>
</body>
</html>
// $Id$
-// JS.
const char* funcversion = "Functions: $Revision$";
#include "functions.h"
// $Id$
-// JS.Wird das File ueberhaupt gebraucht ?
#include <base/function.h>
#include <base/quadrature_lib.h>
#include <numerics/matrices.h>
#include <numerics/vectors.h>
-
-#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/fullmatrix.h>
#include <cmath>
#include <fstream>
#include <iomanip>
#define PRIMEL FELinear<2>
-#define DUEL FEQuadraticSub<2>
// Finite Elements
-static PRIMEL fe_primal;
-// JS.static DUEL fe_dual;
+static PRIMEL fe;
// Quadrature formulae
-static QGauss2<2> qc_primal;
-static QGauss2<1> qf_primal;
-// static QGauss3<2> qc_dual;
-// static QGauss3<1> qf_dual;
+static QGauss2<2> qc;
+static QGauss2<1> qf;
static QGauss5<2> qc_integrate;
static QGauss5<1> qf_integrate;
StraightBoundary<2> stb;
-// JS.ist im Moment noch PureTransportSolution...
Laplace::Laplace()
- : dof_primal(&tr)
+ : dof(&tr)
{
- // JS.Triangulation generieren. Zellränder werden numeriert.
+ // Generate the initial triangulation
tr.create_hypercube(-1.,1.);
- // JS.Freiheitsgrade verteilen. D.h. Zellen numerieren.
- dof_primal.distribute_dofs(fe_primal);
+ // Distribute the degrees of freedom
+ dof.distribute_dofs(fe);
}
Laplace::~Laplace()
{}
-// JS.Gitter verfeinern.
+// Remesh the grid
void
Laplace::remesh(unsigned int steps)
{
if (tr.n_levels() <= 1)
{
- tr.refine_global(1); //JS.Lokal ist execute_coarsening_etc...
+ tr.refine_global(1); // refine globally
}
if (steps)
tr.refine_global(steps);
else
- tr.execute_coarsening_and_refinement();
-
- // JS. Freiheitsgrade neu verteilen.
- dof_primal.distribute_dofs(fe_primal);
- // JS. und dem Problem angemessener nochmal numerieren.
- dof_primal.renumber_dofs(Cuthill_McKee);
- deallog << "Cells " << tr.n_active_cells()
- << " PrimalDoFs " << dof_primal.n_dofs()
- << endl;
+ tr.execute_coarsening_and_refinement(); // refine locally
+
+ // redistribute the degrees of freedom...
+ dof.distribute_dofs(fe);
+ // ...and renumber them so they can be used more efficiently
+ dof.renumber_dofs(Cuthill_McKee);
}
// JS.Primales Problem zusammenstellen.
void
-Laplace::assemble_primal()
+Laplace::assemble()
{
- deallog << "Assembling primal problem" << endl;
- // JS. Platz für neue Matrix mit (?) (2x = quadratisch) Anzahl der Zellen,
- // Anzahl der Kopplungen (Matrix dünn besetzt, für effizientes Speichern)
- matrix_structure.reinit(dof_primal.n_dofs(),dof_primal.n_dofs(),
- dof_primal.max_couplings_between_dofs());
- // JS.Hängende Noden in die Matrix einbauen; d.h. Matrix generieren.
- dof_primal.make_sparsity_pattern(matrix_structure);
+ // Initialize the problem matrix, i.e. reserve storage space
+ matrix_structure.reinit(dof.n_dofs(),dof.n_dofs(),
+ dof.max_couplings_between_dofs());
+ // Generate the matrix structure
+ dof.make_sparsity_pattern(matrix_structure);
hanging_nodes.clear();
- dof_primal.make_constraint_matrix(hanging_nodes);
+ dof.make_constraint_matrix(hanging_nodes);
hanging_nodes.condense(matrix_structure);
- //JS.Problem der Form Au=f.
+ // The problem has the form Au=f.
A.reinit(matrix_structure);
- f.reinit(dof_primal.n_dofs());
+ f.reinit(dof.n_dofs());
- // JS.Ansatzfunktionen auf Zellrändern im Voraus berechnen aus
- // Effizienzgründen.
- FEValues<2> fevalues(fe_primal, qc_primal, UpdateFlags(update_gradients |
+ // Calculate the trial functions on the cell faces.
+ FEValues<2> fevalues(fe, qc, UpdateFlags(update_gradients |
update_JxW_values));
- FEFaceValues<2> ffvalues(fe_primal, qf_primal,
+ FEFaceValues<2> ffvalues(fe, qf,
UpdateFlags(update_JxW_values | update_q_points));
- //JS.Ab hier lokales Problem el... = Finites Element...
- //JS. Index für eine Zelle, für späteren Einbau in globale Matrix.
- vector<int> indices(fe_primal.total_dofs);
- dVector elvec(fe_primal.total_dofs);
+
+ // Integrate the problem locally...
+ vector<int> indices(fe.total_dofs);
+ Vector<float> elvec(fe.total_dofs);
- dFMatrix elmat(fe_primal.total_dofs);
+ FullMatrix<float> elmat(fe.total_dofs);
- // JS.Einmal alle Zellen durchlaufen
- for (DoFHandler<2>::active_cell_iterator c = dof_primal.begin_active()
- ; c != dof_primal.end() ; ++c)
+ 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);
- // JS.Integration des Problems. Diese Schleifenfolge für Effizienz.
- for (unsigned k=0;k<qc_primal.n_quadrature_points;++k)
+ for (unsigned k=0;k<qc.n_quadrature_points;++k)
{
- for (unsigned i=0;i<fe_primal.total_dofs;++i)
+ for (unsigned i=0;i<fe.total_dofs;++i)
{
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.total_dofs;++j)
{
const Point<2> du = fevalues.shape_grad(j,k);
}
}
}
- // JS.Lokale Matrix in globale einbauen.
- for (unsigned i=0;i<fe_primal.total_dofs;++i)
+ // ...and insert the local matrix into the global one.
+ for (unsigned i=0;i<fe.total_dofs;++i)
{
- // f(indices[i]) += elvec(i);
- f(indices[i]) = 0;
-
- for (unsigned j=0;j<fe_primal.total_dofs;++j)
+ f(indices[i]) += elvec(i);
+
+ for (unsigned j=0;j<fe.total_dofs;++j)
{
A.add(indices[i], indices[j], elmat(i,j));
}
}
}
- // JS. Randwerte.
+ // Insert the boundary conditions into the matrix.
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,
+ VectorTools<2>::interpolate_boundary_values(dof,dirichlet_bc,
+ fe,boundary,
boundary_values);
u.reinit(f);
MatrixTools<2>::apply_boundary_values(boundary_values,A,u,f);
-
-
- // JS.Hängende Noden einbauen.
- // hanging_nodes.condense(A);
- //hanging_nodes.condense(f);
}
-// JS. Primales Problem lösen.
+// Solve the primal problem.
void
-Laplace::solve_primal()
+Laplace::solve()
{
- deallog.push("Solve");
-
- // JS.Empfindlichkeit des Lösers einstellen.
+ // Solver control: max 1000 iterations, threshold 1e-10
SolverControl control(1000, 1.e-10);
- // JS. Löser definieren. modifiziertes cg-Verfahren,
- SolverCG<AdvMatrix, dVector> solver(control, mem);
-
- // JS.???
- // u.reinit(f);
+ // Define the solver.
+ SolverCG<AdvMatrix, Vector<float>> solver(control, mem);
- // JS.lösen.
solver.solve(A,u,f);
- // JS.???
hanging_nodes.distribute(u);
-
- deallog.pop();
}
-// JS. Datenausgabe im Gnuplot-Format.
+// Data output for gnuplot.
void Laplace::write_data(const char* name)
{
- deallog << "Writing gnuplot" << endl;
DataOut<2> out;
char fname[100];
- sprintf(fname,"P_%s",name);
-
{
ofstream gnuplot(fname);
out.clear_data_vectors();
- out.attach_dof_handler(dof_primal);
+ out.attach_dof_handler(dof);
out.add_data_vector(u,"solution","kg");
out.write_gnuplot (gnuplot, 1);
#include <grid/dof_constraints.h>
#include <grid/tria_boundary.h>
#include <lac/vector_memory.h>
-#include <lac/dvector.h>
-#include <lac/dsmatrix.h>
+#include <lac/vector.h>
+#include <lac/sparsematrix.h>
-class AdvMatrix :
- public dSMatrix
+class AdvMatrix :
+ public SparseMatrix<double>
{
public:
- void precondition(dVector& dst, const dVector& src) const
+ void precondition(Vector<double>& dst, const Vector<double>& src) const
{
- dSMatrix::precondition_SSOR(dst, src);
+ SparseMatrix::precondition_SSOR(dst, src);
}
};
protected:
Point<2> direction;
Triangulation<2> tr;
- DoFHandler<2> dof_primal;
+ DoFHandler<2> dof;
- dSMatrixStruct matrix_structure;
+ SparseMatrixStruct matrix_structure;
AdvMatrix A;
- dVector u;
- dVector z;
- dVector f;
+ Vector<float> u;
+ Vector<float> z;
+ Vector<float> f;
- PrimitiveVectorMemory<dVector> mem;
+ PrimitiveVectorMemory<float> mem;
ConstraintMatrix hanging_nodes;
~Laplace();
void remesh(unsigned int global_refine = 0);
- void assemble_primal();
- void solve_primal();
+ void assemble();
+ void solve();
void write_data(const char* name);
};
#include <fstream.h>
#include <stdlib.h>
-#include <base/logstream.h>
-#include <base/jobidentifier.h>
-
-ZeroFunction<2> zero;
-
char fname[30];
main(int argc, char** argv)
{
- // JS.Logfile erzeugen, als Stream konzipiert
-
- ofstream logfile("T");
- deallog.attach(logfile);
-
if (argc==1)
cerr << "Usage: " << argv[0] << "firstgrid\nUsing 3"
<< endl;
int firstgrid = 3;
-
- if (argc>=2) firstgrid = atoi(argv[1]);
- deallog << "Firstgrid " << firstgrid << endl;
+ if (argc>=2) firstgrid = atoi(argv[1]);
- // JS.Benötigte Funktionen zur Lösung des Problems
Laplace lap;
- // JS.Logstream ist ein stack, ab hier wird "Adaptive:" vor jede Zeile
- // gestellt
- deallog.push("Adaptive");
- deallog.depth_console(2);
-
for (unsigned step = 0; step < 3 ; ++step)
{
- deallog << "Step " << step << endl;
- // JS.Beim ersten Mal ein verfeinertes Grid erzeugen, firstgrid=Anzahl
- // der Verfeinerungen
+ // Generate a refined grid.
if (!step)
lap.remesh(firstgrid);
else
lap.remesh(1);
}
- deallog.push("Primal");
-
- // JS.exakte Lösung mit 0 auf der rechten Seite
- lap.assemble_primal();
- lap.solve_primal();
-
- // JS.ab hier kein "Adaptive:" mehr
- deallog.pop();
- sprintf(fname,"T%02d",step);
- // JS.Daten zurückschreiben.
+ // Assemble and solve the problem.
+ lap.assemble();
+ lap.solve();
+
+ // Data output.
lap.write_data(fname);
}
- // JS.Logfile sauber abschließen und Schluß
- deallog.pop();
- deallog.detach();
}
<html>
<head>
<title>The Laplace Problem</title>
- <link href="../dealtut.css" rel="StyleSheet" title="DEAL Tutorial">
-i <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de
->">
- <meta name="keywords" content="DEAL,DEAL tutorial">
+ <link href="../dealtut.css" rel="StyleSheet" title="deal.II Tutorial">
+ <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de>">
+ <meta name="keywords" content="deal.II,deal.II tutorial,deal II">
</head>
<body lang="en">
<h1>The Laplace Problem</h1>
<p>
We will start with a simple example for a differential equation to be solved
-with <strong>DEAL</strong>: the Laplace problem. We will try to solve the
+with <strong>deal.II</strong>: the Laplace problem. We will try to solve the
Laplace equation on a square where one of the four boundaries has a constant
and the other three zero potential.</p>
<p>
assemble matrices describing the boundary conditions
</li>
<li>
-let DEAL solve the problem
+let deal.II solve the problem
</li>
<li>
take care of the data output
<li>
<strong><a href="main.html">The main program</a></strong>
<p>
-where we take a look at how DEAL is used</p>
+where we take a look at how deal.II is used</p>
</li>
<li>
<strong><a href="laplace.html">The class Laplace</a></strong>
where some of the details of the problem generation and solution are
discussed</p>
</li>
-<li><p>
-<strong><a href="triangulation.html">Generating a
-triangulation</a></strong></p>
-<p>
-where a triangulation is generated and degrees of freedom are discussed</p>
+<li>
+<strong><a href="triangulation.html">Generating a triangulation</a></strong>
+<p>where a triangulation is generated and degrees of freedom are discussed</p>
</li>
<li>
<strong><a href="assemble.html">Assembling the problem matrix</a></strong>
</li>
<li>
<strong><a href="solution.html">Solving the problem</a></strong>
+<p>
where the problem is solved
+</p>
+</li>
+<li>
+<strong><a href="source.html">The source code</a></strong>
+<p>
+contains link to all the source files.
+</p>
</li>
</ol>
<hr>
<address>
<a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
<p>
-Last modified: Mon 15 Feb 1999
+Last modified: Tue 9 Mar 1999
</p>
</body>
</html>
<code>
class Laplace
{
- Function<2>& exact;
-protected:
+
Triangulation<2> tr;
- DoFHandler<2> dof_primal;
+ DoFHandler<2> dof;
dSMatrixStruct matrix_structure;
LapMatrix A;
<code>
</pre>
<p>
-These few lines define several important elements: The right hand side of the equation
-<em>exact</em>, the triangulation <em>tr</em>, i.e. the grid, and a handler for the degrees
-of freedom for the finite elements <em>dof_primal</em>, all for the two-dimensional case.
-In addition three matrices are defined (the matrix <em>A</em> defining our problem). Note that
-in order to solve any problem at all with DEAL the definitions above are paramount.
+These few lines define several important elements: The triangulation <code>tr</code>, i.e. the grid, and a handler for the degrees
+of freedom for the finite elements <code>dof</code>, all for the two-dimensional case.
+In addition one matrix structure and two matrices are defined (the matrix <code>A</code> defining our problem). Further explanation can be found in the
+chapter <a href="assemble.html#matrix">Assembling the problem</a>. Note that
+in order to solve any problem at all with <acronym>deal.II</acronym> the definitions above are paramount.
</p>
<p>
<pre>
<code>
public:
- Laplace(Function<2>& solution);
+ Laplace();
~Laplace();
</code>
</pre>
<p>
-The next few functions refine the grid - non-adaptively - assemble the primal problem and
+The next few functions refine the grid - non-adaptively - assemble the problem and
call the appropriate solver.
</p>
<pre>
<code>
void remesh(unsigned int global_refine = 0);
- void assemble_primal(const Function<2>& boundary, const Function<2>& rhs);
- void solve_primal();
-</code>
-</pre>
-<!--
- double result(const Function<2>& interior, const Function<2>& boundary);
- double estimate(const Function<2>& boundary, const Function<2>& rhs);
-
- void adapt();
-
- void write_data(const char* name);
-
- void fill_vector(dVector& v, const Function<2>& f) const;
--->
+ void assemble();
+ void solve();
};
</code>
</pre>
<h1>The Laplace Problem</h1>
-<h2>Main Program: <a href="../../../examples/01.laplace/main.cc">main.cc</a></h2>
+<h2>Main Program: <a href="code/main.cc">main.cc</a></h2>
<p>
The main program has several functions:
</p>
<pre>
<code>
-#include "<a href="../../../examples/01.laplace/laplace.h">laplace.h</a>"
-#include "<a href="../../../examples/01.laplace/functions.h">functions.h</a>"
+#include "<a href="code/laplace.h">laplace.h</a>"
+#include "<a href="code/functions.h">functions.h</a>"
</code>
</pre>
<p>
In addition we need to do some refinement (the command line argument was
-previously stored in the variable <var>firstgrid</var>.
+previously stored in the variable <var>firstgrid</var>. During the first
+execution
+of the loop global refinement is done <var>firstgrid</var> times, every
+following time it is done once more.
</p>
<pre><code>
for (unsigned step = 0; step < 3 ; ++step)
<h4><a name="laplace">Problem assemblage and solution</a></h4>
<p>
-Our class assembles and solves the primal problem; the solution is exact (as defined above)
-and the right hand side of the equation is zero (as defined above). If the right
-hand side were not zero we would solve the Poisson equation instead.
+Our class assembles and solves the primal problem.
</p>
<pre><code>
- lap.assemble_primal(boundary, zero);
- lap.solve_primal();
+ lap.assemble();
+ lap.solve();
</code></pre>
<h4><a name="output">Data output</a></h4>
Finally the solution is written to a file.</p>
<pre>
<code>
- sprintf(fname,"T%02d",step);
lap.write_data(fname);
}
}
<body>
<h1>Solving the problem</h1>
+<p>
+Now that the problem matrix is assembled we need to define an appropriate
+solver. In this case we use a <tt>CG</tt>-solver with a maximum of 1000
+iterations and a threshold of 1e-10.
+</p>
+
+<pre><code>
+void
+Laplace::solve_primal()
+{
+
+ SolverControl control(1000, 1.e-10);
+
+ SolverCG<AdvMatrix, dVector> solver(control, mem);
+
+ solver.solve(A,u,f);
+
+ hanging_nodes.distribute(u);
+
+}
+</code></pre>
+
<hr>
<p>
<a href="../index.html">Back to the tutorial index</a>
--- /dev/null
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN"
+ "http://www.w3.org/TR/REC-html40/strict.dtd">
+<html>
+<head>
+<title>The Laplace Problem</title>
+ <link href="../dealtut.css" rel="StyleSheet" title="DEAL Tutorial">
+i <meta name="author" content="Jan Schrage <schrage@gaia.iwr.uni-heidelberg.de
+>">
+ <meta name="keywords" content="DEAL,DEAL tutorial">
+</head>
+<body lang="en">
+<h1>The Source Code for the Solution of the Laplace Problem</h1>
+
+You can have a look at the complete source code for the solution of the
+Laplace problem by following the links below:
+<ul>
+<li><a href="code/main.cc">main.cc</a> containing the main program which
+ sets some parameters and starts the solution process.
+</li>
+<li><a href="code/laplace.h">laplace.h</a> and
+ <a href="code/laplace.cc">laplace.cc</a> which define the class
+ <code>Laplace</code> that assembles the laplace problem.
+</li>
+<li><a href="code/functions.h">functions.h</a> and
+ <a href="func.cc">func.cc</a> which define the function that sets
+ the boundary conditions.
+</li>
+<li><a href="code/Makefile">Makefile</a>. Last but not least. You will
+ probably need to edit the variable <code>root</code> defined at the
+ top of the Makefile. This defines the location of your
+ <acronym>DEAL</acronym> directory.
+</li>
+</ul>
+
+<hr>
+<p>
+<a href="../index.html">Back to the tutorial index</a>
+</p>
+<hr>
+<address>
+<a href="mailto:schrage@gaia.iwr.uni-heidelberg.de">Jan Schrage</a></address>
+<p>
+Last modified: Mon 8 Mar 1999
+</p>
+</body>
+</html>