]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
New class CompressedSparsityPattern. Use it. New test.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 22 Jun 2001 10:36:46 +0000 (10:36 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 22 Jun 2001 10:36:46 +0000 (10:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@4749 0785d39b-7218-0410-832d-ea1e28bc413d

12 files changed:
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/doc/news/2001/c-3-1.html
deal.II/lac/include/lac/compressed_sparsity_pattern.h [new file with mode: 0644]
deal.II/lac/include/lac/sparsity_pattern.h
deal.II/lac/source/compressed_sparsity_pattern.cc [new file with mode: 0644]
deal.II/lac/source/sparsity_pattern.cc
tests/deal.II/Makefile.in
tests/deal.II/matrices.cc
tests/deal.II/sparsity_pattern.cc [new file with mode: 0644]

index ea0a522b89b7ebf124c2b81cdb02abb1e378bf1f..9cca92fac6993fd2e231d4ac51c27f4318ab99b8 100644 (file)
@@ -20,6 +20,7 @@
 #include <base/subscriptor.h>
 
 class SparsityPattern;
+class CompressedSparsityPattern;
 class BlockSparsityPattern;
 template <typename number> class SparseMatrix;
 template <typename number> class BlockSparseMatrix;
@@ -282,6 +283,13 @@ class ConstraintMatrix : public Subscriptor
                                      * patterns.
                                      */
     void condense (BlockSparsityPattern &sparsity) const;
+
+                                    /**
+                                     * Same function as above, but
+                                     * condenses square compressed
+                                     * sparsity patterns.
+                                     */
+    void condense (CompressedSparsityPattern &sparsity) const;
     
                                     /**
                                      * Condense a given matrix. The associated
index 7d537ca86937cb39109123af51ce4cca8c129ddb..aca12523be09938974344e5fa9d53b4493a0bee9 100644 (file)
@@ -192,33 +192,44 @@ class DoFTools
                                      * The actual type of the
                                      * sparsity pattern may be
                                      * @ref{SparsityPattern},
-                                     * @ref{BlockSparsityPattern}, or any
-                                     * other class that satisfies
-                                     * similar requirements.
+                                     * @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.
                                      */
     template <int dim, class SparsityPattern>
-    static void make_sparsity_pattern (const DoFHandler<dim> &dof,
-                                      SparsityPattern       &sparsity_pattern);
+    static
+    void
+    make_sparsity_pattern (const DoFHandler<dim> &dof,
+                          SparsityPattern       &sparsity_pattern);
 
                                     /**
                                      * Locate non-zero entries for
-                                     * mixed methods.  This function
-                                     * does mostly the same as the
-                                     * other @p{make_sparsity_pattern},
-                                     * but it is specialized for
-                                     * mixed finite elements and
-                                     * allows to specify which
-                                     * variables couple in which
-                                     * equation. For example, if
-                                     * wanted to solve the Stokes
-                                     * equations,
-                                     *
-                                     *
+                                     * ector valued finite elements.
+                                     * This function does mostly the
+                                     * same as the previous
+                                     * @p{make_sparsity_pattern}, but
+                                     * it is specialized for vector
+                                     * finite elements and allows to
+                                     * specify which variables couple
+                                     * in which equation. For
+                                     * example, if wanted to solve
+                                     * the Stokes equations,
                                      * @begin{verbatim}
                                      * -\Delta \vec u + \nabla p = 0,
                                      * \div u                    = 0
                                      * @end{verbatim}
-                                     *
                                      * in two space dimensions,
                                      * using stable Q2/Q1 mixed
                                      * elements (using the @ref{FESystem}
@@ -228,7 +239,6 @@ class DoFTools
                                      * rather may want to give the
                                      * following pattern of
                                      * couplings:
-                                     *
                                      * @begin{verbatim}
                                      *   1 0 1
                                      *   0 1 1
@@ -263,15 +273,18 @@ class DoFTools
                                      * The actual type of the
                                      * sparsity pattern may be
                                      * @ref{SparsityPattern},
-                                     * @ref{BlockSparsityPattern}, or any
-                                     * other class that satisfies
+                                     * @ref{CompressedSparsityPattern},
+                                     * @ref{BlockSparsityPattern}, or
+                                     * any other class that satisfies
                                      * similar requirements.
                                      */
-    template<int dim, class SparsityPattern>
-    static void make_sparsity_pattern (const DoFHandler<dim>       &dof,
-                                      const std::vector<std::vector<bool> > &mask,
-                                      SparsityPattern             &sparsity_pattern);
-
+    template <int dim, class SparsityPattern>
+    static
+    void
+    make_sparsity_pattern (const DoFHandler<dim>                 &dof,
+                          const std::vector<std::vector<bool> > &mask,
+                          SparsityPattern                       &sparsity_pattern);
+    
                                     /**
                                      * Write the sparsity structure
                                      * of the matrix composed of the
@@ -293,12 +306,26 @@ class DoFTools
                                      * matrix yourself, using
                                      * @ref{SparsityPattern}@p{::compress()}.
                                      *
-                                     * Since this function is
-                                     * obviously useless in one
-                                     * spatial dimension, it is not
-                                     * implemented.
+                                     * The actual type of the
+                                     * 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.
                                      */
-    template<int dim>
+    template <int dim, class SparsityPattern>
     static void
     make_boundary_sparsity_pattern (const DoFHandler<dim>           &dof,
                                    const std::vector<unsigned int> &dof_to_boundary_mapping,
@@ -306,8 +333,13 @@ class DoFTools
 
                                     /**
                                      * Declaration of same function
-                                     * for different space dimension.
+                                     * for one space dimension. The
+                                     * resulting sparsity pattern
+                                     * couples all dofs of the
+                                     * specified end points to each
+                                     * other.
                                      */
+    template <class SparsityPattern>
     static void
     make_boundary_sparsity_pattern (const DoFHandler<1>             &dof,
                                    const std::vector<unsigned int> &dof_to_boundary_mapping,
@@ -340,22 +372,27 @@ class DoFTools
                                      * the function pointers to null
                                      * pointers.
                                      *
-                                     * Since this function is
-                                     * obviously useless in one
-                                     * spatial dimension, it is not
-                                     * implemented.
+                                     * For the type of the sparsity
+                                     * pattern, the same holds as
+                                     * said above.
                                      */
-    template<int dim>
+    template <int dim, class SparsityPattern>
     static void
     make_boundary_sparsity_pattern (const DoFHandler<dim> &dof,
                                    const typename FunctionMap<dim>::type &boundary_indicators,
                                    const std::vector<unsigned int>  &dof_to_boundary_mapping,
                                    SparsityPattern    &sparsity); 
 
+
                                     /**
                                      * Declaration of same function
-                                     * for different space dimension.
+                                     * for one space dimension. The
+                                     * resulting sparsity pattern
+                                     * couples all dofs of the
+                                     * specified end points to each
+                                     * other.
                                      */
+    template <class SparsityPattern>
     static void
     make_boundary_sparsity_pattern (const DoFHandler<1>             &dof,
                                    const FunctionMap<1>::type      &boundary_indicators,
index 9538c3407c4c95fad66425183105cde6c8371932..e3cc3f26bd8afd7b5b7b6e9ba6aab17ec44db2e1 100644 (file)
@@ -17,6 +17,7 @@
 
 #include <base/memory_consumption.h>
 #include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
 #include <lac/vector.h>
 #include <lac/block_vector.h>
 #include <lac/sparse_matrix.h>
@@ -423,60 +424,227 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const
   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];
-           
+       {
+                                          // regular line. loop over
+                                          // cols. note that this
+                                          // changes the line we are
+                                          // presently working on: we
+                                          // add additional
+                                          // entries. these are put to
+                                          // the end of the
+                                          // row. however, as
+                                          // constrained nodes cannot
+                                          // be constrained to other
+                                          // constrained nodes, nothing
+                                          // will happen if we run into
+                                          // these added nodes, as they
+                                          // can't be distributed
+                                          // further. we might store
+                                          // the position of the last
+                                          // old entry and stop work
+                                          // there, but since operating
+                                          // on the newly added ones
+                                          // only takes two comparisons
+                                          // (column index valid,
+                                          // distribute[column]
+                                          // necessarily ==-1), it is
+                                          // cheaper to not do so and
+                                          // run right until the end of
+                                          // the line
+         const unsigned int row_end = sparsity.get_rowstart_indices()[row+1];
+         for (unsigned int j=sparsity.get_rowstart_indices()[row]; j<row_end; ++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
+                                                    // @p{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. note that
+                                        // here the present row is
+                                        // not touched (unlike above)
+       {
+         const unsigned int row_end = sparsity.get_rowstart_indices()[row+1];
+         for (unsigned int j=sparsity.get_rowstart_indices()[row]; j<row_end; ++j)
                                             // end of row reached?
-           if (column == SparsityPattern::invalid_entry)
+           if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
              break;
            else
-             if (distribute[column] != -1)
-               {
-                                                  // distribute entry at regular
-                                                  // row @p{row} and irregular column
+             {
+               if (distribute[sparsity.get_column_numbers()[j]] == -1)
+                                                  // distribute entry at irregular
+                                                  // row @p{row} and regular 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);
-               };
+                      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 @p{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();
+};
+
+
+
+void ConstraintMatrix::condense (CompressedSparsityPattern &sparsity) const
+{
+  Assert (sorted == true, ExcMatrixNotClosed());
+  Assert (sparsity.n_rows() == sparsity.n_cols(),
+         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
+  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)
+    {
+      if (distribute[row] == -1)
+                                        // regular line. loop over
+                                        // cols. note that as we
+                                        // proceed to distribute
+                                        // cols, the loop may get
+                                        // longer
+       for (unsigned int j=0; j<sparsity.row_length(row); ++j)
+         {
+           const unsigned int column = sparsity.column_number(row,j);
+
+           if (distribute[column] != -1)
+             {
+                                                // distribute entry
+                                                // at regular row
+                                                // @p{row} and
+                                                // irregular column
+                                                // column. note that
+                                                // this changes the
+                                                // line we are
+                                                // presently working
+                                                // on: we add
+                                                // additional
+                                                // entries. if we add
+                                                // another entry at a
+                                                // column behind the
+                                                // present one, we
+                                                // will encounter it
+                                                // later on (but
+                                                // since it can't be
+                                                // further
+                                                // constrained, won't
+                                                // have to do
+                                                // anything about
+                                                // it). if we add it
+                                                // up front of the
+                                                // present column, we
+                                                // will find the
+                                                // present column
+                                                // later on again as
+                                                // it was shifted
+                                                // back (again
+                                                // nothing happens,
+                                                // in particular no
+                                                // endless loop, as
+                                                // when we encounter
+                                                // it the second time
+                                                // we won't be able
+                                                // to add more
+                                                // entries as they
+                                                // all already exist,
+                                                // but we do the same
+                                                // work more often
+                                                // than necessary,
+                                                // and the loop gets
+                                                // longer), so move
+                                                // the cursor one to
+                                                // the right in the
+                                                // case that we add
+                                                // an entry up front
+                                                // that did not exist
+                                                // before. check
+                                                // whether it existed
+                                                // before by tracking
+                                                // the length of this
+                                                // row
+               unsigned int old_rowlength = sparsity.row_length(row);
+               for (unsigned int q=0;
+                    q!=lines[distribute[column]].entries.size();
+                    ++q) 
+                 {
+                   const unsigned int
+                     new_col = lines[distribute[column]].entries[q].first;
+                   
+                   sparsity.add (row, new_col);
+
+                   const unsigned int new_rowlength = sparsity.row_length(row);
+                   if ((new_col < column) && (old_rowlength != new_rowlength))
+                     ++j;
+                   old_rowlength = new_rowlength;
+                 };
+             };
          }
       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 @p{row} and regular column
-                                                // sparsity.colnums[j]
+       for (unsigned int j=0; j<sparsity.row_length(row); ++j)
+         {
+           const unsigned int column = sparsity.column_number(row,j);
+
+           if (distribute[column] == -1)
+                                              // distribute entry at irregular
+                                              // row @p{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,
+                             column);
+           else
+                                              // distribute entry at irregular
+                                              // row @p{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[row]].entries.size(); ++q) 
-                 sparsity.add (lines[distribute[row]].entries[q].first,
-                               sparsity.get_column_numbers()[j]);
-             else
-                                                // distribute entry at irregular
-                                                // row @p{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);
-           };
+                    q!=lines[distribute[sparsity.column_number(row,j)]]
+                                   .entries.size(); ++q)
+                 sparsity.add (lines[distribute[row]].entries[p].first,
+                               lines[distribute[sparsity.column_number(row,j)]]
+                               .entries[q].first);
+         };
     };
-  
-  sparsity.compress();
 };
 
 
index 71268fd57a132aa01d215ac47023f9f057a5b99e..4ace8c073d5947988661895d0ac86932ed537fb5 100644 (file)
@@ -27,6 +27,7 @@
 #include <fe/fe_values.h>
 #include <dofs/dof_tools.h>
 #include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
 #include <lac/block_sparsity_pattern.h>
 #include <lac/vector.h>
 
@@ -68,6 +69,7 @@ DoFTools::make_sparsity_pattern (const DoFHandler<dim> &dof,
 }
 
 
+
 template <int dim, class SparsityPattern>
 void
 DoFTools::make_sparsity_pattern (const DoFHandler<dim>                 &dof,
@@ -115,8 +117,10 @@ DoFTools::make_sparsity_pattern (const DoFHandler<dim>                 &dof,
 }
 
 
+
 #if deal_II_dimension == 1
 
+template <class SparsityPattern>
 void
 DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1>             &dof_handler,
                                          const FunctionMap<1>::type      &function_map,
@@ -160,6 +164,7 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1>             &dof_h
 
 
 
+template <class SparsityPattern>
 void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1> &dof_handler,
                                               const std::vector<unsigned int>  &dof_to_boundary_mapping,
                                               SparsityPattern    &sparsity)
@@ -179,7 +184,7 @@ void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1> &dof_handler,
 #endif
 
 
-template <int dim>
+template <int dim, class SparsityPattern>
 void
 DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
                                          const std::vector<unsigned int>  &dof_to_boundary_mapping,
@@ -192,9 +197,19 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
          ExcDimensionMismatch (sparsity.n_rows(), dof.n_boundary_dofs()));
   Assert (sparsity.n_cols() == dof.n_boundary_dofs(),
          ExcDimensionMismatch (sparsity.n_cols(), dof.n_boundary_dofs()));
-  Assert (*max_element(dof_to_boundary_mapping.begin(),
-                      dof_to_boundary_mapping.end()) == sparsity.n_rows()-1,
-         ExcInternalError());
+#ifdef DEBUG
+  if (true)
+    {
+      unsigned int max_element = 0;
+      for (std::vector<unsigned int>::const_iterator i=dof_to_boundary_mapping.begin();
+          i!=dof_to_boundary_mapping.end(); ++i)
+       if ((*i != DoFHandler<dim>::invalid_dof_index) &&
+           (*i > max_element))
+         max_element = *i;
+      Assert (max_element  == sparsity.n_rows()-1,
+             ExcInternalError());
+    };
+#endif
 
   const unsigned int dofs_per_face = dof.get_fe().dofs_per_face;
   std::vector<unsigned int> dofs_on_this_face(dofs_per_face);
@@ -216,12 +231,6 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
       {
        face->get_dof_indices (dofs_on_this_face);
 
-                                        // make sure all dof indices have a
-                                        // boundary index
-       Assert (*min_element(dofs_on_this_face.begin(),
-                            dofs_on_this_face.end()) != DoFHandler<dim>::invalid_dof_index,
-               ExcInternalError());
-       
                                         // make sparsity pattern for this cell
        for (unsigned int i=0; i<dofs_per_face; ++i)
          for (unsigned int j=0; j<dofs_per_face; ++j) 
@@ -231,7 +240,8 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
 };
 
 
-template <int dim>
+
+template <int dim, class SparsityPattern>
 void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
                                               const typename FunctionMap<dim>::type  &boundary_indicators,
                                               const std::vector<unsigned int>  &dof_to_boundary_mapping,
@@ -270,11 +280,6 @@ void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
       {
        face->get_dof_indices (dofs_on_this_face);
 
-                                        // make sure all dof indices have a
-                                        // boundary index
-       Assert (*min_element(dofs_on_this_face.begin(),
-                            dofs_on_this_face.end()) != DoFHandler<dim>::invalid_dof_index,
-               ExcInternalError());
                                         // make sparsity pattern for this cell
        for (unsigned int i=0; i<dofs_per_face; ++i)
          for (unsigned int j=0; j<dofs_per_face; ++j)
@@ -2163,35 +2168,64 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
                                      SparsityPattern    &,
                                      const FullMatrix<double>&,
                                      const FullMatrix<double>&);
+#endif
+
 template void
 DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
                                          const std::vector<unsigned int>  &,
                                          SparsityPattern    &);
 template void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+                                         const std::vector<unsigned int>  &,
+                                         CompressedSparsityPattern    &);
+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 FunctionMap<deal_II_dimension>::type  &boundary_indicators,
                                          const std::vector<unsigned int>  &dof_to_boundary_mapping,
                                          SparsityPattern    &sparsity);
-#endif
+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,
+                                         CompressedSparsityPattern    &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,
+                                         BlockSparsityPattern    &sparsity);
 
 template void
-DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>dof,
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
                                 SparsityPattern    &sparsity);
 
+template void
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
+                                CompressedSparsityPattern    &sparsity);
+
 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,
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
                                 const std::vector<std::vector<bool> > &mask,
                                 SparsityPattern             &sparsity);
 
 template void 
-DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
+                                const std::vector<std::vector<bool> > &mask,
+                                CompressedSparsityPattern             &sparsity);
+
+template void 
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
                                 const std::vector<std::vector<bool> > &mask,
                                 BlockSparsityPattern        &sparsity);
 
+
 template
 void
 DoFTools::distribute_cell_to_dof_vector (const DoFHandler<deal_II_dimension> &dof_handler,
index 4f3d503521015abc29f70eefe34c2b219d6cf333..2a35c76891bcd6822a9e844d5e0ff52c3e0dce4e 100644 (file)
@@ -305,6 +305,16 @@ documentation, etc</a>.
 <h3>lac</h3>
 
 <ol>
+  <li> <p>
+       New: Class <code class="class">CompressedSparsityPattern</code>
+       may be used as an intermediate form of the <code
+       class="class">SparsityPattern</code> class if memory
+       requirements are tight during construction of the sparsity
+       pattern.
+       <br>
+       (WB 2001/06/22)
+       </p>
+
   <li> <p>
        New: There are now functions 
        <code class="member">SparsityPattern::copy_from</code> and
@@ -441,6 +451,17 @@ documentation, etc</a>.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       New: Implement
+       <code class="member">DoFTools::make_sparsity_pattern</class>,
+       <code class="member">DoFTools::make_boundary_sparsity_pattern</class>,
+       and 
+       <code class="member">ConstraintMatrix::condense</class> to work on
+       the <code class="class">CompressedSparsityPattern</code> class.
+       <br>
+       (WB 2001/06/22)
+       </p>
+
   <li> <p>
        New: <code class="member">FE*Values::get_quadrature</class>
        returns a reference to the quadrature formula used by a
diff --git a/deal.II/lac/include/lac/compressed_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_sparsity_pattern.h
new file mode 100644 (file)
index 0000000..c521288
--- /dev/null
@@ -0,0 +1,366 @@
+//----------------------------  compressed_sparsity_pattern.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  compressed_sparsity_pattern.h  ---------------------------
+#ifndef __deal2__compressed_sparsity_pattern_h
+#define __deal2__compressed_sparsity_pattern_h
+
+
+#include <base/exceptions.h>
+#include <base/subscriptor.h>
+
+template <typename number> class SparseMatrix;
+
+#include <vector>
+#include <set>
+
+
+/**
+ * This class acts as an intermediate form of the
+ * @ref{SparsityPattern} class. From the interface it mostly
+ * represents a @ref{SparsityPattern} object that is kept compressed
+ * at all times. However, since the final sparsity pattern is not
+ * known while constructing it, keeping the pattern compressed at all
+ * times can only be achieved at the expense of either increased
+ * memory or run time consumption upon use. The main purpose of this
+ * class is to avoid some memory bottlenecks, so we chose to implement
+ * it memory conservative, but the chosen data format is too unsuited
+ * to be used for actual matrices. It is therefore necessary to first
+ * copy the data of this object over to an object of type
+ * @ref{SparsityPattern} before using it in actual matrices.
+ *
+ * Another viewpoint is that this class does not need up front
+ * allocation of a certain amount of memory, but grows as necessary.
+ *
+ *
+ * @sect3{Rationale}
+ *
+ * When constructing the sparsity pattern of a matrix, you usually
+ * first have to provide an empty sparsity pattern object with a fixed
+ * maximal number of entries per row. To find out about this maximal
+ * row length, one usually calls the function
+ * @ref{DoFHandler}@p{::max_couplings_per_dof} which returns an
+ * estimate for that quantity. While this estimate is usually quite
+ * good in 2d and exact in 1d, it is often significantly too large in
+ * 3d and especially for higher order elements. Furthermore, normally
+ * only a small fraction of the rows of a matrix will end up having
+ * the maximal number of nonzero entries per row (usually those nodes
+ * adjacent to hanging nodes), most have much less. In effect, the
+ * empty @ref{SparsityPattern} object has allocated much too much
+ * memory. Although this unnecessarily allocated memory is later freed
+ * when @ref{SparsityPattern}@p{::compress} is called, this
+ * overallocation has, with higher order elements and in 3d, sometimes
+ * been so large that the program aborted due to lack of memory.
+ *
+ * This class therefore provides an alternative representation of a
+ * sparsity pattern: we don't specify a maximal row length initially,
+ * but store a set of column indices indicating possible nonzero
+ * entries in the sparsity pattern for each row. This is very much
+ * like the final "compressed" format used in the
+ * @ref{SparsityPattern} object after compression, but uses a less
+ * compact memory storage format, since the exact number of entries
+ * per row is only known a posteriori and since it may change (for the
+ * @ref{SparsityPattern} class, no more changes are allowed after
+ * compressing it). We can therefore not store all the column indices
+ * in a big array, but have to use a vector of sets. This can later be
+ * used to actually initialize a @ref{SparsityPattern} object with the
+ * then final set of necessary indices.
+ *
+ *
+ * @sect3{Interface}
+ *
+ * Since this class is intended as an intermediate replacement of the
+ * @ref{SparsityPattern} class, it has mostly the same interface, with
+ * small changes where necessary. In particular, the @ref{add}
+ * function, and the functions inquiring properties of the sparsity
+ * pattern are the same.
+ *
+ *
+ * @sect3{Usage}
+ *
+ * Use this class as follows:
+ * @begin{verbatim}
+ * CompressedSparsityPattern compressed_pattern (dof_handler.n_dofs());
+ * DoFTools::make_sparsity_pattern (dof_handler,
+ *                                  compressed_pattern);
+ * constraints.condense (compressed_pattern);
+ *
+ * SparsityPattern sp;
+ * sp.copy_from (compressed_pattern);
+ * @end{verbatim}
+ *
+ *
+ * @author Wolfgang Bangerth, 2001
+ */
+class CompressedSparsityPattern : public Subscriptor
+{
+  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.
+                                     */
+    CompressedSparsityPattern ();
+    
+                                    /**
+                                     * Copy constructor. This constructor is
+                                     * only allowed to be called if the matrix
+                                     * structure to be copied is empty. This is
+                                     * so in order to prevent involuntary
+                                     * copies of objects for temporaries, which
+                                     * can use large amounts of computing time.
+                                     * However, copy constructors are needed
+                                     * if yo want to use the STL data types
+                                     * on classes like this, e.g. to write
+                                     * such statements like
+                                     * @p{v.push_back (CompressedSparsityPattern());},
+                                     * with @p{v} a vector of @p{CompressedSparsityPattern}
+                                     * objects.
+                                     *
+                                     * Usually, it is sufficient to
+                                     * use the explicit keyword to
+                                     * disallow unwanted temporaries,
+                                     * but for the STL vectors, this
+                                     * does not work. Since copying a
+                                     * structure like this is not
+                                     * useful anyway because multiple
+                                     * matrices can use the same
+                                     * sparsity structure, copies are
+                                     * only allowed for empty
+                                     * objects, as described above.
+                                     */
+    CompressedSparsityPattern (const CompressedSparsityPattern &);
+
+                                    /**
+                                     * Initialize a rectangular
+                                     * matrix with @p{m} rows and
+                                     * @p{n} columns.
+                                     */
+    CompressedSparsityPattern (const unsigned int m,
+                              const unsigned int n);
+    
+                                    /**
+                                     * Initialize a square matrix of
+                                     * dimension @p{n}.
+                                     */
+    CompressedSparsityPattern (const unsigned int n);
+
+                                    /**
+                                     * Copy operator. For this the
+                                     * same holds as for the copy
+                                     * constructor: it is declared,
+                                     * defined and fine to be called,
+                                     * but the latter only for empty
+                                     * objects.
+                                     */
+    CompressedSparsityPattern & operator = (const CompressedSparsityPattern &);
+    
+                                    /**
+                                     * Reallocate memory and set up
+                                     * data structures for a new
+                                     * matrix with @p{m} rows and
+                                     * @p{n} columns, with at most
+                                     * @p{max_per_row} nonzero
+                                     * entries per row.
+                                     */
+    void reinit (const unsigned int m,
+                const unsigned int n);
+    
+                                    /**
+                                     * Since this object is kept
+                                     * compressed at all times anway,
+                                     * this function does nothing,
+                                     * but is declared to make the
+                                     * interface of this class as
+                                     * much alike as that of the
+                                     * @ref{SparsityPattern} class.
+                                     */
+    void compress ();
+    
+                                    /**
+                                     * Return whether the object is
+                                     * empty. It is empty if no
+                                     * memory is allocated, which is
+                                     * the same as that both
+                                     * dimensions are zero.
+                                     */
+    bool empty () const;
+
+                                    /**
+                                     * Return the maximum number of
+                                     * entries per row. Note that
+                                     * this number may change as
+                                     * entries are added.
+                                     */
+    unsigned int max_entries_per_row () const;
+
+                                    /**
+                                     * Add a nonzero entry to the
+                                     * matrix. If the entry already
+                                     * exists, nothing bad happens.
+                                     */
+    void add (const unsigned int i, 
+             const unsigned int j);
+    
+                                    /**
+                                     * Make the sparsity pattern
+                                     * symmetric by adding the
+                                     * sparsity pattern of the
+                                     * transpose object.
+                                     *
+                                     * This function throws an
+                                     * exception if the sparsity
+                                     * pattern does not represent a
+                                     * square matrix.
+                                     */
+    void symmetrize ();
+    
+                                    /**
+                                     * Print the sparsity of the matrix
+                                     * in a format that @p{gnuplot} understands
+                                     * and which can be used to plot the
+                                     * sparsity pattern in a graphical
+                                     * way. The format consists of pairs
+                                     * @p{i j} of nonzero elements, each
+                                     * representing one entry of this
+                                     * matrix, one per line of the output
+                                     * file. Indices are counted from
+                                     * zero on, as usual. Since sparsity
+                                     * patterns are printed in the same
+                                     * way as matrices are displayed, we
+                                     * print the negative of the column
+                                     * index, which means that the
+                                     * @p{(0,0)} element is in the top left
+                                     * rather than in the bottom left
+                                     * corner.
+                                     *
+                                     * Print the sparsity pattern in
+                                     * gnuplot by setting the data style
+                                     * to dots or points and use the
+                                     * @p{plot} command.
+                                     */
+    void print_gnuplot (std::ostream &out) const;
+
+                                    /**
+                                     * Return number of rows of this
+                                     * matrix, which equals the dimension
+                                     * of the image space.
+                                     */
+    unsigned int n_rows () const;
+
+                                    /**
+                                     * Return number of columns of this
+                                     * matrix, which equals the dimension
+                                     * of the range space.
+                                     */
+    unsigned int n_cols () const;
+
+                                    /**
+                                     * Number of entries in a specific row.
+                                     */
+    unsigned int row_length (const unsigned int row) const;
+
+                                    /**
+                                     * Access to column number field.
+                                     * Return the column number of
+                                     * the @p{index}th entry in @p{row}.
+                                     */
+    unsigned int column_number (const unsigned int row,
+                               const unsigned int index) const;
+
+                                    /**
+                                     * Compute the bandwidth of the matrix
+                                     * represented by this structure. The
+                                     * bandwidth is the maximum of
+                                     * $|i-j|$ for which the index pair
+                                     * $(i,j)$ represents a nonzero entry
+                                     * of the matrix.
+                                     */
+    unsigned int bandwidth () const;
+
+                                    /**
+                                     * Return the number of nonzero elements of
+                                     * this matrix. Actually, it returns the
+                                     * number of entries in the sparsity
+                                     * pattern; if any of the entries should
+                                     * happen to be zero, it is counted
+                                     * anyway.
+                                     *
+                                     * This function may only be called if the
+                                     * matrix struct is compressed. It does not
+                                     * make too much sense otherwise anyway.
+                                     */
+    unsigned int n_nonzero_elements () const;
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidIndex,
+                   int, int,
+                   << "The given index " << arg1
+                   << " should be less than " << arg2 << ".");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcInvalidConstructorCall);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNotSquare);
+    
+  private:
+                                    /**
+                                     * Number of rows that this sparsity
+                                     * structure shall represent.
+                                     */
+    unsigned int rows;
+
+                                    /**
+                                     * Number of columns that this sparsity
+                                     * structure shall represent.
+                                     */
+    unsigned int cols;
+
+                                    /**
+                                     * Actual data: store for each
+                                     * row the set of nonzero
+                                     * entries.
+                                     */
+    std::vector<std::set<unsigned int> > column_indices;
+
+    friend class SparsityPattern;
+};
+
+
+/*---------------------- Inline functions -----------------------------------*/
+
+
+inline
+unsigned int
+CompressedSparsityPattern::n_rows () const
+{
+  return rows;
+};
+
+
+inline
+unsigned int
+CompressedSparsityPattern::n_cols () const
+{
+  return cols;
+};
+
+#endif
index d4baa0cdf3c21d890a37714a9e7f3852318f1a69..2940dbf3fcb907aca0fcce41bfda5f7e5fdaf486 100644 (file)
@@ -18,6 +18,7 @@
 #include <base/subscriptor.h>
 
 template <typename number> class SparseMatrix;
