values(*cache) = local_values(i);
}
+
+
template <int dim, int spacedim, class OutputVector, typename number>
static
void
static
void
set_active_fe_index (DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
- const unsigned int i)
+ const unsigned int i)
{
typedef
dealii::DoFAccessor<dim,DoFHandler<dim,spacedim> >
accessor.dof_handler->levels[accessor.level()]
->active_fe_indices[accessor.present_index] = i;
}
+
+
+
+ template <int dim, int spacedim, typename number, class OutputVector>
+ static
+ void
+ distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
+ const dealii::Vector<number> &local_source,
+ 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,
+ 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<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();
+
+ 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<dealii::Vector<number>*>(&local_source))(0));
+ }
+
+
+
+ template <int dim, int spacedim, typename number, 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)
+ {
+ 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,
+ 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<unsigned int> dofs (n_dofs);
+ accessor.get_dof_indices (dofs);
+
+ // distribute cell vector
+ global_destination.add (dofs, local_source);
+ }
+
+
+
+ template <int dim, int spacedim, typename number, class OutputMatrix>
+ static
+ void
+ distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim> > &accessor,
+ const dealii::FullMatrix<number> &local_source,
+ OutputMatrix &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.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<dim>::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<n_dofs; ++i)
+ global_destination.add(dofs[i], n_dofs, dofs,
+ &local_source(i,0));
+ }
+
+
+
+ template <int dim, int spacedim, typename number, class OutputMatrix>
+ static
+ void
+ distribute_local_to_global (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim> > &accessor,
+ const dealii::FullMatrix<number> &local_source,
+ OutputMatrix &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.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<unsigned int> dofs (n_dofs);
+ accessor.get_dof_indices (dofs);
+
+ // distribute cell vector
+ global_destination.add(dofs,local_source);
+ }
};
}
}
distribute_local_to_global (const Vector<number> &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<unsigned int> dofs (n_dofs);
- this->get_dof_indices (dofs);
-
- // distribute cell vector
- for (unsigned int j=0; j<n_dofs; ++j)
- global_destination(dofs[j]) += local_source(j);
+ internal::DoFCellAccessor::Implementation::
+ distribute_local_to_global (*this,local_source,global_destination);
}
distribute_local_to_global (const FullMatrix<number> &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<unsigned int> dofs (n_dofs);
- this->get_dof_indices (dofs);
-
- // distribute cell matrix
- for (unsigned int i=0; i<n_dofs; ++i)
- for (unsigned int j=0; j<n_dofs; ++j)
- global_destination.add(dofs[i], dofs[j], local_source(i,j));
+ internal::DoFCellAccessor::Implementation::
+ distribute_local_to_global (*this,local_source,global_destination);
}
BlockVectorBase &
operator -= (const BlockVectorBase &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 <typename Number>
+ void add (const std::vector<unsigned int> &indices,
+ const std::vector<Number> &values);
+
+ /**
+ * This is a second collective
+ * add operation. As a
+ * difference, this function
+ * takes a deal.II vector of
+ * values.
+ */
+ template <typename Number>
+ void add (const std::vector<unsigned int> &indices,
+ const Vector<Number> &values);
+
+ /**
+ * Take an address where
+ * <tt>n_elements</tt> are stored
+ * contiguously and add them into
+ * the vector. Handles all cases
+ * which are not covered by the
+ * other two <tt>add()</tt>
+ * functions above.
+ */
+ template <typename Number>
+ void add (const unsigned int n_elements,
+ const unsigned int *indices,
+ const Number *values);
+
/**
* $U(0-DIM)+=s$. Addition of <tt>s</tt>
* to all components. Note that
}
return *this;
}
+
+
+
+template <class VectorType>
+template <typename Number>
+inline
+void
+BlockVectorBase<VectorType> ::add (const std::vector<unsigned int> &indices,
+ const std::vector<Number> &values)
+{
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(), values.size()));
+ add (indices.size(), &indices[0], &values[0]);
+}
template <class VectorType>
-void BlockVectorBase<VectorType>::add (const value_type a)
+template <typename Number>
+inline
+void
+BlockVectorBase<VectorType> ::add (const std::vector<unsigned int> &indices,
+ const Vector<Number> &values)
+{
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(), values.size()));
+ const unsigned int n_indices = indices.size();
+ for (unsigned int i=0; i<n_indices; ++i)
+ (*this)(indices[i]) += values(i);
+}
+
+
+
+template <class VectorType>
+template <typename Number>
+inline
+void
+BlockVectorBase<VectorType> ::add (const unsigned int n_indices,
+ const unsigned int *indices,
+ const Number *values)
{
+ for (unsigned int i=0; i<n_indices; ++i)
+ (*this)(indices[i]) += values[i];
+}
+
+
+template <class VectorType>
+void BlockVectorBase<VectorType>::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)"));
}
const unsigned int n_local_dofs = local_vector.size();
- for (unsigned int i=0; i<n_local_dofs; ++i)
- {
+ if (lines.size() == 0)
+ global_vector.add(local_dof_indices, local_vector);
+ else
+ for (unsigned int i=0; i<n_local_dofs; ++i)
+ {
// let's see if we can use the bool
// 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]);
+ const std::vector<ConstraintLine>::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<n_local_dofs; ++j)
- {
- const std::vector<ConstraintLine>::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<n_local_dofs; ++j)
+ {
+ const std::vector<ConstraintLine>::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; 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;
- }
- }
- }
- }
- 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; 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;
+ }
+ }
+ }
+ }
+ 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; 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(i) * position->entries[j].second;
- }
- }
+ 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(i) * position->entries[j].second;
+ }
+ }
}
#include <boost/lambda/lambda.hpp>
#include <cstdio>
+#include <vector>
DEAL_II_NAMESPACE_OPEN
*/
Vector<Number> & operator -= (const Vector<Number> &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 <typename OtherNumber>
+ void add (const std::vector<unsigned int> &indices,
+ const std::vector<OtherNumber> &values);
+
+ /**
+ * This is a second collective
+ * add operation. As a
+ * difference, this function
+ * takes a deal.II vector of
+ * values.
+ */
+ template <typename OtherNumber>
+ void add (const std::vector<unsigned int> &indices,
+ const Vector<OtherNumber> &values);
+
+ /**
+ * Take an address where
+ * <tt>n_elements</tt> are stored
+ * contiguously and add them into
+ * the vector. Handles all cases
+ * which are not covered by the
+ * other two <tt>add()</tt>
+ * functions above.
+ */
+ template <typename OtherNumber>
+ 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
(factor*boost::lambda::_1),
internal::Vector::minimum_parallel_grain_size);
}
+
+
+
+template <typename Number>
+template <typename Number2>
+inline
+void
+Vector<Number>::add (const std::vector<unsigned int> &indices,
+ const std::vector<Number2> &values)
+{
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(), values.size()));
+ add (indices.size(), &indices[0], &values[0]);
+}
+
+
+
+template <typename Number>
+template <typename Number2>
+inline
+void
+Vector<Number>::add (const std::vector<unsigned int> &indices,
+ const Vector<Number2> &values)
+{
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(), values.size()));
+ add (indices.size(), &indices[0], values.val);
+}
+
+
+
+template <typename Number>
+template <typename Number2>
+inline
+void
+Vector<Number>::add (const unsigned int n_indices,
+ const unsigned int *indices,
+ const Number2 *values)
+{
+ for (unsigned int i=0; i<n_indices; ++i)
+ {
+ Assert (indices[i] < vec_size, ExcIndexRange(indices[i],0,vec_size));
+ Assert (numbers::is_finite(values[i]),
+ ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+ val[indices[i]] += values[i];
+ }
+}
+
template <typename Number>