]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Tutorial/Laplace revisited. Code broken after last changes to deal.II.
authorschrage <schrage@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Mar 1999 14:37:39 +0000 (14:37 +0000)
committerschrage <schrage@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Mar 1999 14:37:39 +0000 (14:37 +0000)
git-svn-id: https://svn.dealii.org/trunk@1035 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/tutorial/chapter-3.laplace/assemble.html
deal.II/doc/tutorial/chapter-3.laplace/code/func.cc
deal.II/doc/tutorial/chapter-3.laplace/code/functions.h
deal.II/doc/tutorial/chapter-3.laplace/code/laplace.cc
deal.II/doc/tutorial/chapter-3.laplace/code/laplace.h
deal.II/doc/tutorial/chapter-3.laplace/code/main.cc
deal.II/doc/tutorial/chapter-3.laplace/index.html
deal.II/doc/tutorial/chapter-3.laplace/laplace.html
deal.II/doc/tutorial/chapter-3.laplace/main.html
deal.II/doc/tutorial/chapter-3.laplace/solution.html
deal.II/doc/tutorial/chapter-3.laplace/source.html [new file with mode: 0644]

index 31abf12950050223a20792d6d17225f5bbfdcbe5..b5ca1ef1c4c0a23d086d99843f909f6bb637a5b6 100644 (file)
@@ -1,13 +1,15 @@
 <!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&lt;2&gt;&amp;exact, const Function&lt;2&gt;&amp;)
-{
-</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&lt;2&gt; fevalues(fe_primal, qc_primal, UpdateFlags(update_gradients |
+  FEValues&lt;dim&gt; fevalues(dof->get_fe(), qc, UpdateFlags(update_gradients |
                                                   update_JxW_values));
-  FEFaceValues&lt;2&gt; ffvalues(fe_primal, qf_primal,
+  FEFaceValues&lt;dim&gt; ffvalues(dof->get_fe(), qf,
                           UpdateFlags(update_JxW_values | update_q_points));
 </code>
 </pre>
@@ -101,39 +99,55 @@ Integration is done locally. Therefore we need appropriate definitions for
 <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&lt;int&gt; indices(fe_primal.total_dofs);
-  dVector elvec(fe_primal.total_dofs);
+  vector&lt;int&gt; 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&lt;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&lt;2&gt;::active_cell_iterator c = dof_primal.begin_active()
-                                       ; c != dof_primal.end() ; ++c)
+  for (DoFHandler&lt;2&gt;::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&lt;qc_primal.n_quadrature_points;++k)
+    for (unsigned k=0;k&lt;qc.n_quadrature_points;++k)
     {
-      for (unsigned i=0;i&lt;fe_primal.total_dofs;++i)
+      for (unsigned i=0;i&lt;fe.total_dofs;++i)
       {
        const Point&lt;2&gt; dv = fevalues.shape_grad(i,k);
        
-       for (unsigned j=0;j&lt;fe_primal.total_dofs;++j)
+       for (unsigned j=0;j&lt;fe.total_dofs;++j)
        {
          const Point&lt;2&gt; du = fevalues.shape_grad(j,k);
          
@@ -146,20 +160,39 @@ discretized Laplace operator. <tt>qc_primal</tt> is a Gaussian quadrature.
     }
 </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&lt;fe.total_dofs;++i)
+    {
+      f(indices[i]) += elvec(i);
+            
+      for (unsigned j=0;j&lt;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>
@@ -167,19 +200,19 @@ This output is used by
 </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&lt;int,double&gt; boundary_values;
   DoFHandler&lt;2&gt;::FunctionMap dirichlet_bc;
-  BoundaryFct bfkt;
-  dirichlet_bc[0]=&amp;bfkt;
-  VectorTools&lt;2&gt;::interpolate_boundary_values(dof_primal,dirichlet_bc,fe_primal,boundary,boundary_values);
+  BoundaryFct bfct;
+  dirichlet_bc[0]=&amp;bfct;
+  VectorTools&lt;2&gt;::interpolate_boundary_values(dof,dirichlet_bc,fe,boundary,boundary_values);
   u.reinit(f);
   MatrixTools&lt;2&gt;::apply_boundary_values(boundary_values,A,u,f);  
 </code></pre>
@@ -188,11 +221,12 @@ 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.
+<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>
@@ -200,15 +234,15 @@ 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
-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>
 
@@ -221,7 +255,7 @@ the equation system which can then be solved.
 <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>
index a5c558ea79063bec02e3449d5679a4e90ff90b7f..8a739c6b2eed19a8384e849288d5b8456ccab65c 100644 (file)
@@ -1,6 +1,5 @@
 // $Id$
 
-// JS. 
 const char* funcversion = "Functions: $Revision$";
 
 #include "functions.h"
index d8b592d352517fda89ab641698ee9f8370432fe5..9b888fd69eb33e304c847406252fe00d42118903 100644 (file)
@@ -1,6 +1,5 @@
 // $Id$
 
-// JS.Wird das File ueberhaupt gebraucht ?
 
 #include <base/function.h>
 
index 0b4d0d5f56f11be061ce3cf7234a0294f619854a..da588e58acfcad9b8044cd3fb8625d36e68dba57 100644 (file)
@@ -20,120 +20,107 @@ const char* laplaceversion = "Laplace: $Revision$";
 #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);
          
@@ -144,74 +131,56 @@ Laplace::assemble_primal()
        }
       }
     }
-    // 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);
index 00cc5a017af1e3f7adcf6c60a706b1f60fcf7d0d..552ea2a1017ccce8cf09b5069ffa30fd2e991ca5 100644 (file)
@@ -7,16 +7,16 @@
 #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);
       }
 };
 
@@ -26,16 +26,16 @@ class Laplace
 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;
 
@@ -46,8 +46,8 @@ public:
   ~Laplace();
 
   void remesh(unsigned int global_refine = 0);
-  void assemble_primal();
-  void solve_primal();
+  void assemble();
+  void solve();
 
   void write_data(const char* name);
 };