+class CompressedSparsityPattern;
 
 #include <vector>
 #include <iterator>
@@ -460,6 +461,15 @@ class SparsityPattern : public Subscriptor
                    const unsigned int    n_cols,
                    const ForwardIterator begin,
                    const ForwardIterator end);
+
+                                    /**
+                                     * Copy data from an object of
+                                     * type
+                                     * @ref{CompressedSparsityPattern}. Previous
+                                     * content of this object is
+                                     * lost.
+                                     */
+    void copy_from (const CompressedSparsityPattern &csp);
     
                                     /**
                                      * Return whether the object is empty. It
diff --git a/deal.II/lac/source/compressed_sparsity_pattern.cc b/deal.II/lac/source/compressed_sparsity_pattern.cc
new file mode 100644 (file)
index 0000000..34f27f8
--- /dev/null
@@ -0,0 +1,219 @@
+//----------------------------  compressed_sparsity_pattern.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  compressed_sparsity_pattern.cc  ---------------------------
+
+
+#include <lac/compressed_sparsity_pattern.h>
+
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+#include <functional>
+
+CompressedSparsityPattern::CompressedSparsityPattern () :
+               rows(0),
+               cols(0)
+{};
+
+
+
+CompressedSparsityPattern::CompressedSparsityPattern (const CompressedSparsityPattern &s) :
+               Subscriptor(),
+               rows(0),
+               cols(0)
+{
+  Assert (s.rows == 0, ExcInvalidConstructorCall());
+  Assert (s.cols == 0, ExcInvalidConstructorCall());
+};
+
+
+
+CompressedSparsityPattern::CompressedSparsityPattern (const unsigned int m,
+                                                     const unsigned int n) 
+               : rows(0),
+                 cols(0)
+{
+  reinit (m,n);
+};
+
+
+
+CompressedSparsityPattern::CompressedSparsityPattern (const unsigned int n)
+               : rows(0),
+                 cols(0)
+{
+  reinit (n,n);
+};
+
+
+
+CompressedSparsityPattern &
+CompressedSparsityPattern::operator = (const CompressedSparsityPattern &s)
+{
+  Assert (s.rows == 0, ExcInvalidConstructorCall());
+  Assert (s.cols == 0, ExcInvalidConstructorCall());
+
+  Assert (rows == 0, ExcInvalidConstructorCall());
+  Assert (cols == 0, ExcInvalidConstructorCall());
+
+  return *this;
+};
+
+
+
+void
+CompressedSparsityPattern::reinit (const unsigned int m,
+                                  const unsigned int n)
+{
+  rows = m;
+  cols = n;
+
+  std::vector<std::set<unsigned int> > new_column_indices (rows);
+  column_indices.swap (new_column_indices);
+}
+
+
+
+void
+CompressedSparsityPattern::compress ()
+{};
+
+
+
+bool
+CompressedSparsityPattern::empty () const
+{
+  return ((rows==0) && (cols==0));
+};
+
+
+
+unsigned int
+CompressedSparsityPattern::max_entries_per_row () const 
+{
+  unsigned int m = 0;
+  for (unsigned int i=1; i<rows; ++i)
+    m = std::max (m, column_indices[i].size());
+
+  return m;
+};
+
+
+
+void
+CompressedSparsityPattern::add (const unsigned int i,
+                               const unsigned int j)
+{
+  Assert (i<rows, ExcInvalidIndex(i,rows));
+  Assert (j<cols, ExcInvalidIndex(j,cols));
+
+                                  // the std::set automatically
+                                  // assures uniqueness and
+                                  // sortedness of the column indices
+  column_indices[i].insert (j);
+};
+
+
+
+void
+CompressedSparsityPattern::symmetrize () 
+{
+  Assert (rows==cols, ExcNotSquare());
+
+                                  // loop over all elements presently
+                                  // in the sparsity pattern and add
+                                  // the transpose element. note:
+                                  //
+                                  // 1. that the sparsity pattern
+                                  // changes which we work on, but
+                                  // not the present row
+                                  //
+                                  // 2. that the @p{add} function can
+                                  // be called on elements that
+                                  // already exist without any harm
+  for (unsigned int row=0; row<rows; ++row)
+    for (std::set<unsigned int>::const_iterator i=column_indices[row].begin();
+        i!=column_indices[row].end(); ++i)
+                                      // add the transpose entry if
+                                      // this is not the diagonal
+      if (row != *i)
+       add (*i, row);
+};
+
+
+
+void
+CompressedSparsityPattern::print_gnuplot (std::ostream &out) const
+{
+  for (unsigned int row=0; row<rows; ++row)
+    for (std::set<unsigned int>::const_iterator i=column_indices[row].begin();
+        i!=column_indices[row].end(); ++i)
+                                      // while matrix entries are
+                                      // usually written (i,j),
+                                      // with i vertical and j
+                                      // horizontal, gnuplot output
+                                      // is x-y, that is we have to
+                                      // exchange the order of
+                                      // output
+      out << *i << " " << -static_cast<signed int>(row) << std::endl;
+
+  AssertThrow (out, ExcIO());
+}
+
+
+
+unsigned int
+CompressedSparsityPattern::row_length (const unsigned int row) const
+{
+  return column_indices[row].size();
+};
+
+
+
+unsigned int
+CompressedSparsityPattern::column_number (const unsigned int row,
+                                         const unsigned int index) const
+{
+  Assert (index < column_indices[row].size(),
+         ExcIndexRange (index, 0, column_indices[row].size()));
+  std::set<unsigned int>::const_iterator p = column_indices[row].begin();
+  std::advance (p, index);
+  return *p;
+};
+
+
+
+unsigned int
+CompressedSparsityPattern::bandwidth () const
+{
+  unsigned int b=0;
+  for (unsigned int row=0; row<rows; ++row)
+    for (std::set<unsigned int>::const_iterator i=column_indices[row].begin();
+        i!=column_indices[row].end(); ++i)
+      if (static_cast<unsigned int>(abs(static_cast<int>(row-*i))) > b)
+       b = abs(static_cast<signed int>(row-*i));
+
+  return b;
+};
+
+
+
+unsigned int
+CompressedSparsityPattern::n_nonzero_elements () const
+{
+  unsigned int n=0;
+  for (unsigned int i=0; i<rows; ++i)
+    n += column_indices[i].size();
+  return n;
+};
index d7fba38f51fec5085f252b34514ec8098006e01e..6a1d1512d2dcd7185fb157f4bb0c892eb99c5dbd 100644 (file)
@@ -12,7 +12,8 @@
 //----------------------------  sparsity_pattern.cc  ---------------------------
 
 
-#include <lac/sparse_matrix.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
 
 #include <iostream>
 #include <iomanip>
@@ -447,6 +448,16 @@ SparsityPattern::compress ()
 
 
 
+void
+SparsityPattern::copy_from (const CompressedSparsityPattern &csp) 
+{
+  copy_from (csp.n_rows(), csp.n_cols(),
+            csp.column_indices.begin(),
+            csp.column_indices.end());
+};
+
+
+
 bool
 SparsityPattern::empty () const
 {
@@ -583,7 +594,8 @@ SparsityPattern::symmetrize ()
                                   // the transpose element. note:
                                   //
                                   // 1. that the sparsity pattern
-                                  // changes which we work on
+                                  // changes which we work on, but
+                                  // not the present row
                                   //
                                   // 2. that the @p{add} function can
                                   // be called on elements that
@@ -599,8 +611,12 @@ SparsityPattern::symmetrize ()
          break;
 
                                         // otherwise add the
-                                        // transpose entry
-       add (colnums[k], row);
+                                        // transpose entry if this is
+                                        // not the diagonal (that
+                                        // would not harm, only take
+                                        // time to check up)
+       if (colnums[k] != row)
+         add (colnums[k], row);
       };
 };
 
index 761b1be9449816b113a277fde2167c901d11c2e5..4cb007ae05be9d0594db696c30f244e957e08ede 100644 (file)
@@ -42,12 +42,13 @@ wave-test-3.exe              : wave-test-3.go                        $(lib-2d)
 support_point_map.exe        : support_point_map.go        $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 filtered_matrix.exe          : filtered_matrix.go          $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 boundaries.exe               : boundaries.go               $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
+sparsity_pattern.exe         : sparsity_pattern.go         $(lib-1d) $(lib-2d) $(lib-1d) $(libraries)
 
 tests = grid_test grid_transform dof_test data_out derivatives gradients constraints mg \
         mglocal block_matrices second_derivatives derivative_approximation \
         matrices error_estimator intergrid_constraints intergrid_map \
         wave-test-3 dof_renumbering support_point_map filtered_matrix \
-       boundaries
+       boundaries sparsity_pattern
 
 ############################################################
 
index 19d526c181570f7edd1255d03702d40a3a5770ec..97e2551f2e9c3da42f2ef79ddc868064f4e174f0 100644 (file)
@@ -1,4 +1,4 @@
-//----------------------------  derivative_approximation.cc  ---------------------------
+//----------------------------  matrices.cc  ---------------------------
 //    $Id$
 //    Version: $Name$
 //
@@ -9,7 +9,7 @@
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  derivative_approximation.cc  ---------------------------
+//----------------------------  matrices.cc  ---------------------------
 
 
 /* Author: Wolfgang Bangerth, University of Heidelberg, 2001 */
