]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Various more updates.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 29 Nov 1999 23:48:01 +0000 (23:48 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 29 Nov 1999 23:48:01 +0000 (23:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@1958 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/tutorial/chapter-1.elements/boundary.html
deal.II/doc/tutorial/chapter-1.elements/matrix_generation.html
deal.II/doc/tutorial/chapter-1.elements/rhs.html

index abd69dfe6203fe71628a6575f5e02d5e649ccba3..21ba9640e80b33a301b99c597194ee1550cc4a5c 100644 (file)
@@ -30,8 +30,10 @@ void VectorTools::interpolate_boundary_values(...)
 </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>
@@ -55,45 +57,44 @@ the handler for the degrees of freedom <code>dof</code> and the corresponding fi
 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&lt;int,double&gt; boundary_values;
-DoFHandler&lt;2&gt;::FunctionMap dirichlet_bc;
-BoundaryFct bfct;
-dirichlet_bc[0]=&amp;bfct;
-VectorTools&lt;2&gt;::interpolate_boundary_values(dof,dirichlet_bc,fe,boundary,boundary_values);
-u.reinit(f.size());
-MatrixTools&lt;2&gt;::apply_boundary_values(boundary_values,A,u,f);  
+                      // and this is some user-defined function
+                      // describing the boundary values
+BoundaryFunction&lt;dim&gt; boundary_function;
+                      // interpolate the DoFs on boundary
+                      // portion '0':
+VectorTools&lt;dim&gt;::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&lt;dim&gt;::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 -->
index a5814d232c54617b1d604bbc12cec51afa38f461..d76fd45fd42e932f373701a09e0141ceccefd1d7 100644 (file)
 
 <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 
@@ -37,10 +37,11 @@ cells has zero value; these zeroes are not stored in order to save memory.
 <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">
@@ -93,20 +94,23 @@ appropriate places. This example initializes a sparse square matrix structure.
 
 
 const unsigned int dim=2; // For example
-SparseMatrixStruct&lt;double&gt; sparsity_pattern;
+SparseMatrixStruct sparsity_pattern;
 SparseMatrix&lt;double&gt; sparse_matrix;
-DoFHandler&lt;dim&gt; dof;
+DoFHandler&lt;dim&gt; 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>
 
@@ -114,10 +118,18 @@ sparse_matrix.reinit(sparsity_pattern);
 <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>
@@ -148,6 +160,75 @@ f.reinit(dof.n_dofs());
 </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&lt;double&gt; 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">      
index d00651a1ffc813f7db2fa08b074bc6d3dfecc6e5..5aa31453a25028994dc22b6f314fdd6cfbbf73e7 100644 (file)
@@ -24,7 +24,8 @@ Having talked a lot about the initialization of matrices and vectors
 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>
@@ -38,26 +39,25 @@ we shall now discuss how to fill them. You have to:
 </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&lt;2&gt; fevalues(fe, qc, UpdateFlags(update_gradients |
-                                          update_JxW_values));
-FEFaceValues&lt;2&gt; ffvalues(fe, qf,
-       UpdateFlags(update_JxW_values | update_q_points));
+FEValues&lt;2&gt; fe_values(fe, quadrature,
+                   UpdateFlags(update_gradients | update_JxW_values));
+FEFaceValues&lt;2&gt; fe_face_values (fe, face_quadrature,
+                              UpdateFlags(update_JxW_values | update_q_points));
 </code>
 </pre>
 
@@ -85,28 +85,30 @@ an index vector that will allow us to reassemble the global matrix later on
 
 <pre class="example">
 <code>
-  vector&lt;int&gt; indices(fe.total_dofs);
-  Vector&lt;double&gt; elvec(fe.total_dofs);
+  vector&lt;int&gt; dof_indices (fe.dofs_per_cell);
+  Vector&lt;double&gt; element_vector (fe.dofs_per_cell);
   
-  FullMatrix&lt;double&gt; elmat(fe.total_dofs);
+  FullMatrix&lt;double&gt; 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&lt;dim&gt;</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.
@@ -115,55 +117,56 @@ gives the discretized Laplace operator.
 <pre class="example">
 <code>
 // Integrate the problem locally...
-vector&lt;int&gt; indices(fe.total_dofs);
-Vector&lt;double&gt; elvec(fe.total_dofs);
-  
-FullMatrix&lt;double&gt; elmat(fe.total_dofs);
-  
-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-&gt;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&lt;qc.n_quadrature_points;++k)
-  {
-    for (unsigned i=0;i&lt;fe.total_dofs;++i)
-    {
-       const Point&lt;2&gt; dv = fevalues.shape_grad(i,k);
+vector&lt;int&gt;        local_dof_indices (fe.dofs_per_cell);
+Vector&lt;double&gt;     element_vector    (fe.dofs_per_cell);
+FullMatrix&lt;double&gt; element_matrix    (fe.dofs_per_cell);
 
-       
-      for (unsigned j=0;j&lt;fe.total_dofs;++j)
+                        // traverse all cells  
+for (DoFHandler&lt;2&gt;::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&lt;quadrature.n_quadrature_points; ++q_point)
+    for (unsigned i=0; i&lt;fe.dofs_per_cell; ++i)
       {
-       const Point&lt;2&gt; du = fevalues.shape_grad(j,k);
+        const Point&lt;dim&gt; dv = fe_values.shape_grad (i, q_point);
+       
+        for (unsigned j=0; j&lt;fe.dofs_per_cell; ++j)
+          {
+           const Point&lt;2&gt; 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&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));
-    }
-  }
-}    
 
+                        // get the global indices of the degrees
+                        // of freedom on the present cell
+  cell-&gt;get_dof_indices (local_dof_indices);
+                        // and insert the local matrix
+                        // into the global one.
+  for (unsigned i=0; i&lt;fe.dofs_per_cell; ++i)
+    for (unsigned j=0; j&lt;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 -->

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.