* Same as the previous function, except that it uses two (possibly) different
* index sets to correctly handle inhomogeneties when the local matrix is
* computed from a combination of two neighboring elements, for example for an
- * edge integral term in dG. Note that in the case that these two elements have
+ * edge integral term in DG. Note that in the case that these two elements have
* different polynomial degree, the local matrix is rectangular.
*
* <tt>local_dof_indices_row<tt> is the set of row indices and
* <tt>local_dof_indices_col<tt> is the set of column indices of the local matrix.
+ * <tt>diagonal=false<tt> says whether the two index sets are equal or not.
+ *
+ * If both index sets are equal, <tt>diagonal<tt> must be set to true or we
+ * simply use the previous function. If both index sets are different (diagonal=false)
+ * the <tt>global_vector<tt> is modified to handle inhomogeneties but no
+ * entries from <tt>local_vector<tt> are added. Note that the edge integrals for inner
+ * edged for DG do not contribute any values to the right hand side.
*/
template <typename VectorType, typename LocalType>
void
const std::vector<size_type> &local_dof_indices_row,
const std::vector<size_type> &local_dof_indices_col,
VectorType &global_vector,
- const FullMatrix<LocalType> &local_matrix) const;
+ const FullMatrix<LocalType> &local_matrix,
+ bool diagonal = false) const;
/**
* Enter a single value into a result vector, obeying constraints.
VectorType &global_vector,
const FullMatrix<LocalType> &local_matrix) const
{
- distribute_local_to_global(local_vector,local_dof_indices,local_dof_indices, global_vector, local_matrix);
+ distribute_local_to_global(local_vector,local_dof_indices,local_dof_indices, global_vector, local_matrix, true);
}
const std::vector<size_type> &local_dof_indices_row,
const std::vector<size_type> &local_dof_indices_col,
VectorType &global_vector,
- const FullMatrix<LocalType> &local_matrix) const
+ const FullMatrix<LocalType> &local_matrix,
+ bool diagonal) const
{
Assert (sorted == true, ExcMatrixNotClosed());
- AssertDimension (local_vector.size(), local_dof_indices_col.size());
+ AssertDimension (local_vector.size(), local_dof_indices_row.size());
AssertDimension (local_matrix.m(), local_dof_indices_row.size());
AssertDimension (local_matrix.n(), local_dof_indices_col.size());
-
- const size_type n_local_dofs = local_vector.size();
+
+ // diagonal checks if we have only one index set (if both are equal
+ // diagonal should be set to true).
+ // If true we do both, assemply of the right hand side (next lines)
+ // and (see further below) modifications of the right hand side
+ // according to the inhomogeneous constraints.
+ // Otherwise we only modify the right hand side according to
+ // local_matrix and the inhomogeneos constraints, and omit the vector add.
+
+ const size_type m_local_dofs = local_dof_indices_row.size();
+ const size_type n_local_dofs = local_dof_indices_col.size();
if (lines.empty())
- global_vector.add(local_dof_indices_col, local_vector);
+ {
+ if (diagonal)
+ global_vector.add(local_dof_indices_row, local_vector);
+ }
else
for (size_type i=0; i<n_local_dofs; ++i)
{
// the constraint
if (is_constrained(local_dof_indices_col[i]) == false)
{
- global_vector(local_dof_indices_col[i]) += local_vector(i);
+ if (diagonal)
+ global_vector(local_dof_indices_row[i]) += local_vector(i);
continue;
}
// constrained. If so, distribute the constraint
const double val = position->inhomogeneity;
if (val != 0)
- for (size_type j=0; j<n_local_dofs; ++j)
+ for (size_type j=0; j<m_local_dofs; ++j)
{
if (is_constrained(local_dof_indices_row[j]) == false)
{
// now distribute the constraint,
// but make sure we don't touch
// the entries of fixed dofs
- for (size_type j=0; j<position->entries.size(); ++j)
+ if (diagonal)
{
- Assert (!(!local_lines.size()
+ for (size_type j=0; j<position->entries.size(); ++j)
+ {
+ Assert (!(!local_lines.size()
|| local_lines.is_element(position->entries[j].first))
- || 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;
+ || 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;
+ }
}
}
}
if (ResidualSimple<VectorType>::constraints == 0)
{
- for (unsigned int i=0; i<i1.size(); ++i)
- (*v)(i1[i]) += vector(i);
-
for (unsigned int j=0; j<i1.size(); ++j)
for (unsigned int k=0; k<i2.size(); ++k)
if (std::fabs(M(j,k)) >= MatrixSimple<MatrixType>::threshold)
}
else
{
- ResidualSimple<VectorType>::constraints->distribute_local_to_global(vector, i1, i2, *v, M);
+ ResidualSimple<VectorType>::constraints->distribute_local_to_global(vector, i1, i2, *v, M, false);
ResidualSimple<VectorType>::constraints->distribute_local_to_global(M, i1, i2, *MatrixSimple<MatrixType>::matrix[index]);
}
}
if (row == column)
assemble(info.matrix(k+m*n*n,false).matrix,
- info.vector(m).block(column), m,
- info.indices_by_block[column]);
+ info.vector(m).block(row), m,
+ info.indices_by_block[row]);
else
assemble(info.matrix(k+m*n*n,false).matrix,
- info.vector(m).block(column), m,
+ info.vector(m).block(row), m,
info.indices_by_block[row],
info.indices_by_block[column]);
}
if (row == column)
{
- assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(column), m,
- info1.indices_by_block[column]);
- assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(column), m,
- info2.indices_by_block[column]);
+ assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(row), m,
+ info1.indices_by_block[row]);
+ assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(row), m,
+ info2.indices_by_block[row]);
}
else
{
- assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(column), m,
+ assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(row), m,
info1.indices_by_block[row], info1.indices_by_block[column]);
- assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(column), m,
+ assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(row), m,
info2.indices_by_block[row], info2.indices_by_block[column]);
}
- assemble(info1.matrix(k+m*n*n,true).matrix, info1.vector(m).block(column), m,
+ assemble(info1.matrix(k+m*n*n,true).matrix, info1.vector(m).block(row), m,
info1.indices_by_block[row], info2.indices_by_block[column]);
- assemble(info2.matrix(k+m*n*n,true).matrix, info2.vector(m).block(column), m,
+ assemble(info2.matrix(k+m*n*n,true).matrix, info2.vector(m).block(row), m,
info2.indices_by_block[row], info1.indices_by_block[column]);
}
}
const std::vector<ConstraintMatrix::size_type> &, \
const std::vector<ConstraintMatrix::size_type> &, \
VectorType &, \
- const FullMatrix<VectorType::value_type> &) const
+ const FullMatrix<VectorType::value_type> &, \
+ bool) const
#define PARALLEL_VECTOR_FUNCTIONS(VectorType) \
template void ConstraintMatrix:: \
const std::vector<ConstraintMatrix::size_type> &, \
const std::vector<ConstraintMatrix::size_type> &, \
VectorType &, \
- const FullMatrix<VectorType::value_type> &) const
-
+ const FullMatrix<VectorType::value_type> &, \
+ bool) const
#ifdef DEAL_II_WITH_PETSC
VECTOR_FUNCTIONS(PETScWrappers::MPI::Vector);
template void ConstraintMatrix::distribute_local_to_global<T<S> > (
const Vector<S>&, const std::vector<types::global_dof_index> &, T<S> &, const FullMatrix<S>&) const;
template void ConstraintMatrix::distribute_local_to_global<T<S> > (
- const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, T<S> &, const FullMatrix<S>&) const;
+ const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, T<S> &, const FullMatrix<S>&, bool) const;
template void ConstraintMatrix::set_zero<T<S> >(T<S> &) const;
}
template void ConstraintMatrix::distribute_local_to_global<parallel::distributed::T<S> > (
const Vector<S>&, const std::vector<types::global_dof_index> &, parallel::distributed::T<S> &, const FullMatrix<S>&) const;
template void ConstraintMatrix::distribute_local_to_global<parallel::distributed::T<S> > (
- const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, parallel::distributed::T<S> &, const FullMatrix<S>&) const;
+ const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, parallel::distributed::T<S> &, const FullMatrix<S>&, bool) const;
template void ConstraintMatrix::set_zero<parallel::distributed::T<S> >(parallel::distributed::T<S> &) const;
}
template void ConstraintMatrix::distribute_local_to_global<V > (
const Vector<V::value_type>&, const std::vector<types::global_dof_index> &, V&, const FullMatrix<V::value_type>&) const;
template void ConstraintMatrix::distribute_local_to_global<V > (
- const Vector<V::value_type>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, V&, const FullMatrix<V::value_type>&) const;
+ const Vector<V::value_type>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, V&, const FullMatrix<V::value_type>&, bool) const;
template void ConstraintMatrix::set_zero<V >(V&) const;
}