]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement CompressedBlockSparsityPattern, use it, test it.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 23 Jun 2001 20:14:03 +0000 (20:14 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 23 Jun 2001 20:14:03 +0000 (20:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@4751 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/include/dofs/dof_tools.h
deal.II/deal.II/source/dofs/dof_constraints.cc
deal.II/deal.II/source/dofs/dof_tools.cc
deal.II/lac/include/lac/block_sparsity_pattern.h
deal.II/lac/source/block_sparsity_pattern.cc
tests/deal.II/sparsity_pattern.cc
tests/deal.II/sparsity_pattern.checked

index 9cca92fac6993fd2e231d4ac51c27f4318ab99b8..6b78a0466070653afeb1a19ec03a6570e7631070 100644 (file)
 class SparsityPattern;
 class CompressedSparsityPattern;
 class BlockSparsityPattern;
+class CompressedBlockSparsityPattern;
 template <typename number> class SparseMatrix;
 template <typename number> class BlockSparseMatrix;
 class BlockIndices;
 
+
 /**
  * This class represents the matrix denoting the distribution of the degrees
  * of freedom of hanging nodes.
@@ -290,6 +292,13 @@ class ConstraintMatrix : public Subscriptor
                                      * sparsity patterns.
                                      */
     void condense (CompressedSparsityPattern &sparsity) const;
+
+                                    /**
+                                     * Same function as above, but
+                                     * condenses square compressed
+                                     * sparsity patterns.
+                                     */
+    void condense (CompressedBlockSparsityPattern &sparsity) const;
     
                                     /**
                                      * Condense a given matrix. The associated
index aca12523be09938974344e5fa9d53b4493a0bee9..8a315a0aa9071d35fbb7eb7cf7be5639ab671688 100644 (file)
@@ -193,20 +193,22 @@ class DoFTools
                                      * sparsity pattern may be
                                      * @ref{SparsityPattern},
                                      * @ref{CompressedSparsityPattern},
-                                     * @ref{BlockSparsityPattern}, or
-                                     * any other class that satisfies
-                                     * similar requirements. It is
-                                     * assumed that the size of the
-                                     * sparsity pattern matches the
-                                     * number of degrees of freedom
-                                     * and that enough unused nonzero
-                                     * entries are left to fill the
-                                     * sparsity pattern. The nonzero
-                                     * entries generated by this
-                                     * function are overlaid to
-                                     * possible previous content of
-                                     * the object, that is previously
-                                     * added entries are not deleted.
+                                     * @ref{BlockSparsityPattern},
+                                     * @ref{CompressedBlockSparsityPattern},
+                                     * or any other class that
+                                     * satisfies similar
+                                     * requirements. It is assumed
+                                     * that the size of the sparsity
+                                     * pattern matches the number of
+                                     * degrees of freedom and that
+                                     * enough unused nonzero entries
+                                     * are left to fill the sparsity
+                                     * pattern. The nonzero entries
+                                     * generated by this function are
+                                     * overlaid to possible previous
+                                     * content of the object, that is
+                                     * previously added entries are
+                                     * not deleted.
                                      */
     template <int dim, class SparsityPattern>
     static
@@ -274,9 +276,11 @@ class DoFTools
                                      * sparsity pattern may be
                                      * @ref{SparsityPattern},
                                      * @ref{CompressedSparsityPattern},
-                                     * @ref{BlockSparsityPattern}, or
-                                     * any other class that satisfies
-                                     * similar requirements.
+                                     * @ref{BlockSparsityPattern},
+                                     * @ref{CompressedBlockSparsityPattern},
+                                     * or any other class that
+                                     * satisfies similar
+                                     * requirements.
                                      */
     template <int dim, class SparsityPattern>
     static
@@ -310,20 +314,22 @@ class DoFTools
                                      * sparsity pattern may be
                                      * @ref{SparsityPattern},
                                      * @ref{CompressedSparsityPattern},
-                                     * @ref{BlockSparsityPattern}, or
-                                     * any other class that satisfies
-                                     * similar requirements. It is
-                                     * assumed that the size of the
-                                     * sparsity pattern matches the
-                                     * number of degrees of freedom
-                                     * and that enough unused nonzero
-                                     * entries are left to fill the
-                                     * sparsity pattern. The nonzero
-                                     * entries generated by this
-                                     * function are overlaid to
-                                     * possible previous content of
-                                     * the object, that is previously
-                                     * added entries are not deleted.
+                                     * @ref{BlockSparsityPattern},
+                                     * @ref{CompressedBlockSparsityPattern},
+                                     * or any other class that
+                                     * satisfies similar
+                                     * requirements. It is assumed
+                                     * that the size of the sparsity
+                                     * pattern matches the number of
+                                     * degrees of freedom and that
+                                     * enough unused nonzero entries
+                                     * are left to fill the sparsity
+                                     * pattern. The nonzero entries
+                                     * generated by this function are
+                                     * overlaid to possible previous
+                                     * content of the object, that is
+                                     * previously added entries are
+                                     * not deleted.
                                      */
     template <int dim, class SparsityPattern>
     static void
index e3cc3f26bd8afd7b5b7b6e9ba6aab17ec44db2e1..96884c75f5885e2a07d5158ae640defffc425e90 100644 (file)
@@ -792,6 +792,138 @@ void ConstraintMatrix::condense (BlockSparsityPattern &sparsity) const
 
 
 
+void ConstraintMatrix::condense (CompressedBlockSparsityPattern &sparsity) const
+{
+  Assert (sorted == true, ExcMatrixNotClosed());
+  Assert (sparsity.n_rows() == sparsity.n_cols(),
+         ExcMatrixNotSquare());
+  Assert (sparsity.n_block_rows() == sparsity.n_block_cols(),
+         ExcMatrixNotSquare());
+  Assert (sparsity.get_column_indices() == sparsity.get_row_indices(),
+         ExcMatrixNotSquare());
+  
+  const BlockIndices &
+    index_mapping = sparsity.get_column_indices();
+
+  const unsigned int n_blocks = sparsity.n_block_rows();
+  
+                                  // store for each index whether it
+                                  // must be distributed or not. If entry
+                                  // is -1, no distribution is necessary.
+                                  // otherwise, the number states which
+                                  // line in the constraint matrix handles
+                                  // this index
+  std::vector<int> distribute (sparsity.n_rows(), -1);
+  
+  for (unsigned int c=0; c<lines.size(); ++c)
+    distribute[lines[c].line] = static_cast<signed int>(c);
+
+  const unsigned int n_rows = sparsity.n_rows();
+  for (unsigned int row=0; row<n_rows; ++row)
+    {
+                                      // get index of this row
+                                      // within the blocks
+      const std::pair<unsigned int,unsigned int>
+       block_index = index_mapping.global_to_local(row);
+      const unsigned int block_row = block_index.first;
+      const unsigned int local_row = block_index.second;
+      
+      if (distribute[row] == -1)
+                                        // regular line. loop over
+                                        // all columns and see
+                                        // whether this column must
+                                        // be distributed. note that
+                                        // as we proceed to
+                                        // distribute cols, the loop
+                                        // over cols may get longer.
+                                        //
+                                        // don't try to be clever
+                                        // here as in the algorithm
+                                        // for the
+                                        // CompressedSparsityPattern,
+                                        // as that would be much more
+                                        // complicated here. after
+                                        // all, we know that
+                                        // compressed patterns are
+                                        // inefficient...
+       {
+
+                                          // to loop over all entries
+                                          // in this row, we have to
+                                          // loop over all blocks in
+                                          // this blockrow and the
+                                          // corresponding row
+                                          // therein
+         for (unsigned int block_col=0; block_col<n_blocks; ++block_col)
+           {
+             const CompressedSparsityPattern &
+               block_sparsity = sparsity.block(block_row, block_col);
+
+             for (unsigned int j=0; j<block_sparsity.row_length(local_row); ++j)
+               {
+                 const unsigned int global_col
+                   = index_mapping.local_to_global(block_col,
+                                                   block_sparsity.column_number(local_row,j));
+                   
+                 if (distribute[global_col] != -1)
+                                                    // distribute entry at regular
+                                                    // row @p{row} and irregular column
+                                                    // global_col
+                   {
+                     for (unsigned int q=0;
+                          q!=lines[distribute[global_col]]
+                                         .entries.size(); ++q)
+                       sparsity.add (row,
+                                     lines[distribute[global_col]].entries[q].first);
+                   };
+               };
+           };
+       }
+      else
+       {
+                                          // row must be
+                                          // distributed. split the
+                                          // whole row into the
+                                          // chunks defined by the
+                                          // blocks
+         for (unsigned int block_col=0; block_col<n_blocks; ++block_col)
+           {
+             const CompressedSparsityPattern &
+               block_sparsity = sparsity.block(block_row,block_col);
+      
+             for (unsigned int j=0; j<block_sparsity.row_length(local_row); ++j)
+               {
+                 const unsigned int global_col
+                   = index_mapping.local_to_global (block_col,
+                                                    block_sparsity.column_number(local_row,j));
+                   
+                 if (distribute[global_col] == -1)
+                                                    // distribute entry at irregular
+                                                    // row @p{row} and regular column
+                                                    // global_col.
+                   {
+                     for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q) 
+                       sparsity.add (lines[distribute[row]].entries[q].first,
+                                     global_col);
+                   }
+                 else
+                                                    // distribute entry at irregular
+                                                    // row @p{row} and irregular column
+                                                    // @p{global_col}
+                   {
+                     for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
+                       for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q)
+                         sparsity.add (lines[distribute[row]].entries[p].first,
+                                       lines[distribute[global_col]].entries[q].first);
+                   };
+               };
+           };
+       };
+    };
+};
+
+
+
 unsigned int ConstraintMatrix::n_constraints () const
 {
   return lines.size();
index 4ace8c073d5947988661895d0ac86932ed537fb5..a9c93175e6e48e93072946b95a496c96a60db721 100644 (file)
@@ -2182,6 +2182,12 @@ template void
 DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
                                          const std::vector<unsigned int>  &,
                                          BlockSparsityPattern    &);
