]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Remove the template argument from the MatrixTools and MatrixCreator classes. This...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 14 Aug 2001 14:19:46 +0000 (14:19 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 14 Aug 2001 14:19:46 +0000 (14:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@4882 0785d39b-7218-0410-832d-ea1e28bc413d

17 files changed:
deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/source/grid/grid_generator.cc
deal.II/deal.II/source/numerics/matrices.all_dimensions.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/matrices.cc
deal.II/deal.II/source/numerics/vectors.cc
deal.II/doc/news/2001/c-3-1.html
deal.II/examples/step-3/step-3.cc
deal.II/examples/step-4/step-4.cc
deal.II/examples/step-5/step-5.cc
deal.II/examples/step-6/step-6.cc
deal.II/examples/step-7/step-7.cc
deal.II/examples/step-8/step-8.cc
tests/deal.II/block_matrices.cc
tests/deal.II/filtered_matrix.cc
tests/deal.II/matrices.cc
tests/deal.II/wave-test-3.cc

index dd0a4d48bc1eb9e67e965587df2f69beba863228..34ae55012cd8e51fd67ed0abc1d4f73917a97284 100644 (file)
@@ -23,7 +23,6 @@ template <int dim> class FEValues;
 template <int dim> class FEFaceValues;
 template <int dim> class FESubfaceValues;
 template <int dim> class FESystem;
-template <int dim> class MatrixCreator;
 
 
 /**
index d723e3011de652f8af05e70fc19f764c4c3047d0..91435db479189e0550a332852a16834bf9e40faf 100644 (file)
@@ -12,6 +12,7 @@
 
 #include <base/exceptions.h>
 #include <base/thread_management.h>
+#include <dofs/function_map.h>
 #include <map>
 
 
@@ -29,7 +30,6 @@ template <int dim> class Mapping;
 template <int dim> class DoFHandler;
 template <int dim> class MGDoFHandler;
 template <int dim> class FEValues;
-template <int dim> class FunctionMap;
 
 
 
@@ -183,8 +183,6 @@ template <int dim> class FunctionMap;
  *
  * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001
  */
-//TODO: [WB] remove template arg and make functions static once FunctionMap is globalized
-template <int dim>
 class MatrixCreator
 {
   public:
@@ -200,6 +198,7 @@ class MatrixCreator
                                      * See the general doc of this class
                                      * for more information.
                                      */
+    template <int dim>
     static void create_mass_matrix (const Mapping<dim>       &mapping,
                                    const DoFHandler<dim>    &dof,
                                    const Quadrature<dim>    &q,
@@ -211,6 +210,7 @@ class MatrixCreator
                                      * function, see above, with
                                      * @p{mapping=MappingQ1<dim>()}.
                                      */
+    template <int dim>
     static void create_mass_matrix (const DoFHandler<dim>    &dof,
                                    const Quadrature<dim>    &q,
                                    SparseMatrix<double>     &matrix,
@@ -229,6 +229,7 @@ class MatrixCreator
                                      * See the general doc of this
                                      * class for more information.
                                      */
+    template <int dim>
     static void create_mass_matrix (const Mapping<dim>       &mapping,
                                    const DoFHandler<dim>    &dof,
                                    const Quadrature<dim>    &q,
@@ -242,6 +243,7 @@ class MatrixCreator
                                      * function, see above, with
                                      * @p{mapping=MappingQ1<dim>()}.
                                      */
+    template <int dim>
     static void create_mass_matrix (const DoFHandler<dim>    &dof,
                                    const Quadrature<dim>    &q,
                                    SparseMatrix<double>     &matrix,
@@ -269,22 +271,38 @@ class MatrixCreator
                                      * See the general doc of this
                                      * class for more information.
                                      */
+    template <int dim>
     static
     void create_boundary_mass_matrix (const Mapping<dim>       &mapping,
                                      const DoFHandler<dim>    &dof,
                                      const Quadrature<dim-1>  &q,
                                      SparseMatrix<double>     &matrix,
-                                     const typename FunctionMap<dim>::type        &boundary_functions,
+                                     const typename FunctionMap<dim>::type &boundary_functions,
                                      Vector<double>           &rhs_vector,
                                      std::vector<unsigned int>&dof_to_boundary_mapping,
                                      const Function<dim> * const a = 0);
 
+                                    /**
+                                     * Same function, but for 1d.
+                                     */
+    static
+    void create_boundary_mass_matrix (const Mapping<1>       &mapping,
+                                     const DoFHandler<1>    &dof,
+                                     const Quadrature<0>    &q,
+                                     SparseMatrix<double>   &matrix,
+                                     const FunctionMap<1>::type &boundary_functions,
+                                     Vector<double>         &rhs_vector,
+                                     std::vector<unsigned int>&dof_to_boundary_mapping,
+                                     const Function<1> * const a = 0);
+
+
                                     /**
                                      * Calls the
                                      * @p{create_boundary_mass_matrix}
                                      * function, see above, with
                                      * @p{mapping=MappingQ1<dim>()}.
                                      */
+    template <int dim>
     static
     void create_boundary_mass_matrix (const DoFHandler<dim>    &dof,
                                      const Quadrature<dim-1>  &q,
@@ -307,6 +325,7 @@ class MatrixCreator
                                      * See the general doc of this
                                      * class for more information.
                                      */
+    template <int dim>
     static void create_laplace_matrix (const Mapping<dim>       &mapping,
                                       const DoFHandler<dim>    &dof,
                                       const Quadrature<dim>    &q,
@@ -318,6 +337,7 @@ class MatrixCreator
                                      * function, see above, with
                                      * @p{mapping=MappingQ1<dim>()}.
                                      */
+    template <int dim>
     static void create_laplace_matrix (const DoFHandler<dim>    &dof,
                                       const Quadrature<dim>    &q,
                                       SparseMatrix<double>     &matrix,
@@ -330,6 +350,7 @@ class MatrixCreator
                                      * See the general doc of this
                                      * class for more information.
                                      */
+    template <int dim>
     static void create_level_laplace_matrix (unsigned int             level,
                                             const MGDoFHandler<dim>& dof,
                                             const Quadrature<dim>&   q,
@@ -350,6 +371,7 @@ class MatrixCreator
                                      * See the general doc of this
                                      * class for more information.
                                      */
+    template <int dim>
     static void create_laplace_matrix (const Mapping<dim>       &mapping,
                                       const DoFHandler<dim>    &dof,
                                       const Quadrature<dim>    &q,
@@ -363,6 +385,7 @@ class MatrixCreator
                                      * function, see above, with
                                      * @p{mapping=MappingQ1<dim>()}.
                                      */
+    template <int dim>
     static void create_laplace_matrix (const DoFHandler<dim>    &dof,
                                       const Quadrature<dim>    &q,
                                       SparseMatrix<double>     &matrix,
@@ -378,17 +401,50 @@ class MatrixCreator
   private:
                                     /**
                                      * Convenience abbreviation for
-                                     * DoF handler cell iterators.
-                                     */
-    typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
-
-                                    /**
-                                     * Pair of iterators denoting a
-                                     * half-open range.
+                                     * pairs of DoF handler cell
+                                     * iterators. This type works
+                                     * just like a
+                                     * @p{std::pair<iterator,iterator>}
+                                     * but is templatized on the
+                                     * space dimension.
                                      */
-    typedef std::pair<active_cell_iterator,active_cell_iterator> IteratorRange;
+    template <int dim>
+    struct IteratorRange 
+    {
+                                        /**
+                                         * Typedef for the iterator type.
+                                         */
+       typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+
+                                        /**
+                                         * Abbreviation for a pair of
+                                         * iterators.
+                                         */
+       typedef typename std::pair<active_cell_iterator,active_cell_iterator> iterator_pair;
+       
+                                        /**
+                                         * Constructor. Initialize
+                                         * the two values by the
+                                         * given values.
+                                         */
+       IteratorRange (const active_cell_iterator &first,
+                      const active_cell_iterator &second);
+
+                                        /**
+                                         * Constructor taking a pair
+                                         * of values for
+                                         * initialization.
+                                         */
+       IteratorRange (const iterator_pair &ip);
+       
+                                        /**
+                                         * Pair of iterators denoting
+                                         * a half-open range.
+                                         */
+       active_cell_iterator first, second;
+    };
+    
     
-
                                     /**
                                      * Version of the same function
                                      * (without suffix @p{_1}) with
@@ -399,13 +455,14 @@ class MatrixCreator
                                      * used to synchronise access to
                                      * the matrix.
                                      */
+    template <int dim>
     static
     void create_mass_matrix_1 (const Mapping<dim>       &mapping,
                               const DoFHandler<dim>    &dof,
                               const Quadrature<dim>    &q,
                               SparseMatrix<double>     &matrix,
                               const Function<dim> * const a,
-                              const IteratorRange      &range,
+                              const IteratorRange<dim>  range,
                               Threads::ThreadMutex     &mutex);
 
                                     /**
@@ -418,6 +475,7 @@ class MatrixCreator
                                      * used to synchronise access to
                                      * the matrix.
                                      */
+    template <int dim>
     static
     void create_mass_matrix_2 (const Mapping<dim>       &mapping,
                               const DoFHandler<dim>    &dof,
@@ -426,7 +484,7 @@ class MatrixCreator
                               const Function<dim>      &rhs,
                               Vector<double>           &rhs_vector,
                               const Function<dim> * const a,
-                              const IteratorRange      &range,
+                              const IteratorRange<dim>  range,
                               Threads::ThreadMutex     &mutex);
 
                                     /**
@@ -439,13 +497,14 @@ class MatrixCreator
                                      * used to synchronise access to
                                      * the matrix.
                                      */
+    template <int dim>
     static
     void create_laplace_matrix_1 (const Mapping<dim>       &mapping,
                                  const DoFHandler<dim>    &dof,
                                  const Quadrature<dim>    &q,
                                  SparseMatrix<double>     &matrix,
                                  const Function<dim> * const a,
-                                 const IteratorRange      &range,
+                                 const IteratorRange<dim>  range,
                                  Threads::ThreadMutex     &mutex);
 
                                     /**
@@ -458,6 +517,7 @@ class MatrixCreator
                                      * used to synchronise access to
                                      * the matrix.
                                      */
+    template <int dim>
     static
     void create_laplace_matrix_2 (const Mapping<dim>       &mapping,
                                  const DoFHandler<dim>    &dof,
@@ -466,7 +526,7 @@ class MatrixCreator
                                  const Function<dim>      &rhs,
                                  Vector<double>           &rhs_vector,
                                  const Function<dim> * const a,
-                                 const IteratorRange      &range,
+                                 const IteratorRange<dim>  range,
                                  Threads::ThreadMutex     &mutex);
 
                                     /**
@@ -479,6 +539,7 @@ class MatrixCreator
                                      * used to synchronise access to
                                      * the matrix.
                                      */
+    template <int dim>
     static
     void create_boundary_mass_matrix_1 (const Mapping<dim>       &mapping,
                                        const DoFHandler<dim>    &dof,
@@ -488,36 +549,11 @@ class MatrixCreator
                                        Vector<double>           &rhs_vector,
                                        std::vector<unsigned int>&dof_to_boundary_mapping,
                                        const Function<dim> * const a,
-                                       const IteratorRange      &range,
+                                       const IteratorRange<dim>  range,
                                        Threads::ThreadMutex     &mutex);
 };
 
 
-/* -------------- declaration of explicit specializations ------------- */
-
-
-template <> void MatrixCreator<1>::create_boundary_mass_matrix (
-  const Mapping<1>          &,
-  const DoFHandler<1>       &,
-  const Quadrature<0>       &,
-  SparseMatrix<double>      &,
-  const FunctionMap<1>::type&,
-  Vector<double>            &,
-  std::vector<unsigned int> &,
-  const Function<1>         * const);
-
-template <> void MatrixCreator<1>::create_boundary_mass_matrix_1 (
-  const Mapping<1>          &,
-  const DoFHandler<1>       &,
-  const Quadrature<0>       &,
-  SparseMatrix<double>      &,
-  const FunctionMap<1>::type&,
-  Vector<double>            &,
-  std::vector<unsigned int> &,
-  const Function<1> * const ,
-  const IteratorRange       &,
-  Threads::ThreadMutex      &);
-
 
 /**
  * Provide a collection of functions operating on matrices. These include
@@ -628,8 +664,7 @@ template <> void MatrixCreator<1>::create_boundary_mass_matrix_1 (
  * 
  * @author Wolfgang Bangerth, 1998, 2000
  */
-template <int dim>
-class MatrixTools : public MatrixCreator<dim>
+class MatrixTools : public MatrixCreator
 {
   public:
                                     /**
index bbe0106d2ac3f508e34a9728ad8a2f585872e62e..1d809861794c6762899f6cb93cd9bdccf07c6008 100644 (file)
@@ -757,7 +757,7 @@ void GridGenerator::laplace_transformation (Triangulation<dim> &tria,
   
   QGauss4<dim> quadrature;
   
-  MatrixCreator<dim>::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S);
+  MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S);
 
                                   // set up the boundary values for
                                   // the laplace problem
diff --git a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc
new file mode 100644 (file)
index 0000000..c712281
--- /dev/null
@@ -0,0 +1,526 @@
+//----------------------------  matrices.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001 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.
+//
+//----------------------------  matrices.cc  ---------------------------
+
+
+#include <numerics/matrices.h>
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
+
+#include <algorithm>
+
+
+
+template <typename number>
+void
+MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+                                   SparseMatrix<number>  &matrix,
+                                   Vector<number>   &solution,
+                                   Vector<number>   &right_hand_side,
+                                   const bool        preserve_symmetry)
+{
+  Assert (matrix.n() == matrix.m(),
+         ExcDimensionsDontMatch(matrix.n(), matrix.m()));
+  Assert (matrix.n() == right_hand_side.size(),
+         ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
+  Assert (matrix.n() == solution.size(),
+         ExcDimensionsDontMatch(matrix.n(), solution.size()));
+                                  // if no boundary values are to be applied
+                                  // simply return
+  if (boundary_values.size() == 0)
+    return;
+
+
+  std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
+                                               endd = boundary_values.end();
+  const unsigned int n_dofs             = matrix.m();
+  const SparsityPattern    &sparsity    = matrix.get_sparsity_pattern();
+  const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
+  const unsigned int *sparsity_colnums  = sparsity.get_column_numbers();
+
+                                  // if a diagonal entry is zero
+                                  // later, then we use another
+                                  // number instead. take it to be
+                                  // the first nonzero diagonal
+                                  // element of the matrix, or 1 if
+                                  // there is no such thing
+  number first_nonzero_diagonal_entry = 1;
+  for (unsigned int i=0; i<n_dofs; ++i)
+    if (matrix.diag_element(i) != 0)
+      {
+       first_nonzero_diagonal_entry = matrix.diag_element(i);
+       break;
+      };
+
+  
+  for (; dof != endd; ++dof)
+    {
+      Assert (dof->first < n_dofs, ExcInternalError());
+      
+      const unsigned int dof_number = dof->first;
+                                      // for each boundary dof:
+      
+                                      // set entries of this line
+                                      // to zero except for the diagonal
+                                      // entry. Note that the diagonal
+                                      // entry is always the first one
+                                      // for square matrices, i.e.
+                                      // we shall not set
+                                      // matrix.global_entry(
+                                      //     sparsity_rowstart[dof.first])
+      const unsigned int last = sparsity_rowstart[dof_number+1];
+      for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+       matrix.global_entry(j) = 0.;
+
+
+                                      // set right hand side to
+                                      // wanted value: if main diagonal
+                                      // entry nonzero, don't touch it
+                                      // and scale rhs accordingly. If
+                                      // zero, take the first main
+                                      // diagonal entry we can find, or
+                                      // one if no nonzero main diagonal
+                                      // element exists. Normally, however,
+                                      // the main diagonal entry should
+                                      // not be zero.
+                                      //
+                                      // store the new rhs entry to make
+                                      // the gauss step more efficient
+      number new_rhs;
+      if (matrix.diag_element(dof_number) != 0.0)
+       new_rhs = right_hand_side(dof_number)
+               = dof->second * matrix.diag_element(dof_number);
+      else
+       {
+                                          // use the SparseMatrix::
+                                          // to work around a bug in
+                                          // egcs
+         matrix.SparseMatrix<number>::set (dof_number, dof_number,
+                                           first_nonzero_diagonal_entry);
+         new_rhs = right_hand_side(dof_number)
+                 = dof->second * first_nonzero_diagonal_entry;
+       };
+
+
+                                      // if the user wants to have
+                                      // the symmetry of the matrix
+                                      // preserved, and if the
+                                      // sparsity pattern is
+                                      // symmetric, then do a Gauss
+                                      // elimination step with the
+                                      // present row
+      if (preserve_symmetry)
+       {
+                                          // store the only nonzero entry
+                                          // of this line for the Gauss
+                                          // elimination step
+         const number diagonal_entry = matrix.diag_element(dof_number);
+         
+                                          // we have to loop over all
+                                          // rows of the matrix which
+                                          // have a nonzero entry in
+                                          // the column which we work
+                                          // in presently. if the
+                                          // sparsity pattern is
+                                          // symmetric, then we can
+                                          // get the positions of
+                                          // these rows cheaply by
+                                          // looking at the nonzero
+                                          // column numbers of the
+                                          // present row. we need not
+                                          // look at the first entry,
+                                          // since that is the
+                                          // diagonal element and
+                                          // thus the present row
+         for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+           {
+             const unsigned int row = sparsity_colnums[j];
+
+                                              // find the position of
+                                              // element
+                                              // (row,dof_number)
+             const unsigned int *
+               p = std::lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
+                                    &sparsity_colnums[sparsity_rowstart[row+1]],
+                                    dof_number);
+
+                                              // check whether this line has
+                                              // an entry in the regarding column
+                                              // (check for ==dof_number and
+                                              // != next_row, since if
+                                              // row==dof_number-1, *p is a
+                                              // past-the-end pointer but points
+                                              // to dof_number anyway...)
+                                              //
+                                              // there should be such an entry!
+             Assert ((*p == dof_number) &&
+                     (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
+                     ExcInternalError());
+
+             const unsigned int global_entry
+               = (p - &sparsity_colnums[sparsity_rowstart[0]]);
+             
+                                              // correct right hand side
+             right_hand_side(row) -= matrix.global_entry(global_entry) /
+                                     diagonal_entry * new_rhs;
+             
+                                              // set matrix entry to zero
+             matrix.global_entry(global_entry) = 0.;
+           };
+       };
+
+                                      // preset solution vector
+      solution(dof_number) = dof->second;
+    };
+};
+
+
+
+
+
+void
+MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+                                   BlockSparseMatrix<double>  &matrix,
+                                   BlockVector<double>   &solution,
+                                   BlockVector<double>   &right_hand_side,
+                                   const bool             preserve_symmetry)
+{
+  const unsigned int blocks = matrix.n_block_rows();
+  
+  Assert (matrix.n() == matrix.m(),
+         ExcDimensionsDontMatch(matrix.n(), matrix.m()));
+  Assert (matrix.n() == right_hand_side.size(),
+         ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
+  Assert (matrix.n() == solution.size(),
+         ExcDimensionsDontMatch(matrix.n(), solution.size()));
+  Assert (matrix.n_block_rows() == matrix.n_block_cols(),
+         ExcMatrixNotBlockSquare());
+  Assert (matrix.get_sparsity_pattern().get_row_indices() == 
+         matrix.get_sparsity_pattern().get_column_indices(),
+         ExcMatrixNotBlockSquare());
+  Assert (matrix.get_sparsity_pattern().get_column_indices() ==
+         solution.get_block_indices (),
+         ExcBlocksDontMatch ());
+  Assert (matrix.get_sparsity_pattern().get_row_indices() ==
+         right_hand_side.get_block_indices (),
+         ExcBlocksDontMatch ());
+  
+  
+                                  // if no boundary values are to be applied
+                                  // simply return
+  if (boundary_values.size() == 0)
+    return;
+
+
+  std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
+                                               endd = boundary_values.end();
+  const unsigned int n_dofs = matrix.m();
+  const BlockSparsityPattern &
+    sparsity_pattern = matrix.get_sparsity_pattern();
+
+                                  // if a diagonal entry is zero
+                                  // later, then we use another
+                                  // number instead. take it to be
+                                  // the first nonzero diagonal
+                                  // element of the matrix, or 1 if
+                                  // there is no such thing
+  double first_nonzero_diagonal_entry = 0;
+  for (unsigned int diag_block=0; diag_block<blocks; ++diag_block)
+    {
+      for (unsigned int i=0; i<matrix.block(diag_block,diag_block).n(); ++i)
+       if (matrix.block(diag_block,diag_block).diag_element(i) != 0)
+         {
+           first_nonzero_diagonal_entry 
+             = matrix.block(diag_block,diag_block).diag_element(i);
+           break;
+         };
+                                      // check whether we have found
+                                      // something in the present
+                                      // block
+      if (first_nonzero_diagonal_entry != 0)
+       break;
+    };
+                                  // nothing found on all diagonal
+                                  // blocks? if so, use 1.0 instead
+  if (first_nonzero_diagonal_entry == 0)
+    first_nonzero_diagonal_entry = 1;
+  
+  
+                                  // pointer to the mapping between
+                                  // global and block indices. since
+                                  // the row and column mappings are
+                                  // equal, store a pointer on only
+                                  // one of them
+  const BlockIndices &
+    index_mapping = sparsity_pattern.get_column_indices();
+  
+                                  // now loop over all boundary dofs
+  for (; dof != endd; ++dof)
+    {
+      Assert (dof->first < n_dofs, ExcInternalError());
+
+                                      // get global index and index
+                                      // in the block in which this
+                                      // dof is located
+      const unsigned int dof_number = dof->first;
+      const std::pair<unsigned int,unsigned int>
+       block_index = index_mapping.global_to_local (dof_number);
+
+                                      // for each boundary dof:
+      
+                                      // set entries of this line
+                                      // to zero except for the diagonal
+                                      // entry. Note that the diagonal
+                                      // entry is always the first one
+                                      // for square matrices, i.e.
+                                      // we shall not set
+                                      // matrix.global_entry(
+                                      //     sparsity_rowstart[dof.first])
+                                      // of the diagonal block
+      for (unsigned int block_col=0; block_col<blocks; ++block_col)
+       {
+         const SparsityPattern &
+           local_sparsity = sparsity_pattern.block(block_index.first,
+                                                   block_col);
+
+                                          // find first and last
+                                          // entry in the present row
+                                          // of the present
+                                          // block. exclude the main
+                                          // diagonal element, which
+                                          // is the diagonal element
+                                          // of a diagonal block,
+                                          // which must be a square
+                                          // matrix so the diagonal
+                                          // element is the first of
+                                          // this row.
+         const unsigned int 
+           last  = local_sparsity.get_rowstart_indices()[block_index.second+1],
+           first = (block_col == block_index.first ?
+                    local_sparsity.get_rowstart_indices()[block_index.second]+1 :
+                    local_sparsity.get_rowstart_indices()[block_index.second]);
+         
+         for (unsigned int j=first; j<last; ++j)
+           matrix.block(block_index.first,block_col).global_entry(j) = 0.;
+       };
+      
+
+                                      // set right hand side to
+                                      // wanted value: if main diagonal
+                                      // entry nonzero, don't touch it
+                                      // and scale rhs accordingly. If
+                                      // zero, take the first main
+                                      // diagonal entry we can find, or
+                                      // one if no nonzero main diagonal
+                                      // element exists. Normally, however,
+                                      // the main diagonal entry should
+                                      // not be zero.
+                                      //
+                                      // store the new rhs entry to make
+                                      // the gauss step more efficient
+      double new_rhs;
+      if (matrix.block(block_index.first, block_index.first)
+         .diag_element(block_index.second) != 0.0)
+       new_rhs = dof->second * 
+                 matrix.block(block_index.first, block_index.first)
+                 .diag_element(block_index.second);
+      else
+       {
+         matrix.block(block_index.first, block_index.first)
+           .diag_element(block_index.second)
+           = first_nonzero_diagonal_entry;
+         new_rhs = dof->second * first_nonzero_diagonal_entry;
+       };
+      right_hand_side.block(block_index.first)(block_index.second)
+       = new_rhs;
+
+
+                                      // if the user wants to have
+                                      // the symmetry of the matrix
+                                      // preserved, and if the
+                                      // sparsity pattern is
+                                      // symmetric, then do a Gauss
+                                      // elimination step with the
+                                      // present row. this is a
+                                      // little more complicated for
+                                      // block matrices.
+      if (preserve_symmetry)
+       {
+                                          // store the only nonzero entry
+                                          // of this line for the Gauss
+                                          // elimination step
+         const double diagonal_entry 
+           = matrix.block(block_index.first,block_index.first)
+           .diag_element(block_index.second);
+         
+                                          // we have to loop over all
+                                          // rows of the matrix which
+                                          // have a nonzero entry in
+                                          // the column which we work
+                                          // in presently. if the
+                                          // sparsity pattern is
+                                          // symmetric, then we can
+                                          // get the positions of
+                                          // these rows cheaply by
+                                          // looking at the nonzero
+                                          // column numbers of the
+                                          // present row.
+                                          //
+                                          // note that if we check
+                                          // whether row @p{row} in
+                                          // block (r,c) is non-zero,
+                                          // then we have to check
+                                          // for the existence of
+                                          // column @p{row} in block
+                                          // (c,r), i.e. of the
+                                          // transpose block
+         for (unsigned int block_row=0; block_row<blocks; ++block_row)
+           {
+                                              // get pointers to the
+                                              // sparsity patterns of
+                                              // this block and of
+                                              // the transpose one
+             const SparsityPattern &this_sparsity
+               = sparsity_pattern.block (block_row, block_index.first);
+             const SparsityPattern &transpose_sparsity
+               = sparsity_pattern.block (block_index.first, block_row);
+             
+                                              // traverse the row of
+                                              // the transpose block
+                                              // to find the
+                                              // interesting rows in
+                                              // the present block.
+                                              // don't use the
+                                              // diagonal element of
+                                              // the diagonal block
+             const unsigned int
+               first = (block_index.first == block_row ?
+                        transpose_sparsity.get_rowstart_indices()[block_index.second]+1 :
+                        transpose_sparsity.get_rowstart_indices()[block_index.second]),
+               last  = transpose_sparsity.get_rowstart_indices()[block_index.second+1];
+             
+             for (unsigned int j=first; j<last; ++j)
+               {
+                                                  // get the number
+                                                  // of the column in
+                                                  // this row in
+                                                  // which a nonzero
+                                                  // entry is. this
+                                                  // is also the row
+                                                  // of the transpose
+                                                  // block which has
+                                                  // an entry in the
+                                                  // interesting row
+                 const unsigned int row = transpose_sparsity.get_column_numbers()[j];
+
+                                                  // find the
+                                                  // position of
+                                                  // element
+                                                  // (row,dof_number)
+                                                  // in this block
+                                                  // (not in the
+                                                  // transpose
+                                                  // one). note that
+                                                  // we have to take
+                                                  // care of special
+                                                  // cases with
+                                                  // square
+                                                  // sub-matrices
+                 const unsigned int *p = 0;
+                 if (this_sparsity.n_rows() == this_sparsity.n_cols())
+                   {
+                     if (this_sparsity.get_column_numbers()
+                         [this_sparsity.get_rowstart_indices()[row]]
+                         ==
+                         block_index.second)
+                       p = &this_sparsity.get_column_numbers()
+                           [this_sparsity.get_rowstart_indices()[row]];
+                     else
+                       p = std::lower_bound(&this_sparsity.get_column_numbers()
+                                            [this_sparsity.get_rowstart_indices()[row]+1],
+                                            &this_sparsity.get_column_numbers()
+                                            [this_sparsity.get_rowstart_indices()[row+1]],
+                                            block_index.second);
+                   }
+                 else
+                   p = std::lower_bound(&this_sparsity.get_column_numbers()
+                                        [this_sparsity.get_rowstart_indices()[row]],
+                                        &this_sparsity.get_column_numbers()
+                                        [this_sparsity.get_rowstart_indices()[row+1]],
+                                        block_index.second);
+
+                                                  // check whether this line has
+                                                  // an entry in the regarding column
+                                                  // (check for ==dof_number and
+                                                  // != next_row, since if
+                                                  // row==dof_number-1, *p is a
+                                                  // past-the-end pointer but points
+                                                  // to dof_number anyway...)
+                                                  //
+                                                  // there should be
+                                                  // such an entry!
+                                                  // note, however,
+                                                  // that this
+                                                  // assertion will
+                                                  // fail sometimes
+                                                  // if the sparsity
+                                                  // pattern is not
+                                                  // symmetric!
+                 Assert ((*p == block_index.second) &&
+                         (p != &this_sparsity.get_column_numbers()
+                          [this_sparsity.get_rowstart_indices()[row+1]]),
+                         ExcInternalError());
+                 
+                 const unsigned int global_entry
+                   = (p
+                      -
+                      &this_sparsity.get_column_numbers()
+                      [this_sparsity.get_rowstart_indices()[0]]);
+
+                                                  // correct right hand side
+                 right_hand_side.block(block_row)(row)
+                   -= matrix.block(block_row,block_index.first).global_entry(global_entry) /
+                   diagonal_entry * new_rhs;
+                 
+                                                  // set matrix entry to zero
+                 matrix.block(block_row,block_index.first).global_entry(global_entry) = 0.;
+               };
+           };
+       };
+
+                                      // preset solution vector
+      solution.block(block_index.first)(block_index.second) = dof->second;
+    };
+};
+
+
+
+
+// explicit instantiations
+
+template
+void
+MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+                                   SparseMatrix<double>  &matrix,
+                                   Vector<double>   &solution,
+                                   Vector<double>   &right_hand_side,
+                                   const bool        preserve_symmetry);
+template
+void
+MatrixTools::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+                                   SparseMatrix<float>  &matrix,
+                                   Vector<float>   &solution,
+                                   Vector<float>   &right_hand_side,
+                                   const bool        preserve_symmetry);
+
index a9ceafe83d716bb7ab189f8c650b522344bc369f..8d59abe53d3f7446fad785a3f26b4658254401b7 100644 (file)
@@ -13,6 +13,7 @@
 
 
 #include <base/function.h>
+#include <base/thread_management.h>
 #include <base/multithread_info.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
@@ -42,20 +43,44 @@ using namespace std;
 
 
 template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
-                                            const DoFHandler<dim>    &dof,
-                                            const Quadrature<dim>    &q,
-                                            SparseMatrix<double>     &matrix,
-                                            const Function<dim> * const coefficient)
+inline
+MatrixCreator::IteratorRange<dim>::
+IteratorRange (const active_cell_iterator &first,
+                                                 const active_cell_iterator &second)
+               :
+               first (first),
+               second (second)
+{};
+
+
+
+template <int dim>
+inline
+MatrixCreator::IteratorRange<dim>::IteratorRange (const iterator_pair &ip)
+               :
+               first (ip.first),
+               second (ip.second)
+{};
+
+
+
+
+template <int dim>
+void MatrixCreator::create_mass_matrix (const Mapping<dim>       &mapping,
+                                       const DoFHandler<dim>    &dof,
+                                       const Quadrature<dim>    &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const Function<dim> * const coefficient)
 {
   const unsigned int n_threads = multithread_info.n_default_threads;
   Threads::ThreadManager thread_manager;
 
                                   // define starting and end point
                                   // for each thread
-  std::vector<IteratorRange> thread_ranges
-    = Threads::split_range<active_cell_iterator> (dof.begin_active(),
-                                                 dof.end(), n_threads);
+  typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;  
+  const std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
+    thread_ranges = Threads::split_range<active_cell_iterator> (dof.begin_active(),
+                                                                dof.end(), n_threads);
 
                                   // mutex to synchronise access to
                                   // the matrix
@@ -64,7 +89,7 @@ void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
                                   // then assemble in parallel
   for (unsigned int thread=0; thread<n_threads; ++thread)
     Threads::spawn (thread_manager,
-                   Threads::encapsulate(&MatrixCreator<dim>::create_mass_matrix_1)
+                   Threads::encapsulate(&MatrixCreator::template create_mass_matrix_1<dim>)
                    .collect_args (mapping, dof, q, matrix, coefficient,
                                   thread_ranges[thread], mutex));
   thread_manager.wait ();  
@@ -74,13 +99,13 @@ void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
 
 // TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
 template <int dim>
-void MatrixCreator<dim>::create_mass_matrix_1 (const Mapping<dim>       &mapping,
-                                              const DoFHandler<dim>    &dof,
-                                              const Quadrature<dim>    &q,
-                                              SparseMatrix<double>     &matrix,
-                                              const Function<dim> * const coefficient,
-                                              const IteratorRange      &range,
-                                              Threads::ThreadMutex     &mutex)
+void MatrixCreator::create_mass_matrix_1 (const Mapping<dim>       &mapping,
+                                         const DoFHandler<dim>    &dof,
+                                         const Quadrature<dim>    &q,
+                                         SparseMatrix<double>     &matrix,
+                                         const Function<dim> * const coefficient,
+                                         const IteratorRange<dim>  range,
+                                         Threads::ThreadMutex     &mutex)
 {
   UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values);
   if (coefficient != 0)
@@ -151,10 +176,10 @@ void MatrixCreator<dim>::create_mass_matrix_1 (const Mapping<dim>       &mapping
 
 
 template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
-                                            const Quadrature<dim>    &q,
-                                            SparseMatrix<double>     &matrix,
-                                            const Function<dim> * const coefficient)
+void MatrixCreator::create_mass_matrix (const DoFHandler<dim>    &dof,
+                                       const Quadrature<dim>    &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const Function<dim> * const coefficient)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   static const MappingQ1<dim> mapping;
@@ -163,20 +188,21 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 
 
 template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
-                                            const DoFHandler<dim>    &dof,
-                                            const Quadrature<dim>    &q,
-                                            SparseMatrix<double>     &matrix,
-                                            const Function<dim>      &rhs,
-                                            Vector<double>           &rhs_vector,
-                                            const Function<dim> * const coefficient)
+void MatrixCreator::create_mass_matrix (const Mapping<dim>       &mapping,
+                                       const DoFHandler<dim>    &dof,
+                                       const Quadrature<dim>    &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const Function<dim>      &rhs,
+                                       Vector<double>           &rhs_vector,
+                                       const Function<dim> * const coefficient)
 {
   const unsigned int n_threads = multithread_info.n_default_threads;
   Threads::ThreadManager thread_manager;
 
                                   // define starting and end point
                                   // for each thread
-  std::vector<IteratorRange> thread_ranges
+  typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;  
+  std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
     = Threads::split_range<active_cell_iterator> (dof.begin_active(),
                                                  dof.end(), n_threads);
 
@@ -187,8 +213,8 @@ void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
                                   // then assemble in parallel
   for (unsigned int thread=0; thread<n_threads; ++thread)
     Threads::spawn (thread_manager,
-                   Threads::encapsulate(&MatrixCreator<dim>::
-                                        create_mass_matrix_2)
+                   Threads::encapsulate(&MatrixCreator::template
+                                        create_mass_matrix_2<dim>)
                    .collect_args (mapping, dof, q, matrix, rhs,
                                   rhs_vector, coefficient,
                                   thread_ranges[thread], mutex));
@@ -200,15 +226,15 @@ void MatrixCreator<dim>::create_mass_matrix (const Mapping<dim>       &mapping,
 // TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
 template <int dim>
 void
-MatrixCreator<dim>::create_mass_matrix_2 (const Mapping<dim>       &mapping,
-                                         const DoFHandler<dim>    &dof,
-                                         const Quadrature<dim>    &q,
-                                         SparseMatrix<double>     &matrix,
-                                         const Function<dim>      &rhs,
-                                         Vector<double>           &rhs_vector,
-                                         const Function<dim> * const coefficient,
-                                         const IteratorRange      &range,
-                                         Threads::ThreadMutex     &mutex)
+MatrixCreator::create_mass_matrix_2 (const Mapping<dim>       &mapping,
+                                    const DoFHandler<dim>    &dof,
+                                    const Quadrature<dim>    &q,
+                                    SparseMatrix<double>     &matrix,
+                                    const Function<dim>      &rhs,
+                                    Vector<double>           &rhs_vector,
+                                    const Function<dim> * const coefficient,
+                                    const IteratorRange<dim>  range,
+                                    Threads::ThreadMutex     &mutex)
 {
   UpdateFlags update_flags = UpdateFlags(update_values |
                                         update_q_points |
@@ -297,12 +323,12 @@ MatrixCreator<dim>::create_mass_matrix_2 (const Mapping<dim>       &mapping,
 
 
 template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
-                                            const Quadrature<dim>    &q,
-                                            SparseMatrix<double>     &matrix,
-                                            const Function<dim>      &rhs,
-                                            Vector<double>           &rhs_vector,
-                                            const Function<dim> * const coefficient)
+void MatrixCreator::create_mass_matrix (const DoFHandler<dim>    &dof,
+                                       const Quadrature<dim>    &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const Function<dim>      &rhs,
+                                       Vector<double>           &rhs_vector,
+                                       const Function<dim> * const coefficient)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   static const MappingQ1<dim> mapping;
@@ -312,18 +338,17 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 
 #if deal_II_dimension == 1
 
-template <>
-void MatrixCreator<1>::create_boundary_mass_matrix (const Mapping<1>          &,
-                                                   const DoFHandler<1>       &,
-                                                   const Quadrature<0>       &,
-                                                   SparseMatrix<double>      &,
-                                                   const FunctionMap<1>::type&,
-                                                   Vector<double>            &,
-                                                   std::vector<unsigned int> &,
-                                                   const Function<1>         * const)
+void MatrixCreator::create_boundary_mass_matrix (const Mapping<1>          &,
+                                                const DoFHandler<1>       &,
+                                                const Quadrature<0>       &,
+                                                SparseMatrix<double>      &,
+                                                const FunctionMap<1>::type&,
+                                                Vector<double>            &,
+                                                std::vector<unsigned int> &,
+                                                const Function<1>         * const)
 {
                                   // what would that be in 1d? the
-                                  // identity matrix on tehe boundary
+                                  // identity matrix on the boundary
                                   // dofs?
   Assert (false, ExcNotImplemented());
 };
@@ -335,21 +360,22 @@ void MatrixCreator<1>::create_boundary_mass_matrix (const Mapping<1>          &,
 // TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
 template <int dim>
 void
-MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>        &mapping,
-                                                const DoFHandler<dim>     &dof,
-                                                const Quadrature<dim-1>   &q,
-                                                SparseMatrix<double>      &matrix,
-                                                const typename FunctionMap<dim>::type         &boundary_functions,
-                                                Vector<double>            &rhs_vector,
-                                                std::vector<unsigned int> &dof_to_boundary_mapping,
-                                                const Function<dim> * const a)
+MatrixCreator::create_boundary_mass_matrix (const Mapping<dim>        &mapping,
+                                           const DoFHandler<dim>     &dof,
+                                           const Quadrature<dim-1>   &q,
+                                           SparseMatrix<double>      &matrix,
+                                           const typename FunctionMap<dim>::type         &boundary_functions,
+                                           Vector<double>            &rhs_vector,
+                                           std::vector<unsigned int> &dof_to_boundary_mapping,
+                                           const Function<dim> * const a)
 {
   const unsigned int n_threads = multithread_info.n_default_threads;
   Threads::ThreadManager thread_manager;
 
                                   // define starting and end point
                                   // for each thread
-  std::vector<IteratorRange> thread_ranges
+  typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;  
+  std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
     = Threads::split_range<active_cell_iterator> (dof.begin_active(),
                                                  dof.end(), n_threads);
 
@@ -360,8 +386,8 @@ MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>        &mapp
                                   // then assemble in parallel
   for (unsigned int thread=0; thread<n_threads; ++thread)
     Threads::spawn (thread_manager,
-                   Threads::encapsulate(&MatrixCreator<dim>::
-                                        create_boundary_mass_matrix_1)
+                   Threads::encapsulate(&MatrixCreator::template 
+                                        create_boundary_mass_matrix_1<dim>)
                    .collect_args (mapping, dof, q, matrix,
                                   boundary_functions, rhs_vector,
                                   dof_to_boundary_mapping, a,
@@ -370,31 +396,10 @@ MatrixCreator<dim>::create_boundary_mass_matrix (const Mapping<dim>        &mapp
 };
 
 
-#if deal_II_dimension == 1
-
-template <>
-void
-MatrixCreator<1>::
-create_boundary_mass_matrix_1 (const Mapping<1>          &,
-                              const DoFHandler<1>       &,
-                              const Quadrature<0>       &,
-                              SparseMatrix<double>      &,
-                              const FunctionMap<1>::type&,
-                              Vector<double>            &,
-                              std::vector<unsigned int> &,
-                              const Function<1> * const ,
-                              const IteratorRange       &,
-                              Threads::ThreadMutex      &)
-{
-  Assert (false, ExcNotImplemented());
-};
-
-#endif
-
 
 template <int dim>
 void
-MatrixCreator<dim>::
+MatrixCreator::
 create_boundary_mass_matrix_1 (const Mapping<dim>        &mapping,
                               const DoFHandler<dim>     &dof,
                               const Quadrature<dim-1>   &q,
@@ -403,7 +408,7 @@ create_boundary_mass_matrix_1 (const Mapping<dim>        &mapping,
                               Vector<double>            &rhs_vector,
                               std::vector<unsigned int> &dof_to_boundary_mapping,
                               const Function<dim> * const a,
-                              const IteratorRange       &range,
+                              const IteratorRange<dim>   range,
                               Threads::ThreadMutex      &mutex)
 {
   const FiniteElement<dim> &fe = dof.get_fe();
@@ -678,8 +683,8 @@ create_boundary_mass_matrix_1 (const Mapping<dim>        &mapping,
              rhs_vector(dof_to_boundary_mapping[dofs[j]]) += cell_vector(j);
            else
              {
-                                                  // compare here for relative
-                                                  // smallness
+                                                // compare here for relative
+                                                // smallness
                Assert (std::fabs(cell_vector(j)) <= 1e-10 * max_diag_entry,
                        ExcInternalError());
              };
@@ -689,13 +694,13 @@ create_boundary_mass_matrix_1 (const Mapping<dim>        &mapping,
 
 
 template <int dim>
-void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>     &dof,
-                                                     const Quadrature<dim-1>   &q,
-                                                     SparseMatrix<double>      &matrix,
-                                                     const typename FunctionMap<dim>::type &rhs,
-                                                     Vector<double>            &rhs_vector,
-                                                     std::vector<unsigned int> &dof_to_boundary_mapping,
-                                                     const Function<dim> * const a)
+void MatrixCreator::create_boundary_mass_matrix (const DoFHandler<dim>     &dof,
+                                                const Quadrature<dim-1>   &q,
+                                                SparseMatrix<double>      &matrix,
+                                                const typename FunctionMap<dim>::type &rhs,
+                                                Vector<double>            &rhs_vector,
+                                                std::vector<unsigned int> &dof_to_boundary_mapping,
+                                                const Function<dim> * const a)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   static const MappingQ1<dim> mapping;
@@ -706,18 +711,19 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>
 
 
 template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mapping,
-                                               const DoFHandler<dim>    &dof,
-                                               const Quadrature<dim>    &q,
-                                               SparseMatrix<double>     &matrix,
-                                               const Function<dim> * const coefficient)
+void MatrixCreator::create_laplace_matrix (const Mapping<dim>       &mapping,
+                                          const DoFHandler<dim>    &dof,
+                                          const Quadrature<dim>    &q,
+                                          SparseMatrix<double>     &matrix,
+                                          const Function<dim> * const coefficient)
 {
   const unsigned int n_threads = multithread_info.n_default_threads;
   Threads::ThreadManager thread_manager;
 
                                   // define starting and end point
                                   // for each thread
-  std::vector<IteratorRange> thread_ranges
+  typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;  
+  std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
     = Threads::split_range<active_cell_iterator> (dof.begin_active(),
                                                  dof.end(), n_threads);
 
@@ -728,7 +734,8 @@ void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mappin
                                   // then assemble in parallel
   for (unsigned int thread=0; thread<n_threads; ++thread)
     Threads::spawn (thread_manager,
-                   Threads::encapsulate(&MatrixCreator<dim>::create_laplace_matrix_1)
+                   Threads::encapsulate(&MatrixCreator::template
+                                        create_laplace_matrix_1<dim>)
                    .collect_args (mapping, dof, q, matrix, coefficient,
                                   thread_ranges[thread], mutex));
   thread_manager.wait ();  
@@ -738,13 +745,13 @@ void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mappin
 
 // TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
 template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix_1 (const Mapping<dim>       &mapping,
-                                                 const DoFHandler<dim>    &dof,
-                                                 const Quadrature<dim>    &q,
-                                                 SparseMatrix<double>     &matrix,
-                                                 const Function<dim> * const coefficient,
-                                                 const IteratorRange      &range,
-                                                 Threads::ThreadMutex     &mutex)
+void MatrixCreator::create_laplace_matrix_1 (const Mapping<dim>       &mapping,
+                                            const DoFHandler<dim>    &dof,
+                                            const Quadrature<dim>    &q,
+                                            SparseMatrix<double>     &matrix,
+                                            const Function<dim> * const coefficient,
+                                            const IteratorRange<dim>  range,
+                                            Threads::ThreadMutex     &mutex)
 {
   UpdateFlags update_flags = UpdateFlags(update_JxW_values |
                                         update_gradients);
@@ -818,10 +825,10 @@ void MatrixCreator<dim>::create_laplace_matrix_1 (const Mapping<dim>       &mapp
 
 
 template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
-                                               const Quadrature<dim>    &q,
-                                               SparseMatrix<double>     &matrix,
-                                               const Function<dim> * const coefficient)
+void MatrixCreator::create_laplace_matrix (const DoFHandler<dim>    &dof,
+                                          const Quadrature<dim>    &q,
+                                          SparseMatrix<double>     &matrix,
+                                          const Function<dim> * const coefficient)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   static const MappingQ1<dim> mapping;
@@ -832,55 +839,56 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
 //TODO:[GK,RH] maybe recreate this function
 /*
 
-template <int dim>
-void MatrixCreator<dim>::create_level_laplace_matrix (unsigned int level,
-                                                     const MGDoFHandler<dim>    &dof,
-                                                     const Quadrature<dim>    &q,
-                                                     SparseMatrix<float>     &matrix,
-                                                     const Function<dim> * const a)
-{
+  template <int dim>
+  void MatrixCreator::create_level_laplace_matrix (unsigned int level,
+  const MGDoFHandler<dim>    &dof,
+  const Quadrature<dim>    &q,
+  SparseMatrix<float>     &matrix,
+  const Function<dim> * const a)
+  {
   Vector<double> dummy;   // no entries, should give an error if accessed
   UpdateFlags update_flags = UpdateFlags(update_gradients |
-                                        update_JxW_values);
+  update_JxW_values);
   if (a != 0)
-    update_flags = UpdateFlags(update_flags | update_q_points);
+  update_flags = UpdateFlags(update_flags | update_q_points);
   const Assembler<dim>::AssemblerData data (dof,
-                                           true, false,  // assemble matrix but not rhs
-                                           matrix, dummy,
-                                           q, update_flags);
+  true, false,  // assemble matrix but not rhs
+  matrix, dummy,
+  q, update_flags);
   TriaIterator<dim, Assembler<dim> >
-    assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
-              dof.get_tria().begin(level)->level(),
-              dof.get_tria().begin(level)->index(),
-              &data);
+  assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
+  dof.get_tria().begin(level)->level(),
+  dof.get_tria().begin(level)->index(),
+  &data);
   LaplaceMatrix<dim> equation (0, a);
   do 
-    {
-      assembler->assemble (equation);
-      ++assembler
-    }
+  {
+  assembler->assemble (equation);
+  ++assembler
+  }
   while ( (assembler.state()==valid) && (assembler->level() == level) );
-};
+  };
 
 */
 
 
 
 template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mapping,
-                                               const DoFHandler<dim>    &dof,
-                                               const Quadrature<dim>    &q,
-                                               SparseMatrix<double>     &matrix,
-                                               const Function<dim>      &rhs,
-                                               Vector<double>           &rhs_vector,
-                                               const Function<dim> * const coefficient)
+void MatrixCreator::create_laplace_matrix (const Mapping<dim>       &mapping,
+                                          const DoFHandler<dim>    &dof,
+                                          const Quadrature<dim>    &q,
+                                          SparseMatrix<double>     &matrix,
+                                          const Function<dim>      &rhs,
+                                          Vector<double>           &rhs_vector,
+                                          const Function<dim> * const coefficient)
 {
   const unsigned int n_threads = multithread_info.n_default_threads;
   Threads::ThreadManager thread_manager;
 
                                   // define starting and end point
                                   // for each thread
-  std::vector<IteratorRange> thread_ranges
+  typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;  
+  std::vector<std::pair<active_cell_iterator,active_cell_iterator> > thread_ranges
     = Threads::split_range<active_cell_iterator> (dof.begin_active(),
                                                  dof.end(), n_threads);
 
@@ -891,8 +899,8 @@ void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mappin
                                   // then assemble in parallel
   for (unsigned int thread=0; thread<n_threads; ++thread)
     Threads::spawn (thread_manager,
-                   Threads::encapsulate(&MatrixCreator<dim>::
-                                        create_laplace_matrix_2)
+                   Threads::encapsulate(&MatrixCreator::template
+                                        create_laplace_matrix_2<dim>)
                    .collect_args (mapping, dof, q, matrix, rhs,
                                   rhs_vector, coefficient,
                                   thread_ranges[thread], mutex));
@@ -904,15 +912,15 @@ void MatrixCreator<dim>::create_laplace_matrix (const Mapping<dim>       &mappin
 // TODO:[RH] extend this function to use vector valued coefficient functions for system elements.
 template <int dim>
 void
-MatrixCreator<dim>::create_laplace_matrix_2 (const Mapping<dim>       &mapping,
-                                            const DoFHandler<dim>    &dof,
-                                            const Quadrature<dim>    &q,
-                                            SparseMatrix<double>     &matrix,
-                                            const Function<dim>      &rhs,
-                                            Vector<double>           &rhs_vector,
-                                            const Function<dim> * const coefficient,
-                                            const IteratorRange      &range,
-                                            Threads::ThreadMutex     &mutex)
+MatrixCreator::create_laplace_matrix_2 (const Mapping<dim>       &mapping,
+                                       const DoFHandler<dim>    &dof,
+                                       const Quadrature<dim>    &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const Function<dim>      &rhs,
+                                       Vector<double>           &rhs_vector,
+                                       const Function<dim> * const coefficient,
+                                       const IteratorRange<dim>  range,
+                                       Threads::ThreadMutex     &mutex)
 {
   UpdateFlags update_flags = UpdateFlags(update_values    |
                                         update_gradients |
@@ -1005,12 +1013,12 @@ MatrixCreator<dim>::create_laplace_matrix_2 (const Mapping<dim>       &mapping,
 
 
 template <int dim>
-void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
-                                               const Quadrature<dim>    &q,
-                                               SparseMatrix<double>     &matrix,
-                                               const Function<dim>      &rhs,
-                                               Vector<double>           &rhs_vector,
-                                               const Function<dim> * const coefficient)
+void MatrixCreator::create_laplace_matrix (const DoFHandler<dim>    &dof,
+                                          const Quadrature<dim>    &q,
+                                          SparseMatrix<double>     &matrix,
+                                          const Function<dim>      &rhs,
+                                          Vector<double>           &rhs_vector,
+                                          const Function<dim> * const coefficient)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   static const MappingQ1<dim> mapping;  
@@ -1019,517 +1027,86 @@ void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
 
 
 
-template <int dim>
-template <typename number>
-void
-MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
-                                        SparseMatrix<number>  &matrix,
-                                        Vector<number>   &solution,
-                                        Vector<number>   &right_hand_side,
-                                        const bool        preserve_symmetry)
-{
-  Assert (matrix.n() == matrix.m(),
-         ExcDimensionsDontMatch(matrix.n(), matrix.m()));
-  Assert (matrix.n() == right_hand_side.size(),
-         ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
-  Assert (matrix.n() == solution.size(),
-         ExcDimensionsDontMatch(matrix.n(), solution.size()));
-                                  // if no boundary values are to be applied
-                                  // simply return
-  if (boundary_values.size() == 0)
-    return;
-
-
-  std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
-                                               endd = boundary_values.end();
-  const unsigned int n_dofs             = matrix.m();
-  const SparsityPattern    &sparsity    = matrix.get_sparsity_pattern();
-  const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
-  const unsigned int *sparsity_colnums  = sparsity.get_column_numbers();
-
-                                  // if a diagonal entry is zero
-                                  // later, then we use another
-                                  // number instead. take it to be
-                                  // the first nonzero diagonal
-                                  // element of the matrix, or 1 if
-                                  // there is no such thing
-  number first_nonzero_diagonal_entry = 1;
-  for (unsigned int i=0; i<n_dofs; ++i)
-    if (matrix.diag_element(i) != 0)
-      {
-       first_nonzero_diagonal_entry = matrix.diag_element(i);
-       break;
-      };
-
-  
-  for (; dof != endd; ++dof)
-    {
-      Assert (dof->first < n_dofs, ExcInternalError());
-      
-      const unsigned int dof_number = dof->first;
-                                      // for each boundary dof:
-      
-                                      // set entries of this line
-                                      // to zero except for the diagonal
-                                      // entry. Note that the diagonal
-                                      // entry is always the first one
-                                      // for square matrices, i.e.
-                                      // we shall not set
-                                      // matrix.global_entry(
-                                      //     sparsity_rowstart[dof.first])
-      const unsigned int last = sparsity_rowstart[dof_number+1];
-      for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
-       matrix.global_entry(j) = 0.;
-
-
-                                      // set right hand side to
-                                      // wanted value: if main diagonal
-                                      // entry nonzero, don't touch it
-                                      // and scale rhs accordingly. If
-                                      // zero, take the first main
-                                      // diagonal entry we can find, or
-                                      // one if no nonzero main diagonal
-                                      // element exists. Normally, however,
-                                      // the main diagonal entry should
-                                      // not be zero.
-                                      //
-                                      // store the new rhs entry to make
-                                      // the gauss step more efficient
-      number new_rhs;
-      if (matrix.diag_element(dof_number) != 0.0)
-       new_rhs = right_hand_side(dof_number)
-               = dof->second * matrix.diag_element(dof_number);
-      else
-       {
-                                          // use the SparseMatrix::
-                                          // to work around a bug in
-                                          // egcs
-         matrix.SparseMatrix<number>::set (dof_number, dof_number,
-                                           first_nonzero_diagonal_entry);
-         new_rhs = right_hand_side(dof_number)
-                 = dof->second * first_nonzero_diagonal_entry;
-       };
-
-
-                                      // if the user wants to have
-                                      // the symmetry of the matrix
-                                      // preserved, and if the
-                                      // sparsity pattern is
-                                      // symmetric, then do a Gauss
-                                      // elimination step with the
-                                      // present row
-      if (preserve_symmetry)
-       {
-                                          // store the only nonzero entry
-                                          // of this line for the Gauss
-                                          // elimination step
-         const number diagonal_entry = matrix.diag_element(dof_number);
-         
-                                          // we have to loop over all
-                                          // rows of the matrix which
-                                          // have a nonzero entry in
-                                          // the column which we work
-                                          // in presently. if the
-                                          // sparsity pattern is
-                                          // symmetric, then we can
-                                          // get the positions of
-                                          // these rows cheaply by
-                                          // looking at the nonzero
-                                          // column numbers of the
-                                          // present row. we need not
-                                          // look at the first entry,
-                                          // since that is the
-                                          // diagonal element and
-                                          // thus the present row
-         for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
-           {
-             const unsigned int row = sparsity_colnums[j];
-
-                                              // find the position of
-                                              // element
-                                              // (row,dof_number)
-             const unsigned int *
-               p = std::lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
-                                    &sparsity_colnums[sparsity_rowstart[row+1]],
-                                    dof_number);
-
-                                              // check whether this line has
-                                              // an entry in the regarding column
-                                              // (check for ==dof_number and
-                                              // != next_row, since if
-                                              // row==dof_number-1, *p is a
-                                              // past-the-end pointer but points
-                                              // to dof_number anyway...)
-                                              //
-                                              // there should be such an entry!
-             Assert ((*p == dof_number) &&
-                     (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
-                     ExcInternalError());
-
-             const unsigned int global_entry
-               = (p - &sparsity_colnums[sparsity_rowstart[0]]);
-             
-                                              // correct right hand side
-             right_hand_side(row) -= matrix.global_entry(global_entry) /
-                                     diagonal_entry * new_rhs;
-             
-                                              // set matrix entry to zero
-             matrix.global_entry(global_entry) = 0.;
-           };
-       };
-
-                                      // preset solution vector
-      solution(dof_number) = dof->second;
-    };
-};
-
-
-
-
-
-template <int dim>
-void
-MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
-                                        BlockSparseMatrix<double>  &matrix,
-                                        BlockVector<double>   &solution,
-                                        BlockVector<double>   &right_hand_side,
-                                        const bool             preserve_symmetry)
-{
-  const unsigned int blocks = matrix.n_block_rows();
-  
-  Assert (matrix.n() == matrix.m(),
-         ExcDimensionsDontMatch(matrix.n(), matrix.m()));
-  Assert (matrix.n() == right_hand_side.size(),
-         ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
-  Assert (matrix.n() == solution.size(),
-         ExcDimensionsDontMatch(matrix.n(), solution.size()));
-  Assert (matrix.n_block_rows() == matrix.n_block_cols(),
-         ExcMatrixNotBlockSquare());
-  Assert (matrix.get_sparsity_pattern().get_row_indices() == 
-         matrix.get_sparsity_pattern().get_column_indices(),
-         ExcMatrixNotBlockSquare());
-  Assert (matrix.get_sparsity_pattern().get_column_indices() ==
-         solution.get_block_indices (),
-         ExcBlocksDontMatch ());
-  Assert (matrix.get_sparsity_pattern().get_row_indices() ==
-         right_hand_side.get_block_indices (),
-         ExcBlocksDontMatch ());
-  
-  
-                                  // if no boundary values are to be applied
-                                  // simply return
-  if (boundary_values.size() == 0)
-    return;
-
-
-  std::map<unsigned int,double>::const_iterator dof  = boundary_values.begin(),
-                                               endd = boundary_values.end();
-  const unsigned int n_dofs = matrix.m();
-  const BlockSparsityPattern &
-    sparsity_pattern = matrix.get_sparsity_pattern();
-
-                                  // if a diagonal entry is zero
-                                  // later, then we use another
-                                  // number instead. take it to be
-                                  // the first nonzero diagonal
-                                  // element of the matrix, or 1 if
-                                  // there is no such thing
-  double first_nonzero_diagonal_entry = 0;
-  for (unsigned int diag_block=0; diag_block<blocks; ++diag_block)
-    {
-      for (unsigned int i=0; i<matrix.block(diag_block,diag_block).n(); ++i)
-       if (matrix.block(diag_block,diag_block).diag_element(i) != 0)
-         {
-           first_nonzero_diagonal_entry 
-             = matrix.block(diag_block,diag_block).diag_element(i);
-           break;
-         };
-                                      // check whether we have found
-                                      // something in the present
-                                      // block
-      if (first_nonzero_diagonal_entry != 0)
-       break;
-    };
-                                  // nothing found on all diagonal
-                                  // blocks? if so, use 1.0 instead
-  if (first_nonzero_diagonal_entry == 0)
-    first_nonzero_diagonal_entry = 1;
-  
-  
-                                  // pointer to the mapping between
-                                  // global and block indices. since
-                                  // the row and column mappings are
-                                  // equal, store a pointer on only
-                                  // one of them
-  const BlockIndices &
-    index_mapping = sparsity_pattern.get_column_indices();
-  
-                                  // now loop over all boundary dofs
-  for (; dof != endd; ++dof)
-    {
-      Assert (dof->first < n_dofs, ExcInternalError());
-
-                                      // get global index and index
-                                      // in the block in which this
-                                      // dof is located
-      const unsigned int dof_number = dof->first;
-      const std::pair<unsigned int,unsigned int>
-       block_index = index_mapping.global_to_local (dof_number);
-
-                                      // for each boundary dof:
-      
-                                      // set entries of this line
-                                      // to zero except for the diagonal
-                                      // entry. Note that the diagonal
-                                      // entry is always the first one
-                                      // for square matrices, i.e.
-                                      // we shall not set
-                                      // matrix.global_entry(
-                                      //     sparsity_rowstart[dof.first])
-                                      // of the diagonal block
-      for (unsigned int block_col=0; block_col<blocks; ++block_col)
-       {
-         const SparsityPattern &
-           local_sparsity = sparsity_pattern.block(block_index.first,
-                                                   block_col);
-
-                                          // find first and last
-                                          // entry in the present row
-                                          // of the present
-                                          // block. exclude the main
-                                          // diagonal element, which
-                                          // is the diagonal element
-                                          // of a diagonal block,
-                                          // which must be a square
-                                          // matrix so the diagonal
-                                          // element is the first of
-                                          // this row.
-         const unsigned int 
-           last  = local_sparsity.get_rowstart_indices()[block_index.second+1],
-           first = (block_col == block_index.first ?
-                    local_sparsity.get_rowstart_indices()[block_index.second]+1 :
-                    local_sparsity.get_rowstart_indices()[block_index.second]);
-         
-         for (unsigned int j=first; j<last; ++j)
-           matrix.block(block_index.first,block_col).global_entry(j) = 0.;
-       };
-      
-
-                                      // set right hand side to
-                                      // wanted value: if main diagonal
-                                      // entry nonzero, don't touch it
-                                      // and scale rhs accordingly. If
-                                      // zero, take the first main
-                                      // diagonal entry we can find, or
-                                      // one if no nonzero main diagonal
-                                      // element exists. Normally, however,
-                                      // the main diagonal entry should
-                                      // not be zero.
-                                      //
-                                      // store the new rhs entry to make
-                                      // the gauss step more efficient
-      double new_rhs;
-      if (matrix.block(block_index.first, block_index.first)
-         .diag_element(block_index.second) != 0.0)
-       new_rhs = dof->second * 
-                 matrix.block(block_index.first, block_index.first)
-                 .diag_element(block_index.second);
-      else
-       {
-         matrix.block(block_index.first, block_index.first)
-           .diag_element(block_index.second)
-           = first_nonzero_diagonal_entry;
-         new_rhs = dof->second * first_nonzero_diagonal_entry;
-       };
-      right_hand_side.block(block_index.first)(block_index.second)
-       = new_rhs;
-
-
-                                      // if the user wants to have
-                                      // the symmetry of the matrix
-                                      // preserved, and if the
-                                      // sparsity pattern is
-                                      // symmetric, then do a Gauss
-                                      // elimination step with the
-                                      // present row. this is a
-                                      // little more complicated for
-                                      // block matrices.
-      if (preserve_symmetry)
-       {
-                                          // store the only nonzero entry
-                                          // of this line for the Gauss
-                                          // elimination step
-         const double diagonal_entry 
-           = matrix.block(block_index.first,block_index.first)
-           .diag_element(block_index.second);
-         
-                                          // we have to loop over all
-                                          // rows of the matrix which
-                                          // have a nonzero entry in
-                                          // the column which we work
-                                          // in presently. if the
-                                          // sparsity pattern is
-                                          // symmetric, then we can
-                                          // get the positions of
-                                          // these rows cheaply by
-                                          // looking at the nonzero
-                                          // column numbers of the
-                                          // present row.
-                                          //
-                                          // note that if we check
-                                          // whether row @p{row} in
-                                          // block (r,c) is non-zero,
-                                          // then we have to check
-                                          // for the existence of
-                                          // column @p{row} in block
-                                          // (c,r), i.e. of the
-                                          // transpose block
-         for (unsigned int block_row=0; block_row<blocks; ++block_row)
-           {
-                                              // get pointers to the
-                                              // sparsity patterns of
-                                              // this block and of
-                                              // the transpose one
-             const SparsityPattern &this_sparsity
-               = sparsity_pattern.block (block_row, block_index.first);
-             const SparsityPattern &transpose_sparsity
-               = sparsity_pattern.block (block_index.first, block_row);
-             
-                                              // traverse the row of
-                                              // the transpose block
-                                              // to find the
-                                              // interesting rows in
-                                              // the present block.
-                                              // don't use the
-                                              // diagonal element of
-                                              // the diagonal block
-             const unsigned int
-               first = (block_index.first == block_row ?
-                        transpose_sparsity.get_rowstart_indices()[block_index.second]+1 :
-                        transpose_sparsity.get_rowstart_indices()[block_index.second]),
-               last  = transpose_sparsity.get_rowstart_indices()[block_index.second+1];
-             
-             for (unsigned int j=first; j<last; ++j)
-               {
-                                                  // get the number
-                                                  // of the column in
-                                                  // this row in
-                                                  // which a nonzero
-                                                  // entry is. this
-                                                  // is also the row
-                                                  // of the transpose
-                                                  // block which has
-                                                  // an entry in the
-                                                  // interesting row
-                 const unsigned int row = transpose_sparsity.get_column_numbers()[j];
-
-                                                  // find the
-                                                  // position of
-                                                  // element
-                                                  // (row,dof_number)
-                                                  // in this block
-                                                  // (not in the
-                                                  // transpose
-                                                  // one). note that
-                                                  // we have to take
-                                                  // care of special
-                                                  // cases with
-                                                  // square
-                                                  // sub-matrices
-                 const unsigned int *p = 0;
-                 if (this_sparsity.n_rows() == this_sparsity.n_cols())
-                   {
-                     if (this_sparsity.get_column_numbers()
-                         [this_sparsity.get_rowstart_indices()[row]]
-                         ==
-                         block_index.second)
-                       p = &this_sparsity.get_column_numbers()
-                           [this_sparsity.get_rowstart_indices()[row]];
-                     else
-                       p = std::lower_bound(&this_sparsity.get_column_numbers()
-                                            [this_sparsity.get_rowstart_indices()[row]+1],
-                                            &this_sparsity.get_column_numbers()
-                                            [this_sparsity.get_rowstart_indices()[row+1]],
-                                            block_index.second);
-                   }
-                 else
-                   p = std::lower_bound(&this_sparsity.get_column_numbers()
-                                        [this_sparsity.get_rowstart_indices()[row]],
-                                        &this_sparsity.get_column_numbers()
-                                        [this_sparsity.get_rowstart_indices()[row+1]],
-                                        block_index.second);
-
-                                                  // check whether this line has
-                                                  // an entry in the regarding column
-                                                  // (check for ==dof_number and
-                                                  // != next_row, since if
-                                                  // row==dof_number-1, *p is a
-                                                  // past-the-end pointer but points
-                                                  // to dof_number anyway...)
-                                                  //
-                                                  // there should be
-                                                  // such an entry!
-                                                  // note, however,
-                                                  // that this
-                                                  // assertion will
-                                                  // fail sometimes
-                                                  // if the sparsity
-                                                  // pattern is not
-                                                  // symmetric!
-                 Assert ((*p == block_index.second) &&
-                         (p != &this_sparsity.get_column_numbers()
-                          [this_sparsity.get_rowstart_indices()[row+1]]),
-                         ExcInternalError());
-                 
-                 const unsigned int global_entry
-                   = (p
-                      -
-                      &this_sparsity.get_column_numbers()
-                      [this_sparsity.get_rowstart_indices()[0]]);
-
-                                                  // correct right hand side
-                 right_hand_side.block(block_row)(row)
-                   -= matrix.block(block_row,block_index.first).global_entry(global_entry) /
-                   diagonal_entry * new_rhs;
-                 
-                                                  // set matrix entry to zero
-                 matrix.block(block_row,block_index.first).global_entry(global_entry) = 0.;
-               };
-           };
-       };
-
-                                      // preset solution vector
-      solution.block(block_index.first)(block_index.second) = dof->second;
-    };
-};
-
-
 
 
 
 // explicit instantiations
 
-template class MatrixCreator<deal_II_dimension>;
-template class MatrixTools<deal_II_dimension>;
-
-
+template
+void MatrixCreator::create_mass_matrix (const Mapping<deal_II_dimension>       &mapping,
+                                       const DoFHandler<deal_II_dimension>    &dof,
+                                       const Quadrature<deal_II_dimension>    &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_mass_matrix (const DoFHandler<deal_II_dimension>    &dof,
+                                       const Quadrature<deal_II_dimension>    &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_mass_matrix (const Mapping<deal_II_dimension>       &mapping,
+                                       const DoFHandler<deal_II_dimension>    &dof,
+                                       const Quadrature<deal_II_dimension>    &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const Function<deal_II_dimension>      &rhs,
+                                       Vector<double>           &rhs_vector,
+                                       const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_mass_matrix (const DoFHandler<deal_II_dimension>    &dof,
+                                       const Quadrature<deal_II_dimension>    &q,
+                                       SparseMatrix<double>     &matrix,
+                                       const Function<deal_II_dimension>      &rhs,
+                                       Vector<double>           &rhs_vector,
+                                       const Function<deal_II_dimension> * const coefficient);
+
+#if deal_II_dimension != 1
 template
 void
-MatrixTools<deal_II_dimension>::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
-                      SparseMatrix<double>  &matrix,
-                      Vector<double>   &solution,
-                      Vector<double>   &right_hand_side,
-                      const bool        preserve_symmetry);
+MatrixCreator::create_boundary_mass_matrix (const Mapping<deal_II_dimension>        &mapping,
+                                           const DoFHandler<deal_II_dimension>     &dof,
+                                           const Quadrature<deal_II_dimension-1>   &q,
+                                           SparseMatrix<double>      &matrix,
+                                           const FunctionMap<deal_II_dimension>::type         &boundary_functions,
+                                           Vector<double>            &rhs_vector,
+                                           std::vector<unsigned int> &dof_to_boundary_mapping,
+                                           const Function<deal_II_dimension> * const a);
+#endif
 
 template
-void
-MatrixTools<deal_II_dimension>::
-apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
-                      SparseMatrix<float>  &matrix,
-                      Vector<float>   &solution,
-                      Vector<float>   &right_hand_side,
-                      const bool       preserve_symmetry);
+void MatrixCreator::create_boundary_mass_matrix (const DoFHandler<deal_II_dimension>     &dof,
+                                                const Quadrature<deal_II_dimension-1>   &q,
+                                                SparseMatrix<double>      &matrix,
+                                                const FunctionMap<deal_II_dimension>::type &rhs,
+                                                Vector<double>            &rhs_vector,
+                                                std::vector<unsigned int> &dof_to_boundary_mapping,
+                                                const Function<deal_II_dimension> * const a);
+template
+void MatrixCreator::create_laplace_matrix (const DoFHandler<deal_II_dimension>    &dof,
+                                          const Quadrature<deal_II_dimension>    &q,
+                                          SparseMatrix<double>     &matrix,
+                                          const Function<deal_II_dimension> * const coefficient);
+
+template
+void MatrixCreator::create_laplace_matrix (const Mapping<deal_II_dimension>       &mapping,
+                                          const DoFHandler<deal_II_dimension>    &dof,
+                                          const Quadrature<deal_II_dimension>    &q,
+                                          SparseMatrix<double>     &matrix,
+                                          const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_laplace_matrix (const Mapping<deal_II_dimension>       &mapping,
+                                          const DoFHandler<deal_II_dimension>    &dof,
+                                          const Quadrature<deal_II_dimension>    &q,
+                                          SparseMatrix<double>     &matrix,
+                                          const Function<deal_II_dimension>      &rhs,
+                                          Vector<double>           &rhs_vector,
+                                          const Function<deal_II_dimension> * const coefficient);
+template
+void MatrixCreator::create_laplace_matrix (const DoFHandler<deal_II_dimension>    &dof,
+                                          const Quadrature<deal_II_dimension>    &q,
+                                          SparseMatrix<double>     &matrix,
+                                          const Function<deal_II_dimension>      &rhs,
+                                          Vector<double>           &rhs_vector,
+                                          const Function<deal_II_dimension> * const coefficient);
 
 
 
index ca905bcf16950787e7e0640c0a03daadd47eee38..6893354929e675f243df420497a4a4706a2afaed 100644 (file)
@@ -400,16 +400,16 @@ void VectorTools::project (const Mapping<dim>       &mapping,
   SparseMatrix<double> mass_matrix (sparsity);
   Vector<double> tmp (mass_matrix.n());
 
-  MatrixCreator<dim>::create_mass_matrix (mapping, dof, quadrature, mass_matrix);
+  MatrixCreator::create_mass_matrix (mapping, dof, quadrature, mass_matrix);
   
   VectorTools::create_right_hand_side (mapping, dof, quadrature, function, tmp);
 
   constraints.condense (mass_matrix);
   constraints.condense (tmp);
   if (boundary_values.size() != 0)
-    MatrixTools<dim>::apply_boundary_values (boundary_values,
-                                            mass_matrix, vec, tmp,
-                                            true);
+    MatrixTools::apply_boundary_values (boundary_values,
+                                       mass_matrix, vec, tmp,
+                                       true);
 
   SolverControl           control(1000,1e-16);
   PrimitiveVectorMemory<> memory;
@@ -876,9 +876,9 @@ VectorTools::project_boundary_values (const Mapping<dim>       &mapping,
   Vector<double>       rhs(sparsity.n_rows());
 
 
-  MatrixCreator<dim>::create_boundary_mass_matrix (mapping, dof, q, 
-                                                  mass_matrix, boundary_functions,
-                                                  rhs, dof_to_boundary_mapping);
+  MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, 
+                                             mass_matrix, boundary_functions,
+                                             rhs, dof_to_boundary_mapping);
 
                                   // same thing as above: if dim>=3 we need
                                   // to consider constraints
index 0ca7cf63deb8c832773fffbdc2d131bcbc62d312..e195ea91b1a881671c65ed429222fc3af7e14e63 100644 (file)
@@ -531,6 +531,16 @@ documentation, etc</a>.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       Changed: The <code class="class">MatrixCreator</code> and <code
+       class="class">MatrixTools</code> class have lost their template
+       arguments. Rather, the individual functions now are templated
+       on their own, so the compiler can pick out the correct space
+       dimension on its own.
+       <br>
+       (WB 2001/08/15)
+       </p>
+
   <li> <p>
        Extended: <code class="member">ConstraintMatrix::merge</code>
        can now handle arguments which further constrain the present object.
index f1f7379396be8a3eac1479d0f0de3ea0df1b99ae..3bbc06d8fe8bf747db84f995e5a0add41c8bfb74 100644 (file)
@@ -669,10 +669,10 @@ void LaplaceProblem::assemble_system ()
                                   // system of equations
                                   // accordingly. This is done by the
                                   // following function call:
-  MatrixTools<2>::apply_boundary_values (boundary_values,
-                                        system_matrix,
-                                        solution,
-                                        system_rhs);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
 };
 
 
index 1998d4dfc2dba6e01a444201997dedcdecf3a5c5..7ab1d213f1abe0326bb38f562680bf2acf14a761 100644 (file)
@@ -468,10 +468,10 @@ void LaplaceProblem<dim>::assemble_system ()
                                            0,
                                            BoundaryValues<dim>(),
                                            boundary_values);
-  MatrixTools<dim>::apply_boundary_values (boundary_values,
-                                          system_matrix,
-                                          solution,
-                                          system_rhs);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
 };
 
 
index dbe1414ac98b91d187bad5f38b90e142f5706476..8397c531a4286e69e864f15ebeefe3eb838f79c7 100644 (file)
@@ -571,10 +571,10 @@ void LaplaceProblem<dim>::assemble_system ()
                                            0,
                                            ZeroFunction<dim>(),
                                            boundary_values);
-  MatrixTools<dim>::apply_boundary_values (boundary_values,
-                                          system_matrix,
-                                          solution,
-                                          system_rhs);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
 };
 
 
index 71bbe914e6d2c035214b841dddb4d6732686dd10..e44207b5b0892906bad5ae60433dd9e7b3fc279f 100644 (file)
@@ -578,10 +578,10 @@ void LaplaceProblem<dim>::assemble_system ()
                                            0,
                                            ZeroFunction<dim>(),
                                            boundary_values);
-  MatrixTools<dim>::apply_boundary_values (boundary_values,
-                                          system_matrix,
-                                          solution,
-                                          system_rhs);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
 };
 
 
index 8059b4162edda35179505f56e82d92ea89b90e23..0f1ea668744be46329798038dcd525b5ceb8e413 100644 (file)
@@ -1020,10 +1020,10 @@ void LaplaceProblem<dim>::assemble_system ()
                                            0,
                                            Solution<dim>(),
                                            boundary_values);
-  MatrixTools<dim>::apply_boundary_values (boundary_values,
-                                          system_matrix,
-                                          solution,
-                                          system_rhs);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
 };
 
 
index 9605840eab0a7903b33f851379963d1159244e03..d2980f13117953c329acd3948edd2787544d3922 100644 (file)
@@ -732,10 +732,10 @@ void ElasticProblem<dim>::assemble_system ()
                                            0,
                                            ZeroFunction<dim>(dim),
                                            boundary_values);
-  MatrixTools<dim>::apply_boundary_values (boundary_values,
-                                          system_matrix,
-                                          solution,
-                                          system_rhs);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
 };
 
 
index 6334b73511df0d49605b2bffe6114fdef4d9a62b..7fea52e91d244de82f5b7cd1115bac29e7df1a98 100644 (file)
@@ -309,10 +309,10 @@ void LaplaceProblem<Vector,Matrix,Sparsity>::assemble_system ()
                                            0,
                                            ZeroFunction<2>(),
                                            boundary_values);
-  MatrixTools<2>::apply_boundary_values (boundary_values,
-                                        system_matrix,
-                                        solution,
-                                        system_rhs);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
 };
 
 
index e5fdd1a6eb5ac73db4c0f2da5636054cf080b7a5..0946cbabb3ba23edae90ae3f418701bf3197442a 100644 (file)
@@ -79,7 +79,7 @@ solve_eliminated (std::map<unsigned int,double> &bv,
                  Vector<double>                &u,
                  Vector<double>                &f)
 {
-  MatrixTools<dim>::apply_boundary_values (bv, A, u, f);
+  MatrixTools::apply_boundary_values (bv, A, u, f);
   
   SolverControl control (1000, 1.e-10, false, false);
   PrimitiveVectorMemory<Vector<double> > mem;
index eb794d44b10e3651c8c1b01293024f1d4f12650f..9816d8eb6aa6e177649bddde153a25180c0750a3 100644 (file)
@@ -82,7 +82,7 @@ check_boundary (const DoFHandler<dim> &dof,
   matrix.reinit (sparsity);
   
   Vector<double> rhs (dof.n_boundary_dofs(function_map));
-  MatrixTools<dim>::
+  MatrixTools::
     create_boundary_mass_matrix (mapping, dof,
                                 face_quadrature, matrix,
                                 function_map, rhs,
@@ -168,12 +168,12 @@ check ()
       switch (test)
        {
          case 0:
-               MatrixTools<dim>::
+               MatrixTools::
                  create_mass_matrix (mapping, dof,
                                      quadrature, matrix, &coefficient);
                break;
          case 1:
-               MatrixTools<dim>::
+               MatrixTools::
                  create_laplace_matrix (mapping, dof,
                                         quadrature, matrix, &coefficient);
                break;
index 2d5d2b75264ffa7815ba84ccae23ecbf176bd54b..7110281a81ee6e8192ac2543fcc77277591eae25 100644 (file)
@@ -5321,9 +5321,9 @@ void TimeStep_Dual<dim>::do_timestep ()
       
       VectorTools::interpolate_boundary_values (*dof_handler, 0, boundary_values,
                                                     boundary_value_list);
-      MatrixTools<dim>::apply_boundary_values (boundary_value_list,
-                                              system_matrix, v,
-                                              right_hand_side1);
+      MatrixTools::apply_boundary_values (boundary_value_list,
+                                         system_matrix, v,
+                                         right_hand_side1);
     };
   
   const unsigned int solver_steps1 = solve (system_matrix, v, right_hand_side1);
@@ -5338,9 +5338,9 @@ void TimeStep_Dual<dim>::do_timestep ()
     };
   constraints.condense (right_hand_side2);
   if (dim != 1)
-    MatrixTools<dim>::apply_boundary_values (boundary_value_list,
-                                            system_matrix, u,
-                                            right_hand_side2);
+    MatrixTools::apply_boundary_values (boundary_value_list,
+                                       system_matrix, u,
+                                       right_hand_side2);
   
   if (parameters.extrapolate_old_solutions)
     {
@@ -7095,9 +7095,9 @@ assemble_vectors (right_hand_side1, right_hand_side2);
       VectorTools::interpolate_boundary_values (*dof_handler, 0,
                                                     *(parameters.boundary_values_u),
                                                     boundary_value_list);
-      MatrixTools<dim>::apply_boundary_values (boundary_value_list,
-                                              system_matrix, u,
-                                              right_hand_side1);
+      MatrixTools::apply_boundary_values (boundary_value_list,
+                                         system_matrix, u,
+                                         right_hand_side1);
     };
 
   const unsigned int solver_steps1 = solve (system_matrix, u, right_hand_side1);
@@ -7119,9 +7119,9 @@ assemble_vectors (right_hand_side1, right_hand_side2);
       VectorTools::interpolate_boundary_values (*dof_handler, 0,
                                                     *(parameters.boundary_values_v),
                                                     boundary_value_list);
-      MatrixTools<dim>::apply_boundary_values (boundary_value_list,
-                                              system_matrix, v,
-                                              right_hand_side2);
+      MatrixTools::apply_boundary_values (boundary_value_list,
+                                         system_matrix, v,
+                                         right_hand_side2);
     };
 
 

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.