template <typename number> class FullMatrix;
template <typename number> class SparseMatrix;
template <typename number> class Vector;
+class ConstraintMatrix;
template <typename Accessor> class TriaRawIterator;
void get_dof_values (const InputVector &values,
Vector<number> &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 <tt>Vector<float></tt>,
+ * Vector<double>, or a
+ * BlockVector<double>, 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 <class InputVector, typename ForwardIterator>
+ 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
+ * <tt>Vector<float></tt>,
+ * Vector<double>, or a
+ * BlockVector<double>, 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 <class InputVector, typename ForwardIterator>
+ 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
distribute_local_to_global (const Vector<number> &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 <em>added</em> up
+ * to the elements in the global
+ * vector, rather than just set, since
+ * this is usually what one wants.
+ */
+ template <typename ForwardIterator, typename OutputVector>
+ 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 <em>added</em> 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 <typename ForwardIterator, typename OutputVector>
+ 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
distribute_local_to_global (const FullMatrix<number> &local_source,
OutputMatrix &global_destination) const;
+ /**
+ * This function does what the two
+ * <tt>distribute_local_to_global</tt>
+ * functions with vector and matrix
+ * argument do, but all at once.
+ */
+ template <typename number, typename OutputMatrix, typename OutputVector>
+ void
+ distribute_local_to_global (const FullMatrix<number> &local_matrix,
+ const Vector<number> &local_vector,
+ OutputMatrix &global_matrix,
+ OutputVector &global_vector) const;
+
/**
* @}
*/
#include <base/config.h>
+#include <lac/constraint_matrix.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_levels.h>
#include <dofs/dof_faces.h>
const int index,
const DH *dof_handler)
:
- internal::DoFAccessor::Inheritance<structdim,DH::dimension,DH::space_dimension>::BaseClass (tria,
- level,
- index),
+ internal::DoFAccessor::Inheritance<structdim,DH::dimension,
+ DH::space_dimension>::BaseClass (tria,
+ level,
+ index),
dof_handler(const_cast<DH*>(dof_handler))
{}
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.
* indices that we hold for each
* cell.
*/
- template <int dim, int spacedim, class InputVector, typename number>
+ template <int dim, int spacedim, class InputVector, typename ForwardIterator>
static
void
get_dof_values (const DoFCellAccessor<DoFHandler<dim,spacedim> > &accessor,
const InputVector &values,
- dealii::Vector<number> &local_values)
+ ForwardIterator local_values_begin,
+ ForwardIterator local_values_end)
{
- typedef
- dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
- BaseClass;
- Assert (local_values.size() == accessor.get_fe().dofs_per_cell,
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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());
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; i<accessor.get_fe().dofs_per_cell; ++i, ++cache)
- local_values(i) = values(*cache);
+ ->cell_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);
}
/**
* not have a cache for the local
* DoF indices.
*/
- template <int dim, int spacedim, class InputVector, typename number>
+ template <int dim, int spacedim, class InputVector, typename ForwardIterator>
static
void
get_dof_values (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
const InputVector &values,
- dealii::Vector<number> &local_values)
+ ForwardIterator local_values_begin,
+ ForwardIterator local_values_end)
{
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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<unsigned int> local_dof_indices (dofs_per_cell);
get_dof_indices (accessor, local_dof_indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- local_values(i) = values(local_dof_indices[i]);
+ for (unsigned int i=0; i<dofs_per_cell; ++i, ++local_values_begin)
+ *local_values_begin = values(local_dof_indices[i]);
+ }
+
+
+ /**
+ * A function that collects the
+ * values of degrees of freedom. This
+ * function works for ::DoFHandler
+ * and all template arguments and
+ * uses the data from the cache of
+ * indices that we hold for each
+ * cell.
+ */
+ template <int dim, int spacedim, class InputVector, typename ForwardIterator>
+ static
+ void
+ get_dof_values (const DoFCellAccessor<DoFHandler<dim,spacedim> > &accessor,
+ const ConstraintMatrix &constraints,
+ const InputVector &values,
+ ForwardIterator local_values_begin,
+ ForwardIterator local_values_end)
+ {
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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<dim>::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 <int dim, int spacedim, class InputVector, typename ForwardIterator>
+ static
+ void
+ get_dof_values (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
+ const ConstraintMatrix &constraints,
+ const InputVector &values,
+ ForwardIterator local_values_begin,
+ ForwardIterator local_values_end)
+ {
+ // no caching for hp::DoFHandler
+ // implemented
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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<unsigned int> 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);
}
const dealii::Vector<number> &local_values,
OutputVector &values)
{
- typedef
- dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
- BaseClass;
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
Assert (local_values.size() == accessor.get_fe().dofs_per_cell,
typename BaseClass::ExcVectorDoesNotMatch());
Assert (values.size() == accessor.get_dof_handler().n_dofs(),
// ::DoFHandler only supports a
// single active fe with index
// zero
- typedef
- dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
- BaseClass;
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
Assert (i == 0, typename BaseClass::ExcInvalidObject());
}
set_active_fe_index (DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
const unsigned int i)
{
- typedef
- dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
- BaseClass;
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
Assert (accessor.dof_handler != 0,
typename BaseClass::ExcInvalidObject());
Assert (static_cast<unsigned int>(accessor.level()) <
- template <int dim, int spacedim, typename number, class OutputVector>
+ template <int dim, int spacedim, typename ForwardIterator, class OutputVector>
static
void
distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
- const dealii::Vector<number> &local_source,
- OutputVector &global_destination)
+ ForwardIterator local_source_begin,
+ ForwardIterator local_source_end,
+ OutputVector &global_destination)
{
- typedef
- dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
- BaseClass;
-
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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());
accessor.get_fe().dofs_per_vertex * GeometryInfo<dim>::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 <int dim, int spacedim, typename number, class OutputVector>
+ template <int dim, int spacedim, typename ForwardIterator, class OutputVector>
static
void
distribute_local_to_global (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
- const dealii::Vector<number> &local_source,
- OutputVector &global_destination)
+ ForwardIterator local_source_begin,
+ ForwardIterator local_source_end,
+ OutputVector &global_destination)
{
- typedef
- dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
- BaseClass;
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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<unsigned int> 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 <int dim, int spacedim, typename ForwardIterator, class OutputVector>
+ static
+ void
+ distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
+ const ConstraintMatrix &constraints,
+ ForwardIterator local_source_begin,
+ ForwardIterator local_source_end,
+ OutputVector &global_destination)
+ {
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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<dim>::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 <int dim, int spacedim, typename ForwardIterator, class OutputVector>
+ static
+ void
+ distribute_local_to_global (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
+ const ConstraintMatrix &constraints,
+ ForwardIterator local_source_begin,
+ ForwardIterator local_source_end,
+ OutputVector &global_destination)
+ {
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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
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);
}
const dealii::FullMatrix<number> &local_source,
OutputMatrix &global_destination)
{
- typedef
- dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
- BaseClass;
-
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
Assert (accessor.dof_handler != 0,
typename BaseClass::ExcInvalidObject());
Assert (&accessor.get_fe() != 0,
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<n_dofs; ++i)
global_destination.add(dofs[i], n_dofs, dofs,
&local_source(i,0));
const dealii::FullMatrix<number> &local_source,
OutputMatrix &global_destination)
{
- typedef
- dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
- BaseClass;
-
- Assert (accessor->dof_handler != 0,
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > BaseClass;
+ Assert (accessor.dof_handler != 0,
typename BaseClass::ExcInvalidObject());
Assert (&accessor.get_fe() != 0,
typename BaseClass::ExcInvalidObject());
std::vector<unsigned int> dofs (n_dofs);
accessor.get_dof_indices (dofs);
- // distribute cell vector
+ // distribute cell matrix
global_destination.add(dofs,local_source);
}
+
+
+
+ template <int dim, int spacedim, typename number,
+ class OutputMatrix, typename OutputVector>
+ static
+ void
+ distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
+ const dealii::FullMatrix<number> &local_matrix,
+ const dealii::Vector<number> &local_vector,
+ OutputMatrix &global_matrix,
+ OutputVector &global_vector)
+ {
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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<dim>::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<n_dofs; ++i)
+ {
+ global_matrix.add(dofs[i], n_dofs, dofs, &local_matrix(i,0));
+ global_vector(dofs[i]) += local_vector(i);
+ }
+ }
+
+
+
+ template <int dim, int spacedim, typename number,
+ class OutputMatrix, typename OutputVector>
+ static
+ void
+ distribute_local_to_global (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
+ const dealii::FullMatrix<number> &local_matrix,
+ const dealii::Vector<number> &local_vector,
+ OutputMatrix &global_matrix,
+ OutputVector &global_vector)
+ {
+ typedef dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> > 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<unsigned int> 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);
+ }
};
}
}
Vector<number> &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 <class DH>
+template <class InputVector, typename ForwardIterator>
+void
+DoFCellAccessor<DH>::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 <class DH>
+template <class InputVector, typename ForwardIterator>
+void
+DoFCellAccessor<DH>::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);
}
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 <class DH>
+template <typename ForwardIterator, typename OutputVector>
+inline
+void
+DoFCellAccessor<DH>::
+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 <class DH>
+template <typename ForwardIterator, typename OutputVector>
+inline
+void
+DoFCellAccessor<DH>::
+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);
}
}
+
+template <class DH>
+template <typename number, typename OutputMatrix, typename OutputVector>
+inline
+void
+DoFCellAccessor<DH>::
+distribute_local_to_global (const FullMatrix<number> &local_matrix,
+ const Vector<number> &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
<ol>
+ <li>
+ <p>
+ 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.
+ </p>
+ <br>
+ (Martin Kronbichler 2009/09/30)
+ </li>
+
<li>
<p>
Fixed: SparsityTools::reorder_Cuthill_McKee would produce an error if the
void condense (BlockSparseMatrix<number> &matrix,
BlockVectorType &vector) const;
+ /**
+ * Delete hanging nodes in a
+ * vector. Sets all hanging node
+ * values to zero. The @p
+ * VectorType may be a
+ * Vector<float>, Vector<double>,
+ * BlockVector<tt><...></tt>, a
+ * PETSc or Trilinos vector
+ * wrapper class, or any other
+ * type having the same
+ * interface.
+ */
+ template <class VectorType>
+ void set_zero (VectorType &vec) const;
+
/**
* @}
*/
* 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 <typename VectorType>
+ void
+ distribute_local_to_global (const Vector<double> &local_vector,
+ const std::vector<unsigned int> &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
* <tt>local_matrix</tt> 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 <typename VectorType>
void
distribute_local_to_global (const Vector<double> &local_vector,
const std::vector<unsigned int> &local_dof_indices,
VectorType &global_vector,
- const FullMatrix<double> &local_matrix = FullMatrix<double>()) const;
+ const FullMatrix<double> &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 <typename ForwardIteratorVec, typename ForwardIteratorInd,
+ class VectorType>
+ 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
add_entries_local_to_global (const std::vector<unsigned int> &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<float>, Vector<double>,
- * BlockVector<tt><...></tt>, a
- * PETSc or Trilinos vector
- * wrapper class, or any other
- * type having the same
- * interface.
- */
- template <class VectorType>
- 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 <typename ForwardIteratorVec, typename ForwardIteratorInd,
+ class VectorType>
+ void
+ get_dof_values (const VectorType &global_vector,
+ ForwardIteratorInd local_indices_begin,
+ ForwardIteratorVec local_vector_begin,
+ ForwardIteratorVec local_vector_end) const;
/**
* @}
/**
* 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
*/
bool operator < (const ConstraintLine &) const;
- /**
+ /**
* This operator is likewise weird:
* it checks whether the line
* indices of the two operands are
* 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<ConstraintLine> 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,
* 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<bool> constraint_line_exists;
+ std::vector<const ConstraintLine*> lines_cache;
/**
* Store whether the arrays are sorted.
*/
static bool check_zero_weight (const std::pair<unsigned int, double> &p);
- /**
- * Dummy table that serves as default
- * argument for function
- * <tt>add_entries_local_to_global()</tt>.
- */
- static const Table<2,bool> default_empty_table;
-
/**
* This function actually implements
* the local_to_global function for
const Table<2,bool> &dof_mask,
internal::bool2type<true>) const;
- /**
- * This function returns a pointer to
- * the respective ConstraintLine object
- * if a dof is constrained, and a zero
- * pointer otherwise.
- */
- std::vector<ConstraintLine>::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<ConstraintLine>::iterator
- find_constraint (const unsigned int line);
-
#ifdef DEAL_II_USE_TRILINOS
//TODO: Make use of the following member thread safe
/**
:
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 ()
// 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<ConstraintLine>::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();
}
Assert (line != column,
ExcMessage ("Can't constrain a degree of freedom to itself"));
- std::vector<ConstraintLine>::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
// 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<ConstraintLine*>(lines_cache[line]);
+ Assert (line_ptr != 0, ExcInternalError());
+ Assert (line_ptr->line == line, ExcInternalError());
for (std::vector<std::pair<unsigned int,double> >::const_iterator
p=line_ptr->entries.begin();
p != line_ptr->entries.end(); ++p)
ConstraintMatrix::set_inhomogeneity (const unsigned int line,
const double value)
{
- find_constraint(line)->inhomogeneity = value;
+ ConstraintLine* line_ptr = const_cast<ConstraintLine*>(lines_cache[line]);
+ line_ptr->inhomogeneity = value;
}
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));
}
bool
ConstraintMatrix::is_inhomogeneously_constrained (const unsigned int index) const
{
- const std::vector<ConstraintLine>::const_iterator position
- = find_constraint(index);
- return position!=lines.end() ? position->inhomogeneity != 0 : false;
+ return (is_constrained(index) &&
+ lines_cache[index]->inhomogeneity != 0);
}
+template <class VectorType>
inline
-std::vector<ConstraintMatrix::ConstraintLine>::const_iterator
-ConstraintMatrix::find_constraint (const unsigned int line) const
+void
+ConstraintMatrix::
+distribute_local_to_global (const Vector<double> &local_vector,
+ const std::vector<unsigned int> &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<ConstraintLine>::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 <typename ForwardIteratorVec, typename ForwardIteratorInd,
+ class VectorType>
+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; j<position->entries.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 <typename ForwardIteratorVec, typename ForwardIteratorInd,
+ class VectorType>
inline
-std::vector<ConstraintMatrix::ConstraintLine>::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<ConstraintLine>::iterator line_ptr;
- const std::vector<ConstraintLine>::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; j<position->entries.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
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)
// vector that tells about whether a
// certain constraint exists. then, we
// simply copy over the data.
- const std::vector<ConstraintLine>::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; j<n_local_dofs; ++j)
+ const ConstraintLine * position =
+ lines_cache.size() <= local_dof_indices[i] ? 0 :
+ lines_cache[local_dof_indices[i]];
+
+ const double val = position->inhomogeneity;
+ if (val != 0)
+ for (unsigned int j=0; j<n_local_dofs; ++j)
+ {
+ const ConstraintLine * position_j =
+ lines_cache.size() <= local_dof_indices[j] ? 0 :
+ lines_cache[local_dof_indices[j]];
+
+ if (position_j == 0)
+ global_vector(local_dof_indices[j]) -= val * local_matrix(j,i);
+ else
{
- const std::vector<ConstraintLine>::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; q<position_j->entries.size(); ++q)
{
- const double matrix_entry = local_matrix(j,i);
- if (matrix_entry == 0)
- continue;
-
- for (unsigned int q=0; q<position_j->entries.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
// resolved in a second step.
for (unsigned int i = 0; i<n_local_dofs; ++i)
{
- const std::vector<ConstraintLine>::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;
}
constraint_lines.push_back (std::make_pair<unsigned int,
- const ConstraintLine *>(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());
unsigned int added_rows = 0;
for (unsigned int i = 0; i<n_local_dofs; ++i)
{
- const std::vector<ConstraintLine>::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;
}
constraint_lines.push_back (std::make_pair<unsigned int,
- const ConstraintLine *>(i,&*position));
+ const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
}
Assert (constraint_lines.size() + added_rows == n_local_dofs,
ExcInternalError());
std::vector<std::pair<unsigned int, const ConstraintLine *> > constraint_lines;
for (unsigned int i = 0; i<n_local_dofs; ++i)
{
- const std::vector<ConstraintLine>::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;
}
constraint_lines.push_back (std::make_pair<unsigned int,
- const ConstraintLine *>(i,&*position));
+ const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
}
Assert (constraint_lines.size() + added_rows == n_local_dofs,
ExcInternalError());
// resolved in a second step.
for (unsigned int i = 0; i<n_local_dofs; ++i)
{
- const std::vector<ConstraintLine>::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;
}
constraint_lines.push_back (std::make_pair<unsigned int,
- const ConstraintLine *>(i,&*position));
+ const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
}
Assert (constraint_lines.size() + added_rows == n_local_dofs,
ExcInternalError());
std::vector<std::pair<unsigned int, const ConstraintLine *> > constraint_lines;
for (unsigned int i = 0; i<n_local_dofs; ++i)
{
- const std::vector<ConstraintLine>::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;
}
constraint_lines.push_back (std::make_pair<unsigned int,
- const ConstraintLine *>(i,&*position));
+ const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
}
Assert (constraint_lines.size() + added_rows == n_local_dofs,
ExcInternalError());
// resolved in a second step.
for (unsigned int i = 0; i<n_local_dofs; ++i)
{
- if (constraint_line_exists.size() <= local_dof_indices[i] ||
- constraint_line_exists[local_dof_indices[i]] == false)
+ 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;
continue;
}
- ConstraintLine index_comparison;
- index_comparison.line = local_dof_indices[i];
-
- const std::vector<ConstraintLine>::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<unsigned int,
- const ConstraintLine *>(i,&*position));
+ const ConstraintLine *>(i,lines_cache[local_dof_indices[i]]));
}
Assert (constraint_lines.size() + added_rows == n_local_dofs,
ExcInternalError());
- // Static member variable
-const Table<2,bool> ConstraintMatrix::default_empty_table = Table<2,bool>();
-
-
-
bool
ConstraintMatrix::check_zero_weight (const std::pair<unsigned int, double> &p)
{
const std::vector<std::pair<unsigned int,double> > &col_val_pairs)
{
Assert (sorted==false, ExcMatrixIsClosed());
+ Assert (is_constrained(line), ExcLineInexistant(line));
- std::vector<ConstraintLine>::iterator line_ptr;
- const std::vector<ConstraintLine>::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<ConstraintLine*>(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
// 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<const ConstraintLine*> new_lines (lines_cache.size());
+ for (std::vector<ConstraintLine>::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; i<lines_cache.size(); ++i)
+ if (lines_cache[i] != 0)
+ Assert (i == lines_cache[i]->line, ExcInternalError());
+
// first, strip zero entries, as we
// have to do that only once
for (std::vector<ConstraintLine>::iterator line = lines.begin();
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<ConstraintLine>::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
// make sure that
// entry->first is not the
// index of a line itself
- ConstraintLine test_line;
- test_line.line = entry->first;
- const std::vector<ConstraintLine>::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
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<other_constraints.constraint_line_exists.size(); ++i)
- if (other_constraints.constraint_line_exists[i] == true)
- constraint_line_exists[i] = true;
+ if (other_constraints.lines_cache.size() > lines_cache.size())
+ lines_cache.resize(other_constraints.lines_cache.size());
// first action is to fold into the
// present object possible
// afterwards as well. otherwise
// leave everything in the unsorted
// state
+ for (std::vector<ConstraintLine>::const_iterator line=lines.begin();
+ line!=lines.end(); ++line)
+ lines_cache[line->line] = &*line;
if (object_was_sorted == true)
close ();
}
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<ConstraintLine>::iterator i = lines.begin();
i != lines.end(); i++)
}
{
- std::vector<bool> tmp;
- constraint_line_exists.swap (tmp);
+ std::vector<const ConstraintLine*> tmp;
+ lines_cache.swap (tmp);
}
#ifdef DEAL_II_USE_TRILINOS
void
ConstraintMatrix::distribute (TrilinosWrappers::MPI::Vector &vec) const
{
- Assert (sorted == true, ExcMatrixNotClosed());
ConstraintLine index_comparison;
index_comparison.line = vec.local_range().first;
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<ConstraintLine>::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<ConstraintLine>::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));
}
ConstraintMatrix::memory_consumption () const
{
return (MemoryConsumption::memory_consumption (lines) +
- MemoryConsumption::memory_consumption (constraint_line_exists) +
+ MemoryConsumption::memory_consumption (lines_cache) +
MemoryConsumption::memory_consumption (sorted));
}