From: schrage Assembling the problem
@@ -16,79 +18,75 @@
In order to assemble the matrices we basically need to:
-
-void
-Laplace::assemble_primal(const Function<2>&exact, const Function<2>&)
-{
-
-
-
-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 max_couplings_between_dofs() returns the maximum
-number of couplings between degrees of freedom and allows DEAL
-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 max_couplings_between_dofs()
returns the maximum
+number of couplings between degrees of freedom, i.e. the number of elements in the
+matrix and allows deal.II to generate the matrix structure more
+efficiently, which in effect is done in the next line of code.
-Afterwards the hanging nodes are copied into the matrix, i.e. -the matrix is generated. +Afterwards a constraint matrix for the hanging nodes is created and they +are copied into the matrix structure.
- 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);
-The problem is of the form Au=f:
+The problem is of the form Au=f, we generate the matrix A with the
+structure given by matrix_structure
:
A.reinit(matrix_structure);
- f.reinit(dof_primal.n_dofs());
+ f.reinit(dof.n_dofs());
-
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 dof
+updating the values of the gradients and of the Jacobi determinant multiplied by a
+weight function given by the quadrature qc
.
- 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));
@@ -101,39 +99,55 @@ Integration is done locally. Therefore we need appropriate definitions for
- 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);
Next we traverse all the cells and integrate the Laplace problem using the
-discretized Laplace operator. qc_primal is a Gaussian quadrature.
+discretized Laplace operator. qc is a
+Quadrature<dim&rt;
.
+The outer loop traverses all the points of the quadrature qc
.
+The inner two loops traverse the degrees of freedom of the finite element
+fe
where du
and dv
are the gradients
+with respect to the quadrature points. fevalues.JxW(k)
gives
+the Jacobi determinant multiplied by the weight of the quadrature point
+k
. Taken together the line
+
+elmat(i,j) += fevalues.JxW(k) * du * dv;
+
+
+gives the discretized Laplace operator.
+
- 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);
@@ -146,20 +160,39 @@ discretized Laplace operator. qc_primal is a Gaussian quadrature.
}
+
+The insertion of the local matrix into the global one
+happens in the following piece of code (f
is the right hand
+vector of the problem, A
the problem matrix):
+
+
+ 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));
+ }
+ }
+
+
-There are two DEAL functions relevant for us at the moment: +There are two deal.II functions relevant for us at the moment:
-static_void interpolate_boundary_values(...)
+void VectorTools::interpolate_boundary_values(...)
-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.
@@ -167,19 +200,19 @@ This output is used by
-static void apply_boundary_values(...)
+void MatrixTools::apply_boundary_values(...)
-that inserts the proper boundary conditions into the equation system: +that inserts the proper boundary conditions into the system of equations:
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);
@@ -188,11 +221,12 @@ First, we need a few definitions:
boundary_values
maps boundary values computed by interpolate_boundary_values
to boundary indicators,i.e. to boundaries.
+boundary_values
maps DoF indices at the boundary computed by interpolate_boundary_values
to their respective values.
dirichlet_bc
maps boundary functions, supplied by us, to boundary indicators. The boundary functions compute the boundary values.
+dirichlet_bc
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.
bfkt
is a function returning sin(x)*sin(y)
+bfct
is a function returning cos(2*PI*x)*sin(2*PI*y)
, thereby supplying boundary values.
@@ -200,15 +234,15 @@ This may seem a bit confusing. What actually happens is the following:
interpolate_boundary_values
takes the boundary functions
-bfkt
, its relation to boundaries dirichlet_bc
and
-the triangulation dof_primal, fe_primal
and returns a
+bfct
, its relation to boundaries dirichlet_bc
and
+the triangulation dof, fe
and returns a
mapping boundary_values
that maps values instead of functions
to our boundaries. The function looks at all the boundaries. All we
ever need to do is specify the initial triangulation.
apply_boundary_values
subsequently takes that mapping and
-our equation system Au=f and inserts the boundary values into
-the equation system which can then be solved.
+our system of equations Au=f and inserts the boundary values into
+the system of equations which can then be solved.
-Last modified: Fri Feb 12, 1999 +Last modified: Tue Mar 9, 1999
diff --git a/deal.II/doc/tutorial/chapter-3.laplace/code/func.cc b/deal.II/doc/tutorial/chapter-3.laplace/code/func.cc index a5c558ea79..8a739c6b2e 100644 --- a/deal.II/doc/tutorial/chapter-3.laplace/code/func.cc +++ b/deal.II/doc/tutorial/chapter-3.laplace/code/func.cc @@ -1,6 +1,5 @@ // $Id$ -// JS. const char* funcversion = "Functions: $Revision$"; #include "functions.h" diff --git a/deal.II/doc/tutorial/chapter-3.laplace/code/functions.h b/deal.II/doc/tutorial/chapter-3.laplace/code/functions.h index d8b592d352..9b888fd69e 100644 --- a/deal.II/doc/tutorial/chapter-3.laplace/code/functions.h +++ b/deal.II/doc/tutorial/chapter-3.laplace/code/functions.h @@ -1,6 +1,5 @@ // $Id$ -// JS.Wird das File ueberhaupt gebraucht ? #includeWe will start with a simple example for a differential equation to be solved -with DEAL: the Laplace problem. We will try to solve the +with deal.II: 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.
@@ -31,7 +30,7 @@ laplace problem assemble matrices describing the boundary conditions
-where we take a look at how DEAL is used
+where we take a look at how deal.II is used-where a triangulation is generated and degrees of freedom are discussed
+where a triangulation is generated and degrees of freedom are discussed
where the problem is solved +
++contains link to all the source files. +
-Last modified: Mon 15 Feb 1999 +Last modified: Tue 9 Mar 1999
diff --git a/deal.II/doc/tutorial/chapter-3.laplace/laplace.html b/deal.II/doc/tutorial/chapter-3.laplace/laplace.html index bbc49f1be2..cce0c9f999 100644 --- a/deal.II/doc/tutorial/chapter-3.laplace/laplace.html +++ b/deal.II/doc/tutorial/chapter-3.laplace/laplace.html @@ -22,10 +22,9 @@ Let's have a look at the class definition:
class Laplace
{
- Function<2>& exact;
-protected:
+
Triangulation<2> tr;
- DoFHandler<2> dof_primal;
+ DoFHandler<2> dof;
dSMatrixStruct matrix_structure;
LapMatrix A;
@@ -35,11 +34,11 @@ protected:
-These few lines define several important elements: The right hand side of the equation
-exact, the triangulation tr, i.e. the grid, and a handler for the degrees
-of freedom for the finite elements dof_primal, all for the two-dimensional case.
-In addition three matrices are defined (the matrix A 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 tr
, i.e. the grid, and a handler for the degrees
+of freedom for the finite elements dof
, all for the two-dimensional case.
+In addition one matrix structure and two matrices are defined (the matrix A
defining our problem). Further explanation can be found in the
+chapter Assembling the problem. Note that
+in order to solve any problem at all with deal.II the definitions above are paramount.
@@ -49,31 +48,19 @@ The constructor has the task of generating a triangulation, too.
public:
- Laplace(Function<2>& solution);
+ Laplace();
~Laplace();
-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.
void remesh(unsigned int global_refine = 0);
- void assemble_primal(const Function<2>& boundary, const Function<2>& rhs);
- void solve_primal();
-
-
-
+ void assemble();
+ void solve();
};
diff --git a/deal.II/doc/tutorial/chapter-3.laplace/main.html b/deal.II/doc/tutorial/chapter-3.laplace/main.html
index 8093c61e65..a53c8dca34 100644
--- a/deal.II/doc/tutorial/chapter-3.laplace/main.html
+++ b/deal.II/doc/tutorial/chapter-3.laplace/main.html
@@ -12,7 +12,7 @@
The Laplace Problem
-Main Program: main.cc
+Main Program: main.cc
The main program has several functions:
@@ -46,8 +46,8 @@ the function definition for our solution function are included:
-#include "laplace.h"
-#include "functions.h"
+#include "laplace.h"
+#include "functions.h"
@@ -78,7 +78,10 @@ Laplace lap;
In addition we need to do some refinement (the command line argument was
-previously stored in the variable firstgrid.
+previously stored in the variable firstgrid. During the first
+execution
+of the loop global refinement is done firstgrid times, every
+following time it is done once more.
for (unsigned step = 0; step < 3 ; ++step)
@@ -94,13 +97,11 @@ for (unsigned step = 0; step < 3 ; ++step)
Problem assemblage and solution
-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.
- lap.assemble_primal(boundary, zero);
- lap.solve_primal();
+ lap.assemble();
+ lap.solve();
Data output
@@ -108,7 +109,6 @@ hand side were not zero we would solve the Poisson equation instead.
Finally the solution is written to a file.
- sprintf(fname,"T%02d",step);
lap.write_data(fname);
}
}
diff --git a/deal.II/doc/tutorial/chapter-3.laplace/solution.html b/deal.II/doc/tutorial/chapter-3.laplace/solution.html
index 5437cd91a8..d728f9fad5 100644
--- a/deal.II/doc/tutorial/chapter-3.laplace/solution.html
+++ b/deal.II/doc/tutorial/chapter-3.laplace/solution.html
@@ -11,6 +11,28 @@
Solving the problem
+
+Now that the problem matrix is assembled we need to define an appropriate
+solver. In this case we use a CG-solver with a maximum of 1000
+iterations and a threshold of 1e-10.
+
+
+
+void
+Laplace::solve_primal()
+{
+
+ SolverControl control(1000, 1.e-10);
+
+ SolverCG solver(control, mem);
+
+ solver.solve(A,u,f);
+
+ hanging_nodes.distribute(u);
+
+}
+
+
Back to the tutorial index
diff --git a/deal.II/doc/tutorial/chapter-3.laplace/source.html b/deal.II/doc/tutorial/chapter-3.laplace/source.html
new file mode 100644
index 0000000000..c55dbbbe81
--- /dev/null
+++ b/deal.II/doc/tutorial/chapter-3.laplace/source.html
@@ -0,0 +1,46 @@
+
+
+
+The Laplace Problem
+
+i
+
+
+
+The Source Code for the Solution of the Laplace Problem
+
+You can have a look at the complete source code for the solution of the
+Laplace problem by following the links below:
+
+- main.cc containing the main program which
+ sets some parameters and starts the solution process.
+
+- laplace.h and
+ laplace.cc which define the class
+
Laplace
that assembles the laplace problem.
+
+- functions.h and
+ func.cc which define the function that sets
+ the boundary conditions.
+
+- Makefile. Last but not least. You will
+ probably need to edit the variable
root
defined at the
+ top of the Makefile. This defines the location of your
+ DEAL directory.
+
+
+
+
+
+
+
+Jan Schrage
+
+Last modified: Mon 8 Mar 1999
+
+
+