diff --git a/tests/deal.II/sparsity_pattern.cc b/tests/deal.II/sparsity_pattern.cc
new file mode 100644 (file)
index 0000000..ec18708
--- /dev/null
@@ -0,0 +1,175 @@
+//----------------------------  sparsity_pattern.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparsity_pattern.cc  ---------------------------
+
+
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2001 */
+
+// check that the direct generation of the sparsity pattern and that
+// via the CompressedSparsityPattern result in the same
+
+
+#include <base/logstream.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/fe_q.h>
+#include <fe/fe_system.h>
+
+#include <fstream>
+
+
+
+bool operator == (const SparsityPattern &sp1,
+                 const SparsityPattern &sp2)
+{
+  if (sp1.n_nonzero_elements() != sp2.n_nonzero_elements())
+    return false;
+  
+  for (unsigned int i=0; i<sp1.n_nonzero_elements(); ++i)
+    if (sp1.get_column_numbers()[i] !=
+       sp2.get_column_numbers()[i])
+      return false;
+
+  for (unsigned int i=0; i<sp1.n_rows(); ++i)
+    if (sp1.get_rowstart_indices()[i] !=
+       sp2.get_rowstart_indices()[i])
+      return false;
+  
+  return true;
+};
+
+
+
+template <int dim>
+void
+check_boundary (const DoFHandler<dim> &dof)
+{
+  std::vector<unsigned int> dof_to_boundary_mapping;
+  DoFTools::map_dof_to_boundary_indices (dof,
+                                        dof_to_boundary_mapping);
+
+                                  // first way: direct generation
+  SparsityPattern sparsity_1(dof.n_boundary_dofs(),
+                            dof.max_couplings_between_boundary_dofs());
+  DoFTools::make_boundary_sparsity_pattern (dof,
+                                           dof_to_boundary_mapping,
+                                           sparsity_1);
+  sparsity_1.compress ();
+
+                                  // second way: via a CompressedSparsityPattern
+  SparsityPattern sparsity_2;
+  CompressedSparsityPattern csp(dof.n_boundary_dofs());
+  DoFTools::make_boundary_sparsity_pattern (dof,
+                                           dof_to_boundary_mapping,
+                                           csp);
+  sparsity_2.copy_from (csp);
+
+                                  // the exact content of sparsity
+                                  // patterns is checked in other
+                                  // tests, so only make sure that
+                                  // sparsity_[12] are equal
+  deallog << __PRETTY_FUNCTION__
+         << " -- "
+         << (sparsity_1 == sparsity_2 ? "ok" : "failed")
+         << std::endl;
+};
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;  
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  tr.begin_active(2)->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FESystem<dim> element(FE_Q<dim>(1), 1,
+                       FE_Q<dim>(2), 1);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  std::vector<std::vector<bool> > mask (2, std::vector<bool>(2, false));
+  mask[0][0] = mask[1][1] = true;
+
+                                  // 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
+  SparsityPattern sparsity_2;
+  CompressedSparsityPattern csp (dof.n_dofs());
+  DoFTools::make_sparsity_pattern (dof, mask, csp);
+  constraints.condense (csp);
+  sparsity_2.copy_from (csp);
+
+
+                                  // the exact content of sparsity
+                                  // patterns is checked in other
+                                  // tests, so only make sure that
+                                  // sparsity_[12] are equal
+  deallog << __PRETTY_FUNCTION__
+         << " -- "
+         << (sparsity_1 == sparsity_2 ? "ok" : "failed")
+         << std::endl;
+  
+  check_boundary (dof);
+};
+
+
+
+int main ()
+{
+  std::ofstream logfile ("sparsity_pattern.output");
+  logfile.precision (2);
+  logfile.setf(std::ios::fixed);  
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}

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.