#include <base/subscriptor.h>
class SparsityPattern;
+class CompressedSparsityPattern;
class BlockSparsityPattern;
template <typename number> class SparseMatrix;
template <typename number> class BlockSparseMatrix;
* 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
* 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}
* rather may want to give the
* following pattern of
* couplings:
- *
* @begin{verbatim}
* 1 0 1
* 0 1 1
* 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
* 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,
/**
* 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,
* 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,
#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>
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();
};
#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>
}
+
template <int dim, class SparsityPattern>
void
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,
+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)
#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,
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);
{
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)
};
-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,
{
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)
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,
<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
<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
--- /dev/null
+//---------------------------- 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
#include <base/subscriptor.h>
template <typename number> class SparseMatrix;
+class CompressedSparsityPattern;
#include <vector>
#include <iterator>
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
--- /dev/null
+//---------------------------- 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;
+};
//---------------------------- sparsity_pattern.cc ---------------------------
-#include <lac/sparse_matrix.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
#include <iostream>
#include <iomanip>
+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
{
// 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
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);
};
};
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
############################################################
-//---------------------------- derivative_approximation.cc ---------------------------
+//---------------------------- matrices.cc ---------------------------
// $Id$
// Version: $Name$
//
// 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 */
--- /dev/null
+//---------------------------- 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 ();
+}
projects / dealii-svn.git / commitdiff