From: Wolfgang Bangerth Date: Fri, 22 Jun 2001 10:36:46 +0000 (+0000) Subject: New class CompressedSparsityPattern. Use it. New test. X-Git-Tag: v8.0.0~19045 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=de4565a0d977a35b79939ba2f7d0a4f0544155c9;p=dealii.git New class CompressedSparsityPattern. Use it. New test. git-svn-id: https://svn.dealii.org/trunk@4749 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index ea0a522b89..9cca92fac6 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -20,6 +20,7 @@ #include class SparsityPattern; +class CompressedSparsityPattern; class BlockSparsityPattern; template class SparseMatrix; template 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 diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index 7d537ca869..aca12523be 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -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 - static void make_sparsity_pattern (const DoFHandler &dof, - SparsityPattern &sparsity_pattern); + static + void + make_sparsity_pattern (const DoFHandler &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 - static void make_sparsity_pattern (const DoFHandler &dof, - const std::vector > &mask, - SparsityPattern &sparsity_pattern); - + template + static + void + make_sparsity_pattern (const DoFHandler &dof, + const std::vector > &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 + template static void make_boundary_sparsity_pattern (const DoFHandler &dof, const std::vector &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 static void make_boundary_sparsity_pattern (const DoFHandler<1> &dof, const std::vector &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 + template static void make_boundary_sparsity_pattern (const DoFHandler &dof, const typename FunctionMap::type &boundary_indicators, const std::vector &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 static void make_boundary_sparsity_pattern (const DoFHandler<1> &dof, const FunctionMap<1>::type &boundary_indicators, diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 9538c3407c..e3cc3f26bd 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -17,6 +17,7 @@ #include #include +#include #include #include #include @@ -423,60 +424,227 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const for (unsigned int row=0; row distribute(sparsity.n_rows(), -1); + + for (unsigned int c=0; c(c); + + const unsigned int n_rows = sparsity.n_rows(); + for (unsigned int row=0; row #include #include +#include #include #include @@ -68,6 +69,7 @@ DoFTools::make_sparsity_pattern (const DoFHandler &dof, } + template void DoFTools::make_sparsity_pattern (const DoFHandler &dof, @@ -115,8 +117,10 @@ DoFTools::make_sparsity_pattern (const DoFHandler &dof, } + #if deal_II_dimension == 1 +template 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 void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1> &dof_handler, const std::vector &dof_to_boundary_mapping, SparsityPattern &sparsity) @@ -179,7 +184,7 @@ void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1> &dof_handler, #endif -template +template void DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, const std::vector &dof_to_boundary_mapping, @@ -192,9 +197,19 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler& 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::const_iterator i=dof_to_boundary_mapping.begin(); + i!=dof_to_boundary_mapping.end(); ++i) + if ((*i != DoFHandler::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 dofs_on_this_face(dofs_per_face); @@ -216,12 +231,6 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler& 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::invalid_dof_index, - ExcInternalError()); - // make sparsity pattern for this cell for (unsigned int i=0; i& dof, }; -template + +template void DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, const typename FunctionMap::type &boundary_indicators, const std::vector &dof_to_boundary_mapping, @@ -270,11 +280,6 @@ void DoFTools::make_boundary_sparsity_pattern (const DoFHandler& 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::invalid_dof_index, - ExcInternalError()); // make sparsity pattern for this cell for (unsigned int i=0; i& dof, SparsityPattern &, const FullMatrix&, const FullMatrix&); +#endif + template void DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, const std::vector &, SparsityPattern &); template void +DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, + const std::vector &, + CompressedSparsityPattern &); +template void +DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, + const std::vector &, + BlockSparsityPattern &); +template void DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, const FunctionMap::type &boundary_indicators, const std::vector &dof_to_boundary_mapping, SparsityPattern &sparsity); -#endif +template void +DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + CompressedSparsityPattern &sparsity); +template void +DoFTools::make_boundary_sparsity_pattern (const DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + BlockSparsityPattern &sparsity); template void -DoFTools::make_sparsity_pattern (const DoFHandler& dof, +DoFTools::make_sparsity_pattern (const DoFHandler &dof, SparsityPattern &sparsity); +template void +DoFTools::make_sparsity_pattern (const DoFHandler &dof, + CompressedSparsityPattern &sparsity); + template void DoFTools::make_sparsity_pattern (const DoFHandler &dof, BlockSparsityPattern &sparsity); template void -DoFTools::make_sparsity_pattern (const DoFHandler& dof, +DoFTools::make_sparsity_pattern (const DoFHandler &dof, const std::vector > &mask, SparsityPattern &sparsity); template void -DoFTools::make_sparsity_pattern (const DoFHandler& dof, +DoFTools::make_sparsity_pattern (const DoFHandler &dof, + const std::vector > &mask, + CompressedSparsityPattern &sparsity); + +template void +DoFTools::make_sparsity_pattern (const DoFHandler &dof, const std::vector > &mask, BlockSparsityPattern &sparsity); + template void DoFTools::distribute_cell_to_dof_vector (const DoFHandler &dof_handler, diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 4f3d503521..2a35c76891 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -305,6 +305,16 @@ documentation, etc.

