#include <lac/block_sparse_matrix.h>
#ifdef DEAL_II_USE_PETSC
+# include <lac/petsc_parallel_sparse_matrix.h>
# include <lac/petsc_sparse_matrix.h>
+# include <lac/petsc_parallel_vector.h>
# include <lac/petsc_vector.h>
#endif
{
first_nonzero_diagonal_entry = matrix.diag_element(i);
break;
- };
+ }
std::map<unsigned int,double>::const_iterator
}
- // 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
+ // set right hand side to wanted value;
+ // also store the new rhs entry to make
+ // the gauss step (when preserving the
// symmetry of the matrix) more
// efficient
//
// preset solution vector
solution(dof_number) = dof->second;
}
+
+ matrix.compress ();
+ solution.compress ();
+ right_hand_side.compress ();
}
+
+
+void
+MatrixTools::
+apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ PETScWrappers::MPI::SparseMatrix &matrix,
+ PETScWrappers::MPI::Vector &solution,
+ PETScWrappers::MPI::Vector &right_hand_side,
+ const bool preserve_symmetry)
+{
+ // this function works almost exactly as
+ // the one above, with the only exception
+ // that we have to make sure that we only
+ // access elements that belong to the slice
+ // of the matrix that is available locally
+ 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;
+
+
+ // we have to read and write from this
+ // matrix (in this order). this will only
+ // work if we compress the matrix first,
+ // done here:
+ matrix.compress ();
+
+ const unsigned int n_dofs = matrix.m();
+
+ // determine the first nonzero diagonal
+ // entry from within the part of the matrix
+ // that we can see. if we can't find such
+ // an entry, take one
+ PetscScalar first_nonzero_diagonal_entry = 1;
+ for (unsigned int i=matrix.local_range().first;
+ i<matrix.local_range().second; ++i)
+ if (matrix.diag_element(i) != 0)
+ {
+ first_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
+ break;
+ }
+
+ // iterate over all fixed degrees of
+ // freedom that are within the local range
+ // of this matrix
+ std::map<unsigned int,double>::const_iterator
+ dof = boundary_values.begin(),
+ endd = boundary_values.end();
+ for (; dof != endd; ++dof)
+ if ((dof->first >= matrix.local_range().first) &&
+ (dof->first < matrix.local_range().second))
+ {
+ 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;
+ // also 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;
+ }
+
+ matrix.compress ();
+ solution.compress ();
+ right_hand_side.compress ();
+}
+
+
#endif