index 73a8900e245bb971c78a263ba156efa2931c20aa..76eb23f5ebda3466fed790ef497efe7a28541eeb 100644 (file)
@@ -7,43 +7,23 @@
 #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
@@ -51,20 +31,12 @@ main(int argc, char** argv)
       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();
 }
index da405438607731bab8cf1e391c1097404fbbe281..39d39b01b6883559918f5dd9b769473819f73597 100644 (file)
@@ -3,17 +3,16 @@
 <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>
@@ -31,7 +30,7 @@ laplace problem
 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
@@ -50,7 +49,7 @@ where the overall design of the program and its classes are discussed</p>
 <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>
@@ -58,11 +57,9 @@ where we take a look at how DEAL is used</p>
 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>
@@ -73,7 +70,15 @@ and the boundary conditions are set
 </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>
@@ -85,7 +90,7 @@ where the problem is solved
 <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>
index bbc49f1be2ca9bc2d5fee4cef186070c7cbf83b6..cce0c9f99941eebf52e674e3086a4cffdf947f10 100644 (file)
@@ -22,10 +22,9 @@ Let's have a look at the class definition:
 <code>
 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:
 <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>
@@ -49,31 +48,19 @@ The constructor has the task of generating a triangulation, too.
 <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>
index 8093c61e6524a88afe5d6dbc577546ee48aeb706..a53c8dca34dbf72530b2b92a5bbd9f379e148fda 100644 (file)
@@ -12,7 +12,7 @@
 
 <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>
@@ -46,8 +46,8 @@ the function definition for our solution function are included:
 
 <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>
     
@@ -78,7 +78,10 @@ Laplace lap;
 
 <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 &lt; 3 ; ++step)   
@@ -94,13 +97,11 @@ for (unsigned step = 0; step &lt; 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>
@@ -108,7 +109,6 @@ hand side were not zero we would solve the Poisson equation instead.
 Finally the solution is written to a file.</p>
 <pre>
 <code>
-  sprintf(fname,"T%02d",step);    
   lap.write_data(fname);   
  }   
 }
index 5437cd91a844e4b73dff7ab352aebb9520dabb98..d728f9fad5f22a61933ddb570ca8b46a655e6c74 100644 (file)
 <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>
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 (file)
index 0000000..c55dbbb
--- /dev/null
@@ -0,0 +1,46 @@
+<!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>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.