From 60775d6106bb70e95fd489dd1147b6cd27ac3fb2 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 30 Sep 2009 15:01:07 +0000 Subject: [PATCH] Change ConstraintMatrix implementation so that it has fast random access to ConstraintLines. git-svn-id: https://svn.dealii.org/trunk@19618 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/dofs/dof_accessor.h | 128 +++++ .../include/dofs/dof_accessor.templates.h | 433 ++++++++++++--- deal.II/doc/news/changes.h | 12 + deal.II/lac/include/lac/constraint_matrix.h | 493 ++++++++++++------ .../include/lac/constraint_matrix.templates.h | 118 ++--- deal.II/lac/source/constraint_matrix.cc | 120 ++--- 6 files changed, 917 insertions(+), 387 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_accessor.h b/deal.II/deal.II/include/dofs/dof_accessor.h index f3151108ea..36cacee7ed 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.h @@ -26,6 +26,7 @@ DEAL_II_NAMESPACE_OPEN template class FullMatrix; template class SparseMatrix; template class Vector; +class ConstraintMatrix; template class TriaRawIterator; @@ -923,6 +924,76 @@ class DoFCellAccessor : public DoFAccessor void get_dof_values (const InputVector &values, Vector &local_values) const; + /** + * Return the values of the given vector + * restricted to the dofs of this + * cell in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, etc, + * dofs on line 0, dofs on line 1, etc, + * dofs on quad 0, etc. + * + * The vector has to have the + * right size before being passed + * to this function. This + * function is only callable for + * active cells. + * + * The input vector may be either + * a Vector, + * Vector, or a + * BlockVector, or a + * PETSc or Trilinos vector if + * deal.II is compiled to support + * these libraries. It is in the + * responsibility of the caller + * to assure that the types of + * the numbers stored in input + * and output vectors are + * compatible and with similar + * accuracy. + */ + template + void get_dof_values (const InputVector &values, + ForwardIterator local_values_begin, + ForwardIterator local_values_end) const; + + /** + * Return the values of the given vector + * restricted to the dofs of this + * cell in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, etc, + * dofs on line 0, dofs on line 1, etc, + * dofs on quad 0, etc. + * + * The vector has to have the + * right size before being passed + * to this function. This + * function is only callable for + * active cells. + * + * The input vector may be either a + * Vector, + * Vector, or a + * BlockVector, or a PETSc or + * Trilinos vector if deal.II is + * compiled to support these + * libraries. It is in the + * responsibility of the caller to + * assure that the types of the numbers + * stored in input and output vectors + * are compatible and with similar + * accuracy. The ConstraintMatrix + * passed as an argument to this + * function makes sure that constraints + * are correctly distributed when the + * dof values are calculated. + */ + template + void get_dof_values (const ConstraintMatrix &constraints, + const InputVector &values, + ForwardIterator local_values_begin, + ForwardIterator local_values_end) const; + /** * This function is the counterpart to * get_dof_values(): it takes a vector @@ -1111,6 +1182,50 @@ class DoFCellAccessor : public DoFAccessor distribute_local_to_global (const Vector &local_source, OutputVector &global_destination) const; + /** + * Distribute a local (cell based) + * vector in iterator format to a + * global one by mapping the local + * numbering of the degrees of freedom + * to the global one and entering the + * local values into the global vector. + * + * The elements are added up + * to the elements in the global + * vector, rather than just set, since + * this is usually what one wants. + */ + template + void + distribute_local_to_global (ForwardIterator local_source_begin, + ForwardIterator local_source_end, + OutputVector &global_destination) const; + + /** + * Distribute a local (cell based) + * vector in iterator format to a + * global one by mapping the local + * numbering of the degrees of freedom + * to the global one and entering the + * local values into the global vector. + * + * The elements are added up + * to the elements in the global + * vector, rather than just set, since + * this is usually what one + * wants. Moreover, the + * ConstraintMatrix passed to this + * function makes sure that also + * constraints are eliminated in this + * process. + */ + template + void + distribute_local_to_global (const ConstraintMatrix &constraints, + ForwardIterator local_source_begin, + ForwardIterator local_source_end, + OutputVector &global_destination) const; + /** * This function does much the * same as the @@ -1127,6 +1242,19 @@ class DoFCellAccessor : public DoFAccessor distribute_local_to_global (const FullMatrix &local_source, OutputMatrix &global_destination) const; + /** + * This function does what the two + * distribute_local_to_global + * functions with vector and matrix + * argument do, but all at once. + */ + template + void + distribute_local_to_global (const FullMatrix &local_matrix, + const Vector &local_vector, + OutputMatrix &global_matrix, + OutputVector &global_vector) const; + /** * @} */ diff --git a/deal.II/deal.II/include/dofs/dof_accessor.templates.h b/deal.II/deal.II/include/dofs/dof_accessor.templates.h index 312a82e64d..c95b59a760 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.templates.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.templates.h @@ -15,6 +15,7 @@ #include +#include #include #include #include @@ -48,9 +49,10 @@ DoFAccessor (const Triangulation *tria, const int index, const DH *dof_handler) : - internal::DoFAccessor::Inheritance::BaseClass (tria, - level, - index), + internal::DoFAccessor::Inheritance::BaseClass (tria, + level, + index), dof_handler(const_cast(dof_handler)) {} @@ -191,10 +193,10 @@ namespace internal static unsigned int get_dof_index (const dealii::DoFHandler<1,spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, + const unsigned int obj_level, + const unsigned int obj_index, + const unsigned int fe_index, + const unsigned int local_index, internal::int2type<1>) { return dof_handler.levels[obj_level]->lines. @@ -2072,17 +2074,16 @@ namespace internal * indices that we hold for each * cell. */ - template + template static void get_dof_values (const DoFCellAccessor > &accessor, const InputVector &values, - dealii::Vector &local_values) + ForwardIterator local_values_begin, + ForwardIterator local_values_end) { - typedef - dealii::DoFAccessor > - BaseClass; - Assert (local_values.size() == accessor.get_fe().dofs_per_cell, + typedef dealii::DoFAccessor > BaseClass; + Assert (local_values_end-local_values_begin == accessor.get_fe().dofs_per_cell, typename BaseClass::ExcVectorDoesNotMatch()); Assert (values.size() == accessor.get_dof_handler().n_dofs(), typename BaseClass::ExcVectorDoesNotMatch()); @@ -2097,9 +2098,10 @@ namespace internal ExcMessage ("Cell must either be active, or all DoFs must be in vertices")); unsigned int *cache = &accessor.dof_handler->levels[accessor.level()] - ->cell_dof_indices_cache[accessor.present_index * accessor.get_fe().dofs_per_cell]; - for (unsigned int i=0; icell_dof_indices_cache[accessor.present_index * + accessor.get_fe().dofs_per_cell]; + for ( ; local_values_begin != local_values_end; ++local_values_begin, ++cache) + *local_values_begin = values(*cache); } /** @@ -2109,22 +2111,98 @@ namespace internal * not have a cache for the local * DoF indices. */ - template + template static void get_dof_values (const DoFCellAccessor > &accessor, const InputVector &values, - dealii::Vector &local_values) + ForwardIterator local_values_begin, + ForwardIterator local_values_end) { + typedef dealii::DoFAccessor > BaseClass; + // no caching for hp::DoFHandler // implemented + Assert (local_values_end-local_values_begin == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell; std::vector local_dof_indices (dofs_per_cell); get_dof_indices (accessor, local_dof_indices); - for (unsigned int i=0; i + static + void + get_dof_values (const DoFCellAccessor > &accessor, + const ConstraintMatrix &constraints, + const InputVector &values, + ForwardIterator local_values_begin, + ForwardIterator local_values_end) + { + typedef dealii::DoFAccessor > BaseClass; + Assert (local_values_end-local_values_begin == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (values.size() == accessor.get_dof_handler().n_dofs(), + typename BaseClass::ExcVectorDoesNotMatch()); + + // check as in documentation that + // cell is either active, or dofs + // are only in vertices + Assert (!accessor.has_children() + || + (accessor.get_fe().dofs_per_cell == + accessor.get_fe().dofs_per_vertex * GeometryInfo::vertices_per_cell), + ExcMessage ("Cell must either be active, or all DoFs must be in vertices")); + + unsigned int *cache = &accessor.dof_handler->levels[accessor.level()] + ->cell_dof_indices_cache[accessor.present_index * + accessor.get_fe().dofs_per_cell]; + constraints.get_dof_values(values, *cache, local_values_begin, + local_values_end); + } + + /** + * Same function as above except + * that it works for + * hp::DoFHandler objects that do + * not have a cache for the local + * DoF indices. + */ + template + static + void + get_dof_values (const DoFCellAccessor > &accessor, + const ConstraintMatrix &constraints, + const InputVector &values, + ForwardIterator local_values_begin, + ForwardIterator local_values_end) + { + // no caching for hp::DoFHandler + // implemented + typedef dealii::DoFAccessor > BaseClass; + Assert (local_values_end-local_values_begin == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell; + + std::vector local_dof_indices (dofs_per_cell); + get_dof_indices (accessor, local_dof_indices); + + constraints.get_dof_values (values, local_dof_indices.begin(), + local_values_begin, local_values_end); } @@ -2140,9 +2218,7 @@ namespace internal const dealii::Vector &local_values, OutputVector &values) { - typedef - dealii::DoFAccessor > - BaseClass; + typedef dealii::DoFAccessor > BaseClass; Assert (local_values.size() == accessor.get_fe().dofs_per_cell, typename BaseClass::ExcVectorDoesNotMatch()); Assert (values.size() == accessor.get_dof_handler().n_dofs(), @@ -2234,9 +2310,7 @@ namespace internal // ::DoFHandler only supports a // single active fe with index // zero - typedef - dealii::DoFAccessor > - BaseClass; + typedef dealii::DoFAccessor > BaseClass; Assert (i == 0, typename BaseClass::ExcInvalidObject()); } @@ -2248,9 +2322,7 @@ namespace internal set_active_fe_index (DoFCellAccessor > &accessor, const unsigned int i) { - typedef - dealii::DoFAccessor > - BaseClass; + typedef dealii::DoFAccessor > BaseClass; Assert (accessor.dof_handler != 0, typename BaseClass::ExcInvalidObject()); Assert (static_cast(accessor.level()) < @@ -2266,22 +2338,20 @@ namespace internal - template + template static void distribute_local_to_global (const DoFCellAccessor > &accessor, - const dealii::Vector &local_source, - OutputVector &global_destination) + ForwardIterator local_source_begin, + ForwardIterator local_source_end, + OutputVector &global_destination) { - typedef - dealii::DoFAccessor > - BaseClass; - + typedef dealii::DoFAccessor > BaseClass; Assert (accessor.dof_handler != 0, typename BaseClass::ExcInvalidObject()); Assert (&accessor.get_fe() != 0, typename BaseClass::ExcInvalidObject()); - Assert (local_source.size() == accessor.get_fe().dofs_per_cell, + Assert (local_source_end-local_source_begin == accessor.get_fe().dofs_per_cell, typename BaseClass::ExcVectorDoesNotMatch()); Assert (accessor.dof_handler->n_dofs() == global_destination.size(), typename BaseClass::ExcVectorDoesNotMatch()); @@ -2295,38 +2365,109 @@ namespace internal accessor.get_fe().dofs_per_vertex * GeometryInfo::vertices_per_cell), ExcMessage ("Cell must either be active, or all DoFs must be in vertices")); - const unsigned int n_dofs = local_source.size(); + const unsigned int n_dofs = local_source_end - local_source_begin; unsigned int * dofs = &accessor.dof_handler->levels[accessor.level()] ->cell_dof_indices_cache[accessor.present_index * n_dofs]; // distribute cell vector - global_destination.add(n_dofs, dofs,local_source.begin()); + global_destination.add(n_dofs, dofs, local_source_begin); } - template + template static void distribute_local_to_global (const DoFCellAccessor > &accessor, - const dealii::Vector &local_source, - OutputVector &global_destination) + ForwardIterator local_source_begin, + ForwardIterator local_source_end, + OutputVector &global_destination) { - typedef - dealii::DoFAccessor > - BaseClass; + typedef dealii::DoFAccessor > BaseClass; + Assert (accessor.dof_handler != 0, + typename BaseClass::ExcInvalidObject()); + Assert (&accessor.get_fe() != 0, + typename BaseClass::ExcInvalidObject()); + Assert (local_source_end-local_source_begin == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_destination.size(), + typename BaseClass::ExcVectorDoesNotMatch()); + + const unsigned int n_dofs = local_source_end - local_source_begin; + +//TODO[WB/MK]: This function could me made more efficient because it allocates memory, which could be avoided by passing in another argument as a scratch array. This should be fixed eventually + + // get indices of dofs + std::vector dofs (n_dofs); + accessor.get_dof_indices (dofs); + + // distribute cell vector + global_destination.add (n_dofs, dofs.begin(), local_source_begin); + } + + - Assert (accessor->dof_handler != 0, + template + static + void + distribute_local_to_global (const DoFCellAccessor > &accessor, + const ConstraintMatrix &constraints, + ForwardIterator local_source_begin, + ForwardIterator local_source_end, + OutputVector &global_destination) + { + typedef dealii::DoFAccessor > BaseClass; + Assert (accessor.dof_handler != 0, typename BaseClass::ExcInvalidObject()); Assert (&accessor.get_fe() != 0, typename BaseClass::ExcInvalidObject()); - Assert (local_source.size() == accessor.get_fe().dofs_per_cell, + Assert (local_source_end-local_source_begin == accessor.get_fe().dofs_per_cell, typename BaseClass::ExcVectorDoesNotMatch()); Assert (accessor.dof_handler->n_dofs() == global_destination.size(), typename BaseClass::ExcVectorDoesNotMatch()); - const unsigned int n_dofs = local_source.size(); + // check as in documentation that + // cell is either active, or dofs + // are only in vertices + Assert (!accessor.has_children() + || + (accessor.get_fe().dofs_per_cell == + accessor.get_fe().dofs_per_vertex * GeometryInfo::vertices_per_cell), + ExcMessage ("Cell must either be active, or all DoFs must be in vertices")); + + const unsigned int n_dofs = local_source_end - local_source_begin; + + unsigned int * dofs = &accessor.dof_handler->levels[accessor.level()] + ->cell_dof_indices_cache[accessor.present_index * n_dofs]; + + // distribute cell vector + constraints.distribute_local_to_global (local_source_begin, local_source_end, + dofs, global_destination); + } + + + + template + static + void + distribute_local_to_global (const DoFCellAccessor > &accessor, + const ConstraintMatrix &constraints, + ForwardIterator local_source_begin, + ForwardIterator local_source_end, + OutputVector &global_destination) + { + typedef dealii::DoFAccessor > BaseClass; + Assert (accessor.dof_handler != 0, + typename BaseClass::ExcInvalidObject()); + Assert (&accessor.get_fe() != 0, + typename BaseClass::ExcInvalidObject()); + Assert (local_source_end-local_source_begin == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_destination.size(), + typename BaseClass::ExcVectorDoesNotMatch()); + + const unsigned int n_dofs = local_source_end - local_source_begin; //TODO[WB/MK]: This function could me made more efficient because it allocates memory, which could be avoided by passing in another argument as a scratch array. This should be fixed eventually @@ -2335,7 +2476,8 @@ namespace internal accessor.get_dof_indices (dofs); // distribute cell vector - global_destination.add (dofs, local_source); + constraints.distribute_local_to_global (local_source_begin, local_source_end, + dofs.begin(), global_destination); } @@ -2347,10 +2489,7 @@ namespace internal const dealii::FullMatrix &local_source, OutputMatrix &global_destination) { - typedef - dealii::DoFAccessor > - BaseClass; - + typedef dealii::DoFAccessor > BaseClass; Assert (accessor.dof_handler != 0, typename BaseClass::ExcInvalidObject()); Assert (&accessor.get_fe() != 0, @@ -2378,7 +2517,7 @@ namespace internal unsigned int * dofs = &accessor.dof_handler->levels[accessor.level()] ->cell_dof_indices_cache[accessor.present_index * n_dofs]; - // distribute cell matrices + // distribute cell matrix for (unsigned int i=0; i &local_source, OutputMatrix &global_destination) { - typedef - dealii::DoFAccessor > - BaseClass; - - Assert (accessor->dof_handler != 0, + typedef dealii::DoFAccessor > BaseClass; + Assert (accessor.dof_handler != 0, typename BaseClass::ExcInvalidObject()); Assert (&accessor.get_fe() != 0, typename BaseClass::ExcInvalidObject()); @@ -2418,9 +2554,103 @@ namespace internal std::vector dofs (n_dofs); accessor.get_dof_indices (dofs); - // distribute cell vector + // distribute cell matrix global_destination.add(dofs,local_source); } + + + + template + static + void + distribute_local_to_global (const DoFCellAccessor > &accessor, + const dealii::FullMatrix &local_matrix, + const dealii::Vector &local_vector, + OutputMatrix &global_matrix, + OutputVector &global_vector) + { + typedef dealii::DoFAccessor > BaseClass; + Assert (accessor.dof_handler != 0, + typename BaseClass::ExcInvalidObject()); + Assert (&accessor.get_fe() != 0, + typename BaseClass::ExcInvalidObject()); + Assert (local_matrix.m() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (local_matrix.n() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_matrix.m(), + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_matrix.n(), + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (local_vector.size() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_vector.size(), + typename BaseClass::ExcVectorDoesNotMatch()); + + // check as in documentation that + // cell is either active, or dofs + // are only in vertices + Assert (!accessor.has_children() + || + (accessor.get_fe().dofs_per_cell == + accessor.get_fe().dofs_per_vertex * GeometryInfo::vertices_per_cell), + ExcMessage ("Cell must either be active, or all DoFs must be in vertices")); + + const unsigned int n_dofs = accessor.get_fe().dofs_per_cell; + unsigned int * dofs = &accessor.dof_handler->levels[accessor.level()] + ->cell_dof_indices_cache[accessor.present_index * n_dofs]; + + // distribute cell matrices + for (unsigned int i=0; i + static + void + distribute_local_to_global (const DoFCellAccessor > &accessor, + const dealii::FullMatrix &local_matrix, + const dealii::Vector &local_vector, + OutputMatrix &global_matrix, + OutputVector &global_vector) + { + typedef dealii::DoFAccessor > BaseClass; + Assert (accessor.dof_handler != 0, + typename BaseClass::ExcInvalidObject()); + Assert (&accessor.get_fe() != 0, + typename BaseClass::ExcInvalidObject()); + Assert (local_matrix.m() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (local_matrix.n() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_matrix.m(), + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_matrix.n(), + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (local_vector.size() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_vector.size(), + typename BaseClass::ExcVectorDoesNotMatch()); + + const unsigned int n_dofs = local_matrix.size(); + +//TODO[WB/MK]: This function could me made more efficient because it allocates memory, which could be avoided by passing in another argument as a scratch array. + + // get indices of dofs + std::vector dofs (n_dofs); + accessor.get_dof_indices (dofs); + + // distribute cell matrix and vector + global_matrix.add(dofs,local_matrix); + global_vector.add(dofs,local_vector); + } }; } } @@ -2554,7 +2784,35 @@ DoFCellAccessor::get_dof_values (const InputVector &values, Vector &local_values) const { internal::DoFCellAccessor::Implementation - ::get_dof_values (*this, values, local_values); + ::get_dof_values (*this, values, local_values.begin(), local_values.end()); +} + + + +template +template +void +DoFCellAccessor::get_dof_values (const InputVector &values, + ForwardIterator local_values_begin, + ForwardIterator local_values_end) const +{ + internal::DoFCellAccessor::Implementation + ::get_dof_values (*this, values, local_values_begin, local_values_end); +} + + + +template +template +void +DoFCellAccessor::get_dof_values (const ConstraintMatrix &constraints, + const InputVector &values, + ForwardIterator local_values_begin, + ForwardIterator local_values_end) const +{ + internal::DoFCellAccessor::Implementation + ::get_dof_values (*this, constraints, values, + local_values_begin, local_values_end); } @@ -2611,7 +2869,41 @@ distribute_local_to_global (const Vector &local_source, OutputVector &global_destination) const { internal::DoFCellAccessor::Implementation:: - distribute_local_to_global (*this,local_source,global_destination); + distribute_local_to_global (*this, local_source.begin(), + local_source.end(), global_destination); +} + + + +template +template +inline +void +DoFCellAccessor:: +distribute_local_to_global (ForwardIterator local_source_begin, + ForwardIterator local_source_end, + OutputVector &global_destination) const +{ + internal::DoFCellAccessor::Implementation:: + distribute_local_to_global (*this, local_source_begin, + local_source_end, global_destination); +} + + + +template +template +inline +void +DoFCellAccessor:: +distribute_local_to_global (const ConstraintMatrix &constraints, + ForwardIterator local_source_begin, + ForwardIterator local_source_end, + OutputVector &global_destination) const +{ + internal::DoFCellAccessor::Implementation:: + distribute_local_to_global (*this, constraints, local_source_begin, + local_source_end, global_destination); } @@ -2629,6 +2921,23 @@ distribute_local_to_global (const FullMatrix &local_source, } + +template +template +inline +void +DoFCellAccessor:: +distribute_local_to_global (const FullMatrix &local_matrix, + const Vector &local_vector, + OutputMatrix &global_matrix, + OutputVector &global_vector) const +{ + internal::DoFCellAccessor::Implementation:: + distribute_local_to_global (*this,local_matrix,local_vector, + global_matrix,global_vector); +} + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index d2e07c4bea..436f7fea2e 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -341,6 +341,18 @@ inconvenience this causes.
    +
  1. +

    + Improved: The ConstraintMatrix class now uses a cache for random access to + the constraint lines. This considerably increases performance of the + *_local_to_global functions, where such an access pattern is usual. Moreover, + the ConstraintMatrix class has now a function get_dof_values that can import + data from a global vector to a cell vector with respecting the constraints. +

    +
    + (Martin Kronbichler 2009/09/30) +
  2. +
  3. Fixed: SparsityTools::reorder_Cuthill_McKee would produce an error if the diff --git a/deal.II/lac/include/lac/constraint_matrix.h b/deal.II/lac/include/lac/constraint_matrix.h index e590cf36d1..393c2fc529 100644 --- a/deal.II/lac/include/lac/constraint_matrix.h +++ b/deal.II/lac/include/lac/constraint_matrix.h @@ -929,6 +929,21 @@ class ConstraintMatrix : public Subscriptor void condense (BlockSparseMatrix &matrix, BlockVectorType &vector) const; + /** + * Delete hanging nodes in a + * vector. Sets all hanging node + * values to zero. The @p + * VectorType may be a + * Vector, Vector, + * BlockVector<...>, a + * PETSc or Trilinos vector + * wrapper class, or any other + * type having the same + * interface. + */ + template + void set_zero (VectorType &vec) const; + /** * @} */ @@ -976,52 +991,144 @@ class ConstraintMatrix : public Subscriptor * vectors and matrices are fully * assembled. * - * Note that, unless an optional - * FullMatrix object is provided, this - * function will apply all constraints - * as if they were homogeneous and - * throw an exception in case it - * encounters inhomogeneous - * constraints. For correctly setting - * inhomogeneous constraints, you - * should provide an additional matrix - * argument or use one of the functions - * with both matrix and vector - * arguments. + * Note that this function will apply all + * constraints as if they were + * homogeneous. For correctly setting + * inhomogeneous constraints, use the + * similar function with a matrix + * argument or the function with both + * matrix and vector arguments. * - * The optional argument + * Note: This function is not + * thread-safe, so you will need to make + * sure that only on process at a time + * calls this function. + */ + template + void + distribute_local_to_global (const Vector &local_vector, + const std::vector &local_dof_indices, + VectorType &global_vector) const; + + /** + * This function takes a vector of + * local contributions (@p + * local_vector) corresponding to the + * degrees of freedom indices given in + * @p local_dof_indices and distributes + * them to the global vector. In most + * cases, these local contributions + * will be the result of an integration + * over a cell or face of a + * cell. However, as long as @p + * local_vector and @p + * local_dof_indices have the same + * number of elements, this function is + * happy with whatever it is + * given. + * + * In contrast to the similar function in + * the DoFAccessor class, this function + * also takes care of constraints, + * i.e. if one of the elements of @p + * local_dof_indices belongs to a + * constrained node, then rather than + * writing the corresponding element of + * @p local_vector into @p global_vector, + * the element is distributed to the + * entries in the global vector to which + * this particular degree of freedom is + * constrained. + * + * Thus, by using this function to + * distribute local contributions to the + * global object, one saves the call to + * the condense function after the + * vectors and matrices are fully + * assembled. + * + * The fourth argument * local_matrix is intended to * be used in case one wants to apply * inhomogeneous constraints on the - * vector only. Such a situation could - * be where one wants to assemble of a - * right hand side vector on a problem - * with inhomogeneous constraints, but - * the global matrix has been assembled - * previously. A typical example of - * this is a time stepping algorithm - * where the stiffness matrix is - * assembled once, and the right hand - * side updated every time step. Note - * that, however, the entries in the - * columns of the local matrix have to - * be exactly the same as those that - * have been written into the global - * matrix. Otherwise, this function - * will not be able to correctly handle - * inhomogeneities. + * vector only. Such a situation could be + * where one wants to assemble of a right + * hand side vector on a problem with + * inhomogeneous constraints, but the + * global matrix has been assembled + * previously. A typical example of this + * is a time stepping algorithm where the + * stiffness matrix is assembled once, + * and the right hand side updated every + * time step. Note that, however, the + * entries in the columns of the local + * matrix have to be exactly the same as + * those that have been written into the + * global matrix. Otherwise, this + * function will not be able to correctly + * handle inhomogeneities. * * Note: This function is not - * thread-safe, so you will need to - * make sure that only on process at a - * time calls this function. + * thread-safe, so you will need to make + * sure that only on process at a time + * calls this function. */ template void distribute_local_to_global (const Vector &local_vector, const std::vector &local_dof_indices, VectorType &global_vector, - const FullMatrix &local_matrix = FullMatrix()) const; + const FullMatrix &local_matrix) const; + + + /** + * This function takes a pointer to a + * vector of local contributions (@p + * local_vector) corresponding to the + * degrees of freedom indices given in + * @p local_dof_indices and distributes + * them to the global vector. In most + * cases, these local contributions + * will be the result of an integration + * over a cell or face of a + * cell. However, as long as the + * entries in @p local_dof_indices + * indicate reasonable global vector + * entries, this function is happy with + * whatever it is given. + * + * If one of the elements of @p + * local_dof_indices belongs to a + * constrained node, then rather than + * writing the corresponding element of + * @p local_vector into @p + * global_vector, the element is + * distributed to the entries in the + * global vector to which this + * particular degree of freedom is + * constrained. + * + * Thus, by using this function to + * distribute local contributions to + * the global object, one saves the + * call to the condense function after + * the vectors and matrices are fully + * assembled. Note that this function + * completely ignores inhomogeneous + * constraints. + * + * Note: This function is not + * thread-safe, so you will need to + * make sure that only on process at a + * time calls this function. + */ + template + void + distribute_local_to_global (ForwardIteratorVec local_vector_begin, + ForwardIteratorVec local_vector_end, + ForwardIteratorInd local_indices_begin, + VectorType &global_vector) const; /** * This function takes a matrix of @@ -1201,23 +1308,49 @@ class ConstraintMatrix : public Subscriptor add_entries_local_to_global (const std::vector &local_dof_indices, SparsityType &sparsity_pattern, const bool keep_constrained_entries = true, - const Table<2,bool> &dof_mask = default_empty_table) const; - - /** - * Delete hanging nodes in a - * vector. Sets all hanging node - * values to zero. The @p - * VectorType may be a - * Vector, Vector, - * BlockVector<...>, a - * PETSc or Trilinos vector - * wrapper class, or any other - * type having the same - * interface. - */ - template - void set_zero (VectorType &vec) const; + const Table<2,bool> &dof_mask = Table<2,bool>()) const; + /** + * This function imports values from a + * global vector (@p global_vector) by + * applying the constraints to a vector + * of local values, expressed in + * iterator format. In most cases, the + * local values will be identified by + * the local dof values on a + * cell. However, as long as the + * entries in @p local_dof_indices + * indicate reasonable global vector + * entries, this function is happy with + * whatever it is given. + * + * If one of the elements of @p + * local_dof_indices belongs to a + * constrained node, then rather than + * writing the corresponding element of + * @p global_vector into @p + * local_vector, the constraints are + * resolved as the respective + * distribute function does, i.e., the + * local entry is constructed from the + * global entries to which this + * particular degree of freedom is + * constrained. + * + * In contrast to the similar function + * get_dof_values in the DoFAccessor + * class, this function does not need + * the constrained values to be + * correctly set (i.e., distribute to + * be called). + */ + template + void + get_dof_values (const VectorType &global_vector, + ForwardIteratorInd local_indices_begin, + ForwardIteratorVec local_vector_begin, + ForwardIteratorVec local_vector_end) const; /** * @} @@ -1363,9 +1496,9 @@ class ConstraintMatrix : public Subscriptor /** * Value of the inhomogeneity. */ - double inhomogeneity; + double inhomogeneity; - /** + /** * This operator is a bit weird and * unintuitive: it compares the * line numbers of two lines. We @@ -1379,7 +1512,7 @@ class ConstraintMatrix : public Subscriptor */ bool operator < (const ConstraintLine &) const; - /** + /** * This operator is likewise weird: * it checks whether the line * indices of the two operands are @@ -1412,17 +1545,21 @@ class ConstraintMatrix : public Subscriptor * to store the lines. This, however, * would mean a much more fractioned * heap since it allocates many small - * objects, ans would additionally make + * objects, and would additionally make * usage of this matrix much slower. */ std::vector lines; /** - * A list of flags that indicate - * whether there is a constraint line - * for a given degree of freedom - * index. Note that this class has no - * notion of how many degrees of + * A list of pointers that contains the + * address of the ConstraintLine of a + * constrained degree of freedom, or NULL + * if the degree of freedom is not + * constrained. The NULL return value + * returns thus whether there is a + * constraint line for a given degree of + * freedom index. Note that this class + * has no notion of how many degrees of * freedom there really are, so if we * check whether there is a constraint * line for a given degree of freedom, @@ -1430,39 +1567,45 @@ class ConstraintMatrix : public Subscriptor * shorter than the index of the DoF we * check for. * - * This field exists since when adding - * a new constraint line we have to - * figure out whether it already + * This field exists since when adding a + * new constraint line we have to figure + * out whether it already * exists. Previously, we would simply * walk the unsorted list of constraint * lines until we either hit the end or - * found it. This algorithm is O(N) if - * N is the number of constraints, - * which makes it O(N^2) when inserting - * all constraints. For large problems - * with many constraints, this could - * easily take 5-10 per cent of the - * total run time. With this field, we - * can at least save this time when - * checking whether a new constraint - * line already exists. + * found it. This algorithm is O(N) if N + * is the number of constraints, which + * makes it O(N^2) when inserting all + * constraints. For large problems with + * many constraints, this could easily + * take 5-10 per cent of the total run + * time. With this field, we can save + * this time since we find any constraint + * in O(1) time or get to know that it a + * certain degree of freedom is not + * constrained. * * To make things worse, traversing the - * list of existing constraints - * requires reads from many different - * places in memory. Thus, in large 3d - * applications, the add_line() - * function showed up very prominently - * in the overall compute time, mainly - * because it generated a lot of cache + * list of existing constraints requires + * reads from many different places in + * memory. Thus, in large 3d + * applications, the add_line() function + * showed up very prominently in the + * overall compute time, mainly because + * it generated a lot of cache * misses. This should also be fixed by - * using the O(1) algorithm to access - * the fields of this array. + * using the O(1) algorithm to access the + * fields of this array. * * The field is useful in a number of - * other contexts as well, though. + * other contexts as well, e.g. when one + * needs random access to the constraints + * as in all the functions that apply + * constraints on the fly while add cell + * contributions into vectors and + * matrices. */ - std::vector constraint_line_exists; + std::vector lines_cache; /** * Store whether the arrays are sorted. @@ -1479,13 +1622,6 @@ class ConstraintMatrix : public Subscriptor */ static bool check_zero_weight (const std::pair &p); - /** - * Dummy table that serves as default - * argument for function - * add_entries_local_to_global(). - */ - static const Table<2,bool> default_empty_table; - /** * This function actually implements * the local_to_global function for @@ -1540,26 +1676,6 @@ class ConstraintMatrix : public Subscriptor const Table<2,bool> &dof_mask, internal::bool2type) const; - /** - * This function returns a pointer to - * the respective ConstraintLine object - * if a dof is constrained, and a zero - * pointer otherwise. - */ - std::vector::const_iterator - find_constraint (const unsigned int line) const; - - /** - * This function returns a pointer to - * the respective ConstraintLine object - * if a dof is constrained, and a zero - * pointer otherwise. This function is - * used if the constraint matrix has not - * been closed yet. - */ - std::vector::iterator - find_constraint (const unsigned int line); - #ifdef DEAL_II_USE_TRILINOS //TODO: Make use of the following member thread safe /** @@ -1588,7 +1704,7 @@ ConstraintMatrix::ConstraintMatrix (const ConstraintMatrix &constraint_matrix) : Subscriptor (), lines (constraint_matrix.lines), - constraint_line_exists (constraint_matrix.constraint_line_exists), + lines_cache (constraint_matrix.lines_cache), sorted (constraint_matrix.sorted) #ifdef DEAL_II_USE_TRILINOS ,vec_distribute () @@ -1606,22 +1722,36 @@ ConstraintMatrix::add_line (const unsigned int line) // check whether line already exists; it // may, in which case we can just quit - if ((line < constraint_line_exists.size()) + if ((line < lines_cache.size()) && - (constraint_line_exists[line] == true)) - return; + (lines_cache[line] != 0)) + { + Assert (lines_cache[line]->line == line, ExcInternalError()); + return; + } // if necessary enlarge vector of - // existing entries - if (line >= constraint_line_exists.size()) - constraint_line_exists.resize (line+1, false); - constraint_line_exists[line] = true; + // existing entries for cache + if (line >= lines_cache.size()) + lines_cache.resize (line+1,0); + // enlarge lines vector if necessary. + // need to reset the pointers to the + // ConstraintLine entries in that case, + // since new memory will be allocated + if (lines.size() == lines.capacity()) + { + lines.reserve (2*lines.size()+8); + std::vector::const_iterator it = lines.begin(); + for ( ; it != lines.end(); ++it) + lines_cache[it->line] = &*it; + } // push a new line to the end of the // list lines.push_back (ConstraintLine()); lines.back().line = line; lines.back().inhomogeneity = 0.; + lines_cache[line] = &lines.back(); } @@ -1636,8 +1766,6 @@ ConstraintMatrix::add_entry (const unsigned int line, Assert (line != column, ExcMessage ("Can't constrain a degree of freedom to itself")); - std::vector::iterator line_ptr = find_constraint(line); - // if in debug mode, check whether an // entry for this column already // exists and if its the same as @@ -1646,6 +1774,9 @@ ConstraintMatrix::add_entry (const unsigned int line, // in any case: exit the function if an // entry for this column already exists, // since we don't want to enter it twice + ConstraintLine* line_ptr = const_cast(lines_cache[line]); + Assert (line_ptr != 0, ExcInternalError()); + Assert (line_ptr->line == line, ExcInternalError()); for (std::vector >::const_iterator p=line_ptr->entries.begin(); p != line_ptr->entries.end(); ++p) @@ -1666,7 +1797,8 @@ void ConstraintMatrix::set_inhomogeneity (const unsigned int line, const double value) { - find_constraint(line)->inhomogeneity = value; + ConstraintLine* line_ptr = const_cast(lines_cache[line]); + line_ptr->inhomogeneity = value; } @@ -1675,9 +1807,9 @@ inline bool ConstraintMatrix::is_constrained (const unsigned int index) const { - return ((index < constraint_line_exists.size()) + return ((index < lines_cache.size()) && - (constraint_line_exists[index] == true)); + (lines_cache[index] != 0)); } @@ -1686,68 +1818,89 @@ inline bool ConstraintMatrix::is_inhomogeneously_constrained (const unsigned int index) const { - const std::vector::const_iterator position - = find_constraint(index); - return position!=lines.end() ? position->inhomogeneity != 0 : false; + return (is_constrained(index) && + lines_cache[index]->inhomogeneity != 0); } +template inline -std::vector::const_iterator -ConstraintMatrix::find_constraint (const unsigned int line) const +void +ConstraintMatrix:: +distribute_local_to_global (const Vector &local_vector, + const std::vector &local_dof_indices, + VectorType &global_vector) const { - if (is_constrained(line) == false) - return lines.end(); - - Assert (sorted==true, ExcMatrixNotClosed()); - ConstraintLine index_comparison; - index_comparison.line = line; + Assert (local_vector.size() == local_dof_indices.size(), + ExcDimensionMismatch(local_vector.size(), local_dof_indices.size())); + distribute_local_to_global (local_vector.begin(), local_vector.end(), + local_dof_indices.begin(), global_vector); +} - const std::vector::const_iterator - position = std::lower_bound (lines.begin(), - lines.end(), - index_comparison); - Assert (position != lines.end(), ExcInternalError()); - // check whether we've really found the - // right constraint. - Assert (position->line == line, - ExcInternalError()); - return position; +template +inline +void ConstraintMatrix:: + distribute_local_to_global (ForwardIteratorVec local_vector_begin, + ForwardIteratorVec local_vector_end, + ForwardIteratorInd local_indices_begin, + VectorType &global_vector) const +{ + Assert (sorted == true, ExcMatrixNotClosed()); + for ( ; local_vector_begin != local_vector_end; + ++local_vector_begin, ++local_indices_begin) + { + if (is_constrained(*local_indices_begin) == false) + global_vector(*local_indices_begin) += *local_vector_begin; + else + { + const ConstraintLine* position = + lines_cache[*local_indices_begin]; + for (unsigned int j=0; jentries.size(); ++j) + { + Assert (is_constrained(position->entries[j].first) == false, + ExcMessage ("Tried to distribute to a fixed dof.")); + global_vector(position->entries[j].first) + += *local_vector_begin * position->entries[j].second; + } + } + } } +template inline -std::vector::iterator -ConstraintMatrix::find_constraint (const unsigned int line) +void ConstraintMatrix::get_dof_values (const VectorType &global_vector, + ForwardIteratorInd local_indices_begin, + ForwardIteratorVec local_vector_begin, + ForwardIteratorVec local_vector_end) const { - Assert (sorted==false, ExcMatrixIsClosed()); - std::vector::iterator line_ptr; - const std::vector::const_iterator start=lines.begin(); - - // the usual case is that the line where - // a value is entered is the one we - // added last, so we search backward - for (line_ptr=(lines.end()-1); line_ptr!=start; --line_ptr) - if (line_ptr->line == line) - break; - - // if the loop didn't break, then - // line_ptr must be begin(). - // we have an error if that doesn't - // point to 'line' then - Assert (line_ptr->line==line, ExcLineInexistant(line)); - // second check: the respective - // flag in the - // constraint_line_exists field - // must exist - Assert (line < constraint_line_exists.size(), ExcInternalError()); - Assert (constraint_line_exists[line] == true, ExcInternalError()); - - return line_ptr; + Assert (sorted == true, ExcMatrixNotClosed()); + for ( ; local_vector_begin != local_vector_end; + ++local_vector_begin, ++local_indices_begin) + { + if (is_constrained(*local_indices_begin) == false) + *local_vector_begin = global_vector(*local_indices_begin); + else + { + const ConstraintLine* position = + lines_cache[*local_indices_begin]; + typename VectorType::value_type value = position->inhomogeneity; + for (unsigned int j=0; jentries.size(); ++j) + { + Assert (is_constrained(position->entries[j].first) == false, + ExcMessage ("Tried to distribute to a fixed dof.")); + value += (global_vector(position->entries[j].first) * + position->entries[j].second); + } + *local_vector_begin = value; + } + } } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/lac/include/lac/constraint_matrix.templates.h b/deal.II/lac/include/lac/constraint_matrix.templates.h index a6c2ca41c6..b5d535b48f 100644 --- a/deal.II/lac/include/lac/constraint_matrix.templates.h +++ b/deal.II/lac/include/lac/constraint_matrix.templates.h @@ -734,14 +734,10 @@ distribute_local_to_global (const Vector &local_vector, Assert (local_vector.size() == local_dof_indices.size(), ExcDimensionMismatch(local_vector.size(), local_dof_indices.size())); Assert (sorted == true, ExcMatrixNotClosed()); - const bool use_matrix = local_matrix.m() != 0 ? true : false; - if (local_matrix.m() != 0) - { - Assert (local_matrix.m() == local_dof_indices.size(), - ExcDimensionMismatch(local_matrix.m(), local_dof_indices.size())); - Assert (local_matrix.n() == local_dof_indices.size(), - ExcDimensionMismatch(local_matrix.n(), local_dof_indices.size())); - } + Assert (local_matrix.m() == local_dof_indices.size(), + ExcDimensionMismatch(local_matrix.m(), local_dof_indices.size())); + Assert (local_matrix.n() == local_dof_indices.size(), + ExcDimensionMismatch(local_matrix.n(), local_dof_indices.size())); const unsigned int n_local_dofs = local_vector.size(); if (lines.size() == 0) @@ -753,50 +749,41 @@ distribute_local_to_global (const Vector &local_vector, // vector that tells about whether a // certain constraint exists. then, we // simply copy over the data. - const std::vector::const_iterator position = - find_constraint(local_dof_indices[i]); - - if (position==lines.end()) + if (is_constrained(local_dof_indices[i]) == false) { global_vector(local_dof_indices[i]) += local_vector(i); continue; } - if (use_matrix) - { - const double val = position->inhomogeneity; - if (val != 0) - for (unsigned int j=0; jinhomogeneity; + if (val != 0) + for (unsigned int j=0; j::const_iterator - position_j = find_constraint(local_dof_indices[j]); + const double matrix_entry = local_matrix(j,i); + if (matrix_entry == 0) + continue; - if (position_j == lines.end()) - global_vector(local_dof_indices[j]) -= val * local_matrix(j,i); - else + for (unsigned int q=0; qentries.size(); ++q) { - const double matrix_entry = local_matrix(j,i); - if (matrix_entry == 0) - continue; - - for (unsigned int q=0; qentries.size(); ++q) - { - Assert (is_constrained(position_j->entries[q].first) == false, - ExcMessage ("Tried to distribute to a fixed dof.")); - global_vector(position_j->entries[q].first) - -= val * position_j->entries[q].second * matrix_entry; - } + Assert (is_constrained(position_j->entries[q].first) == false, + ExcMessage ("Tried to distribute to a fixed dof.")); + global_vector(position_j->entries[q].first) + -= val * position_j->entries[q].second * matrix_entry; } } - } - else - // in case the constraint is - // inhomogeneous and we have no matrix - // available, this function is not - // appropriate. Throw an exception. - Assert (position->inhomogeneity == 0., - ExcMessage ("Inhomogeneous constraint cannot be condensed " - "without any matrix specified.")); + } // now distribute the constraint, // but make sure we don't touch @@ -1259,9 +1246,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, // resolved in a second step. for (unsigned int i = 0; i::const_iterator position = - find_constraint(local_dof_indices[i]); - if (position == lines.end()) + if (is_constrained(local_dof_indices[i]) == false) { my_indices[added_rows].global_row = local_dof_indices[i]; my_indices[added_rows].local_row = i; @@ -1270,7 +1255,9 @@ distribute_local_to_global (const FullMatrix &local_matrix, } constraint_lines.push_back (std::make_pair(i,&*position)); + const ConstraintLine *>(i,lines_cache[local_dof_indices[i]])); + Assert (lines_cache[local_dof_indices[i]]->line == local_dof_indices[i], + ExcInternalError()); } Assert (constraint_lines.size() + added_rows == n_local_dofs, ExcInternalError()); @@ -1577,9 +1564,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, unsigned int added_rows = 0; for (unsigned int i = 0; i::const_iterator position - = find_constraint (local_dof_indices[i]); - if (position == lines.end()) + if (is_constrained(local_dof_indices[i]) == false) { my_indices[added_rows].global_row = local_dof_indices[i]; my_indices[added_rows].local_row = i; @@ -1588,7 +1573,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, } constraint_lines.push_back (std::make_pair(i,&*position)); + const ConstraintLine *>(i,lines_cache[local_dof_indices[i]])); } Assert (constraint_lines.size() + added_rows == n_local_dofs, ExcInternalError()); @@ -1827,9 +1812,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, std::vector > constraint_lines; for (unsigned int i = 0; i::const_iterator position - = find_constraint(local_dof_indices[i]); - if (position == lines.end()) + if (is_constrained(local_dof_indices[i]) == false) { actual_dof_indices[added_rows] = local_dof_indices[i]; ++added_rows; @@ -1837,7 +1820,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, } constraint_lines.push_back (std::make_pair(i,&*position)); + const ConstraintLine *>(i,lines_cache[local_dof_indices[i]])); } Assert (constraint_lines.size() + added_rows == n_local_dofs, ExcInternalError()); @@ -1918,9 +1901,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, // resolved in a second step. for (unsigned int i = 0; i::const_iterator position - = find_constraint(local_dof_indices[i]); - if (position == lines.end()) + if (is_constrained(local_dof_indices[i]) == false) { my_indices[added_rows].global_row = local_dof_indices[i]; my_indices[added_rows].local_row = i; @@ -1929,7 +1910,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, } constraint_lines.push_back (std::make_pair(i,&*position)); + const ConstraintLine *>(i,lines_cache[local_dof_indices[i]])); } Assert (constraint_lines.size() + added_rows == n_local_dofs, ExcInternalError()); @@ -2144,9 +2125,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, std::vector > constraint_lines; for (unsigned int i = 0; i::const_iterator position - = find_constraint(local_dof_indices[i]); - if (position == lines.end()) + if (is_constrained(local_dof_indices[i]) == false) { actual_dof_indices[added_rows] = local_dof_indices[i]; ++added_rows; @@ -2154,7 +2133,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, } constraint_lines.push_back (std::make_pair(i,&*position)); + const ConstraintLine *>(i,lines_cache[local_dof_indices[i]])); } Assert (constraint_lines.size() + added_rows == n_local_dofs, ExcInternalError()); @@ -2254,8 +2233,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, // resolved in a second step. for (unsigned int i = 0; i &local_dof_indices, continue; } - ConstraintLine index_comparison; - index_comparison.line = local_dof_indices[i]; - - const std::vector::const_iterator - position = std::lower_bound (lines.begin(), - lines.end(), - index_comparison); - Assert (position->line == local_dof_indices[i], - ExcInternalError()); - constraint_lines.push_back (std::make_pair(i,&*position)); + const ConstraintLine *>(i,lines_cache[local_dof_indices[i]])); } Assert (constraint_lines.size() + added_rows == n_local_dofs, ExcInternalError()); diff --git a/deal.II/lac/source/constraint_matrix.cc b/deal.II/lac/source/constraint_matrix.cc index 6d806f72a5..e0bc0c15b6 100644 --- a/deal.II/lac/source/constraint_matrix.cc +++ b/deal.II/lac/source/constraint_matrix.cc @@ -55,11 +55,6 @@ DEAL_II_NAMESPACE_OPEN - // Static member variable -const Table<2,bool> ConstraintMatrix::default_empty_table = Table<2,bool>(); - - - bool ConstraintMatrix::check_zero_weight (const std::pair &p) { @@ -98,21 +93,15 @@ ConstraintMatrix::add_entries (const unsigned int line, const std::vector > &col_val_pairs) { Assert (sorted==false, ExcMatrixIsClosed()); + Assert (is_constrained(line), ExcLineInexistant(line)); - std::vector::iterator line_ptr; - const std::vector::const_iterator start=lines.begin(); - // the usual case is that the line where - // a value is entered is the one we - // added last, so we search backward - for (line_ptr=(lines.end()-1); line_ptr!=start; --line_ptr) - if (line_ptr->line == line) - break; + ConstraintLine * line_ptr = const_cast(lines_cache[line]); + Assert (line_ptr->line == line, ExcInternalError()); // if the loop didn't break, then // line_ptr must be begin(). // we have an error if that doesn't // point to 'line' then - Assert (line_ptr->line==line, ExcLineInexistant(line)); // if in debug mode, check whether an // entry for this column already @@ -157,6 +146,24 @@ void ConstraintMatrix::close () // sort the lines std::sort (lines.begin(), lines.end()); + // update list of pointers and give the + // vector a sharp size since we won't + // modify the size any more after this + // point. + { + std::vector new_lines (lines_cache.size()); + for (std::vector::const_iterator line=lines.begin(); + line!=lines.end(); ++line) + new_lines[line->line] = &*line; + std::swap (lines_cache, new_lines); + } + + // in debug mode: check whether we really + // set the pointers correctly. + for (unsigned int i=0; iline, ExcInternalError()); + // first, strip zero entries, as we // have to do that only once for (std::vector::iterator line = lines.begin(); @@ -229,19 +236,9 @@ void ConstraintMatrix::close () Assert (dof_index != line->line, ExcMessage ("Cycle in constraints detected!")); - // find the line - // corresponding to - // this entry. note - // that we have - // already sorted - // them to make this - // process faster - ConstraintLine test_line; - test_line.line = dof_index; - const std::vector::const_iterator - constrained_line = std::lower_bound (lines.begin(), - lines.end(), - test_line); + const ConstraintLine * constrained_line = lines_cache[dof_index]; + Assert (constrained_line->line == dof_index, + ExcInternalError()); // now we have to // replace an entry @@ -468,15 +465,9 @@ void ConstraintMatrix::close () // make sure that // entry->first is not the // index of a line itself - ConstraintLine test_line; - test_line.line = entry->first; - const std::vector::const_iterator - test_line_position = std::lower_bound (lines.begin(), - lines.end(), - test_line); - Assert ((test_line_position == lines.end()) - || - (test_line_position->line != entry->first), + const ConstraintLine * it = entry->first < lines_cache.size() ? + lines_cache[entry->first] : 0; + Assert (it == 0, ExcDoFConstrainedToConstrainedDoF(line->line, entry->first)); }; #endif @@ -543,15 +534,8 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints) const bool object_was_sorted = sorted; sorted = false; - // before we even start: merge the - // two flag arrays - if (other_constraints.constraint_line_exists.size() > - constraint_line_exists.size()) - constraint_line_exists.resize (other_constraints.constraint_line_exists.size(), - false); - for (unsigned int i=0; i lines_cache.size()) + lines_cache.resize(other_constraints.lines_cache.size()); // first action is to fold into the // present object possible @@ -705,6 +689,9 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints) // afterwards as well. otherwise // leave everything in the unsorted // state + for (std::vector::const_iterator line=lines.begin(); + line!=lines.end(); ++line) + lines_cache[line->line] = &*line; if (object_was_sorted == true) close (); } @@ -713,8 +700,7 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints) void ConstraintMatrix::shift (const unsigned int offset) { - constraint_line_exists.insert (constraint_line_exists.begin(), offset, - false); + lines_cache.insert (lines_cache.begin(), offset, 0); for (std::vector::iterator i = lines.begin(); i != lines.end(); i++) @@ -737,8 +723,8 @@ void ConstraintMatrix::clear () } { - std::vector tmp; - constraint_line_exists.swap (tmp); + std::vector tmp; + lines_cache.swap (tmp); } #ifdef DEAL_II_USE_TRILINOS @@ -1871,7 +1857,6 @@ template<> void ConstraintMatrix::distribute (TrilinosWrappers::MPI::Vector &vec) const { - Assert (sorted == true, ExcMatrixNotClosed()); ConstraintLine index_comparison; index_comparison.line = vec.local_range().first; @@ -1955,39 +1940,14 @@ bool ConstraintMatrix::is_identity_constrained (const unsigned int index) const if (is_constrained(index) == false) return false; - if (sorted == true) - { - ConstraintLine index_comparison; - index_comparison.line = index; + const ConstraintLine * p = lines_cache[index]; + Assert (p->line == index, ExcInternalError()); - const std::vector::const_iterator - p = std::lower_bound (lines.begin (), - lines.end (), - index_comparison); // return if an entry for this // line was found and if it has // only one entry equal to 1.0 - // - // note that lower_bound only - // returns a valid iterator if - // 'index' is less than the - // largest line index in out - // constraints list - return ((p != lines.end()) && - (p->line == index) && - (p->entries.size() == 1) && - (p->entries[0].second == 1.0)); - } - else - { - for (std::vector::const_iterator i=lines.begin(); - i!=lines.end(); ++i) - if (i->line == index) - return ((i->entries.size() == 1) && - (i->entries[0].second == 1.0)); - - return false; - } + return ((p->entries.size() == 1) && + (p->entries[0].second == 1.0)); } @@ -2087,7 +2047,7 @@ unsigned int ConstraintMatrix::memory_consumption () const { return (MemoryConsumption::memory_consumption (lines) + - MemoryConsumption::memory_consumption (constraint_line_exists) + + MemoryConsumption::memory_consumption (lines_cache) + MemoryConsumption::memory_consumption (sorted)); } -- 2.39.5