]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Remove unfinished tutorial program.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 8 Oct 2011 15:42:30 +0000 (15:42 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 8 Oct 2011 15:42:30 +0000 (15:42 +0000)
git-svn-id: https://svn.dealii.org/branches/releases/Branch-7-1@24577 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-37/Makefile [deleted file]
deal.II/examples/step-37/doc/builds-on [deleted file]
deal.II/examples/step-37/doc/intro.dox [deleted file]
deal.II/examples/step-37/doc/kind [deleted file]
deal.II/examples/step-37/doc/results.dox [deleted file]
deal.II/examples/step-37/doc/tooltip [deleted file]
deal.II/examples/step-37/step-37.cc [deleted file]

diff --git a/deal.II/examples/step-37/Makefile b/deal.II/examples/step-37/Makefile
deleted file mode 100644 (file)
index 2c39a27..0000000
+++ /dev/null
@@ -1,143 +0,0 @@
-# $Id$
-
-
-# For the small projects Makefile, you basically need to fill in only
-# four fields.
-#
-# The first is the name of the application. It is assumed that the
-# application name is the same as the base file name of the single C++
-# file from which the application is generated.
-target = step-37
-
-# The second field determines whether you want to run your program in
-# debug or optimized mode. The latter is significantly faster, but no
-# run-time checking of parameters and internal states is performed, so
-# you should set this value to `on' while you develop your program,
-# and to `off' when running production computations.
-debug-mode = on
-
-
-# As third field, we need to give the path to the top-level deal.II
-# directory. You need to adjust this to your needs. Since this path is
-# probably the most often needed one in the Makefile internals, it is
-# designated by a single-character variable, since that can be
-# reference using $D only, i.e. without the parentheses that are
-# required for most other parameters, as e.g. in $(target).
-D = ../../
-
-
-# The last field specifies the names of data and other files that
-# shall be deleted when calling `make clean'. Object and backup files,
-# executables and the like are removed anyway. Here, we give a list of
-# files in the various output formats that deal.II supports.
-clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2
-
-
-
-
-#
-#
-# Usually, you will not need to change anything beyond this point.
-#
-#
-# The next statement tells the `make' program where to find the
-# deal.II top level directory and to include the file with the global
-# settings
-include $D/common/Make.global_options
-
-
-# Since the whole project consists of only one file, we need not
-# consider difficult dependencies. We only have to declare the
-# libraries which we want to link to the object file. deal.II has two
-# libraries: one for the debug mode version of the
-# application and one for optimized mode.
-libs.g   := $(lib-deal2.g)
-libs.o   := $(lib-deal2.o)
-
-
-# We now use the variable defined above to switch between debug and
-# optimized mode to select the set of libraries to link with. Included
-# in the list of libraries is the name of the object file which we
-# will produce from the single C++ file. Note that by default we use
-# the extension .g.o for object files compiled in debug mode and .o for
-# object files in optimized mode (or whatever local default on your
-# system is instead of .o)
-ifeq ($(debug-mode),on)
-  libraries = $(target).g.$(OBJEXT) $(libs.g)
-else
-  libraries = $(target).$(OBJEXT) $(libs.o)
-endif
-
-
-# Now comes the first production rule: how to link the single object
-# file produced from the single C++ file into the executable. Since
-# this is the first rule in the Makefile, it is the one `make' selects
-# if you call it without arguments.
-$(target)$(EXEEXT) : $(libraries)
-       @echo ============================ Linking $@
-       @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS)
-
-
-# To make running the application somewhat independent of the actual
-# program name, we usually declare a rule `run' which simply runs the
-# program. You can then run it by typing `make run'. This is also
-# useful if you want to call the executable with arguments which do
-# not change frequently. You may then want to add them to the
-# following rule:
-run: $(target)$(EXEEXT)
-       @echo ============================ Running $<
-       @./$(target)$(EXEEXT)
-
-
-# As a last rule to the `make' program, we define what to do when
-# cleaning up a directory. This usually involves deleting object files
-# and other automatically created files such as the executable itself,
-# backup files, and data files. Since the latter are not usually quite
-# diverse, you needed to declare them at the top of this file.
-clean:
-       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
-
-
-# Since we have not yet stated how to make an object file from a C++
-# file, we should do so now. Since the many flags passed to the
-# compiler are usually not of much interest, we suppress the actual
-# command line using the `at' sign in the first column of the rules
-# and write the string indicating what we do instead.
-./%.g.$(OBJEXT) :
-       @echo "==============debug========= $(<F)  ->  $@"
-       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
-./%.$(OBJEXT) :
-       @echo "==============optimized===== $(<F)  ->  $@"
-       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
-
-
-# The following statement tells make that the rules `run' and `clean'
-# are not expected to produce files of the same name as Makefile rules
-# usually do.
-.PHONY: run clean
-
-
-# Finally there is a rule which you normally need not care much about:
-# since the executable depends on some include files from the library,
-# besides the C++ application file of course, it is necessary to
-# re-generate the executable when one of the files it depends on has
-# changed. The following rule creates a dependency file
-# `Makefile.dep', which `make' uses to determine when to regenerate
-# the executable. This file is automagically remade whenever needed,
-# i.e. whenever one of the cc-/h-files changed. Make detects whether
-# to remake this file upon inclusion at the bottom of this file.
-#
-# If the creation of Makefile.dep fails, blow it away and fail
-Makefile.dep: $(target).cc Makefile \
-              $(shell echo $D/include/deal.II/*/*.h)
-       @echo ============================ Remaking $@
-       @$D/common/scripts/make_dependencies  $(INCLUDE) -B. $(target).cc \
-               > $@ \
-         || (rm -f $@ ; false)
-       @if test -s $@ ; then : else rm $@ ; fi
-
-# To make the dependencies known to `make', we finally have to include
-# them:
-include Makefile.dep
-
-
diff --git a/deal.II/examples/step-37/doc/builds-on b/deal.II/examples/step-37/doc/builds-on
deleted file mode 100644 (file)
index 42c2846..0000000
+++ /dev/null
@@ -1 +0,0 @@
-step-16
diff --git a/deal.II/examples/step-37/doc/intro.dox b/deal.II/examples/step-37/doc/intro.dox
deleted file mode 100644 (file)
index 9d488ec..0000000
+++ /dev/null
@@ -1,406 +0,0 @@
-<br>
-
-<i>
-This program was contributed by Katharina Kormann and Martin
-Kronbichler.
-
-This program is currently under construction.
-
-The algorithm for the matrix-vector product is built upon the report &quot;MPI
-parallelization of a cell-based matrix-vector product for finite elements. An
-application from quantum dynamics&quot; by Katharina Kormann, Uppsala
-University, June 2009.  
-</i>
-
-
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
-This example shows how to implement a matrix-free method, that is, a method
-that does not explicitly store the matrix elements, for a
-second-order Poisson equation with variable coefficients on a fairly 
-unstructured mesh representing a circle.
-
-<h3>Matrix-vector product implementation</h3>
-
-In order to find out how we can write a code that performs a matrix-vector
-product, but does not need to store the matrix elements, let us start at
-looking how some finite-element related matrix <i>A</i> is assembled:
-@f{eqnarray*}
-A = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T A_\mathrm{cell}
-P_\mathrm{cell,{loc-glob}}.
-@f}
-In this formula, the matrix <i>P</i><sub>cell,loc-glob</sub> is a rectangular
-matrix that defines the index mapping from local degrees of freedom in the
-current cell to the global degrees of freedom. The information from which this
-operator can be built is usually encoded in the <code>local_dof_indices</code>
-variable we have always used in the assembly of matrices.
-
-If we are to perform a matrix-vector product, we can hence use that
-@f{eqnarray*}
-y &=& A\cdot x = \left(\sum_{\text{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
-A_\mathrm{cell} P_\mathrm{cell,{loc-glob}}\right) \cdot x
-\\
-&=& \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
-A_\mathrm{cell} x_\mathrm{cell}
-\\
-&=& \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
-y_\mathrm{cell},
-@f}
-where <i>x</i><sub>cell</sub> are the values of <i>x</i> at the degrees of freedom
-of the respective cell, and  <i>x</i><sub>cell</sub> correspondingly for the result.
-A naive attempt to implement the local action of the Laplacian would hence be
-to use the following code:
-@code
-MatrixFree::vmult (Vector<double>       &dst,
-                  const Vector<double> &src) const
-{
-  dst = 0;
-
-  QGauss<dim>  quadrature_formula(fe.degree+1);
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                           update_gradients | update_JxW_values);
-
-  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-  const unsigned int   n_q_points    = quadrature_formula.size();
-
-  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-  Vector<double>       cell_src (dofs_per_cell),
-                      cell_dst (dofs_per_cell);
-
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
-  typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    {
-      cell_matrix = 0;
-      fe_values.reinit (cell);
-
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-            cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
-                                 fe_values.shape_grad(j,q_point) *
-                                 fe_values.JxW(q_point));
-
-      cell->get_dof_indices (local_dof_indices);
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        cell_src(i) = src(local_dof_indices(i));
-
-      cell_matrix.vmult (cell_dst, cell_src);
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        dst(local_dof_indices(i)) += cell_dst;
-    }
-}
-@endcode
-
-Here we neglected boundary conditions as well as any hanging nodes we may
-have, though neither would be very difficult to include using the
-ConstraintMatrix class. Note how we first generate the local matrix in the
-usual way. To form the actual product as expressed in the above formula, we
-read in the values of <code>src</code> of the cell-related degrees of freedom
-(the action of <i>P</i><sub>cell,loc-glob</sub>), multiply by the local matrix
-(the action of <i>A</i><sub>cell</sub>), and finally add the result to the
-destination vector <code>dst</code> (the action of
-<i>P</i><sub>cell,loc-glob</sub><sup>T</sup>, added over all the elements). It
-is not more difficult than that, in principle.
-
-While this code is completely correct, it is very slow. For every cell, we
-generate a local matrix, which takes three nested loops with as many
-elements as there are degrees of freedom on the actual cell to compute. The
-multiplication itself is then done by two nested loops, which means that it
-is much cheaper.
-
-One way to improve this is to realize that conceptually the local
-matrix can be thought of as the product of three matrices,
-@f{eqnarray*}
-A_\mathrm{cell} = B_\mathrm{cell}^T D_\mathrm{cell} B_\mathrm{cell},
-@f}
-where for the example of the Laplace operator the (<i>q</i>*dim+<i>d,i</i>)-th
-element of <i>B</i><sub>cell</sub> is given by
-<code>fe_values.shape_grad(i,q)[d]</code>. The matrix consists of
-<code>dim*n_q_points</code> rows and @p dofs_per_cell columns). The matrix
-<i>D</i><sub>cell</sub> is diagonal and contains the values
-<code>fe_values.JxW(q)</code> (or, rather, @p dim copies of it).
-
-Every numerical analyst learns in one of her first classes that for
-forming a product of the form
-@f{eqnarray*}
-A_\mathrm{cell}\cdot x_\mathrm{cell} = B_\mathrm{cell} D_\mathrm{cell} 
-                    B_\mathrm{cell}^T \cdot x_\mathrm{cell},
-@f}
-one should never form the matrix-matrix products, but rather multiply with the
-vector from right to left so that only three successive matrix-vector products
-are formed. To put this into code, we can write:
-@code
-...
-  for (; cell!=endc; ++cell)
-    {
-      fe_values.reinit (cell);
-
-      cell->get_dof_indices (local_dof_indices);
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        cell_src(i) = src(local_dof_indices(i));
-
-      temp_vector = 0;
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-        for (unsigned int d=0; d<dim; ++d)
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-           temp_vector(q_point*dim+d) += fe_values.shape_grad(i,q_point)[d] *
-                                         cell_src(i);
-
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-        for (unsigned int d=0; d<dim; ++d)
-         temp_vector(q_point*dim+d) *= fe_values.JxW(q_point);
-
-      cell_dst = 0;
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-          for (unsigned int d=0; d<dim; ++d)
-           cell_dst(i) += fe_values.shape_grad(i,q_point)[d] *
-                          temp_vector(q_point*dim+d);
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        dst(local_dof_indices(i)) += cell_dst(i);
-    }
-}
-@endcode
-
-This removed the three nested loops in the calculation of the local matrix
-(here the loop over <i>d</i> is a not really a loop, rather two or three
-operations). What happens is as follows: We first transform the vector of
-values on the local dofs to a vector of gradients on the quadrature
-points. In the second loop, we multiply these gradients by the integration
-weight. The third loop applies the second gradient (in transposed form), so
-that we get back to a vector of (Laplacian) values on the cell dofs.
-
-This improves the situation a lot and reduced the complexity of the product
-from something like $\mathcal {O}(\mathrm{dofs\_per\_cell}^3)$ to $\mathcal
-{O}(\mathrm{dofs\_per\_cell}^2)$. In fact, all the remainder is just to make
-a slightly more clever use of data in order to gain some extra speed. It does
-not change the code structure, though.
-
-The bottleneck in the above code is the operations done by the call
-<code>fe_values.reinit(cell)</code>, which take about as much time as the
-other steps together (at least if the mesh is unstructured; deal.II can
-recognize that the gradients are often unchanged on structured meshes). That
-is certainly not ideal and we would like to do better than this. What the
-reinit function does is to calculate the gradient in real space by
-transforming the gradient on the reference cell using the Jacobian of the
-transformation from real to reference cell. This is done for each basis
-function on the cell, for each quadrature point. The Jacobian does not depend on
-the basis function, but it is different on different quadrature points in
-general. The trick is now to factor out the Jacobian transformation and first
-apply the operation that leads us to <code>temp_vector</code> only with the
-gradient on the reference cell. That transforms the vector of values on the
-local dofs to a vector of gradients on the quadrature points. There, we first
-apply the Jacobian that we factored out from the gradient, then we apply the
-weights of the quadrature, and we apply with the transposed Jacobian for
-preparing the third loop which again uses the gradients on the unit cell.
-
-Let us again write this in terms of matrices. Let the matrix
-<i>B</i><sub>cell</sub> denote the cell-related gradient matrix, with each row
-containing the values of the quadrature points. It is constructed by a
-matrix-matrix product as
-@f{eqnarray*}
-B_\mathrm{cell} = J_\mathrm{cell} B_\mathrm{ref\_cell},
-@f}
-where <i>B</i><sub>ref_cell</sub> denotes the gradient on the reference cell
-and <i>J</i><sub>cell</sub> denotes the Jacobian
-transformation. <i>J</i><sub>cell</sub> is block-diagonal, and the blocks size
-is equal to the dimension of the problem. Each diagonal block is the Jacobian
-transformation that goes from the reference cell to the real cell.
-
-Putting things together, we find that
-@f{eqnarray*}
-A_\mathrm{cell} = B_\mathrm{cell}^T D B_\mathrm{cell} 
-               = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T 
-                 D_\mathrm{cell} 
-                 J_\mathrm{cell} B_\mathrm{ref\_cell},
-@f}
-so we calculate the product (starting the local product from the right)
-@f{eqnarray*}
-y_\mathrm{cell} = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T D J_\mathrm{cell} 
-B_\mathrm{ref\_cell} x_\mathrm{cell}, \quad
-y = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
-y_\mathrm{cell}.
-@f}
-@code
-...
-  FEValues<dim> fe_values_reference (fe, quadrature_formula,
-                                    update_gradients);
-  Triangulation<dim> reference_cell;
-  GridGenerator::hyper_cube(reference_cell, 0., 1.);
-  fe_values_reference.reinit (reference_cell.begin());
-
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                           update_inverse_jacobians | update_JxW_values);
-
-  for (; cell!=endc; ++cell)
-    {
-      fe_values.reinit (cell);
-
-      cell->get_dof_indices (local_dof_indices);
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        cell_src(i) = src(local_dof_indices(i));
-
-      temp_vector = 0;
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-        for (unsigned int d=0; d<dim; ++d)
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-           temp_vector(q_point*dim+d) +=
-             fe_values_reference.shape_grad(i,q_point)[d] * cell_src(i);
-
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-        {
-          // apply the Jacobian of the mapping from unit to real cell
-         Tensor<1,dim> temp;
-         for (unsigned int d=0; d<dim; ++d)
-           temp[d] = temp_vector(q_point*dim+d);
-         for (unsigned int d=0; d<dim; ++d)
-           {
-             double sum = 0;
-             for (unsigned int e=0; e<dim; ++e)
-               sum += fe_values.inverse_jacobian(q_point)[d][e] *
-                      temp[e];
-             temp_vector(q_point*dim+d) = sum;
-           }
-
-          // multiply with integration weight
-         for (unsigned int d=0; d<dim; ++d)
-           temp_vector(q_point*dim+d) *= fe_values.JxW(q_point);
-
-          // apply the transpose of the Jacobian of the mapping from unit
-         // to real cell
-         for (unsigned int d=0; d<dim; ++d)
-           temp[d] = temp_vector(q_point*dim+d);
-         for (unsigned int d=0; d<dim; ++d)
-           {
-             double sum = 0;
-             for (unsigned int e=0; e<dim; ++e)
-               sum += fe_values.inverse_jacobian(q_point)[e][d] *
-                      temp[e];
-             temp_vector(q_point*dim+d) = sum;
-           }
-        }
-
-      cell_dst = 0;
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-          for (unsigned int d=0; d<dim; ++d)
-           cell_dst(i) += fe_values_reference.shape_grad(i,q_point)[d] *
-                          temp_vector(q_point*dim+d);
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        dst(local_dof_indices(i)) += cell_dst(i);
-    }
-}
-@endcode
-
-Note how we create an additional FEValues object for the reference cell
-gradients and how we initialize it to the reference cell. The actual
-derivative data is then applied by the Jacobians (deal.II calls the Jacobian
-matrix from unit to real cell inverse_jacobian, because the transformation
-direction in deal.II is from real to unit cell).
-
-To sum up, we want to look at the additional costs introduced by this
-written-out matrix-matrix-vector product compared to a sparse matrix. To first
-approximation, we have increased the number of operations for the local
-matrix-vector product by a factor of 4 in 2D and 6 in 3D, i.e., two
-matrix-vector products with matrices that have <code>dim</code> as many
-columns than before (here we assume that the number of quadrature points is
-the same as the number of degrees of freedom per cell, which is usual for
-scalar problems). Then, we also need to keep in mind that we touch some
-degrees of freedom several times because they belong to several cells. This
-also increases computational costs. A realistic value compared to a sparse
-matrix is that we now have to perform about 10 times as many operations (a bit
-less in 2D, a bit more in 3D).
-
-The above is, in essence, what happens in the code below and if you have
-difficulties in understanding the implementation, you should try to first
-understand what happens in the code above. In the actual implementation
-there are a few more points done to be even more efficient, namely:
-<ul>
-  <li>We pre-compute the inverse of the Jacobian of the transformation and
-      store it in an extra array. This allows us to fuse the three
-      operations <i>J</i><sub>cell</sub><sup>T</sup><i> D</i><sub>cell</sub>
-      <i>J</i><sub>cell</sub> (apply Jacobian, multiply by 
-      weights, apply transposed
-      Jacobian) into one second-rank tensor that is also symmetric (so we
-      only need to store half the tensor).
-  <li>We work on several cells at once when we apply the gradients of the
-      unit cell (it is always the same matrix with the reference cell
-      data). This allows us to replace the matrix-vector product by a
-      matrix-matrix product (several vectors of cell-data form a matrix),
-      which enables a faster implementation. Obviously, we need some adapted
-      data structures for that, but it isn't too hard to provide that. What
-      is nice is that dense matrix-matrix products are close to today's
-      processors' peak performance if the matrices are neither too small nor
-      too large &mdash; and these operations are the most expensive part in
-      the implementation shown here. 
-</ul>
-
-The implementation of the matrix-free matrix-vector product shown in this
-tutorial is slower than a matrix-vector product using a sparse matrix for
-linear and quadratic elements, but on par with third order elements and faster
-for even higher order elements. An additional gain with this implementation is
-that we do not have to build the sparse matrix itself, which can also be quite
-expensive depending on the underlying differential equation.
-
-
-<h3>Combination with multigrid</h3>
-
-Above, we have gone to significant lengths to implement a matrix-vector
-product that does not actually store the matrix elements. In many user codes,
-however, one wants more than just performing some uncertain number of
-matrix-vector products &mdash; one wants to do as little of these operations
-as possible when solving linear equation systems. In theory, we could use the
-CG method without preconditioning; however, that would not be very
-efficient. Rather, one uses preconditioners for improving speed. On the other
-hand, most of the more frequently used preconditioners such as SSOR, ILU or
-algebraic multigrid (AMG) can now no longer be used here because their
-implementation requires knowledge of the elements of the system matrix.
-
-One solution is to use multigrid methods as shown in
-step-16. They are known to be very fast, and they are suitable for our
-purpose since they can be designed based purely on matrix-vector products. All
-one needs to do is to find a smoother that works with matrix-vector products
-only (our choice requires knowledge of the diagonal entries of the matrix,
-though). One such candidate would be a damped Jacobi iteration, but that is
-often not sufficiently good in damping high-frequency errors.
-A Chebyshev preconditioner, eventually, is what we use here. It can be
-seen as an extension of the Jacobi method by using Chebyshev polynomials. With
-degree zero, the Jacobi method with optimal damping parameter is retrieved,
-whereas higher order corrections improve the smoothing properties if some
-parameters are suitably chosen. The effectiveness of Chebyshev smoothing in
-multigrid has been demonstrated, e.g., in the article <i>M. Adams, M. Brezina,
-J. Hu, R. Tuminaro. Parallel multigrid smoothers: polynomial versus
-Gauss&ndash;Seidel, J. Comput. Phys. 188:593&ndash;610, 2003</i>. This
-publication also identifies one more advantage of Chebyshev smoothers that we
-exploit here, namely that they are easy to parallelize, whereas
-SOR/Gauss&ndash;Seidel smoothing relies on substitutions, which can often only
-be parallelized by working on diagonal sub-blocks of the matrix, which
-decreases efficiency.
-
-The implementation into the multigrid framework is then straightforward. This
-program is based on an earlier version of step-16 that demonstrated multigrid
-on uniformly refined grids. However, the present matrix-free techniques would
-obviously also apply to the adaptive meshes the current step-16 uses.
-
-
-<h3>The test case</h3>
-
-In order to demonstrate the capabilities of the method, we work on a rather
-general Poisson problem, based on a more or less unstructured mesh (where
-the Jacobians are different from cell to cell), higher order mappings to a
-curved boundary, and a non-constant coefficient in the equation. If we
-worked on a constant-coefficient case with structured mesh, we could
-decrease the operation count by a factor of 4 in 2D and 6 in 3D by building
-a local matrix (which is then the same for all cells), and doing the
-products as in the first developing step of the above code pieces.
diff --git a/deal.II/examples/step-37/doc/kind b/deal.II/examples/step-37/doc/kind
deleted file mode 100644 (file)
index c1d9154..0000000
+++ /dev/null
@@ -1 +0,0 @@
-techniques
diff --git a/deal.II/examples/step-37/doc/results.dox b/deal.II/examples/step-37/doc/results.dox
deleted file mode 100644 (file)
index 6e00f0b..0000000
+++ /dev/null
@@ -1,311 +0,0 @@
-<h1>Results</h1>
-
-<h3>Program output</h3>
-
-Since this example solves the same problem as step-5 (except for
-a different coefficient), we refer to the graphical output there. Here, we
-evaluate some aspects of the multigrid solver.
-
-When we run this program in 2D for quadratic ($Q_2$) elements, we get the
-following output:
-@code
-Cycle 0
-Number of degrees of freedom: 337
-System matrix memory consumption: 0.02573 MiB.
-Multigrid objects memory consumption: 0.05083 MiB.
-Convergence in 10 CG iterations.
-
-Cycle 1
-Number of degrees of freedom: 1313
-System matrix memory consumption: 0.09257 MiB.
-Multigrid objects memory consumption: 0.1794 MiB.
-Convergence in 10 CG iterations.
-
-Cycle 2
-Number of degrees of freedom: 5185
-System matrix memory consumption: 0.3553 MiB.
-Multigrid objects memory consumption: 0.6779 MiB.
-Convergence in 10 CG iterations.
-
-Cycle 3
-Number of degrees of freedom: 20609
-System matrix memory consumption: 1.397 MiB.
-Multigrid objects memory consumption: 2.645 MiB.
-Convergence in 10 CG iterations.
-
-Cycle 4
-Number of degrees of freedom: 82177
-System matrix memory consumption: 5.546 MiB.
-Multigrid objects memory consumption: 10.46 MiB.
-Convergence in 10 CG iterations.
-
-Cycle 5
-Number of degrees of freedom: 328193
-System matrix memory consumption: 22.11 MiB.
-Multigrid objects memory consumption: 41.65 MiB.
-Convergence in 10 CG iterations.
-@endcode
-
-As in step-16, we see that the number of CG iterations remains constant with
-increasing number of degrees of freedom. We can also see that the various
-objects we have to store for the multigrid method on the individual levels of
-our mesh together make up about twice as much as the matrix on the finest
-level.
-
-Not much changes if we run the
-program in three spatial dimensions, with the exception that the multilevel
-objects now take up comparatively less space (because in 3d, each level has
-only one eighth the number of cells of the next finer one, whereas in 2d this
-factor if one quarter):
-
-@code
-Cycle 0
-Number of degrees of freedom: 517
-System matrix memory consumption: 0.1001 MiB.
-Multigrid objects memory consumption: 0.1463 MiB.
-Convergence in 9 CG iterations.
-
-Cycle 1
-Number of degrees of freedom: 3817
-System matrix memory consumption: 0.6613 MiB.
-Multigrid objects memory consumption: 0.8896 MiB.
-Convergence in 10 CG iterations.
-
-Cycle 2
-Number of degrees of freedom: 29521
-System matrix memory consumption: 5.1 MiB.
-Multigrid objects memory consumption: 6.653 MiB.
-Convergence in 10 CG iterations.
-
-Cycle 3
-Number of degrees of freedom: 232609
-System matrix memory consumption: 40.4 MiB.
-Multigrid objects memory consumption: 52.24 MiB.
-Convergence in 11 CG iterations.
-
-Cycle 4
-Number of degrees of freedom: 1847617
-System matrix memory consumption: 322 MiB.
-Multigrid objects memory consumption: 415.1 MiB.
-Convergence in 11 CG iterations.
-@endcode
-
-
-<h3>Comparison with a sparse matrix</h3>
-
-In order to understand the capabilities of this class, we compare the memory
-consumption and execution (wallclock) time for assembly and 50 matrix-vector
-products (MV) on a <b>3D problem with one million unknowns</b> the classical
-sparse matrix implementation (SpM) and the MatrixFree implementation shown
-here (M-F). Both matrices are based on @p double %numbers. The program is run
-on a 2.8 GHz Opteron processor with the <a href="http://www.amd.com/acml">ACML
-BLAS</a>. We present results running on one core core and four cores,
-respectively. Moreover, we measure the time it takes to construct the
-individual matrices and filling them with data (@p setup and @p assemble
-functions). The sparse matrix is initialized using a
-CompressedSimpleSparsityPattern for calling the
-DoFTools::make_sparsity_pattern function, and then copied to a SparsityPattern
-object. The boundary nodes are eliminated using the ConstraintMatrix class, so
-that only elements that are actually nonzero are stored in the matrix.
-
-<table align="center" border="1">
-  <tr>
-    <th>&nbsp;</th>
-    <th colspan="2">Memory consumption</th>
-    <th colspan="2">Time assembly</th>
-    <th colspan="2">Time 50 MV, 1 CPU</th>
-    <th colspan="2">Time 50 MV, 4 CPUs</th>
-  </tr>
-  <tr>
-    <th align="center">element order</th>
-    <th align="center">SpM</th>
-    <th align="center">M-F</th>
-    <th align="center">SpM</th>
-    <th align="center">M-F</th>
-    <th align="center">SpM</th>
-    <th align="center">M-F</th>
-    <th align="center">SpM</th>
-    <th align="center">M-F</th>
-  </tr>
-  <tr>
-    <td align="center">1</td>
-    <td align="center">299 MiB</td>
-    <td align="center">394 MiB</td>
-    <td align="center">8.09 s</td>
-    <td align="center">3.43 s</td>
-    <td align="center">5.50 s</td>
-    <td align="center">22.4 s</td>
-    <td align="center">4.30 s</td>
-    <td align="center">11.0 s</td>
-  </tr>
-  <tr>
-    <td align="center">2</td>
-    <td align="center">698 MiB</td>
-    <td align="center">177 MiB</td>
-    <td align="center">12.43 s</td>
-    <td align="center">1.32 s</td>
-    <td align="center">12.0 s</td>
-    <td align="center">18.6 s</td>
-    <td align="center">9.10 s</td>
-    <td align="center">6.31 s</td>
-  </tr>
-  <tr>
-    <td align="center">3</td>
-    <td align="center">1295 MiB</td>
-    <td align="center">124 MiB</td>
-    <td align="center">41.1 s</td>
-    <td align="center">1.31 s</td>
-    <td align="center">21.2 s</td>
-    <td align="center">23.7 s</td>
-    <td align="center">16.0 s</td>
-    <td align="center">7.43 s</td>
-  </tr>
-  <tr>
-    <td align="center">4</td>
-    <td align="center">2282 MiB</td>
-    <td align="center">107 MiB</td>
-    <td align="center">117 s</td>
-    <td align="center">1.97 s</td>
-    <td align="center">40.8 s</td>
-    <td align="center">36.3 s</td>
-    <td align="center">19.7 s</td>
-    <td align="center">10.9 s</td>
-  </tr>
-  <tr>
-    <td align="center">5</td>
-    <td align="center">3597 MiB</td>
-    <td align="center">96.4 MiB</td>
-    <td align="center">510 s</td>
-    <td align="center">5.52 s</td>
-    <td align="center">75.7 s</td>
-    <td align="center">53.9 s</td>
-    <td align="center">29.3 s</td>
-    <td align="center">15.7 s</td>
-  </tr>
-  <tr>
-    <td align="center">6</td>
-    <td align="center">5679 MiB</td>
-    <td align="center">96.3 MiB</td>
-    <td align="center">2389 s</td>
-    <td align="center">26.1 s</td>
-    <td align="center">135 s</td>
-    <td align="center">79.1 s</td>
-    <td align="center">45.8 s</td>
-    <td align="center">24.3 s</td>
-  </tr>
-</table>
-
-There are a few interesting things with the %numbers in this table. 
-
-Firstly, we see the disappointing fact that for <b>linear elements</b> the
-MatrixFree class does actually consume more memory than a SparseMatrix with
-its SparsityPattern, despite the efforts made in this program. As mentioned
-earlier, this is mostly because the Transformation data is stored for every
-quadrature point. For each quadrature point, the transformation consists of
-six doubles, and there are about eight times as many quadrature points as
-there are degrees of freedom. In first approximation, this means that the
-matrix consumes 384 (= @p sizeof(double) * 6 * 8) bytes for each degree of
-freedom. On the other hand, the sparse matrix has a bandwidth of 27 or less,
-so each dof gives rise to at most 324 (= 27 * 12) bytes. A more clever
-implementation would try to compress the Jacobian transformation data, by
-exploiting similarities between the mappings within the cells, as well as from
-one cell to the next. This could dramatically reduce the memory requirements,
-and hence, increase the speed for lower-order implementations.
-
-Secondly, we observe that the memory requirements for a SparseMatrix grow
-quickly as the order of the elements increases. This is because there are
-increasingly many entries in each row, which exist due to more degrees of
-freedom that couple to each other. The matrix-free implementation does not
-suffer from this drawback. Here, the memory consumption decreases instead,
-since there are less DoFs that are shared among elements, which decreases the
-relative amount of quadrature points. Regarding the execution speed, we see
-that the matrix-free variant gets more competitive with higher order, and it
-does scale better when run on multiple processors (3.5 speedup with four
-processors compared to the serial case, compared to 2-2.5 speedup for the
-SparseMatrix). The advantage in %parallel scaling was expected, because the
-matrix-free variant is less memory-bound for higher order implementations, so
-that the additional computing power from many cores can better be exploited.
-
-A third thing, which is unrelated to this tutorial program, is the fact that
-standard matrix assembly gets really slow for high order elements. The %numbers
-shown here are based on the usual routines that many other tutorial programs
-make use of. A closer analysis of this shows that the cell data does not fit
-into cache anymore. One could circumvent this problem by writing the assembly
-as a matrix-matrix product, and using (cache-aware) BLAS implementations.
-
-For completeness, here comes a similar table for a <b>2D problem with 5.7
-million unknowns</b>. Since the excess in work for the matrix-free
-implementation is less compared to 3D, the implementation is more competitive
-for lower-order elements.
-
-<table align="center" border="1">
-  <tr>
-    <th>&nbsp;</th>
-    <th colspan="2">Memory consumption</th>
-    <th colspan="2">Time assembly</th>
-    <th colspan="2">Time 50 MV, 4 CPUs</th>
-  </tr>
-  <tr>
-    <th align="center">element order</th>
-    <th align="center">SpM</th>
-    <th align="center">M-F</th>
-    <th align="center">SpM</th>
-    <th align="center">M-F</th>
-    <th align="center">SpM</th>
-    <th align="center">M-F</th>
-  </tr>
-  <tr>
-    <td align="center">1</td>
-    <td align="center">659 MiB</td>
-    <td align="center">661 MiB</td>
-    <td align="center">18.8 s</td>
-    <td align="center">6.45 s</td>
-    <td align="center">11.0 s</td>
-    <td align="center">28.8 s</td>
-  </tr>
-  <tr>
-    <td align="center">2</td>
-    <td align="center">1119 MiB</td>
-    <td align="center">391 MiB</td>
-    <td align="center">15.6 s</td>
-    <td align="center">2.46 s</td>
-    <td align="center">17.1 s</td>
-    <td align="center">16.2 s</td>
-  </tr>
-  <tr>
-    <td align="center">3</td>
-    <td align="center">1711 MiB</td>
-    <td align="center">318 MiB</td>
-    <td align="center">17.4 s</td>
-    <td align="center">1.82 s</td>
-    <td align="center">23.1 s</td>
-    <td align="center">13.7 s</td>
-  </tr>
-  <tr>
-    <td align="center">4</td>
-    <td align="center">2434 MiB</td>
-    <td align="center">285 MiB</td>
-    <td align="center">24.2 s</td>
-    <td align="center">1.34 s</td>
-    <td align="center">31.1 s</td>
-    <td align="center">14.6 s</td>
-  </tr>
-  <tr>
-    <td align="center">5</td>
-    <td align="center">3289 MiB</td>
-    <td align="center">266 MiB</td>
-    <td align="center">35.9 s</td>
-    <td align="center">1.26 s</td>
-    <td align="center">29.6 s</td>
-    <td align="center">16.7 s</td>
-  </tr>
-  <tr>
-    <td align="center">6</td>
-    <td align="center">4274 MiB</td>
-    <td align="center">254 MiB</td>
-    <td align="center">58.0 s</td>
-    <td align="center">1.12 s</td>
-    <td align="center">35.9 s</td>
-    <td align="center">19.4 s</td>
-  </tr>
-</table>
diff --git a/deal.II/examples/step-37/doc/tooltip b/deal.II/examples/step-37/doc/tooltip
deleted file mode 100644 (file)
index d469776..0000000
+++ /dev/null
@@ -1 +0,0 @@
-Matrix-free methods. Multigrid.
diff --git a/deal.II/examples/step-37/step-37.cc b/deal.II/examples/step-37/step-37.cc
deleted file mode 100644 (file)
index bbc6b85..0000000
+++ /dev/null
@@ -1,1639 +0,0 @@
-/* $Id$ */
-/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2009 */
-
-/*    $Id$       */
-/*                                                                */
-/*    Copyright (C) 2009, 2010 by the deal.II authors                   */
-/*                                                                */
-/*    This file is subject to QPL and may not be  distributed     */
-/*    without copyright and license information. Please refer     */
-/*    to the file deal.II/doc/license.html for the  text  and     */
-/*    further information on this license.                        */
-
-
-                                // To start with the include files are more
-                                // or less the same as in step-16:
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/work_stream.h>
-
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
-
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/grid/grid_generator.h>
-
-#include <deal.II/multigrid/multigrid.h>
-#include <deal.II/multigrid/mg_dof_handler.h>
-#include <deal.II/multigrid/mg_dof_accessor.h>
-#include <deal.II/multigrid/mg_transfer.h>
-#include <deal.II/multigrid/mg_tools.h>
-#include <deal.II/multigrid/mg_coarse.h>
-#include <deal.II/multigrid/mg_smoother.h>
-#include <deal.II/multigrid/mg_matrix.h>
-
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/numerics/vectors.h>
-
-#include <fstream>
-#include <sstream>
-
-using namespace dealii;
-
-
-
-                                // @sect3{Equation data}
-
-                                // We define a variable coefficient function
-                                // for the Poisson problem. It is similar to
-                                // the function in step-5 but we use the form
-                                // $a(\mathbf x)=\frac{1}{0.1 + \|\bf x\|^2}$
-                                // instead of a discontinuous one. It is
-                                // merely to demonstrate the possibilities of
-                                // this implementation, rather than making
-                                // much sense physically.
-template <int dim>
-class Coefficient : public Function<dim>
-{
-  public:
-    Coefficient ()  : Function<dim>() {}
-
-    virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
-
-    virtual void value_list (const std::vector<Point<dim> > &points,
-                            std::vector<double>            &values,
-                            const unsigned int              component = 0) const;
-};
-
-
-
-template <int dim>
-double Coefficient<dim>::value (const Point<dim> &p,
-                               const unsigned int /*component*/) const
-{
-  return 1./(0.1+p.square());
-}
-
-
-
-template <int dim>
-void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
-                                  std::vector<double>            &values,
-                                  const unsigned int              component) const
-{
-  Assert (values.size() == points.size(),
-         ExcDimensionMismatch (values.size(), points.size()));
-  Assert (component == 0,
-         ExcIndexRange (component, 0, 1));
-
-  const unsigned int n_points = points.size();
-
-  for (unsigned int i=0; i<n_points; ++i)
-    values[i] = 1./(0.1+points[i].square());
-}
-
-
-
-                                // @sect3{Matrix-free implementation}
-
-                                // In this program, we want to make
-                                // use of the ability of deal.II to
-                                // runs things in %parallel if compute
-                                // resources are available. We will
-                                // follow the general framework laid
-                                // out in the @ref threads module and
-                                // use the WorkStream class to do
-                                // operations on the range of all
-                                // cells.
-                                //
-                                // To this end, we first have to have
-                                // a few declarations that we use for
-                                // defining the %parallel layout of
-                                // the vector multiplication function
-                                // with the WorkStream concept in the
-                                // Matrix-free class. These comprise
-                                // so-called scratch data that we use
-                                // for calculating cell-related
-                                // information, and copy data that is
-                                // eventually used in a separate
-                                // function for writing local data
-                                // into the global vector. The reason
-                                // for this split-up definition is
-                                // that many threads at a time can
-                                // execute the local multiplications
-                                // (and filling up the copy data),
-                                // but than that copy data needs to
-                                // be worked on by one process at a
-                                // time.
-namespace WorkStreamData
-{
-  template <typename number>
-  struct ScratchData
-  {
-      ScratchData ();
-      ScratchData (const ScratchData &scratch);
-      FullMatrix<number> solutions;
-  };
-
-  template<typename number>
-  ScratchData<number>::ScratchData ()
-                 :
-                 solutions ()
-  {}
-
-  template<typename number>
-  ScratchData<number>::ScratchData (const ScratchData &)
-                 :
-                 solutions ()
-  {}
-
-  template <typename number>
-  struct CopyData : public ScratchData<number>
-  {
-      CopyData ();
-      CopyData (const CopyData &scratch);
-      unsigned int first_cell;
-      unsigned int n_dofs;
-  };
-
-  template <typename number>
-  CopyData<number>::CopyData ()
-                 :
-                 ScratchData<number> ()
-  {}
-
-  template <typename number>
-  CopyData<number>::CopyData (const CopyData &)
-                 :
-                 ScratchData<number> ()
-  {}
-
-}
-
-
-
-                                // Next comes the implementation of the
-                                // matrix-free class. It provides some
-                                // standard information we expect for
-                                // matrices (like returning the dimensions
-                                // of the matrix), it implements
-                                // matrix-vector multiplications in several
-                                // forms, and it provides functions for
-                                // filling the matrix with data.
-                                //
-                                // We choose to make this class generic,
-                                // i.e., we do not implement the actual
-                                // differential operator (here: Laplace
-                                // operator) directly in this class.  We
-                                // instead let the actual transformation
-                                // (which happens on the level of quadrature
-                                // points, see the discussion in the
-                                // introduction) be a template parameter that
-                                // is implemented by another class. We then
-                                // only have to store a list of these objects
-                                // for each quadrature point on each cell in
-                                // a big list &ndash; we choose a
-                                // <code>Table<2,Transformation></code> data
-                                // format) &ndash; and call a transform
-                                // command of the @p Transformation
-                                // class. This template magic makes it easy
-                                // to reuse this MatrixFree class for other
-                                // problems that are based on a symmetric
-                                // operation without the need for substantial
-                                // changes.
-template <typename number, class Transformation>
-class MatrixFree : public Subscriptor
-{
-  public:
-    MatrixFree ();
-
-    void reinit (const unsigned int        n_dofs,
-                const unsigned int        n_cells,
-                const FullMatrix<double> &cell_matrix,
-                const unsigned int        n_points_per_cell);
-    void clear();
-
-    unsigned int m () const;
-    unsigned int n () const;
-    ConstraintMatrix & get_constraints ();
-
-    void set_local_dof_indices (const unsigned int               cell_no,
-                               const std::vector<unsigned int> &local_dof_indices);
-    void set_derivative_data (const unsigned int    cell_no,
-                             const unsigned int    quad_point,
-                             const Transformation &trans_in);
-
-    template <typename number2>
-    void vmult (Vector<number2> &dst,
-               const Vector<number2> &src) const;
-    template <typename number2>
-    void Tvmult (Vector<number2> &dst,
-                const Vector<number2> &src) const;
-    template <typename number2>
-    void vmult_add (Vector<number2> &dst,
-                   const Vector<number2> &src) const;
-    template <typename number2>
-    void Tvmult_add (Vector<number2> &dst,
-                    const Vector<number2> &src) const;
-
-    number el (const unsigned int row,
-              const unsigned int col) const;
-    void calculate_diagonal () const;
-
-    std::size_t memory_consumption () const;
-
-                                    // The private member variables of the
-                                    // @p MatrixFree class are a
-                                    // small matrix that does the
-                                    // transformation from solution values to
-                                    // quadrature points, a list with the
-                                    // mapping between local degrees of freedom
-                                    // and global degrees of freedom for each
-                                    // cell (stored as a two-dimensional array,
-                                    // where each row corresponds to one
-                                    // cell, and the columns within individual
-                                    // cells are the local degrees of freedom),
-                                    // the transformation variable for
-                                    // implementing derivatives, a constraint
-                                    // matrix for handling boundary conditions
-                                    // as well as a few other variables that
-                                    // store matrix properties.
-  private:
-    typedef std::vector<std::pair<unsigned int,unsigned int> >::const_iterator
-    CellChunkIterator;
-    template <typename number2>
-    void local_vmult (CellChunkIterator                    cell_range,
-                     WorkStreamData::ScratchData<number> &scratch,
-                     WorkStreamData::CopyData<number>    &copy,
-                     const Vector<number2>               &src) const;
-
-    template <typename number2>
-    void
-    copy_local_to_global (const WorkStreamData::CopyData<number> &copy,
-                         Vector<number2>                        &dst) const;
-
-    FullMatrix<number>      B_ref_cell;
-    Table<2,unsigned int>   indices_local_to_global;
-    Table<2,Transformation> derivatives;
-
-    ConstraintMatrix        constraints;
-
-    mutable Vector<number>  diagonal_values;
-    mutable bool            diagonal_is_calculated;
-
-    struct MatrixSizes
-    {
-       unsigned int n_dofs, n_cells;
-       unsigned int m, n;
-       unsigned int n_points, n_comp;
-       std::vector<std::pair<unsigned int,unsigned int> > chunks;
-    }  matrix_sizes;
-};
-
-
-
-                                // This is the constructor of the @p
-                                // MatrixFree class. All it does is to
-                                // subscribe to the general deal.II @p
-                                // Subscriptor scheme that makes sure that we
-                                // do not delete an object of this class as
-                                // long as it used somewhere else, e.g. in a
-                                // preconditioner.
-template <typename number, class Transformation>
-MatrixFree<number,Transformation>::MatrixFree ()
-               :
-               Subscriptor()
-{}
-
-
-
-                                // The next functions return the
-                                // number of rows and columns of the
-                                // global matrix (i.e. the dimensions
-                                // of the operator this class
-                                // represents, the point of this
-                                // tutorial program was, after all,
-                                // that we don't actually store the
-                                // elements of the rows and columns
-                                // of this operator). Since the
-                                // matrix is square, the returned
-                                // numbers are the same.
-template <typename number, class Transformation>
-unsigned int
-MatrixFree<number,Transformation>::m () const
-{
-  return matrix_sizes.n_dofs;
-}
-
-
-
-template <typename number, class Transformation>
-unsigned int
-MatrixFree<number,Transformation>::n () const
-{
-  return matrix_sizes.n_dofs;
-}
-
-
-
-                                // One more function that just returns an
-                                // %internal variable. Note that the user
-                                // will need to change this variable, so it
-                                // returns a non-constant reference to the
-                                // ConstraintMatrix.
-template <typename number, class Transformation>
-ConstraintMatrix &
-MatrixFree<number,Transformation>::get_constraints ()
-{
-  return constraints;
-}
-
-
-
-                                // The following function takes a vector of
-                                // local dof indices on cell level and writes
-                                // the data into the
-                                // @p indices_local_to_global field
-                                // in order to have fast access to it. It
-                                // performs a few sanity checks like whether
-                                // the sizes in the matrix are set
-                                // correctly. One tiny thing: Whenever we
-                                // enter this function, we probably make some
-                                // modification to the matrix. This means
-                                // that the diagonal of the matrix, which we
-                                // might have computed to have fast access to
-                                // those elements, is invalidated. We set the
-                                // respective flag to @p false.
-template <typename number, class Transformation>
-void MatrixFree<number,Transformation>::
-set_local_dof_indices (const unsigned int               cell_no,
-                      const std::vector<unsigned int> &local_dof_indices)
-{
-  Assert (local_dof_indices.size() == matrix_sizes.m,
-         ExcDimensionMismatch(local_dof_indices.size(),
-                              matrix_sizes.m));
-  for (unsigned int i=0; i<matrix_sizes.m; ++i)
-    {
-      Assert (local_dof_indices[i] < matrix_sizes.n_dofs, ExcInternalError());
-      indices_local_to_global(cell_no,i) = local_dof_indices[i];
-    }
-  diagonal_is_calculated = false;
-}
-
-
-
-                                // Next a function that writes the derivative
-                                // data on a certain cell and a certain
-                                // quadrature point to the array that keeps
-                                // the data around. Even though the array @p
-                                // derivatives stands for the majority of the
-                                // matrix memory consumption, it still pays
-                                // off to have that data around since it
-                                // would be quite expensive to manually
-                                // compute it every time we make a
-                                // matrix-vector product.
-template <typename number, class Transformation>
-void MatrixFree<number,Transformation>::
-set_derivative_data (const unsigned int cell_no,
-                    const unsigned int quad_point,
-                    const Transformation &trans_in)
-{
-  Assert (quad_point < matrix_sizes.n_points, ExcInternalError());
-  derivatives(cell_no,quad_point) = trans_in;
-  diagonal_is_calculated = false;
-}
-
-
-
-                                // Now finally to the central function of the
-                                // matrix-free class, implementing the
-                                // multiplication of the matrix with a
-                                // vector. This function does not actually
-                                // work on all cells of a mesh, but only the
-                                // subset of cells specified by the first
-                                // argument @p cell_range. Since this
-                                // function operates similarly irrespective
-                                // on which cell chunk we are sitting, we can
-                                // call it simultaneously on many processors,
-                                // but with different cell range data.
-                                //
-                                // The goal of this function is to provide
-                                // the multiplication of a vector with the
-                                // local contributions of a set of cells. As
-                                // mentioned in the introduction, if we were
-                                // to deal with a single cell, this would
-                                // amount to performing the product
-                                // @f{eqnarray*}
-                                // P^T_\mathrm{cell,local-global} A_\mathrm{cell}
-                                // P_\mathrm{cell,local-global} x
-                                // @f}
-                                // where
-                                // @f{eqnarray*}
-                                // A_\mathrm{cell} =
-                                // B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T
-                                // D_\mathrm{cell}
-                                // J_\mathrm{cell} B_\mathrm{ref\_cell}
-                                // @f}
-                                // and <i>P</i><sub>cell,local-global</sub>
-                                // is the transformation from local to global
-                                // indices.
-                                //
-                                // To do this, we would have to do the
-                                // following steps:
-                                // <ol>
-                                //   <li> Form $x_\mathrm{cell} =
-                                //   P_\mathrm{cell,local-global} x$. This is
-                                //   done by using the command
-                                //   ConstraintMatrix::get_dof_values.
-                                //   <li> Form $x_1 = B_\mathrm{ref\_cell}
-                                //   x_\mathrm{cell}$. The vector
-                                //   <i>x</i><sub>1</sub> contains the
-                                //   reference cell gradient to the local
-                                //   cell vector.
-                                //   <li> Form $x_2 = J_\mathrm{cell}^T
-                                //   D_\mathrm{cell} J_\mathrm{cell}
-                                //   x_1$. This is a block-diagonal
-                                //   operation, with the block size equal to
-                                //   @p dim. The blocks just
-                                //   correspond to the individual quadrature
-                                //   points. The operation on each quadrature
-                                //   point is implemented by the
-                                //   Transformation class object that this
-                                //   class is equipped with. Compared to the
-                                //   introduction, the matrix
-                                //   <i>D</i><sub>cell</sub> now contains the
-                                //   @p JxW values and the
-                                //   inhomogeneous coefficient.
-                                //   <li> Form $y_\mathrm{cell} =
-                                //   B_\mathrm{ref\_cell}^T x_2$. This gives
-                                //   the local result of the matrix-vector
-                                //   product.
-                                //   <li> Form $y \leftarrow y +
-                                //   P_\mathrm{cell,local-global}^T
-                                //   y_\mathrm{cell}$. This adds the local
-                                //   result to the global vector, which is
-                                //   realized using the method
-                                //   ConstraintMatrix::distribute_local_to_global.
-                                //   Note that we do this in an extra
-                                //   function called
-                                //   @p copy_local_to_global
-                                //   because that operation must not be done
-                                //   in %parallel, in order to avoid two or
-                                //   more processes trying to add to the same
-                                //   positions in the result vector <i>y</i>.
-                                //   </ol>
-                                // The steps 1 to 4 can be done in %parallel
-                                // by multiple processes.
-
-                                // Now, it turns out that the most expensive
-                                // part of the above is the multiplication
-                                // <i>B</i><sub>ref_cell</sub>
-                                // <i>x</i><sub>cell</sub> in the second step
-                                // and the transpose operation in step
-                                // 4. Note that the matrix
-                                // <i>J</i><sup>T</sup><i> D J</i> is
-                                // block-diagonal, and hence, its application
-                                // is cheaper. Since the matrix
-                                // <i>B</i><sub>ref_cell</sub> is the same
-                                // for all cells, all that changes is the
-                                // vector <i>x</i><sub>cell</sub>. Hence,
-                                // nothing prevents us from collecting
-                                // several cell vectors to a (rectangular)
-                                // matrix, and then perform a matrix-matrix
-                                // product. These matrices are both full, but
-                                // not very large, having of the order @p
-                                // dofs_per_cell rows and columns. This is an
-                                // operation that can be much better
-                                // optimized than matrix-vector products. The
-                                // functions @p FullMatrix<number>::mmult and
-                                // @p FullMatrix<number>::mTmult use the BLAS
-                                // dgemm function (as long as BLAS has been
-                                // detected in deal.II configuration), which
-                                // provides optimized kernels for doing this
-                                // product. In our case, a matrix-matrix
-                                // product is between three and five times
-                                // faster than doing the matrix-vector
-                                // product on one cell after the other. The
-                                // variables that hold the solution on the
-                                // respective cell's support points and the
-                                // quadrature points are thus full matrices,
-                                // which we set to the correct size as a
-                                // first action in this function. The number
-                                // of rows in the two matrices @p
-                                // scratch.solutions and @p copy.solutions is
-                                // given by the number of cells they work on,
-                                // and the number of columns is the number of
-                                // degrees of freedom per cell for the first
-                                // and the number of quadrature points times
-                                // the number of components per point for the
-                                // latter.
-template <typename number, class Transformation>
-template <typename number2>
-void
-MatrixFree<number,Transformation>::
-local_vmult (CellChunkIterator                    cell_range,
-            WorkStreamData::ScratchData<number> &scratch,
-            WorkStreamData::CopyData<number>    &copy,
-            const Vector<number2>               &src) const
-{
-  const unsigned int chunk_size = cell_range->second - cell_range->first;
-
-  scratch.solutions.reinit (chunk_size, matrix_sizes.n, true);
-  copy.solutions.reinit    (chunk_size, matrix_sizes.m, true);
-  copy.first_cell         = cell_range->first;
-  copy.n_dofs             = chunk_size*matrix_sizes.m;
-
-  constraints.get_dof_values(src, &indices_local_to_global(copy.first_cell,0),
-                            &copy.solutions(0,0),
-                            &copy.solutions(0,0)+copy.n_dofs);
-
-  copy.solutions.mmult (scratch.solutions, B_ref_cell);
-
-  for (unsigned int i=0, k = copy.first_cell; i<chunk_size; ++i, ++k)
-    for (unsigned int j=0; j<matrix_sizes.n_points; ++j)
-      derivatives(k,j).transform(&scratch.solutions(i, j*matrix_sizes.n_comp));
-
-  scratch.solutions.mTmult (copy.solutions, B_ref_cell);
-}
-
-
-
-template <typename number, class Transformation>
-template <typename number2>
-void
-MatrixFree<number,Transformation>::
-copy_local_to_global (const WorkStreamData::CopyData<number> &copy,
-                     Vector<number2>                        &dst) const
-{
-  constraints.distribute_local_to_global (&copy.solutions(0,0),
-                                         &copy.solutions(0,0)+copy.n_dofs,
-                                         &indices_local_to_global(copy.first_cell,0),
-                                         dst);
-}
-
-
-
-                                // Now to the @p vmult function that is
-                                // called externally: In addition to what we
-                                // do in a @p vmult_add function, we set the
-                                // destination to zero first.
-template <typename number, class Transformation>
-template <typename number2>
-void
-MatrixFree<number,Transformation>::vmult (Vector<number2>       &dst,
-                                         const Vector<number2> &src) const
-{
-  dst = 0;
-  vmult_add (dst, src);
-}
-
-
-
-                                // Transposed matrix-vector products (needed
-                                // for the multigrid operations to be
-                                // well-defined): do the same. Since we
-                                // implement a symmetric operation, we can
-                                // refer to the @p vmult_add operation.
-template <typename number, class Transformation>
-template <typename number2>
-void
-MatrixFree<number,Transformation>::Tvmult (Vector<number2>       &dst,
-                                          const Vector<number2> &src) const
-{
-  dst = 0;
-  Tvmult_add (dst,src);
-}
-
-
-
-template <typename number, class Transformation>
-template <typename number2>
-void
-MatrixFree<number,Transformation>::Tvmult_add (Vector<number2>       &dst,
-                                              const Vector<number2> &src) const
-{
-  vmult_add (dst,src);
-}
-
-
-
-                                // This is the @p vmult_add function that
-                                // multiplies the matrix with vector @p src
-                                // and adds the result to vector @p dst.  We
-                                // include a few sanity checks to make sure
-                                // that the size of the vectors is the same
-                                // as the dimension of the matrix. We call a
-                                // %parallel function that applies the
-                                // multiplication on a chunk of cells at once
-                                // using the WorkStream module (cf. also the
-                                // @ref threads module). The subdivision into
-                                // chunks will be performed in the reinit
-                                // function and is stored in the field @p
-                                // matrix_sizes.chunks. What the rather
-                                // cryptic command to @p std_cxx1x::bind does
-                                // is to transform a function that has
-                                // several arguments (source vector, chunk
-                                // information) into a function which has
-                                // three arguments (in the first case) or one
-                                // argument (in the second), which is what
-                                // the WorkStream::run function expects. The
-                                // placeholders <code>_1, std_cxx1x::_2, _3</code> in
-                                // the local vmult specify variable input
-                                // values, given by the chunk information,
-                                // scratch data and copy data that the
-                                // WorkStream::run function will provide,
-                                // whereas the other arguments to the @p
-                                // local_vmult function are bound: to @p this
-                                // and a constant reference to the @p src in
-                                // the first case, and @p this and a
-                                // reference to the output vector in the
-                                // second. Similarly, the placeholder
-                                // @p _1 argument in the
-                                // @p copy_local_to_global function
-                                // sets the first explicit argument of that
-                                // function, which is of class
-                                // @p CopyData. We need to
-                                // abstractly specify these arguments because
-                                // the tasks defined by different cell chunks
-                                // will be scheduled by the WorkStream class,
-                                // and we will reuse available scratch and
-                                // copy data.
-template <typename number, class Transformation>
-template <typename number2>
-void
-MatrixFree<number,Transformation>::vmult_add (Vector<number2>       &dst,
-                                             const Vector<number2> &src) const
-{
-  Assert (src.size() == n(), ExcDimensionMismatch(src.size(), n()));
-  Assert (dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
-
-  WorkStream::run (matrix_sizes.chunks.begin(), matrix_sizes.chunks.end(),
-                  std_cxx1x::bind(&MatrixFree<number,Transformation>::
-                                  template local_vmult<number2>,
-                                  this, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3, boost::cref(src)),
-                  std_cxx1x::bind(&MatrixFree<number,Transformation>::
-                                  template copy_local_to_global<number2>,
-                                  this, std_cxx1x::_1, boost::ref(dst)),
-                  WorkStreamData::ScratchData<number>(),
-                  WorkStreamData::CopyData<number>(),
-                  2*multithread_info.n_default_threads,1);
-
-                                  // One thing to be cautious about:
-                                  // The deal.II classes expect that
-                                  // the matrix still contains a
-                                  // diagonal entry for constrained
-                                  // dofs (otherwise, the matrix
-                                  // would be singular, which is not
-                                  // what we want). Since the
-                                  // <code>distribute_local_to_global</code>
-                                  // command of the constraint matrix
-                                  // which we used for adding the
-                                  // local elements into the global
-                                  // vector does not do anything with
-                                  // constrained elements, we have to
-                                  // circumvent that problem by
-                                  // artificially setting the
-                                  // diagonal to some non-zero value
-                                  // and adding the source values. We
-                                  // simply set it to one, which
-                                  // corresponds to copying the
-                                  // respective elements of the
-                                  // source vector into the matching
-                                  // entry of the destination vector.
-  for (unsigned int i=0; i<matrix_sizes.n_dofs; ++i)
-    if (constraints.is_constrained(i) == true)
-      dst(i) += 1.0 * src(i);
-}
-
-
-
-                                // The next function initializes the
-                                // structures of the matrix. It writes the
-                                // number of total degrees of freedom in the
-                                // problem as well as the number of cells to
-                                // the MatrixSizes struct and copies the
-                                // small matrix that transforms the solution
-                                // from support points to quadrature
-                                // points. It uses the small matrix for
-                                // determining the number of degrees of
-                                // freedom per cell (number of rows in @p
-                                // B_ref_cell). The number of quadrature
-                                // points needs to be passed through the last
-                                // variable @p n_points_per_cell, since the
-                                // number of columns in the small matrix is
-                                // @p dim*n_points_per_cell for the Laplace
-                                // problem (the Laplacian is a tensor and has
-                                // @p dim components). In this function, we
-                                // also give the fields containing the
-                                // derivative information and the local dof
-                                // indices the correct sizes. They will be
-                                // filled by calling the respective set
-                                // function defined above.
-template <typename number, class Transformation>
-void MatrixFree<number,Transformation>::
-reinit (const unsigned int        n_dofs_in,
-       const unsigned int        n_cells_in,
-       const FullMatrix<double> &B_ref_cell_in,
-       const unsigned int        n_points_per_cell)
-{
-  B_ref_cell = B_ref_cell_in;
-
-  derivatives.reinit (n_cells_in, n_points_per_cell);
-  indices_local_to_global.reinit (n_cells_in, B_ref_cell.m());
-
-  diagonal_is_calculated = false;
-
-  matrix_sizes.n_dofs = n_dofs_in;
-  matrix_sizes.n_cells = n_cells_in;
-  matrix_sizes.m = B_ref_cell.m();
-  matrix_sizes.n = B_ref_cell.n();
-  matrix_sizes.n_points = n_points_per_cell;
-  matrix_sizes.n_comp   = B_ref_cell.n()/matrix_sizes.n_points;
-  Assert(matrix_sizes.n_comp * n_points_per_cell == B_ref_cell.n(),
-        ExcInternalError());
-
-                                  // One thing to make the matrix-vector
-                                  // product with this class efficient is to
-                                  // decide how many cells should be combined
-                                  // to one chunk, which will determine the
-                                  // size of the full matrix that we work
-                                  // on. If we choose too few cells, then the
-                                  // gains from using the matrix-matrix
-                                  // product will not be fully utilized
-                                  // (dgemm tends to provide more efficiency
-                                  // the larger the matrix dimensions get),
-                                  // so we choose at least 60 cells for one
-                                  // chunk (except when there are very few
-                                  // cells, like on the coarse levels of the
-                                  // multigrid scheme). If we choose too
-                                  // many, we will degrade parallelization
-                                  // (we need to have sufficiently
-                                  // independent tasks). We need to also
-                                  // think about the fact that most high
-                                  // performance BLAS implementations
-                                  // internally work with square
-                                  // sub-matrices. Choosing as many cells in
-                                  // a chunk as there are degrees of freedom
-                                  // on each cell (coded in @p
-                                  // matrix_sizes.m) respects the BLAS GEMM
-                                  // design, whenever we exceed 60. Clearly,
-                                  // the chunk size is an
-                                  // architecture-dependent value and the
-                                  // interested user can squeeze out some
-                                  // extra performance by hand-tuning this
-                                  // parameter. Once we have chosen the
-                                  // number of cells we collect in one chunk,
-                                  // we determine how many chunks we have on
-                                  // the given cell range and recalculate the
-                                  // actual chunk size in order to evenly
-                                  // distribute the chunks.
-  const unsigned int divisor = std::max(60U, matrix_sizes.m);
-  const unsigned int n_chunks = std::max (matrix_sizes.n_cells/divisor + 1,
-                                         2*multithread_info.n_default_threads);
-
-  const unsigned int chunk_size = (matrix_sizes.n_cells/n_chunks +
-                                  (matrix_sizes.n_cells%n_chunks>0));
-
-  std::pair<unsigned int, unsigned int> chunk;
-  for (unsigned int i=0; i<n_chunks; ++i)
-    {
-      chunk.first = i*chunk_size;
-      if ((i+1)*chunk_size > matrix_sizes.n_cells)
-       chunk.second = matrix_sizes.n_cells;
-      else
-       chunk.second = (i+1)*chunk_size;
-
-      if (chunk.second > chunk.first)
-       matrix_sizes.chunks.push_back(chunk);
-      else
-       break;
-    }
-}
-
-
-
-                                // Then we need a function if we want to
-                                // delete the content of the matrix,
-                                // e.g. when we are finished with one grid
-                                // level and continue to the next one. Just
-                                // set all the field sizes to 0.
-template <typename number, class Transformation>
-void
-MatrixFree<number,Transformation>::clear ()
-{
-  B_ref_cell.reinit(0,0);
-  derivatives.reinit (0,0);
-  indices_local_to_global.reinit(0,0);
-
-  constraints.clear();
-
-  diagonal_values.reinit (0);
-  diagonal_is_calculated = false;
-
-  matrix_sizes.n_dofs = 0;
-  matrix_sizes.n_cells = 0;
-  matrix_sizes.chunks.clear();
-}
-
-
-
-                                // The next function returns the entries of the
-                                // matrix. Since this class is intended not
-                                // to store the matrix entries, it would make
-                                // no sense to provide all those
-                                // elements. However, diagonal entries are
-                                // explicitly needed for the implementation
-                                // of the Chebyshev smoother that we intend
-                                // to use in the multigrid
-                                // preconditioner. This matrix is equipped
-                                // with a vector that stores the diagonal,
-                                // and we compute it when this function is
-                                // called for the first time.
-template <typename number, class Transformation>
-number
-MatrixFree<number,Transformation>::el (const unsigned int row,
-                                      const unsigned int col) const
-{
-  Assert (row == col, ExcNotImplemented());
-  if (diagonal_is_calculated == false)
-    calculate_diagonal();
-
-  return diagonal_values(row);
-}
-
-
-
-                                // Regarding the calculation of the diagonal,
-                                // remember that this is as simple (or
-                                // complicated) as assembling a right hand
-                                // side in deal.II. Well, it is a bit easier
-                                // to do this within this class since we have
-                                // all the derivative information
-                                // available. What we do is to go through all
-                                // the cells (now in serial, since this
-                                // function should not be called very often
-                                // anyway), then all the degrees of
-                                // freedom. At this place, we first copy the
-                                // first basis functions in all the
-                                // quadrature points to a temporary array,
-                                // apply the derivatives from the Jacobian
-                                // matrix, and finally multiply with the
-                                // second basis function. This is exactly the
-                                // value that would be written into the
-                                // diagonal of a sparse matrix. Note that we
-                                // need to condense hanging node constraints
-                                // and set the constrained diagonals to one.
-template <typename number, class Transformation>
-void
-MatrixFree<number,Transformation>::calculate_diagonal() const
-{
-  diagonal_values.reinit (matrix_sizes.n_dofs);
-  std::vector<number> calculation (matrix_sizes.n);
-  for (unsigned int cell=0; cell<matrix_sizes.n_cells; ++cell)
-    for (unsigned int dof=0; dof<matrix_sizes.m; ++dof)
-      {
-       memcpy (&calculation[0],&B_ref_cell(dof,0),
-               matrix_sizes.n*sizeof(number));
-       for (unsigned int q=0; q<matrix_sizes.n_points; ++q)
-         derivatives(cell,q).transform(&calculation[q*matrix_sizes.n_comp]);
-       double diag_value = 0;
-       for (unsigned int q=0; q<matrix_sizes.n; ++q)
-         diag_value += calculation[q] * B_ref_cell(dof,q);
-       diagonal_values(indices_local_to_global(cell,dof)) += diag_value;
-      }
-  constraints.condense (diagonal_values);
-  for (unsigned int i=0; i<matrix_sizes.n_dofs; ++i)
-    if (constraints.is_constrained(i) == true)
-      diagonal_values(i) = 1.0;
-
-  diagonal_is_calculated = true;
-}
-
-
-
-                                // Eventually, we provide a function that
-                                // calculates how much memory this class
-                                // uses. We just need to sum up the memory
-                                // consumption of the arrays, the
-                                // constraints, the small matrix and of the
-                                // local variables. Just as a remark: In 2D
-                                // and with data type @p double,
-                                // about 80 per cent of the memory
-                                // consumption is due to the
-                                // @p derivatives array, while in 3D
-                                // this number is even 85 per cent.
-template <typename number, class Transformation>
-std::size_t MatrixFree<number,Transformation>::memory_consumption () const
-{
-  std::size_t glob_size = derivatives.memory_consumption() +
-                         indices_local_to_global.memory_consumption() +
-                         constraints.memory_consumption() +
-                         B_ref_cell.memory_consumption() +
-                         diagonal_values.memory_consumption() +
-                         matrix_sizes.chunks.size()*2*sizeof(unsigned int) +
-                         sizeof(*this);
-  return glob_size;
-}
-
-
-
-                                // @sect3{Laplace operator implementation}
-
-                                // This class implements the local action of
-                                // a Laplace operator on a quadrature
-                                // point. This is a very basic class
-                                // implementation, providing functions for
-                                // initialization with a Tensor of rank 2 and
-                                // implementing the @p transform operation
-                                // needed by the @p MatrixFree class. There
-                                // is one point worth noting: The
-                                // quadrature-point related action of the
-                                // Laplace operator is a tensor of rank
-                                // two. It is symmetric since it is the
-                                // product of the inverse Jacobian
-                                // transformation between unit and real cell
-                                // with its transpose (times quadrature
-                                // weights and a coefficient, which are
-                                // scalar), so we can just save the diagonal
-                                // and upper diagonal part. We could use the
-                                // SymmetricTensor<2,dim> class for doing
-                                // this, however, that class is only based on
-                                // @p double %numbers. Since we also want to
-                                // use @p float %numbers for the multigrid
-                                // preconditioner (in order to save memory
-                                // and computing time), we manually implement
-                                // this operator. Note that @p dim is a
-                                // template argument and hence known at
-                                // compile-time, so the compiler knows that
-                                // this symmetric rank-2 tensor has 3 entries
-                                // if used in 2D and 6 entries if used in 3D.
-template <int dim,typename number>
-class LaplaceOperator
-{
-  public:
-    LaplaceOperator ();
-
-    LaplaceOperator (const Tensor<2,dim> &tensor);
-
-    void transform (number * result) const;
-
-    LaplaceOperator<dim,number>&
-    operator = (const Tensor<2,dim> &tensor);
-
-    unsigned int memory_consumption () const;
-
-  private:
-    number transformation[dim*(dim+1)/2];
-};
-
-template<int dim,typename number>
-LaplaceOperator<dim,number>::LaplaceOperator()
-{}
-
-
-template<int dim,typename number>
-LaplaceOperator<dim,number>::LaplaceOperator(const Tensor<2,dim> &tensor)
-{
-  *this = tensor;
-}
-
-                                // Now implement the transformation, which is
-                                // just a so-called contraction
-                                // operation between a tensor of rank two and a
-                                // tensor of rank one. Unfortunately, we
-                                // need to implement this by hand, since we
-                                // chose not to use the
-                                // SymmetricTensor<2,dim> class (note that
-                                // the resulting values are entries in a full
-                                // matrix that consists of doubles or
-                                // floats). It feels a bit unsafe to operate
-                                // on a pointer to the data, but that is the
-                                // only possibility if we do not want to copy
-                                // data back and forth, which is expensive
-                                // since this is the innermost position of
-                                // the loop in the @p vmult
-                                // operation of the MatrixFree class. We need
-                                // to pay attention to the fact that we only
-                                // saved half of the (symmetric) rank-two
-                                // tensor.
-                                //
-                                // At first sight, it seems inefficient that
-                                // we have an @p if clause at this position
-                                // in the code at the innermost loop, but
-                                // note once again that @p dim is known when
-                                // this piece of code is compiled, so the
-                                // compiler can optimize away the @p if
-                                // statement (and actually even inline these
-                                // few lines of code into the @p MatrixFree
-                                // class).
-template <int dim, typename number>
-void LaplaceOperator<dim,number>::transform (number* result) const
-{
-  if (dim == 2)
-    {
-      const number temp = result[0];
-      result[0] = transformation[0] * temp + transformation[1] * result[1];
-      result[1] = transformation[1] * temp + transformation[2] * result[1];
-    }
-  else if (dim == 3)
-    {
-      const number temp1 = result[0];
-      const number temp2 = result[1];
-      result[0] = transformation[0] * temp1 + transformation[1] * temp2 +
-                 transformation[2] * result[2];
-      result[1] = transformation[1] * temp1 + transformation[3] * temp2 +
-                 transformation[4] * result[2];
-      result[2] = transformation[2] * temp1 + transformation[4] * temp2 +
-                 transformation[5] * result[2];
-    }
-  else
-    ExcNotImplemented();
-}
-
-                                // The final function in this group
-                                // takes the content of a rank-2
-                                // tensor and writes it to the field
-                                // @p transformation of
-                                // this class. We save the upper part
-                                // of the symmetric tensor row-wise:
-                                // we first take the (0,0)-entry,
-                                // then the (0,1)-entry, and so
-                                // on. We only implement this for
-                                // dimensions two and three, which
-                                // for the moment should do just
-                                // fine:
-template <int dim, typename number>
-LaplaceOperator<dim,number>&
-LaplaceOperator<dim,number>::operator=(const Tensor<2,dim> &tensor)
-{
-  if (dim == 2)
-    {
-      transformation[0] = tensor[0][0];
-      transformation[1] = tensor[0][1];
-      transformation[2] = tensor[1][1];
-      Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15,
-             ExcInternalError());
-    }
-  else if (dim == 3)
-    {
-      transformation[0] = tensor[0][0];
-      transformation[1] = tensor[0][1];
-      transformation[2] = tensor[0][2];
-      transformation[3] = tensor[1][1];
-      transformation[4] = tensor[1][2];
-      transformation[5] = tensor[2][2];
-      Assert (std::fabs(tensor[1][0]-tensor[0][1])<1e-15,
-             ExcInternalError());
-      Assert (std::fabs(tensor[2][0]-tensor[0][2])<1e-15,
-             ExcInternalError());
-      Assert (std::fabs(tensor[2][1]-tensor[1][2])<1e-15,
-             ExcInternalError());
-    }
-  else
-    ExcNotImplemented();
-  return *this;
-}
-
-
-template<int dim,typename number>
-unsigned int
-LaplaceOperator<dim,number>::memory_consumption () const
-{
-  return sizeof(*this);
-}
-
-
-
-                                // @sect3{LaplaceProblem class}
-
-                                // This class is based on the same
-                                // class in step-16. However, we
-                                // replaced the SparseMatrix<double>
-                                // class by our matrix-free
-                                // implementation, which means that
-                                // we can also skip the sparsity
-                                // patterns.
-template <int dim>
-class LaplaceProblem
-{
-  public:
-    LaplaceProblem (const unsigned int degree);
-    void run ();
-
-  private:
-    void setup_system ();
-    void assemble_system ();
-    void assemble_multigrid ();
-    void solve ();
-    void output_results (const unsigned int cycle) const;
-
-    Triangulation<dim>   triangulation;
-    FE_Q<dim>            fe;
-    MGDoFHandler<dim>    mg_dof_handler;
-
-    MatrixFree<double,LaplaceOperator<dim,double> > system_matrix;
-    typedef MatrixFree<float,LaplaceOperator<dim,float> > MatrixFreeType;
-    MGLevelObject<MatrixFreeType> mg_matrices;
-    FullMatrix<float>             coarse_matrix;
-
-    Vector<double>       solution;
-    Vector<double>       system_rhs;
-};
-
-
-
-template <int dim>
-LaplaceProblem<dim>::LaplaceProblem (const unsigned int degree)
-               :
-                fe (degree),
-               mg_dof_handler (triangulation)
-{}
-
-
-
-                                // @sect4{LaplaceProblem::setup_system}
-
-                                // This is the function of step-16 with
-                                // relevant changes due to the MatrixFree
-                                // class. What we need to do is to somehow
-                                // create a local gradient matrix that does
-                                // not contain any cell-related data
-                                // (gradient on the reference cell). The
-                                // way to get to this matrix is to create
-                                // an FEValues object with gradient
-                                // information on a cell that corresponds
-                                // to the reference cell, which is a cube
-                                // with side length 1. So we create a
-                                // pseudo triangulation, initialize the
-                                // FEValues to the only cell of that
-                                // triangulation, and read off the
-                                // gradients (which we put in a
-                                // FullMatrix). That full matrix is then
-                                // passed to the reinit function of the
-                                // MatrixFree class used as a system matrix
-                                // and, further down, as multigrid matrices
-                                // on the individual levels. We need to
-                                // implement Dirichlet boundary conditions
-                                // here, which is done with the
-                                // ConstraintMatrix function as shown,
-                                // e.g., in step-22.
-template <int dim>
-void LaplaceProblem<dim>::setup_system ()
-{
-  system_matrix.clear();
-  mg_matrices.clear();
-
-  mg_dof_handler.distribute_dofs (fe);
-
-  std::cout << "Number of degrees of freedom: "
-           << mg_dof_handler.n_dofs()
-           << std::endl;
-
-  const unsigned int nlevels = triangulation.n_levels();
-  mg_matrices.resize(0, nlevels-1);
-
-  QGauss<dim>  quadrature_formula(fe.degree+1);
-  FEValues<dim> fe_values_reference (fe, quadrature_formula,
-                                    update_gradients);
-  Triangulation<dim> reference_cell;
-  GridGenerator::hyper_cube (reference_cell, 0, 1);
-  fe_values_reference.reinit (reference_cell.begin());
-  FullMatrix<double> ref_cell_gradients (fe.dofs_per_cell,
-                                        quadrature_formula.size()*dim);
-  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-    {
-      for (unsigned int j=0; j<quadrature_formula.size(); ++j)
-       {
-         for (unsigned int d=0; d<dim; ++d)
-           ref_cell_gradients(i,j*dim+d) = fe_values_reference.shape_grad(i,j)[d];
-       }
-    }
-  system_matrix.reinit (mg_dof_handler.n_dofs(), triangulation.n_active_cells(),
-                       ref_cell_gradients, quadrature_formula.size());
-  VectorTools::interpolate_boundary_values (mg_dof_handler,
-                                           0,
-                                           ZeroFunction<dim>(),
-                                           system_matrix.get_constraints());
-  system_matrix.get_constraints().close();
-  std::cout.precision(4);
-  std::cout << "System matrix memory consumption: "
-           << system_matrix.memory_consumption()/double(1<<20)
-           << " MiB."
-           << std::endl;
-
-  solution.reinit (mg_dof_handler.n_dofs());
-  system_rhs.reinit (mg_dof_handler.n_dofs());
-
-                                  // Next, initialize the matrices for the
-                                  // multigrid method on all the
-                                  // levels. Unfortunately, the function
-                                  // MGTools::make_boundary_list cannot write
-                                  // Dirichlet boundary conditions into a
-                                  // ConstraintMatrix object directly, so we
-                                  // first have to make the boundary list and
-                                  // then manually fill the boundary
-                                  // conditions using the command
-                                  // ConstraintMatrix::add_line. Once this is
-                                  // done, we close the ConstraintMatrix so
-                                  // it can be used for matrix-vector
-                                  // products.
-  typename FunctionMap<dim>::type dirichlet_boundary;
-  ZeroFunction<dim>               homogeneous_dirichlet_bc (1);
-  dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
-  std::vector<std::set<unsigned int> > boundary_indices(triangulation.n_levels());
-  MGTools::make_boundary_list (mg_dof_handler,
-                              dirichlet_boundary,
-                              boundary_indices);
-  for (unsigned int level=0;level<nlevels;++level)
-    {
-      mg_matrices[level].reinit(mg_dof_handler.n_dofs(level),
-                               triangulation.n_cells(level),
-                               ref_cell_gradients,
-                               quadrature_formula.size());
-      std::set<unsigned int>::iterator bc_it = boundary_indices[level].begin();
-      for ( ; bc_it != boundary_indices[level].end(); ++bc_it)
-       mg_matrices[level].get_constraints().add_line(*bc_it);
-      mg_matrices[level].get_constraints().close();
-    }
-  coarse_matrix.reinit (mg_dof_handler.n_dofs(0),
-                       mg_dof_handler.n_dofs(0));
-}
-
-
-
-                                // @sect4{LaplaceProblem::assemble_system}
-
-                                // The assemble function is significantly
-                                // reduced compared to step-16. All we need
-                                // to do is to assemble the right hand side
-                                // and to calculate the cell-dependent part
-                                // of the Laplace operator. The first task is
-                                // standard. The second is also not too hard
-                                // given the discussion in the introduction:
-                                // We need to take the inverse of the
-                                // Jacobian of the transformation from unit
-                                // to real cell, multiply it with its
-                                // transpose and multiply the resulting
-                                // rank-2 tensor with the quadrature weights
-                                // and the coefficient values at the
-                                // quadrature points. To make this work, we
-                                // add the update flag @p
-                                // update_inverse_jacobians to the FEValues
-                                // constructor, and query the inverse of the
-                                // Jacobian in a loop over the quadrature
-                                // points (note that the Jacobian is not
-                                // related to any kind of degrees of freedom
-                                // directly). In the end, we condense the
-                                // constraints from Dirichlet boundary
-                                // conditions away from the right hand side.
-template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
-{
-  QGauss<dim>  quadrature_formula(fe.degree+1);
-  MappingQ<dim> mapping (fe.degree);
-  FEValues<dim> fe_values (mapping, fe, quadrature_formula,
-                          update_values   | update_inverse_jacobians |
-                           update_quadrature_points | update_JxW_values);
-
-  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-  const unsigned int   n_q_points    = quadrature_formula.size();
-
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-  const Coefficient<dim> coefficient;
-  std::vector<double>    coefficient_values (n_q_points);
-
-  unsigned int cell_no = 0;
-
-  typename DoFHandler<dim>::active_cell_iterator cell = mg_dof_handler.begin_active(),
-                                                endc = mg_dof_handler.end();
-  for (; cell!=endc; ++cell, ++cell_no)
-    {
-      cell->get_dof_indices (local_dof_indices);
-      fe_values.reinit (cell);
-      coefficient.value_list (fe_values.get_quadrature_points(),
-                             coefficient_values);
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       {
-         double rhs_val = 0;
-         for (unsigned int q=0; q<n_q_points; ++q)
-           rhs_val += (fe_values.shape_value(i,q) * 1.0 *
-                       fe_values.JxW(q));
-         system_rhs(local_dof_indices[i]) += rhs_val;
-       }
-
-      system_matrix.set_local_dof_indices (cell_no, local_dof_indices);
-      for (unsigned int q=0; q<n_q_points; ++q)
-       system_matrix.set_derivative_data (cell_no, q,
-                                          (transpose
-                                           (fe_values.inverse_jacobian(q)) *
-                                           fe_values.inverse_jacobian(q)) *
-                                          fe_values.JxW(q) *
-                                          coefficient_values[q]);
-    }
-  system_matrix.get_constraints().condense(system_rhs);
-}
-
-
-                                // @sect4{LaplaceProblem::assemble_multigrid}
-
-                                // Here is another assemble
-                                // function. The integration core is
-                                // the same as above. Only the loop
-                                // goes over all existing cells now
-                                // and the results must be entered
-                                // into the correct matrix.
-
-                                // Since we only do multilevel
-                                // preconditioning, no right-hand side is
-                                // assembled here. Compared to step-16, there
-                                // is one new thing here: we manually
-                                // calculate the matrix on the coarsest
-                                // level. In step-16, we could simply copy
-                                // the entries from the respective sparse
-                                // matrix, but this is obviously not possible
-                                // here. We could have integrated this into the
-                                // MatrixFree class as well, but it is simple
-                                // enough, so calculate it here instead.
-template <int dim>
-void LaplaceProblem<dim>::assemble_multigrid ()
-{
-  coarse_matrix = 0;
-  QGauss<dim>  quadrature_formula(fe.degree+1);
-  MappingQ<dim> mapping (fe.degree);
-  FEValues<dim> fe_values (mapping, fe, quadrature_formula,
-                          update_gradients  | update_inverse_jacobians |
-                           update_quadrature_points | update_JxW_values);
-
-  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-  const unsigned int   n_q_points    = quadrature_formula.size();
-
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-  const Coefficient<dim> coefficient;
-  std::vector<double>    coefficient_values (n_q_points);
-
-  std::vector<unsigned int> cell_no(triangulation.n_levels());
-  typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
-                                           endc = mg_dof_handler.end();
-  for (; cell!=endc; ++cell)
-    {
-      const unsigned int level = cell->level();
-      cell->get_mg_dof_indices (local_dof_indices);
-      fe_values.reinit (cell);
-      coefficient.value_list (fe_values.get_quadrature_points(),
-                             coefficient_values);
-
-      mg_matrices[level].set_local_dof_indices (cell_no[level],
-                                               local_dof_indices);
-      for (unsigned int q=0; q<n_q_points; ++q)
-       mg_matrices[level].set_derivative_data (cell_no[level], q,
-                                               (transpose
-                                                (fe_values.inverse_jacobian(q)) *
-                                                fe_values.inverse_jacobian(q)) *
-                                               fe_values.JxW(q) *
-                                               coefficient_values[q]);
-
-      ++cell_no[level];
-      if (level == 0)
-       {
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           for (unsigned int j=0; j<dofs_per_cell; ++j)
-             {
-               double add_value = 0;
-               for (unsigned int q=0; q<n_q_points; ++q)
-                 add_value += (fe_values.shape_grad(i,q) *
-                               fe_values.shape_grad(j,q) *
-                               coefficient_values[q] *
-                               fe_values.JxW(q));
-               coarse_matrix(local_dof_indices[i],
-                             local_dof_indices[j]) += add_value;
-             }
-       }
-    }
-
-                                  // In a final step, we need to
-                                  // condense the boundary conditions
-                                  // on the coarse matrix. There is
-                                  // no built-in function for doing
-                                  // this on a full matrix, so
-                                  // manually delete the rows and
-                                  // columns of the matrix that are
-                                  // constrained.
-  for (unsigned int i=0; i<coarse_matrix.m(); ++i)
-    if (mg_matrices[0].get_constraints().is_constrained(i))
-      for (unsigned int j=0; j<coarse_matrix.n(); ++j)
-       if (i!=j)
-         {
-           coarse_matrix(i,j) = 0;
-           coarse_matrix(j,i) = 0;
-         }
-}
-
-
-
-                                // @sect4{LaplaceProblem::solve}
-
-                                // The solution process again looks like
-                                // step-16. We now use a Chebyshev smoother
-                                // instead of SOR (SOR would be very
-                                // difficult to implement because we do not
-                                // have the matrix elements available
-                                // explicitly, and it is difficult to make it
-                                // work efficiently in %parallel). The
-                                // multigrid classes provide a simple
-                                // interface for using the Chebyshev smoother
-                                // which is defined in a preconditioner
-                                // class: MGSmootherPrecondition.
-template <int dim>
-void LaplaceProblem<dim>::solve ()
-{
-  GrowingVectorMemory<>   vector_memory;
-
-  MGTransferPrebuilt<Vector<double> > mg_transfer;
-  mg_transfer.build_matrices(mg_dof_handler);
-
-  MGCoarseGridHouseholder<float, Vector<double> > mg_coarse;
-  mg_coarse.initialize(coarse_matrix);
-
-  typedef PreconditionChebyshev<MatrixFreeType,Vector<double> > SMOOTHER;
-  MGSmootherPrecondition<MatrixFreeType, SMOOTHER, Vector<double> >
-    mg_smoother(vector_memory);
-
-                                  // Then, we initialize the smoother
-                                  // with our level matrices and the
-                                  // required, additional data for
-                                  // the Chebyshev smoother. In
-                                  // particular, we use a higher
-                                  // polynomial degree for higher
-                                  // order elements, since smoothing
-                                  // gets more difficult for
-                                  // these. Smooth out a range of
-                                  // $[\lambda_{\max}/10,\lambda_{\max}]$. In
-                                  // order to compute the maximum
-                                  // eigenvalue of the corresponding
-                                  // matrix, the Chebyshev
-                                  // initializations performs a few
-                                  // steps of a CG algorithm. Since
-                                  // all we need is a rough estimate,
-                                  // we choose some eight iterations
-                                  // (more if the finite element
-                                  // polynomial degree is larger,
-                                  // less if it is smaller than
-                                  // quadratic).
-  typename SMOOTHER::AdditionalData smoother_data;
-  smoother_data.smoothing_range = 10.;
-  smoother_data.degree = fe.degree;
-  smoother_data.eig_cg_n_iterations = 4+2*fe.degree;
-  mg_smoother.initialize(mg_matrices, smoother_data);
-
-  MGMatrix<MatrixFreeType, Vector<double> >
-    mg_matrix(&mg_matrices);
-
-  Multigrid<Vector<double> > mg(mg_dof_handler,
-                               mg_matrix,
-                               mg_coarse,
-                               mg_transfer,
-                               mg_smoother,
-                               mg_smoother);
-  PreconditionMG<dim, Vector<double>,
-    MGTransferPrebuilt<Vector<double> > >
-  preconditioner(mg_dof_handler, mg, mg_transfer);
-
-                                // Finally, write out the memory
-                                // consumption of the Multigrid object
-                                // (or rather, of its most significant
-                                // components, since there is no built-in
-                                // function for the total multigrid
-                                // object), then create the solver object
-                                // and solve the system. This is very
-                                // easy, and we didn't even see any
-                                // difference in the solve process
-                                // compared to step-16. The magic is all
-                                // hidden behind the implementation of
-                                // the MatrixFree::vmult operation.
-  const unsigned int multigrid_memory
-    = (mg_matrices.memory_consumption() +
-       mg_transfer.memory_consumption() +
-       coarse_matrix.memory_consumption());
-  std::cout << "Multigrid objects memory consumption: "
-           << multigrid_memory/double(1<<20)
-           << " MiB."
-           << std::endl;
-
-  SolverControl           solver_control (1000, 1e-12);
-  SolverCG<>              cg (solver_control);
-
-  cg.solve (system_matrix, solution, system_rhs,
-           preconditioner);
-
-  std::cout << "Convergence in " << solver_control.last_step()
-           << " CG iterations." << std::endl;
-}
-
-
-
-                                // @sect4{LaplaceProblem::output_results}
-
-                                // Here is the data output, which is a
-                                // simplified version of step-5. We use the
-                                // standard VTK output for each grid
-                                // produced in the refinement process.
-template <int dim>
-void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
-{
-  DataOut<dim> data_out;
-
-  data_out.attach_dof_handler (mg_dof_handler);
-  data_out.add_data_vector (solution, "solution");
-  data_out.build_patches ();
-
-  std::ostringstream filename;
-  filename << "solution-"
-          << cycle
-          << ".vtk";
-
-  std::ofstream output (filename.str().c_str());
-  data_out.write_vtk (output);
-}
-
-
-
-                                // @sect4{LaplaceProblem::run}
-
-                                // The function that runs the program is
-                                // very similar to the one in step-16. We
-                                // make less refinement steps in 3D
-                                // compared to 2D, but that's it.
-template <int dim>
-void LaplaceProblem<dim>::run ()
-{
-  for (unsigned int cycle=0; cycle<8-dim; ++cycle)
-    {
-      std::cout << "Cycle " << cycle << std::endl;
-
-      if (cycle == 0)
-       {
-         GridGenerator::hyper_ball(triangulation);
-         static const HyperBallBoundary<dim> boundary;
-         triangulation.set_boundary (0, boundary);
-         triangulation.refine_global (3-dim);
-       }
-      triangulation.refine_global (1);
-      setup_system ();
-      assemble_system ();
-      assemble_multigrid ();
-      solve ();
-      output_results (cycle);
-      std::cout << std::endl;
-    };
-}
-
-
-
-                                // @sect3{The <code>main</code> function}
-
-                                // This is as in all other programs:
-int main ()
-{
-  try
-    {
-      deallog.depth_console (0);
-      LaplaceProblem<2> laplace_problem (2);
-      laplace_problem.run ();
-    }
-  catch (std::exception &exc)
-    {
-      std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      std::cerr << "Exception on processing: " << std::endl
-               << exc.what() << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      return 1;
-    }
-  catch (...)
-    {
-      std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      std::cerr << "Unknown exception!" << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      return 1;
-    }
-
-  return 0;
-}

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.