From f8d972d6e7d2831c365ceb7fa8183f34c933cea3 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 17 May 2000 14:30:39 +0000 Subject: [PATCH] Add apply_boundary_values for block matrices and vectors. At present not implemented, but I check it in to work at it at home. git-svn-id: https://svn.dealii.org/trunk@2880 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/numerics/matrices.cc | 169 ++++++++++++++++++++ 1 file changed, 169 insertions(+) diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index 4d26689880..168e7f2ee7 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -675,6 +675,175 @@ MatrixTools::apply_boundary_values (const map &boundar }; + +/* + +template +template +void +MatrixTools::apply_boundary_values (const map &boundary_values, + BlockSparseMatrix &matrix, + BlockVector &solution, + BlockVector &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::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; ifirst < 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; jsecond * 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; jsecond; + }; +}; + + +*/ + + + template MassMatrix::MassMatrix (const Function * const rhs, const Function * const a) : -- 2.39.5