*
* \subsection{Boundary conditions}
*
- * The #apply_boundar_values# function inserts boundary conditions of
+ * The #apply_boundary_values# function inserts boundary conditions of
* into a system of equations. To actually do this you have to specify
- * a list of degree of freedom indices along with the value this degree of
+ * a list of degree of freedom indices along with the values these degrees of
* freedom shall assume. To see how to get such a list, see the discussion
* of the #VectorTools::interpolate_boundary_values# function.
*
* eliminating all coupling between this degree of freedom and others. Now
* also the column consists only of zeroes, apart from the main diagonal entry.
*
+ * Finding which rows contain an entry in the column for which we are
+ * presently performing a Gauss elimination step is either difficult
+ * or very simple, depending on the circumstances. If the sparsity
+ * pattern is symmetric (whether the matrix is symmetric is irrelevant
+ * here), then we can infer the rows which have a nonzero entry in the
+ * present column by looking at which columns in the present row are
+ * nonempty. In this case, we only need to look into a fixed number of
+ * rows and need not search all rows. On the other hand, if the
+ * sparsity pattern is nonsymmetric, then we need to use an iterative
+ * solver which can handle nonsymmetric matrices in any case, so there
+ * may be no need to do the Gauss elimination anyway. In fact, this is
+ * the way the function works: it takes a parameter
+ * (#elininate_columns#) that specifies whether the sparsity pattern
+ * is symmetric; if so, then the column is eliminated and the right
+ * hand side is also modified accordingly. If not, then only the row
+ * is deleted and the column is not touched at all, and all right hand
+ * side values apart from the one corresponding to the present row
+ * remain unchanged.
+ *
+ * If the sparsity pattern for your matrix is non-symmetric, you must
+ * set the value of this parameter to #false# in any case, since then
+ * we can't eliminate the column without searching all rows, which
+ * would be too expensive (if #N# be the number of rows, and #m# the
+ * number of nonzero elements per row, then eliminating one column is
+ * an #O(N*log(m))# operation, since searching in each row takes
+ * #log(m)# operations). If your sparsity pattern is symmetric, but
+ * your matrix is not, then you might specify #false# as well. If your
+ * sparsity pattern and matrix are both symmetric, you might want to
+ * specify #true# (the complexity of eliminating one row is then
+ * #O(m*log(m))#, since we only have to search #m# rows for the
+ * respective element of the column). Given the fact that #m# is
+ * roughly constant, irrespective of the discretization, and that the
+ * number of boundary nodes is #sqrt(N)# in 2d, the algorithm for
+ * symmetric sparsity patterns is #O(sqrt(N)*m*log(m))#, while it
+ * would be #O(N*sqrt(N)*m*log(m))# for the general case; the latter
+ * is too expensive to be performed.
+ *
* It seems as if we had to make clear not to overwrite the lines of other
* boundary nodes when doing the Gauss elimination step. However, since we
* reset the right hand side when passing such a node, it is not a problem
* the right hand side. We need therefore not take special care of other
* boundary nodes.
*
- * To make solving faster, we preset the solution vector with the right boundary
- * values. Since boundary nodes can never be hanging nodes, and since all other
- * entries of the solution vector are zero, we need not condense the solution
- * vector if the condensation process is done in-place. If done by copying
- * matrix and vectors to smaller ones, it would also be necessary to condense
- * the solution vector to preserve the preset boundary values.
- *
- * It it not clear whether the deletion of coupling between the boundary degree
- * of freedom and other dofs really forces the corresponding entry in the
- * solution vector to have the right value when using iterative solvers,
- * since their search directions may contain components in the direction
- * of the boundary node. For this reason, we perform a very simple line
- * balancing by not setting the main diagonal entry to unity, but rather
- * to the value it had before deleting this line, or to the first nonzero
- * main diagonal entry if it is zero from a previous Gauss elimination
- * step. Of course we have to change
- * the right hand side appropriately. This is not a very good
- * strategy, but it at least should give the main diagonal entry a value
- * in the right order of dimension, which makes the solving process a bit
- * more stable. A refined algorithm would set the entry to the mean of the
- * other diagonal entries, but this seems to be too expensive.
- *
- * Because of the mentioned question, whether or not a preset solution value
- * which does not couple with other degrees of freedom remains its value or
- * not during solving iteratively, it may or may not be necessary to set
- * the correct value after solving again. This question is an open one as of
- * now and may be answered by future experience.
+ * To make solving faster, we preset the solution vector with the
+ * right boundary values. It it not clear whether the deletion of
+ * coupling between the boundary degree of freedom and other dofs
+ * really forces the corresponding entry in the solution vector to
+ * have the right value when using iterative solvers, since their
+ * search directions may contain components in the direction of the
+ * boundary node. For this reason, we perform a very simple line
+ * balancing by not setting the main diagonal entry to unity, but
+ * rather to the value it had before deleting this line, or to the
+ * first nonzero main diagonal entry if it is zero for some reason.
+ * Of course we have to change the right hand side appropriately. This
+ * is not a very good strategy, but it at least should give the main
+ * diagonal entry a value in the right order of dimension, which makes
+ * the solvution process a bit more stable. A refined algorithm would
+ * set the entry to the mean of the other diagonal entries, but this
+ * seems to be too expensive.
*
*
* @author Wolfgang Bangerth, 1998
static void apply_boundary_values (const map<unsigned int,double> &boundary_values,
SparseMatrix<double> &matrix,
Vector<double> &solution,
- Vector<double> &right_hand_side);
+ Vector<double> &right_hand_side,
+ const bool eliminate_columns = true);
/**
* Exception
void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
- const Function<dim> * const a) {
+ const Function<dim> * const a)
+{
Vector<double> dummy; // no entries, should give an error if accessed
UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values);
if (a != 0)
};
+
template <int dim>
void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
const Function<dim> &rhs,
Vector<double> &rhs_vector,
- const Function<dim> * const a) {
+ const Function<dim> * const a)
+{
UpdateFlags update_flags = UpdateFlags(update_values |
update_q_points |
update_JxW_values);
};
+
template <int dim>
void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim> &dof,
- SparseMatrix<double> &matrix) {
+ SparseMatrix<double> &matrix)
+{
const FiniteElement<dim> &fe = dof.get_fe();
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const FunctionMap &,
Vector<double> &,
vector<unsigned int> &,
- const Function<1> *) {
+ const Function<1> *)
+{
Assert (false, ExcNotImplemented());
};
const FunctionMap &rhs,
Vector<double> &rhs_vector,
vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> *a) {
+ const Function<dim> *a)
+{
const FiniteElement<dim> &fe = dof.get_fe();
const unsigned int n_components = fe.n_components();
const bool fe_is_system = (n_components != 1);
};
+
template <int dim>
void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
- const Function<dim> * const a) {
+ const Function<dim> * const a)
+{
const unsigned int n_components = dof.get_fe().n_components();
Assert ((n_components==1) || (a==0), ExcNotImplemented());
while ((++assembler).state() == valid);
};
+
+
/*
template <int dim>
*/
+
+
template <int dim>
void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
};
+
template <int dim>
-void MatrixTools<dim>::apply_boundary_values (const map<unsigned int,double> &boundary_values,
- SparseMatrix<double> &matrix,
- Vector<double> &solution,
- Vector<double> &right_hand_side) {
+void
+MatrixTools<dim>::apply_boundary_values (const map<unsigned int,double> &boundary_values,
+ SparseMatrix<double> &matrix,
+ Vector<double> &solution,
+ Vector<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(),
new_rhs = right_hand_side(dof_number)
= dof->second * first_nonzero_diagonal_entry;
};
-
- // store the only nonzero entry
- // of this line for the Gauss
- // elimination step
- const double diagonal_entry = matrix.diag_element(dof_number);
- // do the Gauss step
- for (unsigned int row=0; row<n_dofs; ++row)
+
+ // 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)
{
- // we need not handle the
- // row we have already cleared
- if (row == dof_number)
- continue;
-
- // check whether the line has
- // an entry in the row corresponding
- // to the dof presently worked with.
- // note again: the first entry is
- // the diagonal entry which we
- // cannot be interested in; following
- // are the other entries in sorted
- // order, so we can use a binary
- // search
- //
- // if this row contains an element
- // for this dof, *p==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...)
- if ((*p == dof_number) &&
- (p != &sparsity_colnums[sparsity_rowstart[row+1]]))
- // this line has an entry
- // in the regarding column
+ // 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]]);
};
};
-
// preset solution vector
solution(dof_number) = dof->second;
};
const Function<dim> * const a) :
Equation<dim> (1),
right_hand_side (rhs),
- coefficient (a) {};
+ coefficient (a)
+{};
+
template <int dim>
void MassMatrix<dim>::assemble (FullMatrix<double> &cell_matrix,
const FEValues<dim> &fe_values,
- const typename DoFHandler<dim>::cell_iterator &) const {
+ const typename DoFHandler<dim>::cell_iterator &) const
+{
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
n_q_points = fe_values.n_quadrature_points;
const FiniteElement<dim> &fe = fe_values.get_fe();
};
+
template <int dim>
void MassMatrix<dim>::assemble (FullMatrix<double> &cell_matrix,
Vector<double> &rhs,
const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &) const {
+ const DoFHandler<dim>::cell_iterator &) const
+{
Assert (right_hand_side != 0, ExcNoRHSSelected());
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
};
+
template <int dim>
void MassMatrix<dim>::assemble (Vector<double> &rhs,
const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &) const {
+ const DoFHandler<dim>::cell_iterator &) const
+{
Assert (right_hand_side != 0, ExcNoRHSSelected());
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
};
+
template <int dim>
LaplaceMatrix<dim>::LaplaceMatrix (const Function<dim> * const rhs,
const Function<dim> * const a) :
coefficient (a) {};
+
template <int dim>
void LaplaceMatrix<dim>::assemble (FullMatrix<double> &cell_matrix,
Vector<double> &rhs,
const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &) const {
+ const DoFHandler<dim>::cell_iterator &) const
+{
Assert (right_hand_side != 0, ExcNoRHSSelected());
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
};
+
template <int dim>
void LaplaceMatrix<dim>::assemble (FullMatrix<double> &cell_matrix,
const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &) const {
+ const DoFHandler<dim>::cell_iterator &) const
+{
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
n_q_points = fe_values.n_quadrature_points;
};
+
template <int dim>
void LaplaceMatrix<dim>::assemble (Vector<double> &rhs,
const FEValues<dim> &fe_values,
- const DoFHandler<dim>::cell_iterator &) const {
+ const DoFHandler<dim>::cell_iterator &) const
+{
Assert (right_hand_side != 0, ExcNoRHSSelected());
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
};
+
template<int dim>
void
MatrixCreator<dim>::create_interpolation_matrix(const FiniteElement<dim> &high,
}
+
+// explicit instantiations
+
template class MatrixCreator<deal_II_dimension>;
template class MatrixTools<deal_II_dimension>;
template class MassMatrix<deal_II_dimension>;