+template void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+                                         const std::vector<unsigned int>  &,
+                                         CompressedBlockSparsityPattern    &);
+
+
 template void
 DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
                                          const FunctionMap<deal_II_dimension>::type  &boundary_indicators,
@@ -2197,6 +2203,13 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& d
                                          const FunctionMap<deal_II_dimension>::type  &boundary_indicators,
                                          const std::vector<unsigned int>  &dof_to_boundary_mapping,
                                          BlockSparsityPattern    &sparsity);
+template void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+                                         const FunctionMap<deal_II_dimension>::type  &boundary_indicators,
+                                         const std::vector<unsigned int>  &dof_to_boundary_mapping,
+                                         CompressedBlockSparsityPattern    &sparsity);
+
+
 
 template void
 DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
@@ -2209,6 +2222,10 @@ DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
 template void
 DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
                                 BlockSparsityPattern                &sparsity);
+template void
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
+                                CompressedBlockSparsityPattern      &sparsity);
+
 
 template void 
 DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
@@ -2225,6 +2242,12 @@ DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
                                 const std::vector<std::vector<bool> > &mask,
                                 BlockSparsityPattern        &sparsity);
 
+template void 
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
+                                const std::vector<std::vector<bool> > &mask,
+                                CompressedBlockSparsityPattern        &sparsity);
+
+
 
 template
 void
