From: wolf Date: Mon, 15 Mar 2004 22:53:17 +0000 (+0000) Subject: Implement boundary dof removal for petsc matrices. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=48ec24547aee4014d459f00a79cc2dab783fb510;p=dealii-svn.git Implement boundary dof removal for petsc matrices. git-svn-id: https://svn.dealii.org/trunk@8754 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index f982a0c8b5..e898ed0ab8 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -36,6 +36,13 @@ template class DoFHandler; template class MGDoFHandler; template class FEValues; +#ifdef DEAL_II_USE_PETSC +namespace PETScWrappers +{ + class MatrixBase; + class VectorBase; +} +#endif /** @@ -654,7 +661,7 @@ class MatrixCreator * modifying the matrix and right hand side vectors for boundary * values. * - * @author Wolfgang Bangerth, 1998, 2000 + * @author Wolfgang Bangerth, 1998, 2000, 2004 */ class MatrixTools : public MatrixCreator { @@ -699,6 +706,28 @@ class MatrixTools : public MatrixCreator BlockVector &right_hand_side, const bool eliminate_columns = true); + /** + * Apply dirichlet boundary conditions to + * the system matrix and vectors as + * described in the general + * documentation. This function works on + * the classes that are used to wrap + * PETSc objects. + * + * For a replacement function, + * see the documentation of the + * @ref{FilteredMatrix} class in + * the @p{LAC} sublibrary. + */ +#ifdef DEAL_II_USE_PETSC + static void + apply_boundary_values (const std::map &boundary_values, + PETScWrappers::MatrixBase &matrix, + PETScWrappers::VectorBase &solution, + PETScWrappers::VectorBase &right_hand_side, + const bool eliminate_columns = true); +#endif + /** * Exception */ diff --git a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc index cf1353c865..073447bc0f 100644 --- a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc @@ -18,6 +18,11 @@ #include #include +#ifdef DEAL_II_USE_PETSC +# include +# include +#endif + #include @@ -42,12 +47,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar return; - std::map::const_iterator dof = boundary_values.begin(), - endd = boundary_values.end(); - const unsigned int n_dofs = matrix.m(); - const SparsityPattern &sparsity = matrix.get_sparsity_pattern(); - const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices(); - const unsigned int *sparsity_colnums = sparsity.get_column_numbers(); + const unsigned int n_dofs = matrix.m(); // if a diagonal entry is zero // later, then we use another @@ -64,6 +64,11 @@ MatrixTools::apply_boundary_values (const std::map &boundar }; + std::map::const_iterator dof = boundary_values.begin(), + endd = boundary_values.end(); + const SparsityPattern &sparsity = matrix.get_sparsity_pattern(); + const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices(); + const unsigned int *sparsity_colnums = sparsity.get_column_numbers(); for (; dof != endd; ++dof) { Assert (dof->first < n_dofs, ExcInternalError()); @@ -187,8 +192,6 @@ MatrixTools::apply_boundary_values (const std::map &boundar - - void MatrixTools::apply_boundary_values (const std::map &boundary_values, BlockSparseMatrix &matrix, @@ -225,11 +228,7 @@ MatrixTools::apply_boundary_values (const std::map &boundar return; - std::map::const_iterator dof = boundary_values.begin(), - endd = boundary_values.end(); const unsigned int n_dofs = matrix.m(); - const BlockSparsityPattern & - sparsity_pattern = matrix.get_sparsity_pattern(); // if a diagonal entry is zero // later, then we use another @@ -258,6 +257,11 @@ MatrixTools::apply_boundary_values (const std::map &boundar if (first_nonzero_diagonal_entry == 0) first_nonzero_diagonal_entry = 1; + + std::map::const_iterator dof = boundary_values.begin(), + endd = boundary_values.end(); + const BlockSparsityPattern & + sparsity_pattern = matrix.get_sparsity_pattern(); // pointer to the mapping between // global and block indices. since @@ -509,6 +513,133 @@ MatrixTools::apply_boundary_values (const std::map &boundar +#ifdef DEAL_II_USE_PETSC + +void +MatrixTools:: +apply_boundary_values (const std::map &boundary_values, + PETScWrappers::MatrixBase &matrix, + PETScWrappers::VectorBase &solution, + PETScWrappers::VectorBase &right_hand_side, + const bool preserve_symmetry) +{ + Assert (matrix.n() == right_hand_side.size(), + ExcDimensionMismatch(matrix.n(), right_hand_side.size())); + Assert (matrix.n() == solution.size(), + ExcDimensionMismatch(matrix.n(), solution.size())); + + // if no boundary values are to be applied + // simply return + if (boundary_values.size() == 0) + return; + + + matrix.compress (); + + 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 + PetscScalar first_nonzero_diagonal_entry = 1; + for (unsigned int i=0; i::const_iterator + dof = boundary_values.begin(), + endd = boundary_values.end(); + for (; dof != endd; ++dof) + { + Assert (dof->first < n_dofs, ExcInternalError()); + + const unsigned int dof_number = dof->first; + + // for each constrained dof: + + // set entries of this line + // to zero except for the diagonal + // entry. + { + PETScWrappers::MatrixBase::const_iterator + p = matrix.begin(dof_number), + e = matrix.end(dof_number); + + // iterate over all elements of this + // row and set elements to zero + // except for the diagonal + // element. note that this is not + // exactly clean programming, since + // we change the matrix underneath, + // while we still keep working with + // the iterators into it + for (; p!=e; ++p) + if (p->column() != dof_number) + matrix.set (dof_number, p->column(), 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 (when preserving the + // symmetry of the matrix) more + // efficient + // + // note that for petsc matrices + // interleaving read with write + // operations is very expensive. thus, + // we here always replace the diagonal + // element, rather than first checking + // whether it is nonzero and in that + // case preserving it. this is + // different from the case of deal.II + // sparse matrices treated in the other + // functions. + PetscScalar new_rhs; + matrix.set (dof_number, dof_number, + first_nonzero_diagonal_entry); + new_rhs = dof->second * first_nonzero_diagonal_entry; + right_hand_side(dof_number) = new_rhs; + + // 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) + { + Assert (false, ExcNotImplemented()); + } + + // preset solution vector + solution(dof_number) = dof->second; + }; +} + +#endif + + + + // explicit instantiations