]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add some more functions for block sparse matrices. Not all yet tested.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 18 May 2000 13:02:53 +0000 (13:02 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 18 May 2000 13:02:53 +0000 (13:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@2889 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/include/dofs/dof_constraints.templates.h
deal.II/deal.II/source/dofs/dof_constraints.cc

index 1c6b88f5d15f409479dd920e371086fa4f175a4c..0dd6b3c9b4e8e50a62566335dc7404dd324edfad 100644 (file)
@@ -234,6 +234,14 @@ class ConstraintMatrix : public Subscriptor
                                      */
     void condense (SparsityPattern &sparsity) const;
 
+                                    /**
+                                     * Same function as above, but
+                                     * condenses square block sparsity
+                                     * patterns.
+                                     */
+    template <int blocks>
+    void condense (BlockSparsityPattern<blocks,blocks> &sparsity) const;
+    
                                     /**
                                      * Condense a given matrix. The associated
                                      * matrix struct should be condensed and
@@ -259,49 +267,94 @@ class ConstraintMatrix : public Subscriptor
     template<typename number>
     void condense (SparseMatrix<number> &matrix) const;
 
+                                    /**
+                                     * Same function as above, but
+                                     * condenses square block sparse
+                                     * matrices.
+                                     */
+    template <typename number, int blocks>
+    void condense (BlockSparseMatrix<number,blocks,blocks> &sparsity) const;
+    
                                     /**
                                      * Condense the given vector #uncondensed#
                                      * into #condensed#. It is the user's
                                      * responsibility to guarantee that all
                                      * entries of #condensed# be zero!
+                                     *
+                                     * The #VectorType# may be a
+                                     * #Vector<float>#,
+                                     * #Vector<double>#,
+                                     * #BlockVector<...>#, or any
+                                     * other type having the same
+                                     * interface.
                                      */
-    template<typename number>
-    void condense (const Vector<number> &uncondensed,
-                  Vector<number>       &condensed) const;
+    template <class VectorType>
+    void condense (const VectorType &uncondensed,
+                  VectorType       &condensed) const;
 
                                     /**
-                                     * Condense the given vector in-place.
+                                     * Condense the given vector
+                                     * in-place. The #VectorType# may
+                                     * be a #Vector<float>#,
+                                     * #Vector<double>#,
+                                     * #BlockVector<...>#, or any
+                                     * other type having the same
+                                     * interface.
                                      */
-    template<typename number>
-    void condense (Vector<number> &vec) const;
+    template <class VectorType>
+    void condense (VectorType &vec) const;
 
                                     /**
-                                     * Re-distribute the elements of the vector
-                                     * #condensed# to #uncondensed#. It is the user's
-                                     * responsibility to guarantee that all
-                                     * entries of #uncondensed# be zero!
+                                     * Re-distribute the elements of
+                                     * the vector #condensed# to
+                                     * #uncondensed#. It is the
+                                     * user's responsibility to
+                                     * guarantee that all entries of
+                                     * #uncondensed# be zero!
+                                     *
+                                     * This function undoes the
+                                     * action of #condense# somehow,
+                                     * but it should be noted that it
+                                     * is not the inverse of
+                                     * #condense#.
                                      *
-                                     * This function undoes the action of
-                                     * #condense# somehow, but it should be noted
-                                     * that it is not the inverse of #condense#-
+                                     * The #VectorType# may be a
+                                     * #Vector<float>#,
+                                     * #Vector<double>#,
+                                     * #BlockVector<...>#, or any
+                                     * other type having the same
+                                     * interface.
                                      */
-    template<typename number>
-    void distribute (const Vector<number> &condensed,
-                    Vector<number>       &uncondensed) const;
+    template <class VectorType>
+    void distribute (const VectorType &condensed,
+                    VectorType       &uncondensed) const;
 
                                     /**
-                                     * Re-distribute the elements of the vector
-                                     * in-place.
+                                     * Re-distribute the elements of
+                                     * the vector in-place. The
+                                     * #VectorType# may be a
+                                     * #Vector<float>#,
+                                     * #Vector<double>#,
+                                     * #BlockVector<...>#, or any
+                                     * other type having the same
+                                     * interface.
                                      */
-    template<typename number>
-    void distribute (Vector<number> &vec) const;
+    template <class VectorType>
+    void distribute (VectorType &vec) const;
     
                                     /**
-                                     * Delete hanging nodes in a vector.
-                                     * Sets all hanging node values to zero.
+                                     * Delete hanging nodes in a
+                                     * vector.  Sets all hanging node
+                                     * values to zero. The
+                                     * #VectorType# may be a
+                                     * #Vector<float>#,
+                                     * #Vector<double>#,
+                                     * #BlockVector<...>#, or any
+                                     * other type having the same
+                                     * interface.
                                      */
-    template<typename number>
-    void set_zero (Vector<number> &vec) const;
+    template <class VectorType>
+    void set_zero (VectorType &vec) const;
 
 
                                     /**
index 31cc4664dc81a30e3eb63a81d77fae5a14469c19..03c0de63d5eb55e9405f7549bb1d456e51e4e68b 100644 (file)
@@ -15,7 +15,6 @@
 
 
 #include <dofs/dof_constraints.h>
-#include <lac/sparse_matrix.h>
 
 
 template<typename number>
@@ -261,10 +260,160 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed) const
 
 
 
-template<typename number>
+template <typename number, int blocks>
 void
-ConstraintMatrix::condense (const Vector<number> &uncondensed,
-                           Vector<number>       &condensed) const
+ConstraintMatrix::condense (BlockSparseMatrix<number,blocks,blocks> &uncondensed) const
+{
+  const BlockSparsityPattern<blocks,blocks> &
+    sparsity = uncondensed.get_sparsity_pattern ();
+
+  Assert (sorted == true, ExcMatrixNotClosed());
+  Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed());
+  Assert (sparsity.n_rows() == sparsity.n_cols(),
+         ExcMatrixNotSquare());
+  Assert (sparsity.get_column_indices() == sparsity.get_row_indices(),
+         ExcMatrixNotSquare());
+
+  const BlockIndices<blocks> &
+    index_mapping = sparsity.get_column_indices();
+  
+                                  // 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
+  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)
+    {
+      if (distribute[row] == -1)
+       {
+                                          // get index of this row
+                                          // within the blocks
+         const pair<unsigned int,unsigned int>
+           block_index = index_mapping.global_to_local(row);
+         
+                                          // regular line. loop over
+                                          // all columns and see
+                                          // whether this column must
+                                          // be distributed
+         for (unsigned int block_col=0; block_col<blocks; ++block_col)
+           {
+             const SparsityPattern &
+               block_sparsity = sparsity.block(block_index.first,block_col);
+             
+             const unsigned int
+               first = block_sparsity.get_rowstart_indices()[block_index.second],
+               last  = block_sparsity.get_rowstart_indices()[block_index.second+1];
+             for (unsigned int j=first; j<last; ++j)
+                                                // end of row reached?
+               if (block_sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+                 {
+                                                    // this should not happen, since
+                                                    // we only operate on compressed
+                                                    // matrices!
+                   Assert (false, ExcMatrixNotClosed());
+                   break;
+                 }
+               else
+                 {
+                   const unsigned int global_col
+                     =index_mapping.local_to_global(block_col,
+                                                    block_sparsity.get_column_numbers()[j]);
+                   
+                   if (distribute[global_col] != -1)
+                                                      // distribute entry at regular
+                                                      // row #row# and irregular column
+                                                      // global_col;
+                                                      // set old entry to zero
+                     {
+                       const double old_value =
+                         uncondensed.block(block_index.first,block_col).global_entry(j);
+                       
+                       for (unsigned int q=0;
+                            q!=lines[distribute[global_col]]
+                                           .entries.size(); ++q)
+                         uncondensed.add (row,
+                                          lines[distribute[global_col]].entries[q].first,
+                                          old_value *
+                                          lines[distribute[global_col]].entries[q].second);
+                       
+                       uncondensed.block(block_index.first,block_col).global_entry(j)
+                         = 0.;
+                     };
+                 };
+           };
+       }
+      else
+       {
+       }
+      
+//                                      // row must be distributed
+//     for (unsigned int j=sparsity.get_rowstart_indices()[row];
+//          j<sparsity.get_rowstart_indices()[row+1]; ++j)
+//                                        // end of row reached?
+//       if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+//         {
+//                                            // this should not happen, since
+//                                            // we only operate on compressed
+//                                            // matrices!
+//           Assert (false, ExcMatrixNotClosed());
+//           break;
+//         }
+//       else
+//         {
+//           if (distribute[sparsity.get_column_numbers()[j]] == -1)
+//                                              // distribute entry at irregular
+//                                              // row #row# and regular column
+//                                              // sparsity.get_column_numbers()[j]. set old
+//                                              // entry to zero
+//             {
+//               for (unsigned int q=0;
+//                    q!=lines[distribute[row]].entries.size(); ++q) 
+//                 uncondensed.add (lines[distribute[row]].entries[q].first,
+//                                  sparsity.get_column_numbers()[j],
+//                                  uncondensed.global_entry(j) *
+//                                  lines[distribute[row]].entries[q].second);
+               
+//               uncondensed.global_entry(j) = 0.;
+//             }
+//           else
+//                                              // distribute entry at irregular
+//                                              // row #row# and irregular column
+//                                              // sparsity.get_column_numbers()[j]
+//                                              // set old entry to one if on main
+//                                              // diagonal, zero otherwise
+//             {
+//               for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
+//                 for (unsigned int q=0;
+//                      q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+//                                     .entries.size(); ++q)
+//                   uncondensed.add (lines[distribute[row]].entries[p].first,
+//                                    lines[distribute[sparsity.get_column_numbers()[j]]]
+//                                    .entries[q].first,
+//                                    uncondensed.global_entry(j) *
+//                                    lines[distribute[row]].entries[p].second *
+//                                    lines[distribute[sparsity.get_column_numbers()[j]]]
+//                                    .entries[q].second);
+               
+//               uncondensed.global_entry(j) = (row == sparsity.get_column_numbers()[j] ?
+//                                              1. : 0. );
+//             };
+//         };
+    };
+};
+
+
+
+template<class VectorType>
+void
+ConstraintMatrix::condense (const VectorType &uncondensed,
+                           VectorType       &condensed) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
   Assert (condensed.size()+n_constraints() == uncondensed.size(),
@@ -335,9 +484,9 @@ ConstraintMatrix::condense (const Vector<number> &uncondensed,
 
 
 
-template<typename number>
+template <class VectorType>
 void
-ConstraintMatrix::condense (Vector<number> &vec) const
+ConstraintMatrix::condense (VectorType &vec) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
 
@@ -366,9 +515,9 @@ ConstraintMatrix::condense (Vector<number> &vec) const
 
 
 
-template<typename number>
+template <class VectorType>
 void
-ConstraintMatrix::set_zero (Vector<number> &vec) const
+ConstraintMatrix::set_zero (VectorType &vec) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
 
@@ -392,10 +541,10 @@ ConstraintMatrix::set_zero (Vector<number> &vec) const
 
 
 
-template<typename number>
+template<class VectorType>
 void
-ConstraintMatrix::distribute (const Vector<number> &condensed,
-                             Vector<number>       &uncondensed) const
+ConstraintMatrix::distribute (const VectorType &condensed,
+                             VectorType       &uncondensed) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
   Assert (condensed.size()+n_constraints() == uncondensed.size(),
@@ -465,9 +614,9 @@ ConstraintMatrix::distribute (const Vector<number> &condensed,
 
 
 
-template<typename number>
+template<class VectorType>
 void
-ConstraintMatrix::distribute (Vector<number> &vec) const
+ConstraintMatrix::distribute (VectorType &vec) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
 
index 1df1b5b143927b4212b1a83ce20d8807562957e3..c6d6aa4ff5b140747bf3c0340e1e1d0b6eea16a8 100644 (file)
 
 
 #include <dofs/dof_constraints.h>
+#include <dofs/dof_constraints.templates.h>
+
 #include <lac/sparsity_pattern.h>
 #include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
 #include <iostream>
 #include <algorithm>
 #include <numeric>
@@ -28,11 +33,15 @@ ConstraintMatrix::ConstraintLine::operator < (const ConstraintLine &a) const
 };
 
 
+
 ConstraintMatrix::ConstraintMatrix () :
-               lines(), sorted(false) {};
+               lines(),
+               sorted(false)
+{};
 
 
-void ConstraintMatrix::add_line (const unsigned int line) {
+void ConstraintMatrix::add_line (const unsigned int line)
+{
   Assert (sorted==false, ExcMatrixIsClosed());
 
                                   // check whether line already exists;
@@ -48,9 +57,11 @@ void ConstraintMatrix::add_line (const unsigned int line) {
 };
 
 
+
 void ConstraintMatrix::add_entry (const unsigned int line,
                                  const unsigned int column,
-                                 const double       value) {
+                                 const double       value)
+{
   Assert (sorted==false, ExcMatrixIsClosed());
 
   vector<ConstraintLine>::iterator line_ptr;
@@ -89,7 +100,9 @@ void ConstraintMatrix::add_entry (const unsigned int line,
 };
 
 
-void ConstraintMatrix::close () {
+
+void ConstraintMatrix::close ()
+{
   Assert (sorted==false, ExcMatrixIsClosed());
 
                                   // sort the entries in the different lines
@@ -148,14 +161,18 @@ void ConstraintMatrix::close () {
 };
 
 
-void ConstraintMatrix::clear () {
+
+void ConstraintMatrix::clear ()
+{
   lines = vector<ConstraintLine>();
   sorted = false;
 };
 
 
+
 void ConstraintMatrix::condense (const SparsityPattern &uncondensed,
-                                SparsityPattern       &condensed) const {
+                                SparsityPattern       &condensed) const
+{
   Assert (sorted == true, ExcMatrixNotClosed());
   Assert (uncondensed.is_compressed() == true, ExcMatrixNotClosed());
   Assert (uncondensed.n_rows() == uncondensed.n_cols(),
@@ -262,7 +279,9 @@ void ConstraintMatrix::condense (const SparsityPattern &uncondensed,
 };
 
 
-void ConstraintMatrix::condense (SparsityPattern &sparsity) const {
+
+void ConstraintMatrix::condense (SparsityPattern &sparsity) const
+{
   Assert (sorted == true, ExcMatrixNotClosed());
   Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed());
   Assert (sparsity.n_rows() == sparsity.n_cols(),
@@ -276,6 +295,91 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const {
                                   // this index
   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);
+
+  unsigned int n_rows = sparsity.n_rows();
+  for (unsigned int row=0; row<n_rows; ++row)
+    {
+      if (distribute[row] == -1)
+                                        // regular line. loop over cols
+       for (unsigned int j=sparsity.get_rowstart_indices()[row];
+            j<sparsity.get_rowstart_indices()[row+1]; ++j)
+         {
+           const unsigned int column = sparsity.get_column_numbers()[j];
+           
+                                            // end of row reached?
+           if (column == SparsityPattern::invalid_entry)
+             break;
+           else
+             if (distribute[column] != -1)
+               {
+                                                  // distribute entry at regular
+                                                  // row #row# and irregular column
+                                                  // sparsity.colnums[j]
+                 for (unsigned int q=0;
+                      q!=lines[distribute[column]].entries.size();
+                      ++q) 
+                   sparsity.add (row,
+                                 lines[distribute[column]].entries[q].first);
+               };
+         }
+      else
+                                        // row must be distributed
+       for (unsigned int j=sparsity.get_rowstart_indices()[row];
+            j<sparsity.get_rowstart_indices()[row+1]; ++j)
+                                          // end of row reached?
+         if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
+           break;
+         else
+           {
+             if (distribute[sparsity.get_column_numbers()[j]] == -1)
+                                                // distribute entry at irregular
+                                                // row #row# and regular column
+                                                // sparsity.colnums[j]
+               for (unsigned int q=0;
+                    q!=lines[distribute[row]].entries.size(); ++q) 
+                 sparsity.add (lines[distribute[row]].entries[q].first,
+                               sparsity.get_column_numbers()[j]);
+             else
+                                                // distribute entry at irregular
+                                                // row #row# and irregular column
+                                                // sparsity.get_column_numbers()[j]
+               for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
+                 for (unsigned int q=0;
+                      q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+                                     .entries.size(); ++q)
+                   sparsity.add (lines[distribute[row]].entries[p].first,
+                                 lines[distribute[sparsity.get_column_numbers()[j]]]
+                                 .entries[q].first);
+           };
+    };
+  
+  sparsity.compress();
+};
+
+
+
+template <int blocks>
+void ConstraintMatrix::condense (BlockSparsityPattern<blocks,blocks> &sparsity) const
+{
+  Assert (false, ExcInternalError());
+/*
+  Assert (sorted == true, ExcMatrixNotClosed());
+  Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed());
+  Assert (sparsity.n_rows() == sparsity.n_cols(),
+         ExcMatrixNotSquare());
+  Assert (sparsity.get_column_indices() == sparsity.get_row_indices(),
+         ExcMatrixNotSquare());
+  
+                                  // 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
+  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);
 
@@ -337,15 +441,18 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const {
     };
   
   sparsity.compress();
+*/  
 };
 
 
+
 unsigned int ConstraintMatrix::n_constraints () const
 {
   return lines.size();
 };
 
 
+
 bool ConstraintMatrix::is_constrained (const unsigned int index) const 
 {
   if (sorted == true)
@@ -397,32 +504,82 @@ void ConstraintMatrix::print (ostream &out) const {
   AssertThrow (out, ExcIO());
 };
 
-#include <dofs/dof_constraints.templates.h>
 
-#define number double
-template void ConstraintMatrix::condense(const SparseMatrix<number> &uncondensed,
-                                 SparseMatrix<number>       &condensed) const;
-template void ConstraintMatrix::condense(SparseMatrix<number> &uncondensed) const;
-template void ConstraintMatrix::condense(const Vector<number> &uncondensed,
-                           Vector<number>       &condensed) const;
-template void ConstraintMatrix::condense(Vector<number> &vec) const;
-template void ConstraintMatrix::set_zero(Vector<number> &vec) const;
-template void ConstraintMatrix::distribute(const Vector<number> &condensed,
-                                     Vector<number>       &uncondensed) const;
-template void ConstraintMatrix::distribute(Vector<number> &vec) const;
-
-#undef number
-#define number float
-
-template void ConstraintMatrix::condense(const SparseMatrix<number> &uncondensed,
-                                 SparseMatrix<number>       &condensed) const;
-template void ConstraintMatrix::condense(SparseMatrix<number> &uncondensed) const;
-template void ConstraintMatrix::condense(const Vector<number> &uncondensed,
-                           Vector<number>       &condensed) const;
-template void ConstraintMatrix::condense(Vector<number> &vec) const;
-template void ConstraintMatrix::set_zero(Vector<number> &vec) const;
-template void ConstraintMatrix::distribute(const Vector<number> &condensed,
-                                     Vector<number>       &uncondensed) const;
-template void ConstraintMatrix::distribute(Vector<number> &vec) const;
-
-#undef number
+
+
+
+// explicit instantiations
+//
+// define a list of functions for vectors and matrices, respectively,
+// where the vector/matrix can be replaced using a preprocessor
+// variable VectorType/MatrixType. note that we cannot do so by using
+// a preprocessor function with one arg, since
+// #vector_functions(BlockVector<2,double>)# is not recognized as one
+// arg, and putting parentheses around the arg yields incorrect
+// syntax...
+
+#define vector_functions \
+  template void ConstraintMatrix::condense(const VectorType &uncondensed,\
+                                          VectorType       &condensed) const;\
+  template void ConstraintMatrix::condense(VectorType &vec) const;\
+  template void ConstraintMatrix::set_zero(VectorType &vec) const;\
+  template void ConstraintMatrix::distribute(const VectorType &condensed,\
+                                            VectorType       &uncondensed) const;\
+  template void ConstraintMatrix::distribute(VectorType &vec) const;
+
+
+
+#define matrix_functions_1 \
+  template void ConstraintMatrix::condense(const MatrixType &uncondensed,\
+                                          MatrixType       &condensed) const;
+#define matrix_functions_2 \
+  template void ConstraintMatrix::condense(MatrixType &uncondensed) const;
+
+
+
+#define VectorType Vector<float>
+vector_functions;
+#undef VectorType
+
+#define VectorType Vector<double>
+vector_functions;
+#undef VectorType
+
+#define VectorType BlockVector<2,double>
+vector_functions;
+#undef VectorType
+
+#define VectorType BlockVector<3,double>
+vector_functions;
+#undef VectorType
+
+
+
+#define MatrixType SparseMatrix<float>
+matrix_functions_1;
+matrix_functions_2;
+#undef MatrixType
+
+#define MatrixType SparseMatrix<double>
+matrix_functions_1;
+matrix_functions_2;
+#undef MatrixType
+
+// block sparse matrices are only implemented for one of the two matrix functions
+#define MatrixType BlockSparseMatrix<double,2,2>
+matrix_functions_2;
+#undef MatrixType
+
+#define MatrixType BlockSparseMatrix<double,3,3>
+matrix_functions_2;
+#undef MatrixType
+
+
+#define MatrixType BlockSparsityPattern<2,2>
+matrix_functions_2;
+#undef MatrixType
+
+#define MatrixType BlockSparsityPattern<3,3>
+matrix_functions_2;
+#undef MatrixType
+

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.