index ca38ead500a743ba7e349be1b25a65116e46966d..d7480cd448986ba2639ac802cde18b4c4512ecd1 100644 (file)
 #include <base/subscriptor.h>
 #include <base/smartpointer.h>
 #include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
 #include <lac/block_indices.h>
 
 
+template <typename number> class BlockSparseMatrix;
+class BlockSparsityPattern;
+class CompressedBlockSparsityPattern;
+
+
 /**
- * This is a block version of the sparsity pattern class. It has not
- * much functionality, but only administrated an array of sparsity
- * pattern objects and delegates work to them. It has mostly the same
- * interface as has the @p{SparsityPattern} class, and simply transforms
+ * This is the base class for block versions of the sparsity pattern
+ * and compressed sparsity pattern classes. It has not much
+ * functionality, but only administrates an array of sparsity pattern
+ * objects and delegates work to them. It has mostly the same
+ * interface as has the @ref{SparsityPattern} and
+ * @ref{CompressedSparsityPattern} classes, and simply transforms
  * calls to its member functions to calls to the respective member
  * functions of the member sparsity patterns.
  *
- * The largest difference between the @p{SparsityPattern} class and
- * this class is that mostly, the matrices have different properties
- * and you will want to work on the blocks making up the matrix rather
- * than the whole matrix. You can access the different blocks using
- * the @p{block(row,col)} function.
+ * The largest difference between the @ref{SparsityPattern} and
+ * @ref{CompressedSparsityPattern} classes and this class is that
+ * mostly, the matrices have different properties and you will want to
+ * work on the blocks making up the matrix rather than the whole
+ * matrix. You can access the different blocks using the
+ * @p{block(row,col)} function.
  *
  * Attention: this object is not automatically notified if the size of
  * one of its subobjects' size is changed. After you initialize the
  * that all sub-matrices in a column have to have the same number of
  * columns.
  *
- * @author Wolfgang Bangerth, 2000
+ * You will in general not want to use this class, but one of the
+ * derived classes.
+ *
+ * @author Wolfgang Bangerth, 2000, 2001
  */
