]> https://gitweb.dealii.org/ - dealii.git/commitdiff
working on constraints for matrix blocks
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 21 Mar 2010 19:22:22 +0000 (19:22 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 21 Mar 2010 19:22:22 +0000 (19:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@20871 0785d39b-7218-0410-832d-ea1e28bc413d

19 files changed:
deal.II/base/include/base/mg_level_object.h [moved from deal.II/deal.II/include/multigrid/mg_level_object.h with 97% similarity]
deal.II/deal.II/include/multigrid/mg_block_smoother.h
deal.II/deal.II/include/multigrid/mg_matrix.h
deal.II/deal.II/include/multigrid/mg_smoother.h
deal.II/deal.II/include/multigrid/mg_transfer.h
deal.II/deal.II/include/multigrid/mg_transfer_block.h
deal.II/deal.II/include/multigrid/mg_transfer_component.h
deal.II/deal.II/include/multigrid/multigrid.h
deal.II/deal.II/include/numerics/mesh_worker.h
deal.II/deal.II/include/numerics/mesh_worker_assembler.h
deal.II/deal.II/include/numerics/mesh_worker_output.h
deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h
deal.II/deal.II/source/multigrid/mg_base.cc
deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc
deal.II/deal.II/source/multigrid/mg_tools.cc
deal.II/examples/step-12/step-12.cc
deal.II/lac/include/lac/constraint_matrix.h
deal.II/lac/include/lac/constraint_matrix.templates.h
deal.II/lac/include/lac/matrix_block.h

similarity index 97%
rename from deal.II/deal.II/include/multigrid/mg_level_object.h
rename to deal.II/base/include/base/mg_level_object.h
index 1b97b476a81b0856e74e02fcaa68412c20da8f9b..7b75ed41547023499c94f4127a6c5cc324415d01 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 2003, 2004, 2005, 2006, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -21,9 +21,6 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-/*!@addtogroup mg */
-/*@{*/
-
 /**
  * An array with an object for each level.  The purpose of this class
  * is mostly to store objects and allow access by level number, even
@@ -34,7 +31,9 @@ DEAL_II_NAMESPACE_OPEN
  * In most cases, the objects which are stored on each levels, are
  * either matrices or vectors.
  *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2005
+ * @ingroup mg
+ * @ingroup data
+ * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2005, 2010
  */
 template<class Object>
 class MGLevelObject : public Subscriptor
@@ -116,7 +115,6 @@ class MGLevelObject : public Subscriptor
     std::vector<std_cxx1x::shared_ptr<Object> > objects;
 };
 
-/*@}*/
 
 /* ------------------------------------------------------------------- */
 
index 3065fc7e90812218afdf5f69b11e425b1e274ec3..b51afecf49537bac74fa78494843f4ad8ed3a5c2 100644 (file)
@@ -20,7 +20,7 @@
 #include <lac/vector_memory.h>
 #include <lac/block_vector.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <vector>
 
 DEAL_II_NAMESPACE_OPEN
index d66cf2c3de47707ebfc9b23e577c8fe919bd3324..274b1eb3eac5e3282f2df23215f10a0c3f1a79f7 100644 (file)
@@ -16,7 +16,7 @@
 #include <lac/vector.h>
 #include <lac/sparse_matrix.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 DEAL_II_NAMESPACE_OPEN
 
index d298470a95f83aae5a6f7ad35bb0b6876079c6d3..e28f7cd486c706b1b251f212c4f619eab14e292b 100644 (file)
@@ -19,7 +19,7 @@
 #include <lac/pointer_matrix.h>
 #include <lac/vector_memory.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <vector>
 
 DEAL_II_NAMESPACE_OPEN
index 8c63dc237a8c3336f8719107f368b673bd9fd07e..7db0779cc55ce5e3ce02241fdafdf74f29d5138a 100644 (file)
@@ -27,7 +27,7 @@
 #include <lac/vector_memory.h>
 
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 
index bc895bc3ecfcd8b0631528f0ead850688955f3ec..a9f7451ef21396692205bcd67fc6f3589b497f1a 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
+//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -27,7 +27,7 @@
 #include <lac/block_matrix_array.h>
 
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 
index 2d1e1de35293c47146eab72403abf14257d4bc9f..3a209635f4d310982a3b87e66806b5be52374606 100644 (file)
@@ -27,7 +27,7 @@
 #include <lac/vector_memory.h>
 
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 
index 8442d467ef6b79a53b4b62891c2f0a7fa51ca5c4..fbcccdce6bcfd69dce042c9d13079665813590e4 100644 (file)
@@ -21,7 +21,7 @@
 #include <lac/sparse_matrix.h>
 #include <lac/vector.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <multigrid/mg_dof_handler.h>
 
 #include <vector>
index 3602b4331f6199be1178bbae8b6b3b42eb62db7a..6f9eda3ae4ab007cbd9c1f081118f38a7206004e 100644 (file)
@@ -18,7 +18,7 @@
 #include <base/geometry_info.h>
 #include <base/named_data.h>
 #include <lac/block_indices.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <numerics/mesh_worker_vector_selector.h>
 #include <numerics/mesh_worker_info.h>
 #include <numerics/mesh_worker_assembler.h>
index 01371224c78a9c1a3158e7a3ff208de63b4ad009..e1af235b0b2c2d86406dbaefc60986608a1f330f 100644 (file)
@@ -17,7 +17,7 @@
 #include <base/named_data.h>
 #include <base/smartpointer.h>
 #include <lac/block_vector.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 DEAL_II_NAMESPACE_OPEN
index 0e557c67afa80760cef46f4b5b7445cd61fc6970..1ae30f1ec81c3063a9958dd3705792d88151d21b 100644 (file)
@@ -17,7 +17,7 @@
 #include <base/named_data.h>
 #include <base/smartpointer.h>
 #include <lac/block_vector.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 DEAL_II_NAMESPACE_OPEN
index 80312e8f881b3932eddea8dfe7642216d32cb558..41f709fd65c0103a6625fd814f70c283fa72cc98 100644 (file)
@@ -16,7 +16,7 @@
 #include <base/named_data.h>
 #include <base/tensor.h>
 #include <base/smartpointer.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 DEAL_II_NAMESPACE_OPEN
 
index 90476e759b2a0630c26fcc916c8411dcf805eba7..e7088cdfd79eafe9f94ac1ad1bd71cb98c0f983a 100644 (file)
@@ -17,7 +17,7 @@
 #include <fe/fe.h>
 #include <multigrid/mg_dof_handler.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <multigrid/mg_tools.h>
 
 
index e85e42f8fde3cd53039de03d1b769d5a1b010396..8a7fbba6f1f4b1497261b86edc4d568d18a650f4 100644 (file)
@@ -18,7 +18,7 @@
 #include <lac/block_vector.h>
 #include <multigrid/mg_tools.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 #include <vector>
 #include <algorithm>
index 92eaf98f415ed182b1d7a4c5723d026f84105cf6..035e80a42c957e7f19fa609dd9923649066a9d3c 100644 (file)
@@ -26,7 +26,7 @@
 #include <multigrid/mg_dof_accessor.h>
 #include <multigrid/mg_tools.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe.h>
 
index 8eff89bddf586ce3be512f2ebccee5d2a13fef5c..64322d22856e472f27bf3cf8885156e15866d880 100644 (file)
@@ -456,14 +456,14 @@ void Step12<dim>::integrate_cell_term (DoFInfo& dinfo, CellInfo& info)
 }
 
                                 // Now the same for the boundary terms. Note
