From: Martin Kronbichler Date: Mon, 21 Sep 2009 06:48:08 +0000 (+0000) Subject: Work on a few TODOs in DoFCellAccessor::distribute_local_to_global. Let these impleme... X-Git-Tag: v8.0.0~7092 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=84fc8c185e1a2e2f4431ed3826a2c4a82236399b;p=dealii.git Work on a few TODOs in DoFCellAccessor::distribute_local_to_global. Let these implementations use the cell_dof_indices_cache for regular DoFHandler functions. Implement a collective add function for all vectors. git-svn-id: https://svn.dealii.org/trunk@19483 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 26618816a3..17ea766c97 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.templates.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.templates.h @@ -2163,6 +2163,8 @@ namespace internal values(*cache) = local_values(i); } + + template static void @@ -2244,7 +2246,7 @@ namespace internal static void set_active_fe_index (DoFCellAccessor > &accessor, - const unsigned int i) + const unsigned int i) { typedef dealii::DoFAccessor > @@ -2261,6 +2263,165 @@ namespace internal accessor.dof_handler->levels[accessor.level()] ->active_fe_indices[accessor.present_index] = i; } + + + + template + static + void + distribute_local_to_global (const DoFCellAccessor > &accessor, + const dealii::Vector &local_source, + 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, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_destination.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 = local_source.size(); + + 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, + &(*const_cast*>(&local_source))(0)); + } + + + + template + static + void + distribute_local_to_global (const DoFCellAccessor > &accessor, + const dealii::Vector &local_source, + 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, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_destination.size(), + typename BaseClass::ExcVectorDoesNotMatch()); + + const unsigned int n_dofs = local_source.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. 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 (dofs, local_source); + } + + + + template + static + void + distribute_local_to_global (const DoFCellAccessor > &accessor, + const dealii::FullMatrix &local_source, + OutputMatrix &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.m() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (local_source.n() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_destination.m(), + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_destination.n(), + typename BaseClass::ExcMatrixDoesNotMatch()); + + // 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.m(); + + 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_source, + OutputMatrix &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.m() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (local_source.n() == accessor.get_fe().dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_destination.m(), + typename BaseClass::ExcMatrixDoesNotMatch()); + Assert (accessor.dof_handler->n_dofs() == global_destination.n(), + typename BaseClass::ExcMatrixDoesNotMatch()); + + const unsigned int n_dofs = local_source.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 vector + global_destination.add(dofs,local_source); + } }; } } @@ -2450,26 +2611,8 @@ DoFCellAccessor:: distribute_local_to_global (const Vector &local_source, OutputVector &global_destination) const { - Assert (this->dof_handler != 0, - typename BaseClass::ExcInvalidObject()); - Assert (&this->get_fe() != 0, - typename BaseClass::ExcInvalidObject()); - Assert (local_source.size() == this->get_fe().dofs_per_cell, - typename BaseClass::ExcVectorDoesNotMatch()); - Assert (this->dof_handler->n_dofs() == global_destination.size(), - typename BaseClass::ExcVectorDoesNotMatch()); - - const unsigned int n_dofs = local_source.size(); - -//TODO[WB]: This function could me made more efficient. First, it allocates memory, which could be avoided by passing in another argument as a scratch array (or by using the dof_indices cache for non-hp DoFHandlers, see the functions in the Implementation class above). second, the elementwise access is really slow if we use PETSc vectors/matrices. This should be fixed eventually - - // get indices of dofs - std::vector dofs (n_dofs); - this->get_dof_indices (dofs); - - // distribute cell vector - for (unsigned int j=0; j:: distribute_local_to_global (const FullMatrix &local_source, OutputMatrix &global_destination) const { - Assert (this->dof_handler != 0, - typename BaseClass::ExcInvalidObject()); - Assert (&this->get_fe() != 0, - typename BaseClass::ExcInvalidObject()); - Assert (local_source.m() == this->get_fe().dofs_per_cell, - typename BaseClass::ExcVectorDoesNotMatch()); - Assert (local_source.m() == local_source.n(), - typename BaseClass::ExcMatrixDoesNotMatch()); - Assert (this->dof_handler->n_dofs() == global_destination.m(), - typename BaseClass::ExcMatrixDoesNotMatch()); - Assert (global_destination.m() == global_destination.n(), - typename BaseClass::ExcMatrixDoesNotMatch()); - - const unsigned int n_dofs = local_source.m(); - -//TODO[WB]: This function could me made more efficient. First, it allocates memory, which could be avoided by passing in another argument as a scratch array (or by using the dof_indices cache for non-hp DoFHandlers, see the functions in the Implementation class above). second, the elementwise access is really slow if we use PETSc vectors/matrices. This should be fixed eventually - - // get indices of dofs - std::vector dofs (n_dofs); - this->get_dof_indices (dofs); - - // distribute cell matrix - for (unsigned int i=0; i + void add (const std::vector &indices, + const std::vector &values); + + /** + * This is a second collective + * add operation. As a + * difference, this function + * takes a deal.II vector of + * values. + */ + template + void add (const std::vector &indices, + const Vector &values); + + /** + * Take an address where + * n_elements are stored + * contiguously and add them into + * the vector. Handles all cases + * which are not covered by the + * other two add() + * functions above. + */ + template + void add (const unsigned int n_elements, + const unsigned int *indices, + const Number *values); + /** * $U(0-DIM)+=s$. Addition of s * to all components. Note that @@ -1865,13 +1903,56 @@ BlockVectorBase::operator -= (const BlockVectorBase& v) } return *this; } + + + +template +template +inline +void +BlockVectorBase ::add (const std::vector &indices, + const std::vector &values) +{ + Assert (indices.size() == values.size(), + ExcDimensionMismatch(indices.size(), values.size())); + add (indices.size(), &indices[0], &values[0]); +} template -void BlockVectorBase::add (const value_type a) +template +inline +void +BlockVectorBase ::add (const std::vector &indices, + const Vector &values) +{ + Assert (indices.size() == values.size(), + ExcDimensionMismatch(indices.size(), values.size())); + const unsigned int n_indices = indices.size(); + for (unsigned int i=0; i +template +inline +void +BlockVectorBase ::add (const unsigned int n_indices, + const unsigned int *indices, + const Number *values) { + for (unsigned int i=0; i +void BlockVectorBase::add (const value_type a) +{ Assert (numbers::is_finite(a), ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); diff --git a/deal.II/lac/include/lac/constraint_matrix.templates.h b/deal.II/lac/include/lac/constraint_matrix.templates.h index 3fd91ec336..a6c2ca41c6 100644 --- a/deal.II/lac/include/lac/constraint_matrix.templates.h +++ b/deal.II/lac/include/lac/constraint_matrix.templates.h @@ -744,68 +744,71 @@ distribute_local_to_global (const Vector &local_vector, } const unsigned int n_local_dofs = local_vector.size(); - for (unsigned int i=0; i::const_iterator position = - find_constraint(local_dof_indices[i]); + const std::vector::const_iterator position = + find_constraint(local_dof_indices[i]); - if (position==lines.end()) - { - global_vector(local_dof_indices[i]) += local_vector(i); - continue; - } + if (position==lines.end()) + { + 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; j::const_iterator - position_j = find_constraint(local_dof_indices[j]); + if (use_matrix) + { + const double val = position->inhomogeneity; + if (val != 0) + for (unsigned int j=0; j::const_iterator + position_j = find_constraint(local_dof_indices[j]); - if (position_j == lines.end()) - global_vector(local_dof_indices[j]) -= val * local_matrix(j,i); - else - { - 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; - } - } - } - } - else + if (position_j == lines.end()) + global_vector(local_dof_indices[j]) -= val * local_matrix(j,i); + else + { + 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; + } + } + } + } + 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.")); + 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 // the entries of fixed dofs - 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(i) * position->entries[j].second; - } - } + 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(i) * position->entries[j].second; + } + } } diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index 5a818e8a18..62fadd1b57 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -468,7 +468,7 @@ namespace SparseMatrixIterators * @; others can be generated in application programs (see the * section on @ref Instantiations in the manual). * - * @author Essentially everyone who has ever worked on deal.II, 1994-2007 + * @author Essentially everyone who has ever worked on deal.II, 1994-2009 */ template class SparseMatrix : public virtual Subscriptor diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h index e9d78239a5..ccec876778 100644 --- a/deal.II/lac/include/lac/vector.h +++ b/deal.II/lac/include/lac/vector.h @@ -21,6 +21,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -653,6 +654,43 @@ class Vector : public Subscriptor */ Vector & operator -= (const Vector &V); + /** + * A collective add operation: + * This funnction adds a whole + * set of values stored in @p + * values to the vector + * components specified by @p + * indices. + */ + template + void add (const std::vector &indices, + const std::vector &values); + + /** + * This is a second collective + * add operation. As a + * difference, this function + * takes a deal.II vector of + * values. + */ + template + void add (const std::vector &indices, + const Vector &values); + + /** + * Take an address where + * n_elements are stored + * contiguously and add them into + * the vector. Handles all cases + * which are not covered by the + * other two add() + * functions above. + */ + template + void add (const unsigned int n_elements, + const unsigned int *indices, + const OtherNumber *values); + /** * Addition of @p s to all * components. Note that @p s is a @@ -1254,7 +1292,55 @@ Vector::scale (const Number factor) (factor*boost::lambda::_1), internal::Vector::minimum_parallel_grain_size); } + + + +template +template +inline +void +Vector::add (const std::vector &indices, + const std::vector &values) +{ + Assert (indices.size() == values.size(), + ExcDimensionMismatch(indices.size(), values.size())); + add (indices.size(), &indices[0], &values[0]); +} + + + +template +template +inline +void +Vector::add (const std::vector &indices, + const Vector &values) +{ + Assert (indices.size() == values.size(), + ExcDimensionMismatch(indices.size(), values.size())); + add (indices.size(), &indices[0], values.val); +} + + + +template +template +inline +void +Vector::add (const unsigned int n_indices, + const unsigned int *indices, + const Number2 *values) +{ + for (unsigned int i=0; i