-class BlockSparsityPattern : public Subscriptor
+template <class SparsityPatternBase>
+class BlockSparsityPatternBase : public Subscriptor
 {
   public:
                                     /**
@@ -72,7 +85,7 @@ class BlockSparsityPattern : public Subscriptor
                                      * structure usable by calling
                                      * the @p{reinit} function.
                                      */
-    BlockSparsityPattern ();
+    BlockSparsityPatternBase ();
 
                                     /**
                                      * Initialize the matrix with the
@@ -82,8 +95,8 @@ class BlockSparsityPattern : public Subscriptor
                                      * to call @p{collect_args} after
                                      * you assign them sizes.
                                      */
-    BlockSparsityPattern (const unsigned int n_rows,
-                         const unsigned int n_columns);
+    BlockSparsityPatternBase (const unsigned int n_rows,
+                             const unsigned int n_columns);
 
                                     /**
                                      * Copy constructor. This
@@ -96,12 +109,12 @@ class BlockSparsityPattern : public Subscriptor
                                      * @p{SparsityPattern}, see there
                                      * for the details.
                                      */
-    BlockSparsityPattern (const BlockSparsityPattern &bsp);
+    BlockSparsityPatternBase (const BlockSparsityPatternBase &bsp);
 
                                     /**
                                      * Destructor.
                                      */
-    ~BlockSparsityPattern ();
+    ~BlockSparsityPatternBase ();
     
                                     /**
                                      * Resize the matrix. This
@@ -140,7 +153,7 @@ class BlockSparsityPattern : public Subscriptor
                                      * but the latter only for empty
                                      * objects.
                                      */
-    BlockSparsityPattern & operator = (const BlockSparsityPattern &);
+    BlockSparsityPatternBase & operator = (const BlockSparsityPatternBase &);
 
                                     /**
                                      * This function collects the
@@ -160,7 +173,7 @@ class BlockSparsityPattern : public Subscriptor
                                      * Access the block with the
                                      * given coordinates.
                                      */
-    SparsityPattern &
+    SparsityPatternBase &
     block (const unsigned int row,
           const unsigned int column);
     
@@ -170,7 +183,7 @@ class BlockSparsityPattern : public Subscriptor
                                      * given coordinates. Version for
                                      * constant objects.
                                      */
-    const SparsityPattern &
+    const SparsityPatternBase &
     block (const unsigned int row,
           const unsigned int column) const;    
 
@@ -290,21 +303,6 @@ class BlockSparsityPattern : public Subscriptor
                                      */
     unsigned int n_nonzero_elements () const;
 
-                                    /**
-                                     * Return whether the structure
-                                     * is compressed or not,
-                                     * i.e. whether all sub-matrices
-                                     * are compressed.
-                                     */
-    bool is_compressed () const;
-
-                                    /**
-                                     * Determine an estimate for the
-                                     * memory consumption (in bytes)
-                                     * of this object.
-                                     */
-    unsigned int memory_consumption () const;
-
                                     /**
                                      * Exception
                                      */
@@ -334,7 +332,7 @@ class BlockSparsityPattern : public Subscriptor
                    << "The number of blocks " << arg1 << " and " << arg2
                    << " are different.");
     
-  private:
+  protected:
 
                                     /**
                                      * Number of block rows.
@@ -349,7 +347,7 @@ class BlockSparsityPattern : public Subscriptor
                                     /**
                                      * Array of sparsity patterns.
                                      */
