};
+
+/*
+
+template <int dim>
+template <int blocks>
+void
+MatrixTools<dim>::apply_boundary_values (const map<unsigned int,double> &boundary_values,
+ BlockSparseMatrix<double,blocks,blocks> &matrix,
+ BlockVector<blocks,double> &solution,
+ BlockVector<blocks,double> &right_hand_side,
+ const bool preserve_symmetry)
+{
+ Assert (matrix.n() == matrix.m(),
+ ExcDimensionsDontMatch(matrix.n(), matrix.m()));
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionsDontMatch(matrix.n(), solution.size()));
+ Assert (matrix.get_row_indices() == matrix.get_column_indices(),
+ ExcMatrixNotBlockSquare());
+
+ // if no boundary values are to be applied
+ // simply return
+ if (boundary_values.size() == 0)
+ return;
+
+
+ map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
+ endd = boundary_values.end();
+ const unsigned int n_dofs = matrix.m();
+
+ // if a diagonal entry is zero
+ // later, then we use another
+ // number instead. take it to be
+ // the first nonzero diagonal
+ // element of the matrix, or 1 if
+ // there is no such thing
+ double first_nonzero_diagonal_entry = 1;
+ for (unsigned int i=0; i<n_dofs; ++i)
+ if (matrix.diag_element(i) != 0)
+ {
+ first_nonzero_diagonal_entry = matrix.diag_element(i);
+ break;
+ };
+
+
+ for (; dof != endd; ++dof)
+ {
+ Assert (dof->first < n_dofs, ExcInternalError());
+
+ const unsigned int dof_number = dof->first;
+ // for each boundary dof:
+
+ // set entries of this line
+ // to zero except for the diagonal
+ // entry. Note that the diagonal
+ // entry is always the first one
+ // for square matrices, i.e.
+ // we shall not set
+ // matrix.global_entry(
+ // sparsity_rowstart[dof.first])
+ const unsigned int last = sparsity_rowstart[dof_number+1];
+ for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+ matrix.global_entry(j) = 0.;
+
+
+ // set right hand side to
+ // wanted value: if main diagonal
+ // entry nonzero, don't touch it
+ // and scale rhs accordingly. If
+ // zero, take the first main
+ // diagonal entry we can find, or
+ // one if no nonzero main diagonal
+ // element exists. Normally, however,
+ // the main diagonal entry should
+ // not be zero.
+ //
+ // store the new rhs entry to make
+ // the gauss step more efficient
+ double new_rhs;
+ if (matrix.diag_element(dof_number) != 0.0)
+ new_rhs = right_hand_side(dof_number)
+ = dof->second * matrix.diag_element(dof_number);
+ else
+ {
+ matrix.set (dof_number, dof_number,
+ first_nonzero_diagonal_entry);
+ new_rhs = right_hand_side(dof_number)
+ = dof->second * first_nonzero_diagonal_entry;
+ };
+
+
+ // if the user wants to have
+ // the symmetry of the matrix
+ // preserved, and if the
+ // sparsity pattern is
+ // symmetric, then do a Gauss
+ // elimination step with the
+ // present row
+ if (preserve_symmetry)
+ {
+ // store the only nonzero entry
+ // of this line for the Gauss
+ // elimination step
+ const double diagonal_entry = matrix.diag_element(dof_number);
+
+ // we have to loop over all
+ // rows of the matrix which
+ // have a nonzero entry in
+ // the column which we work
+ // in presently. if the
+ // sparsity pattern is
+ // symmetric, then we can
+ // get the positions of
+ // these rows cheaply by
+ // looking at the nonzero
+ // column numbers of the
+ // present row. we need not
+ // look at the first entry,
+ // since that is the
+ // diagonal element and
+ // thus the present row
+ for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+ {
+ const unsigned int row = sparsity_colnums[j];
+
+ // find the position of
+ // element
+ // (row,dof_number)
+ const unsigned int *
+ p = lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
+ &sparsity_colnums[sparsity_rowstart[row+1]],
+ dof_number);
+
+ // check whether this line has
+ // an entry in the regarding column
+ // (check for ==dof_number and
+ // != next_row, since if
+ // row==dof_number-1, *p is a
+ // past-the-end pointer but points
+ // to dof_number anyway...)
+ //
+ // there should be such an entry!
+ Assert ((*p == dof_number) &&
+ (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
+ ExcInternalError());
+
+ const unsigned int global_entry
+ = (p - &sparsity_colnums[sparsity_rowstart[0]]);
+
+ // correct right hand side
+ right_hand_side(row) -= matrix.global_entry(global_entry) /
+ diagonal_entry * new_rhs;
+
+ // set matrix entry to zero
+ matrix.global_entry(global_entry) = 0.;
+ };
+ };
+
+ // preset solution vector
+ solution(dof_number) = dof->second;
+ };
+};
+
+
+*/
+
+
+
template <int dim>
MassMatrix<dim>::MassMatrix (const Function<dim> * const rhs,
const Function<dim> * const a) :