-                                // that now we use FEFaceValuesBase, the base
+                                // that now we use FEValuesBase, the base
                                 // class for both FEFaceValues and
                                 // FESubfaceValues, in order to get access to
                                 // normal vectors.
 template <int dim>
 void Step12<dim>::integrate_boundary_term (DoFInfo& dinfo, FaceInfo& info)
 {
-  const FEFaceValuesBase<dim>& fe_v = info.fe_values();
+  const FEValuesBase<dim>& fe_v = info.fe_values();
   FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;
   Vector<double>& local_vector = dinfo.vector(0).block(0);
 
@@ -512,14 +512,14 @@ void Step12<dim>::integrate_face_term (DoFInfo& dinfo1, DoFInfo& dinfo2,
 {
                                   // For quadrature points, weights,
                                   // etc., we use the
-                                  // FEFaceValuesBase object of the
+                                  // FEValuesBase object of the
                                   // first argument.
-  const FEFaceValuesBase<dim>& fe_v = info1.fe_values();
+  const FEValuesBase<dim>& fe_v = info1.fe_values();
 
                                   // For additional shape functions,
                                   // we have to ask the neighbors
-                                  // FEFaceValuesBase.
-  const FEFaceValuesBase<dim>& fe_v_neighbor = info2.fe_values();
+                                  // FEValuesBase.
+  const FEValuesBase<dim>& fe_v_neighbor = info2.fe_values();
 
                                   // Then we get references to the
                                   // four local matrices. The letters
index c4c21c47276671ab40c79d402f786b858a8dc03f..79ffcfe2b3f2a7af88156501a410e0b7e67be4cf 100644 (file)
@@ -1825,6 +1825,18 @@ class ConstraintMatrix : public Subscriptor
                                 const Table<2,bool>             &dof_mask,
                                 internal::bool2type<true>) const;
 
+                                    /**
+                                     * Internal helper function for
+                                     * distribute_local_to_global
+                                     * function.
+                                     *
+                                     * Creates a list of affected
+                                     * rows for sistribution.
+                                     */
+    void
+    make_sorted_row_list (const std::vector<unsigned int> &local_dof_indices,
+                         internals::GlobalRowsFromLocal  &global_rows) const;
+
                                     /**
                                      * Internal helper function for
                                      * distribute_local_to_global function.
index 5f20f9465f7e5c7e8201452211df53a26197ba4b..d9a2516ad3ef7233b394184fa76e9ef41494be1b 100644 (file)
@@ -1722,6 +1722,85 @@ namespace internals
 } // end of namespace internals
 
 
+void
+ConstraintMatrix::
+make_sorted_row_list (const std::vector<unsigned int> &local_dof_indices,
+                     internals::GlobalRowsFromLocal  &global_rows) const
+{
+  const unsigned int n_local_dofs = local_dof_indices.size();
+                                  // when distributing the local data to
+                                  // the global matrix, we can quite
+                                  // cheaply sort the indices (obviously,
+                                  // this introduces the need for
+                                  // allocating some memory on the way, but
+                                  // we need to do this only for rows,
+                                  // whereas the distribution process
+                                  // itself goes over rows and
+                                  // columns). This has the advantage that
+                                  // when writing into the global matrix,
+                                  // we can make use of the sortedness.
+
+                                  // so the first step is to create a
+                                  // sorted list of all row values that are
+                                  // possible. these values are either the
+                                  // rows from unconstrained dofs, or some
+                                  // indices introduced by dofs constrained
+                                  // to a combination of some other
+                                  // dofs. regarding the data type, choose
+                                  // an STL vector of a pair of unsigned
+                                  // ints (for global columns) and internal
+                                  // data (containing local columns +
+                                  // possible jumps from
+                                  // constraints). Choosing an STL map or
+                                  // anything else M.K. knows of would be
+                                  // much more expensive here!
+
+                                  // cache whether we have to resolve any
+                                  // indirect rows generated from resolving
+                                  // constrained dofs.
+  unsigned int added_rows = 0;
+
+                                  // first add the indices in an unsorted
+                                  // way and only keep track of the
+                                  // constraints that appear. They are
+                                  // resolved in a second step.
+  for (unsigned int i = 0; i<n_local_dofs; ++i)
+    {
+      if (is_constrained(local_dof_indices[i]) == false)
+       {
+         global_rows.global_row(added_rows)  = local_dof_indices[i];
+         global_rows.local_row(added_rows++) = i;
+         continue;
+       }
+      global_rows.insert_constraint(i);
+    }
+  global_rows.sort();
+
+  const unsigned int n_constrained_rows = n_local_dofs-added_rows;
+  for (unsigned int i=0; i<n_constrained_rows; ++i)
+    {
+      const unsigned int local_row = global_rows.last_constrained_local_row();
+      Assert (local_row < n_local_dofs, ExcIndexRange(local_row, 0, n_local_dofs));
+      const unsigned int global_row = local_dof_indices[local_row];
+      Assert (is_constrained(global_row), ExcInternalError());
+      const ConstraintLine & position =
+       lines[lines_cache[calculate_line_index(global_row)]];
+      if (position.inhomogeneity != 0)
+       global_rows.set_last_row_inhomogeneous();
+      else
+       global_rows.pop_back();
+      for (unsigned int q=0; q<position.entries.size(); ++q)
+       global_rows.insert_index (position.entries[q].first,
+                                 local_row,
+                                 position.entries[q].second);
+    }
+}
+
+
+
+
+//TODO: This function is DANGEROUS, because it does more than it claims!
+
 
                                // Basic idea of setting up a list of
                                // all global dofs: first find all rows and columns
index 7d251af375b8864fa74da5e63be8193505bd0ed4..668fa620cd42f3eef339fa31d1251aaf8d864cdf 100644 (file)
@@ -27,30 +27,52 @@ DEAL_II_NAMESPACE_OPEN
  * A wrapper around a matrix object, storing the coordinates in a
  * block matrix as well.
  *
- * @note They are the position in the global matrix, where we don't use
- * BlockSparseMatrix. Matrices may have the same coordinates for the
- * following reason: A preconditioner for the Oseen system can be
- * built as a block system, where the pressure block is of the form
- * M^-1 F A^-1 with m the mass matrix (u,v), A the Laplacian and F the
- * advection diffusion operator. So, I need to build three matrices
- * for a single block in my system + in some cases a pressure
- * stabilization.
+ * This class is an alternative to BlockMatrixBase, if you only want
+ * to generate a single block of the system, not the whole
+ * system. Using the add() functions of this class, it is possible to
+ * use the standard assembling functions used for block matrices, but
+ * only enter in one of the blocks and still avoiding the index
+ * computations innvolved.
+
+ * The reason for this class is, that we may need a different number
+ * of matrices for different blocks in a block system. For example, a
+ * preconditioner for the Oseen system can be built as a block system,
+ * where the pressure block is of the form
+ * <b>M</b><sup>-1</sup><b>FA</b><sup>-1</sup> with <b>M</b> the
+ * pressure mass matrix, <b>A</b> the pressure Laplacian and <b>F</b>
+ * the advection diffusion operator applied to the pressure
+ * space. Since only a single matrix is needed for the other blocks,
+ * using BlockSparseMatrix or similar would be a waste of memory.
+ *
+ * While the add() functions make a MatrixBlock appear like a block
+ * matrix for assembling, the functions vmult(), Tvmult(),
+ * vmult_add(), and Tvmult_add() make it behave like a MATRIX, when it
+ * comes to applying it to a vector. This behavior allows us to store
+ * MatrixBlock objects in vectors, for instance in MGLevelObject
+ * without extracting the #matrix first. 
  *
  * MatrixBlock comes handy when using BlockMatrixArray. Once the
  * MatrixBlock has been properly initalized and filled, it can be used
- * as:
+ * in the simplest case as:
  * @code
- * std::vector<MatrixBlock<SparseMatrix<double> > > blocks;
+ * MatrixBlockVector<SparseMatrix<double> > > blocks;
  *
  * ...
  *
  * BlockMatrixArray matrix (n_blocks, n_blocks);
  *
  * for (unsigned int i=0;i<blocks.size;++i)
- *   matrix.enter(blocks[i].row, blocks[i].column, blocks[i].matrix);
+ *   matrix.enter(blocks.block(i).row, blocks.block(i).column, blocks.matrix(i));
  * @endcode
  *
+ * Here, we have not gained very much, except that we do not need to
+ * set up empty blocks in the block system.
+ *
+ * @todo Example for the product preconditioner of the pressure Schur
+ * complement.
  *
+ * @ingroup Matrix2
+ * @ingroup vector_valued
  * @author Guido Kanschat, 2006
  */
 template <class MATRIX>
@@ -79,33 +101,34 @@ class MatrixBlock
 
                                     /**
                                      * Add <tt>value</tt> to the
-                                     * element (<i>i,j</i>).  Throws
+                                     * element (<i>i,j</i>). Throws
                                      * an error if the entry does not
-                                     * exist or if <tt>value</tt> is
-                                     * not a finite number. Still, it
-                                     * is allowed to store zero
-                                     * values in non-existent fields.
+                                     * exist or if it is in a
+                                     * different block.
                                      */
     void add (const unsigned int i,
               const unsigned int j,
              const typename MATRIX::value_type value);
 
                                        /**
-                                        * Add all elements given in a
-                                        * FullMatrix<double> into sparse
+                                        * Add all elements in a
+                                        * FullMatrix into sparse
                                         * matrix locations given by
-                                        * <tt>indices</tt>. In other words,
-                                        * this function adds the elements in
-                                        * <tt>full_matrix</tt> to the
-                                        * respective entries in calling
-                                        * matrix, using the local-to-global
-                                        * indexing specified by
-                                        * <tt>indices</tt> for both the rows
-                                        * and the columns of the
-                                        * matrix. This function assumes a
-                                        * quadratic sparse matrix and a
-                                        * quadratic full_matrix, the usual
-                                        * situation in FE calculations.
+                                        * <tt>indices</tt>. This function
+                                        * assumes a quadratic sparse
+                                        * matrix and a quadratic
+                                        * full_matrix.  The global
+                                        * locations are translated
+                                        * into locations in this block
+                                        * and ExcBlockIndexMismatch is
+                                        * thrown, if the global index
+                                        * does not point into the
+                                        * block refered to by #row and
+                                        * #column.
+                                       *
+                                       * @todo
+                                       * <tt>elide_zero_values<tt> is
+                                       * currently ignored.
                                        *
                                        * The optional parameter
                                        * <tt>elide_zero_values</tt> can be
@@ -123,11 +146,33 @@ class MatrixBlock
              const bool                       elide_zero_values = true);
 
                                        /**
-                                        * Same function as before, but now
-                                        * including the possibility to use
-                                        * rectangular full_matrices and
-                                        * different local-to-global indexing
-                                        * on rows and columns, respectively.
+                                       * Add all elements in a
+                                        * FullMatrix into global
+                                        * locations given by
+                                        * <tt>row_indices</tt> and
+                                        * <tt>col_indices</tt>,
+                                        * respectively. The global
+                                        * locations are translated
+                                        * into locations in this block
+                                        * and ExcBlockIndexMismatch is
+                                        * thrown, if the global index
+                                        * does not point into the
+                                        * block refered to by #row and
+                                        * #column.
+                                       *
+                                       * @todo
+                                       * <tt>elide_zero_values<tt> is
+                                       * currently ignored.
+                                       *
+                                       * The optional parameter
+                                       * <tt>elide_zero_values</tt> can be
+                                       * used to specify whether zero
+                                       * values should be added anyway or
+                                       * these should be filtered away and
+                                       * only non-zero data is added. The
+                                       * default value is <tt>true</tt>,
+                                       * i.e., zero values won't be added
+                                       * into the matrix.
                                        */
     template <typename number>
     void add (const std::vector<unsigned int>& row_indices,
@@ -137,10 +182,27 @@ class MatrixBlock
 
                                        /**
                                         * Set several elements in the
-                                        * specified row of the matrix with
-                                        * column indices as given by
-                                        * <tt>col_indices</tt> to the
-                                        * respective value.
+                                        * specified row of the matrix
+                                        * with column indices as given
+                                        * by <tt>col_indices</tt> to
+                                        * the respective value. This
+                                        * is the function doing thye
+                                        * actual work for the ones
+                                        * adding full matrices. The
+                                        * global locations
+                                        * <tt>row_index</tt> and
+                                        * <tt>col_indices</tt> are
+                                        * translated into locations in
+                                        * this block and
+                                        * ExcBlockIndexMismatch is
+                                        * thrown, if the global index
+                                        * does not point into the
+                                        * block refered to by #row and
+                                        * #column.
+                                       *
+                                       * @todo
+                                       * <tt>elide_zero_values<tt> is
+                                       * currently ignored.
                                        *
                                        * The optional parameter
                                        * <tt>elide_zero_values</tt> can be
@@ -153,7 +215,7 @@ class MatrixBlock
                                        * into the matrix.
                                        */
     template <typename number>
-    void add (const unsigned int               row,
+    void add (const unsigned int               row_index,
              const std::vector<unsigned int>& col_indices,
              const std::vector<number>&       values,
              const bool                       elide_zero_values = true);
@@ -183,6 +245,50 @@ class MatrixBlock
              const bool          elide_zero_values = true,
              const bool          col_indices_are_sorted = false);
 
+                                    /**
+                                     * Matrix-vector-multiplication,
+                                     * forwarding to the same
+                                     * function in MATRIX. No index
+                                     * computations are done, thus,
+                                     * the vectors need to have sizes
+                                     * matching #matrix.
+                                     */
+    template<class VECTOR>
+    void vmult (VECTOR& w, const VECTOR& v) const;
+    
+                                    /**
+                                     * Matrix-vector-multiplication,
+                                     * forwarding to the same
+                                     * function in MATRIX. No index
+                                     * computations are done, thus,
+                                     * the vectors need to have sizes
+                                     * matching #matrix.
+                                     */
+    template<class VECTOR>
+    void vmult_add (VECTOR& w, const VECTOR& v) const;
+    
+                                    /**
+                                     * Matrix-vector-multiplication,
+                                     * forwarding to the same
+                                     * function in MATRIX. No index
+                                     * computations are done, thus,
+                                     * the vectors need to have sizes
+                                     * matching #matrix.
+                                     */
+    template<class VECTOR>
+    void Tvmult (VECTOR& w, const VECTOR& v) const;
+
+                                    /**
+                                     * Matrix-vector-multiplication,
+                                     * forwarding to the same
+                                     * function in MATRIX. No index
+                                     * computations are done, thus,
+                                     * the vectors need to have sizes
+                                     * matching #matrix.
+                                     */
+    template<class VECTOR>
+    void Tvmult_add (VECTOR& w, const VECTOR& v) const;
+
                                     /**
                                      * The block number computed from
                                      * an index by using
@@ -229,6 +335,8 @@ class MatrixBlock
  * pointers, in order to allow for copying and rearranging. Each
  * matrix block can be identified by name.
  *
+ * @relates MatrixBlock
+ * @ingroup vector_valued
  * @author Baerbel Janssen, Guido Kanschat, 2010
  */
 template <class MATRIX>
@@ -264,6 +372,48 @@ class MatrixBlockVector : public NamedData<boost::shared_ptr<MatrixBlock<MATRIX>
 };
 
 
+/**
+ * A vector of MGLevelObject<MatrixBlock>, which is implemented using shared
+ * pointers, in order to allow for copying and rearranging. Each
+ * matrix block can be identified by name.
+ *
+ * @relates MatrixBlock
+ * @ingroup vector_valued
+ * @author Baerbel Janssen, Guido Kanschat, 2010
+ */
+template <class MATRIX>
+class MGMatrixBlockVector : public NamedData<boost::shared_ptr<MatrixBlock<MATRIX> > >
+{
+  public:
+                                    /**
+                                     * Add a new matrix block at the
+                                     * position <tt(row,column)</tt>
+                                     * in the block system.
+                                     */
+    void add(unsigned int row, unsigned int column,
+            const std::string& name,
+            const BlockIndices* block_indices);
+                                    /**
+                                     * Access a constant reference to
+                                     * the block at position <i>i</i>.
+                                     */
+    const MatrixBlock<MATRIX>& block(unsigned int i) const;
+                                    /**
+                                     * Access a reference to
+                                     * the block at position <i>i</i>.
+                                     */
+    MatrixBlock<MATRIX>& block(unsigned int i);
+                                    /**
+                                     * Access the matrix at position
+                                     * <i>i</i> for read and write
+                                     * access.
+                                     */
+    MATRIX& matrix(unsigned int i);
+    
+    
+};
+
+
 //----------------------------------------------------------------------//
 
 template <class MATRIX>
@@ -336,6 +486,43 @@ MatrixBlock<MATRIX>::add (const std::vector<unsigned int>&         row_indices,
 }
 
 
+template <class MATRIX>
+template <typename number>
+inline
+void
+MatrixBlock<MATRIX>::add (const unsigned int   b_row,
+                         const unsigned int   n_cols,
+                         const unsigned int  *col_indices,
+                         const number        *values,
+                         const bool,
+                         const bool)
+{
+  Assert(block_indices != 0, ExcNotInitialized());
+  const std::pair<unsigned int, unsigned int> bi
+    = block_indices->global_to_local(b_row);
+
+                                  // In debug mode, we check whether
+                                  // all indices are in the correct
+                                  // block.
+
+                                  // Actually, for the time being, we
+                                  // leave it at this. While it may
+                                  // not be the most efficient way,
+                                  // it is at least thread safe.
+//#ifdef DEBUG
+  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
+
+  for (unsigned int j=0;j<n_cols;++j)
+    {
+      const std::pair<unsigned int, unsigned int> bj
+       = block_indices->global_to_local(col_indices[j]);
+      Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
+      
+      matrix.add(bi.second, bj.second, values[j]);
+    }
+//#endif
+}
+
 
 template <class MATRIX>
 template <typename number>
@@ -374,45 +561,46 @@ MatrixBlock<MATRIX>::add (const unsigned int               row,
 
 
 template <class MATRIX>
-template <typename number>
+template <class VECTOR>
 inline
 void
-MatrixBlock<MATRIX>::add (const unsigned int   b_row,
-                         const unsigned int   n_cols,
-                         const unsigned int  *col_indices,
-                         const number        *values,
-                         const bool,
-                         const bool)
+MatrixBlock<MATRIX>::vmult (VECTOR& w, const VECTOR& v) const
 {
-  Assert(block_indices != 0, ExcNotInitialized());
-  const std::pair<unsigned int, unsigned int> bi
-    = block_indices->global_to_local(b_row);
+  matrix.vmult(w,v);
+}
 
-                                  // In debug mode, we check whether
-                                  // all indices are in the correct
-                                  // block.
 
-                                  // Actually, for the time being, we
-                                  // leave it at this. While it may
-                                  // not be the most efficient way,
-                                  // it is at least thread safe.
-//#ifdef DEBUG
-  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
+template <class MATRIX>
+template <class VECTOR>
+inline
+void
+MatrixBlock<MATRIX>::vmult_add (VECTOR& w, const VECTOR& v) const
+{
+  matrix.vmult_add(w,v);
+}
 
-  for (unsigned int j=0;j<n_cols;++j)
-    {
-      const std::pair<unsigned int, unsigned int> bj
-       = block_indices->global_to_local(col_indices[j]);
-      Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
-      
-      matrix.add(bi.second, bj.second, values[j]);
-    }
-  
-//#endif
 
+template <class MATRIX>
+template <class VECTOR>
+inline
+void
+MatrixBlock<MATRIX>::Tvmult (VECTOR& w, const VECTOR& v) const
+{
+  matrix.Tvmult(w,v);
+}
+
+
+template <class MATRIX>
+template <class VECTOR>
+inline
+void
+MatrixBlock<MATRIX>::Tvmult_add (VECTOR& w, const VECTOR& v) const
+{
+  matrix.Tvmult_add(w,v);
 }
 
 
+
 //----------------------------------------------------------------------//
 
 template <class MATRIX>

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.