-    std::vector<std::vector<SmartPointer<SparsityPattern> > > sub_objects;
+    std::vector<std::vector<SmartPointer<SparsityPatternBase> > > sub_objects;
 
                                     /**
                                      * Object storing and managing
@@ -379,23 +377,132 @@ class BlockSparsityPattern : public Subscriptor
 
 
 
+/**
+ * This class extends the base class to implement an array of sparsity
+ * patterns that can be used by block sparse matrix objects. It only
+ * adds a few additional member functions, but the main interface
+ * stems from the base class, see there for more information.
+ *
+ * @author Wolfgang Bangerth, 2000, 2001
+ */
+class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
+{
+  public:
+    
+                                    /**
+                                     * Initialize the matrix empty,
+                                     * that is with no memory
+                                     * allocated. This is useful if
+                                     * you want such objects as
+                                     * member variables in other
+                                     * classes. You can make the
+                                     * structure usable by calling
+                                     * the @p{reinit} function.
+                                     */
+    BlockSparsityPattern ();
+
+                                    /**
+                                     * Initialize the matrix with the
+                                     * given number of block rows and
+                                     * columns. The blocks themselves
+                                     * are still empty, and you have
+                                     * to call @p{collect_args} after
+                                     * you assign them sizes.
+                                     */
+    BlockSparsityPattern (const unsigned int n_rows,
+                         const unsigned int n_columns);
+
+                                    /**
+                                     * Return whether the structure
+                                     * is compressed or not,
+                                     * i.e. whether all sub-matrices
+                                     * are compressed.
+                                     */
+    bool is_compressed () const;    
+
+                                    /**
+                                     * Determine an estimate for the
+                                     * memory consumption (in bytes)
+                                     * of this object.
+                                     */
+    unsigned int memory_consumption () const;
+
+                                    /**
+                                     * Copy data from an object of
+                                     * type
+                                     * @ref{CompressedBlockSparsityPattern},
+                                     * i.e. resize this object to the
+                                     * size of the given argument,
+                                     * and copy over the contents of
+                                     * each of the
+                                     * subobjects. Previous content
+                                     * of this object is lost.
+                                     */
+    void copy_from (const CompressedBlockSparsityPattern &csp);
+};
+
+
+
+/**
+ * This class extends the base class to implement an array of
+ * compressed sparsity patterns that can be used to initialize objects
+ * of type @ref{BlockSparsityPattern}. It does not add additional
+ * member functions, but rather acts as a @p{typedef} to introduce the
+ * name of this class, without requiring the user to specify the
+ * templated name of the base class. For information on the interface
+ * of this class refer to the base class.
+ *
+ * @author Wolfgang Bangerth, 2000, 2001
+ */
+class CompressedBlockSparsityPattern : public BlockSparsityPatternBase<CompressedSparsityPattern>
+{
+  public:
+    
+                                    /**
+                                     * Initialize the matrix empty,
+                                     * that is with no memory
+                                     * allocated. This is useful if
+                                     * you want such objects as
+                                     * member variables in other
+                                     * classes. You can make the
+                                     * structure usable by calling
+                                     * the @p{reinit} function.
+                                     */
+    CompressedBlockSparsityPattern ();
+
+                                    /**
+                                     * Initialize the matrix with the
+                                     * given number of block rows and
+                                     * columns. The blocks themselves
+                                     * are still empty, and you have
+                                     * to call @p{collect_args} after
+                                     * you assign them sizes.
+                                     */
+    CompressedBlockSparsityPattern (const unsigned int n_rows,
+                                   const unsigned int n_columns);
+};
+
+
+
 /*---------------------- Template functions -----------------------------------*/
 
 
 
+template <class SparsityPatternBase>
 inline
-SparsityPattern &
-BlockSparsityPattern::block (const unsigned int row,
-                            const unsigned int column)
+SparsityPatternBase &
+BlockSparsityPatternBase<SparsityPatternBase>::block (const unsigned int row,
+                                                 const unsigned int column)
 {
   return *sub_objects[row][column];
 };
 
 
 
+template <class SparsityPatternBase>
 inline