lac

    +
  1. + New: Class CompressedSparsityPattern + may be used as an intermediate form of the SparsityPattern class if memory + requirements are tight during construction of the sparsity + pattern. +
    + (WB 2001/06/22) +

    +
  2. New: There are now functions SparsityPattern::copy_from and @@ -441,6 +451,17 @@ documentation, etc.

    deal.II

      +
    1. + New: Implement + DoFTools::make_sparsity_pattern, + DoFTools::make_boundary_sparsity_pattern, + and + ConstraintMatrix::condense to work on + the CompressedSparsityPattern class. +
      + (WB 2001/06/22) +

      +
    2. New: FE*Values::get_quadrature 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 index 0000000000..c5212885e5 --- /dev/null +++ b/deal.II/lac/include/lac/compressed_sparsity_pattern.h @@ -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 +#include + +template class SparseMatrix; + +#include +#include + + +/** + * 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 > 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 diff --git a/deal.II/lac/include/lac/sparsity_pattern.h b/deal.II/lac/include/lac/sparsity_pattern.h index d4baa0cdf3..2940dbf3fc 100644 --- a/deal.II/lac/include/lac/sparsity_pattern.h +++ b/deal.II/lac/include/lac/sparsity_pattern.h @@ -18,6 +18,7 @@ #include template class SparseMatrix; +class CompressedSparsityPattern; #include #include @@ -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 index 0000000000..34f27f8f8f --- /dev/null +++ b/deal.II/lac/source/compressed_sparsity_pattern.cc @@ -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 + +#include +#include +#include +#include +#include +#include + +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 > 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::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::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(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::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::const_iterator i=column_indices[row].begin(); + i!=column_indices[row].end(); ++i) + if (static_cast(abs(static_cast(row-*i))) > b) + b = abs(static_cast(row-*i)); + + return b; +}; + + + +unsigned int +CompressedSparsityPattern::n_nonzero_elements () const +{ + unsigned int n=0; + for (unsigned int i=0; i +#include +#include #include #include @@ -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); }; }; diff --git a/tests/deal.II/Makefile.in b/tests/deal.II/Makefile.in index 761b1be944..4cb007ae05 100644 --- a/tests/deal.II/Makefile.in +++ b/tests/deal.II/Makefile.in @@ -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 ############################################################ diff --git a/tests/deal.II/matrices.cc b/tests/deal.II/matrices.cc index 19d526c181..97e2551f2e 100644 --- a/tests/deal.II/matrices.cc +++ b/tests/deal.II/matrices.cc @@ -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 index 0000000000..ec18708351 --- /dev/null +++ b/tests/deal.II/sparsity_pattern.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + +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 +void +check_boundary (const DoFHandler &dof) +{ + std::vector 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 +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 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 element(FE_Q(1), 1, + FE_Q(2), 1); + DoFHandler 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 > mask (2, std::vector(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 (); +}