-const SparsityPattern &
-BlockSparsityPattern::block (const unsigned int row,
+const SparsityPatternBase &
+BlockSparsityPatternBase<SparsityPatternBase>::block (const unsigned int row,
                             const unsigned int column) const
 {
   return *sub_objects[row][column];
@@ -403,28 +510,31 @@ BlockSparsityPattern::block (const unsigned int row,
 
 
 
+template <class SparsityPatternBase>
 inline
 const BlockIndices &
-BlockSparsityPattern::get_row_indices () const
+BlockSparsityPatternBase<SparsityPatternBase>::get_row_indices () const
 {
   return row_indices;
 };
 
 
 
+template <class SparsityPatternBase>
 inline
 const BlockIndices &
-BlockSparsityPattern::get_column_indices () const
+BlockSparsityPatternBase<SparsityPatternBase>::get_column_indices () const
 {
   return column_indices;
 };
 
 
 
+template <class SparsityPatternBase>
 inline
 void
-BlockSparsityPattern::add (const unsigned int i,
-                          const unsigned int j)
+BlockSparsityPatternBase<SparsityPatternBase>::add (const unsigned int i,
+                                               const unsigned int j)
 {
                                   // if you get an error here, are
                                   // you sure you called
@@ -438,18 +548,20 @@ BlockSparsityPattern::add (const unsigned int i,
 
 
 
+template <class SparsityPatternBase>
 inline
 unsigned int
-BlockSparsityPattern::n_block_cols () const
+BlockSparsityPatternBase<SparsityPatternBase>::n_block_cols () const
 {
   return columns;
 }
 
 
 
+template <class SparsityPatternBase>
 inline
 unsigned int
-BlockSparsityPattern::n_block_rows () const
+BlockSparsityPatternBase<SparsityPatternBase>::n_block_rows () const
 {
   return rows;
 }
index 3ffa3bb07a4d7ec8724a92d3bec6613f5b91371b..f6ef2959aa2e432c88f5cfb973d4ca7cd6ead410 100644 (file)
@@ -16,7 +16,8 @@
 
 
 
-BlockSparsityPattern::BlockSparsityPattern ()
+template <class SparsityPatternBase>
+BlockSparsityPatternBase<SparsityPatternBase>::BlockSparsityPatternBase ()
                :
                rows (0),
                columns (0)
@@ -24,8 +25,10 @@ BlockSparsityPattern::BlockSparsityPattern ()
 
 
 
-BlockSparsityPattern::BlockSparsityPattern (const unsigned int r,
-                                           const unsigned int c)
+template <class SparsityPatternBase>
+BlockSparsityPatternBase<SparsityPatternBase>::
+BlockSparsityPatternBase (const unsigned int r,
+                     const unsigned int c)
                :
                rows (0),
                columns (0)
@@ -35,7 +38,8 @@ BlockSparsityPattern::BlockSparsityPattern (const unsigned int r,
 
 
 
-BlockSparsityPattern::~BlockSparsityPattern ()
+template <class SparsityPatternBase>
+BlockSparsityPatternBase<SparsityPatternBase>::~BlockSparsityPatternBase ()
 {
                                   // clear all memory
   reinit (0,0);
@@ -43,15 +47,16 @@ BlockSparsityPattern::~BlockSparsityPattern ()
 
 
 
+template <class SparsityPatternBase>
 void
-BlockSparsityPattern::reinit (const unsigned int r,
-                             const unsigned int c)
+BlockSparsityPatternBase<SparsityPatternBase>::reinit (const unsigned int r,
+                                                      const unsigned int c)
 {
                                   // delete previous content
   for (unsigned int i=0; i<rows; ++i)
     for (unsigned int j=0; j<columns; ++j)
       {
-       SparsityPattern * sp = sub_objects[i][j];
+       SparsityPatternBase * sp = sub_objects[i][j];
        sub_objects[i][j] = 0;
        delete sp;
       };
@@ -60,18 +65,20 @@ BlockSparsityPattern::reinit (const unsigned int r,
                                   // set new sizes
   rows = r;
   columns = c;
-  sub_objects.resize (rows, std::vector<SmartPointer<SparsityPattern> > (columns));
+  sub_objects.resize (rows, std::vector<SmartPointer<SparsityPatternBase> > (columns));
 
                                   // allocate new objects
   for (unsigned int i=0; i<rows; ++i)
     for (unsigned int j=0; j<columns; ++j)
-      sub_objects[i][j] = new SparsityPattern;
+      sub_objects[i][j] = new SparsityPatternBase;
 };
 
 
 
-BlockSparsityPattern &
-BlockSparsityPattern::operator = (const BlockSparsityPattern &bsp)
+template <class SparsityPatternBase>
+BlockSparsityPatternBase<SparsityPatternBase> &
+BlockSparsityPatternBase<SparsityPatternBase>::
+operator = (const BlockSparsityPatternBase<SparsityPatternBase> &bsp)
 {
   Assert (rows == bsp.rows, ExcIncompatibleSizes(rows, bsp.rows));
   Assert (columns == bsp.columns, ExcIncompatibleSizes(columns, bsp.columns));
@@ -84,12 +91,12 @@ BlockSparsityPattern::operator = (const BlockSparsityPattern &bsp)
 
   return *this;
 };
-
   
 
 
+template <class SparsityPatternBase>
 void
-BlockSparsityPattern::collect_sizes ()
+BlockSparsityPatternBase<SparsityPatternBase>::collect_sizes ()
 {
   std::vector<unsigned int> row_sizes (rows);
   std::vector<unsigned int> col_sizes (columns);
@@ -126,8 +133,9 @@ BlockSparsityPattern::collect_sizes ()
 
 
 
+template <class SparsityPatternBase>
 void
-BlockSparsityPattern::compress ()
+BlockSparsityPatternBase<SparsityPatternBase>::compress ()
 {
   for (unsigned int i=0; i<rows; ++i)
     for (unsigned int j=0; j<columns; ++j)
@@ -136,8 +144,9 @@ BlockSparsityPattern::compress ()
 
 
 
+template <class SparsityPatternBase>
 bool
-BlockSparsityPattern::empty () const
+BlockSparsityPatternBase<SparsityPatternBase>::empty () const
 {
   for (unsigned int i=0; i<rows; ++i)
     for (unsigned int j=0; j<columns; ++j)
@@ -148,8 +157,9 @@ BlockSparsityPattern::empty () const
 
 
 
+template <class SparsityPatternBase>
 unsigned int
-BlockSparsityPattern::max_entries_per_row () const
+BlockSparsityPatternBase<SparsityPatternBase>::max_entries_per_row () const
 {
   unsigned int max_entries = 0;
   for (unsigned int block_row=0; block_row<rows; ++block_row)
@@ -165,8 +175,10 @@ BlockSparsityPattern::max_entries_per_row () const
 };
 
 
+
+template <class SparsityPatternBase>
 unsigned int
-BlockSparsityPattern::n_rows () const
+BlockSparsityPatternBase<SparsityPatternBase>::n_rows () const
 {
                                   // only count in first column, since
                                   // all rows should be equivalent
@@ -178,8 +190,9 @@ BlockSparsityPattern::n_rows () const
 
 
 
+template <class SparsityPatternBase>
 unsigned int
-BlockSparsityPattern::n_cols () const
+BlockSparsityPatternBase<SparsityPatternBase>::n_cols () const
 {
                                   // only count in first row, since
                                   // all rows should be equivalent
@@ -191,10 +204,9 @@ BlockSparsityPattern::n_cols () const
 
 
 
-
-
+template <class SparsityPatternBase>
 unsigned int
-BlockSparsityPattern::n_nonzero_elements () const
+BlockSparsityPatternBase<SparsityPatternBase>::n_nonzero_elements () const
 {
   unsigned int count = 0;
   for (unsigned int i=0; i<rows; ++i)
@@ -205,6 +217,20 @@ BlockSparsityPattern::n_nonzero_elements () const
 
 
 
+BlockSparsityPattern::BlockSparsityPattern ()
+{};
+
+
+
+BlockSparsityPattern::BlockSparsityPattern (const unsigned int n_rows,
+                                           const unsigned int n_columns)
+               :
+               BlockSparsityPatternBase<SparsityPattern>(n_rows,
+                                                         n_columns)
+{};
+
+
+
 bool
 BlockSparsityPattern::is_compressed () const
 {
@@ -216,6 +242,7 @@ BlockSparsityPattern::is_compressed () const
 };
 
 
+
 unsigned int
 BlockSparsityPattern::memory_consumption () const
 {
@@ -231,3 +258,43 @@ BlockSparsityPattern::memory_consumption () const
   
   return mem;
 };
+
+
+
+void
+BlockSparsityPattern::copy_from  (const CompressedBlockSparsityPattern &csp)
+{
+                                  // delete old content, set block
+                                  // sizes anew
+  reinit (csp.n_block_rows(), csp.n_block_cols());
+
+                                  // copy over blocks
+  for (unsigned int i=0; i<rows; ++i)
+    for (unsigned int j=0; j<rows; ++j)
+      block(i,j).copy_from (csp.block(i,j));
+
+                                  // and finally enquire their new
+                                  // sizes
+  collect_sizes();
+};
+
+
+
+CompressedBlockSparsityPattern::CompressedBlockSparsityPattern ()
+{};
+
+
+
+CompressedBlockSparsityPattern::
+CompressedBlockSparsityPattern (const unsigned int n_rows,
+                               const unsigned int n_columns)
+               :
+               BlockSparsityPatternBase<CompressedSparsityPattern>(n_rows,
+                                                                   n_columns)
+{};
+
+
+
+// explicit instantiations
+template class BlockSparsityPatternBase<SparsityPattern>;
+template class BlockSparsityPatternBase<CompressedSparsityPattern>;
index ec1870835160f893ebda978d1d1519f794f36f0b..de8c17e3769d81039b72cd503c67e7fb766c4347 100644 (file)
@@ -20,6 +20,7 @@
 
 #include <base/logstream.h>
 #include <lac/sparsity_pattern.h>
+#include <lac/block_sparsity_pattern.h>
 #include <lac/compressed_sparsity_pattern.h>
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
@@ -56,6 +57,25 @@ bool operator == (const SparsityPattern &sp1,
 
 
 
+bool operator == (const BlockSparsityPattern &sp1,
+                 const BlockSparsityPattern &sp2)
+{
+  if (sp1.n_block_rows() != sp2.n_block_rows())
+    return false;
+  
+  if (sp1.n_block_cols() != sp2.n_block_cols())
+    return false;
+  
+  for (unsigned int i=0; i<sp1.n_block_rows(); ++i)
+    for (unsigned int j=0; j<sp1.n_block_cols(); ++j)
+      if (!(sp1.block(i,j) == sp2.block(i,j)))
+       return false;
+  
+  return true;
+};
+
+
+
 template <int dim>
 void
 check_boundary (const DoFHandler<dim> &dof)
@@ -84,7 +104,7 @@ check_boundary (const DoFHandler<dim> &dof)
                                   // patterns is checked in other
                                   // tests, so only make sure that
                                   // sparsity_[12] are equal
-  deallog << __PRETTY_FUNCTION__
+  deallog << "Check boundary"
          << " -- "
          << (sparsity_1 == sparsity_2 ? "ok" : "failed")
          << std::endl;
@@ -127,29 +147,76 @@ check ()
   std::vector<std::vector<bool> > mask (2, std::vector<bool>(2, false));
   mask[0][0] = mask[1][1] = true;
 
+  
+//--------------- Regular sparsity pattern checks -----------------
+  
                                   // first way: directly
   SparsityPattern sparsity_1 (dof.n_dofs(), dof.n_dofs());
   DoFTools::make_sparsity_pattern (dof, mask, sparsity_1);
   constraints.condense (sparsity_1);
   sparsity_1.compress ();
 
-                                  // second way: via CompressedSparsityPattern
+                                  // second way: via
+                                  // CompressedSparsityPattern
   SparsityPattern sparsity_2;
-  CompressedSparsityPattern csp (dof.n_dofs());
-  DoFTools::make_sparsity_pattern (dof, mask, csp);
-  constraints.condense (csp);
-  sparsity_2.copy_from (csp);
+  CompressedSparsityPattern csp_2 (dof.n_dofs());
+  DoFTools::make_sparsity_pattern (dof, mask, csp_2);
+  constraints.condense (csp_2);
+  sparsity_2.copy_from (csp_2);
 
 
                                   // the exact content of sparsity
                                   // patterns is checked in other
                                   // tests, so only make sure that
                                   // sparsity_[12] are equal
-  deallog << __PRETTY_FUNCTION__
+  deallog << "Check 1:"
          << " -- "
          << (sparsity_1 == sparsity_2 ? "ok" : "failed")
          << std::endl;
+
+
+  
+//--------------- Block sparsity pattern checks -----------------
+
+  const unsigned int n  = dof.n_dofs();
+  const unsigned int n1 = n/3;
+  const unsigned int n2 = n - n1;
+
+  BlockSparsityPattern sparsity_3(2,2);
+  sparsity_3.block(0,0).reinit (n1,n1,n);
+  sparsity_3.block(1,0).reinit (n2,n1,n);
+  sparsity_3.block(0,1).reinit (n1,n2,n);
+  sparsity_3.block(1,1).reinit (n2,n2,n);
+  sparsity_3.collect_sizes ();
+
+  DoFTools::make_sparsity_pattern (dof, sparsity_3);
+  constraints.condense (sparsity_3);
+  sparsity_3.compress ();
+
+  BlockSparsityPattern sparsity_4;
+  CompressedBlockSparsityPattern csp_4(2,2);
+  csp_4.block(0,0).reinit (n1,n1);
+  csp_4.block(1,0).reinit (n2,n1);
+  csp_4.block(0,1).reinit (n1,n2);
+  csp_4.block(1,1).reinit (n2,n2);
+  csp_4.collect_sizes ();
+
+  DoFTools::make_sparsity_pattern (dof, csp_4);
+  constraints.condense (csp_4);
+  csp_4.compress ();
+
+  sparsity_4.copy_from (csp_4);
+
+  deallog << "Check 2:"
+         << " -- "
+         << (sparsity_3 == sparsity_4 ? "ok" : "failed")
+         << std::endl;
+
   
+//--------------- Sparsity pattern checks for
+//                boundary sparsity generators -----------------
+
+                                  // check boundary matrices
   check_boundary (dof);
 };
 
index 5caca2e11949bdd71d4074f92430069a683ffbdf..76b2767418423000624d425a21dddf7d8815f975 100644 (file)
@@ -1,7 +1,10 @@
 
-DEAL:1d::void check<1>() -- ok
-DEAL:1d::void check_boundary<1>(const DoFHandler<1> &) -- ok
-DEAL:2d::void check<2>() -- ok
-DEAL:2d::void check_boundary<2>(const DoFHandler<2> &) -- ok
-DEAL:3d::void check<3>() -- ok
-DEAL:3d::void check_boundary<3>(const DoFHandler<3> &) -- ok
+DEAL:1d::Check 1: -- ok
+DEAL:1d::Check 2: -- ok
+DEAL:1d::Check boundary -- ok
+DEAL:2d::Check 1: -- ok
+DEAL:2d::Check 2: -- ok
+DEAL:2d::Check boundary -- ok
+DEAL:3d::Check 1: -- ok
+DEAL:3d::Check 2: -- ok
+DEAL:3d::Check boundary -- ok

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.