#ifdef DEAL_II_WITH_PETSC
- namespace internal
+namespace internal
+{
+ namespace PETScWrappers
{
- namespace PETScWrappers
+ template <typename PETScMatrix, typename PETScVector>
+ void
+ apply_boundary_values (const std::map<types::global_dof_index,PetscScalar> &boundary_values,
+ PETScMatrix &matrix,
+ PETScVector &solution,
+ PETScVector &right_hand_side,
+ const bool eliminate_columns)
{
- template <typename PETScMatrix, typename PETScVector>
- void
- apply_boundary_values (const std::map<types::global_dof_index,PetscScalar> &boundary_values,
- PETScMatrix &matrix,
- PETScVector &solution,
- PETScVector &right_hand_side,
- const bool eliminate_columns)
- {
- (void)eliminate_columns;
- Assert (eliminate_columns == false, ExcNotImplemented());
-
- 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, then
- // jump straight to the compress() calls that we still have
- // to perform because they are collective operations
- if (boundary_values.size() > 0)
- {
- const std::pair<types::global_dof_index, types::global_dof_index> local_range
- = matrix.local_range();
- Assert (local_range == right_hand_side.local_range(),
- ExcInternalError());
- Assert (local_range == solution.local_range(),
- ExcInternalError());
-
- // 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 average_nonzero_diagonal_entry = 1;
- for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
- if (matrix.diag_element(i) != PetscScalar ())
- {
- average_nonzero_diagonal_entry = std::abs(matrix.diag_element(i));
- break;
- }
-
- // figure out which rows of the matrix we
- // have to eliminate on this processor
- std::vector<types::global_dof_index> constrained_rows;
- for (std::map<types::global_dof_index,PetscScalar>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- constrained_rows.push_back (dof->first);
-
- // then eliminate these rows and set
- // their diagonal entry to what we have
- // determined above. 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.
- matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
-
- std::vector<types::global_dof_index> indices;
- std::vector<PetscScalar> solution_values;
- for (std::map<types::global_dof_index,PetscScalar>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- {
- indices.push_back (dof->first);
- solution_values.push_back (dof->second);
- }
- solution.set (indices, solution_values);
-
- // now also set appropriate values for
- // the rhs
- for (unsigned int i=0; i<solution_values.size(); ++i)
- solution_values[i] *= average_nonzero_diagonal_entry;
-
- right_hand_side.set (indices, solution_values);
- }
- else
- {
- // clear_rows() is a collective operation so we still have to call
- // it:
- std::vector<types::global_dof_index> constrained_rows;
- matrix.clear_rows (constrained_rows, 1.);
- }
+ (void)eliminate_columns;
+ Assert (eliminate_columns == false, ExcNotImplemented());
+
+ 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, then
+ // jump straight to the compress() calls that we still have
+ // to perform because they are collective operations
+ if (boundary_values.size() > 0)
+ {
+ const std::pair<types::global_dof_index, types::global_dof_index> local_range
+ = matrix.local_range();
+ Assert (local_range == right_hand_side.local_range(),
+ ExcInternalError());
+ Assert (local_range == solution.local_range(),
+ ExcInternalError());
+
+ // 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 average_nonzero_diagonal_entry = 1;
+ for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
+ if (matrix.diag_element(i) != PetscScalar ())
+ {
+ average_nonzero_diagonal_entry = std::abs(matrix.diag_element(i));
+ break;
+ }
+
+ // figure out which rows of the matrix we
+ // have to eliminate on this processor
+ std::vector<types::global_dof_index> constrained_rows;
+ for (std::map<types::global_dof_index,PetscScalar>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ constrained_rows.push_back (dof->first);
+
+ // then eliminate these rows and set
+ // their diagonal entry to what we have
+ // determined above. 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.
+ matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
+
+ std::vector<types::global_dof_index> indices;
+ std::vector<PetscScalar> solution_values;
+ for (std::map<types::global_dof_index,PetscScalar>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ {
+ indices.push_back (dof->first);
+ solution_values.push_back (dof->second);
+ }
+ solution.set (indices, solution_values);
+
+ // now also set appropriate values for
+ // the rhs
+ for (unsigned int i=0; i<solution_values.size(); ++i)
+ solution_values[i] *= average_nonzero_diagonal_entry;
+
+ right_hand_side.set (indices, solution_values);
+ }
+ else
+ {
+ // clear_rows() is a collective operation so we still have to call
+ // it:
+ std::vector<types::global_dof_index> constrained_rows;
+ matrix.clear_rows (constrained_rows, 1.);
+ }
- // clean up
- solution.compress (VectorOperation::insert);
- right_hand_side.compress (VectorOperation::insert);
- }
+ // clean up
+ solution.compress (VectorOperation::insert);
+ right_hand_side.compress (VectorOperation::insert);
}
}
-
-
-
- void
- apply_boundary_values (const std::map<types::global_dof_index,PetscScalar> &boundary_values,
- PETScWrappers::SparseMatrix &matrix,
- PETScWrappers::Vector &solution,
- PETScWrappers::Vector &right_hand_side,
- const bool eliminate_columns)
- {
- // simply redirect to the generic function
- // used for both petsc matrix types
- internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
- right_hand_side, eliminate_columns);
- }
-
-
-
- void
- apply_boundary_values (const std::map<types::global_dof_index,PetscScalar> &boundary_values,
- PETScWrappers::MPI::SparseMatrix &matrix,
- PETScWrappers::MPI::Vector &solution,
- PETScWrappers::MPI::Vector &right_hand_side,
- const bool eliminate_columns)
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,PetscScalar> &boundary_values,
+ PETScWrappers::SparseMatrix &matrix,
+ PETScWrappers::Vector &solution,
+ PETScWrappers::Vector &right_hand_side,
+ const bool eliminate_columns)
+{
+ // simply redirect to the generic function
+ // used for both petsc matrix types
+ internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
+ right_hand_side, eliminate_columns);
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,PetscScalar> &boundary_values,
+ PETScWrappers::MPI::SparseMatrix &matrix,
+ PETScWrappers::MPI::Vector &solution,
+ PETScWrappers::MPI::Vector &right_hand_side,
+ const bool eliminate_columns)
+{
+ // simply redirect to the generic function
+ // used for both petsc matrix types
+ internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
+ right_hand_side, eliminate_columns);
+}
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,PetscScalar> &boundary_values,
+ PETScWrappers::MPI::BlockSparseMatrix &matrix,
+ PETScWrappers::MPI::BlockVector &solution,
+ PETScWrappers::MPI::BlockVector &right_hand_side,
+ const bool eliminate_columns)
+{
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionMismatch(matrix.n(), solution.size()));
+ Assert (matrix.n_block_rows() == matrix.n_block_cols(),
+ ExcNotQuadratic());
+
+ const unsigned int n_blocks = matrix.n_block_rows();
+
+ // We need to find the subdivision
+ // into blocks for the boundary values.
+ // To this end, generate a vector of
+ // maps with the respective indices.
+ std::vector<std::map<dealii::types::global_dof_index,PetscScalar> > block_boundary_values(n_blocks);
{
- // simply redirect to the generic function
- // used for both petsc matrix types
- internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
- right_hand_side, eliminate_columns);
+ int block = 0;
+ dealii::types::global_dof_index offset = 0;
+ for (std::map<types::global_dof_index,PetscScalar>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ {
+ if (dof->first >= matrix.block(block,0).m() + offset)
+ {
+ offset += matrix.block(block,0).m();
+ block++;
+ }
+ const types::global_dof_index index = dof->first - offset;
+ block_boundary_values[block].insert(std::pair<types::global_dof_index, PetscScalar> (index,dof->second));
+ }
}
-
- void
- apply_boundary_values (const std::map<types::global_dof_index,PetscScalar> &boundary_values,
- PETScWrappers::MPI::BlockSparseMatrix &matrix,
- PETScWrappers::MPI::BlockVector &solution,
- PETScWrappers::MPI::BlockVector &right_hand_side,
- const bool eliminate_columns)
- {
- Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
- Assert (matrix.n() == solution.size(),
- ExcDimensionMismatch(matrix.n(), solution.size()));
- Assert (matrix.n_block_rows() == matrix.n_block_cols(),
- ExcNotQuadratic());
-
- const unsigned int n_blocks = matrix.n_block_rows();
-
- // We need to find the subdivision
- // into blocks for the boundary values.
- // To this end, generate a vector of
- // maps with the respective indices.
- std::vector<std::map<dealii::types::global_dof_index,PetscScalar> > block_boundary_values(n_blocks);
+ // Now call the non-block variants on
+ // the diagonal subblocks and the
+ // solution/rhs.
+ for (unsigned int block=0; block<n_blocks; ++block)
+ internal::PETScWrappers::apply_boundary_values(block_boundary_values[block],
+ matrix.block(block,block),
+ solution.block(block),
+ right_hand_side.block(block),
+ eliminate_columns);
+
+ // Finally, we need to do something
+ // about the off-diagonal matrices. This
+ // is luckily not difficult. Just clear
+ // the whole row.
+ for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
{
- int block = 0;
- dealii::types::global_dof_index offset = 0;
+ const std::pair<types::global_dof_index, types::global_dof_index> local_range
+ = matrix.block(block_m,0).local_range();
+
+ std::vector<types::global_dof_index> constrained_rows;
for (std::map<types::global_dof_index,PetscScalar>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
+ dof = block_boundary_values[block_m].begin();
+ dof != block_boundary_values[block_m].end();
++dof)
- {
- if (dof->first >= matrix.block(block,0).m() + offset)
- {
- offset += matrix.block(block,0).m();
- block++;
- }
- const types::global_dof_index index = dof->first - offset;
- block_boundary_values[block].insert(std::pair<types::global_dof_index, PetscScalar> (index,dof->second));
- }
- }
-
- // Now call the non-block variants on
- // the diagonal subblocks and the
- // solution/rhs.
- for (unsigned int block=0; block<n_blocks; ++block)
- internal::PETScWrappers::apply_boundary_values(block_boundary_values[block],
- matrix.block(block,block),
- solution.block(block),
- right_hand_side.block(block),
- eliminate_columns);
-
- // Finally, we need to do something
- // about the off-diagonal matrices. This
- // is luckily not difficult. Just clear
- // the whole row.
- for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
- {
- const std::pair<types::global_dof_index, types::global_dof_index> local_range
- = matrix.block(block_m,0).local_range();
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ constrained_rows.push_back (dof->first);
- std::vector<types::global_dof_index> constrained_rows;
- for (std::map<types::global_dof_index,PetscScalar>::const_iterator
- dof = block_boundary_values[block_m].begin();
- dof != block_boundary_values[block_m].end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- constrained_rows.push_back (dof->first);
-
- for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
- if (block_m != block_n)
- matrix.block(block_m,block_n).clear_rows(constrained_rows);
- }
- }
+ for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
+ if (block_m != block_n)
+ matrix.block(block_m,block_n).clear_rows(constrained_rows);
+ }
+}
#endif
#ifdef DEAL_II_WITH_TRILINOS
- namespace internal
+namespace internal
+{
+ namespace TrilinosWrappers
{
- namespace TrilinosWrappers
+ template <typename TrilinosMatrix, typename TrilinosVector>
+ void
+ apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
+ TrilinosMatrix &matrix,
+ TrilinosVector &solution,
+ TrilinosVector &right_hand_side,
+ const bool eliminate_columns)
{
- template <typename TrilinosMatrix, typename TrilinosVector>
- void
- apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
- TrilinosMatrix &matrix,
- TrilinosVector &solution,
- TrilinosVector &right_hand_side,
- const bool eliminate_columns)
- {
- Assert (eliminate_columns == false, ExcNotImplemented());
- (void)eliminate_columns;
-
- Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
- Assert (matrix.n() == solution.size(),
- ExcDimensionMismatch(matrix.m(), solution.size()));
-
- // if no boundary values are to be applied, then
- // jump straight to the compress() calls that we still have
- // to perform because they are collective operations
- if (boundary_values.size() > 0)
- {
- const std::pair<types::global_dof_index, types::global_dof_index> local_range
- = matrix.local_range();
- Assert (local_range == right_hand_side.local_range(),
- ExcInternalError());
- Assert (local_range == solution.local_range(),
- ExcInternalError());
-
- // 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
- TrilinosScalar average_nonzero_diagonal_entry = 1;
- for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
- if (matrix.diag_element(i) != 0)
- {
- average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
- break;
- }
-
- // figure out which rows of the matrix we
- // have to eliminate on this processor
- std::vector<types::global_dof_index> constrained_rows;
- for (std::map<types::global_dof_index,TrilinosScalar>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- constrained_rows.push_back (dof->first);
-
- // then eliminate these rows and
- // set their diagonal entry to
- // what we have determined
- // above. if the value already is
- // nonzero, it will be preserved,
- // in accordance with the basic
- // matrix classes in deal.II.
- matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
-
- std::vector<types::global_dof_index> indices;
- std::vector<TrilinosScalar> solution_values;
- for (std::map<types::global_dof_index,TrilinosScalar>::const_iterator
- dof = boundary_values.begin();
- dof != boundary_values.end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- {
- indices.push_back (dof->first);
- solution_values.push_back (dof->second);
- }
- solution.set (indices, solution_values);
-
- // now also set appropriate
- // values for the rhs
- for (unsigned int i=0; i<solution_values.size(); ++i)
- solution_values[i] *= matrix.diag_element(indices[i]);
-
- right_hand_side.set (indices, solution_values);
- }
- else
- {
- // clear_rows() is a collective operation so we still have to call
- // it:
- std::vector<types::global_dof_index> constrained_rows;
- matrix.clear_rows (constrained_rows, 1.);
- }
-
- // clean up
- matrix.compress (VectorOperation::insert);
- solution.compress (VectorOperation::insert);
- right_hand_side.compress (VectorOperation::insert);
- }
-
-
-
- template <typename TrilinosMatrix, typename TrilinosBlockVector>
- void
- apply_block_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
- TrilinosMatrix &matrix,
- TrilinosBlockVector &solution,
- TrilinosBlockVector &right_hand_side,
- const bool eliminate_columns)
- {
- Assert (eliminate_columns == false, ExcNotImplemented());
-
- Assert (matrix.n() == right_hand_side.size(),
- ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
- Assert (matrix.n() == solution.size(),
- ExcDimensionMismatch(matrix.n(), solution.size()));
- Assert (matrix.n_block_rows() == matrix.n_block_cols(),
- ExcNotQuadratic());
-
- const unsigned int n_blocks = matrix.n_block_rows();
-
- // We need to find the subdivision
- // into blocks for the boundary values.
- // To this end, generate a vector of
- // maps with the respective indices.
- std::vector<std::map<types::global_dof_index,TrilinosScalar> > block_boundary_values(n_blocks);
+ Assert (eliminate_columns == false, ExcNotImplemented());
+ (void)eliminate_columns;
+
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionMismatch(matrix.m(), solution.size()));
+
+ // if no boundary values are to be applied, then
+ // jump straight to the compress() calls that we still have
+ // to perform because they are collective operations
+ if (boundary_values.size() > 0)
{
- int block=0;
- types::global_dof_index offset = 0;
+ const std::pair<types::global_dof_index, types::global_dof_index> local_range
+ = matrix.local_range();
+ Assert (local_range == right_hand_side.local_range(),
+ ExcInternalError());
+ Assert (local_range == solution.local_range(),
+ ExcInternalError());
+
+ // 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
+ TrilinosScalar average_nonzero_diagonal_entry = 1;
+ for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
+ if (matrix.diag_element(i) != 0)
+ {
+ average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
+ break;
+ }
+
+ // figure out which rows of the matrix we
+ // have to eliminate on this processor
+ std::vector<types::global_dof_index> constrained_rows;
+ for (std::map<types::global_dof_index,TrilinosScalar>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ constrained_rows.push_back (dof->first);
+
+ // then eliminate these rows and
+ // set their diagonal entry to
+ // what we have determined
+ // above. if the value already is
+ // nonzero, it will be preserved,
+ // in accordance with the basic
+ // matrix classes in deal.II.
+ matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
+
+ std::vector<types::global_dof_index> indices;
+ std::vector<TrilinosScalar> solution_values;
for (std::map<types::global_dof_index,TrilinosScalar>::const_iterator
dof = boundary_values.begin();
dof != boundary_values.end();
++dof)
- {
- if (dof->first >= matrix.block(block,0).m() + offset)
- {
- offset += matrix.block(block,0).m();
- block++;
- }
- const types::global_dof_index index = dof->first - offset;
- block_boundary_values[block].insert(
- std::pair<types::global_dof_index, TrilinosScalar> (index,dof->second));
- }
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ {
+ indices.push_back (dof->first);
+ solution_values.push_back (dof->second);
+ }
+ solution.set (indices, solution_values);
+
+ // now also set appropriate
+ // values for the rhs
+ for (unsigned int i=0; i<solution_values.size(); ++i)
+ solution_values[i] *= matrix.diag_element(indices[i]);
+
+ right_hand_side.set (indices, solution_values);
+ }
+ else
+ {
+ // clear_rows() is a collective operation so we still have to call
+ // it:
+ std::vector<types::global_dof_index> constrained_rows;
+ matrix.clear_rows (constrained_rows, 1.);
}
- // Now call the non-block variants on
- // the diagonal subblocks and the
- // solution/rhs.
- for (unsigned int block=0; block<n_blocks; ++block)
- TrilinosWrappers::apply_boundary_values(block_boundary_values[block],
- matrix.block(block,block),
- solution.block(block),
- right_hand_side.block(block),
- eliminate_columns);
-
- // Finally, we need to do something
- // about the off-diagonal matrices. This
- // is luckily not difficult. Just clear
- // the whole row.
- for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
- {
- const std::pair<types::global_dof_index, types::global_dof_index> local_range
- = matrix.block(block_m,0).local_range();
-
- std::vector<types::global_dof_index> constrained_rows;
- for (std::map<types::global_dof_index,TrilinosScalar>::const_iterator
- dof = block_boundary_values[block_m].begin();
- dof != block_boundary_values[block_m].end();
- ++dof)
- if ((dof->first >= local_range.first) &&
- (dof->first < local_range.second))
- constrained_rows.push_back (dof->first);
-
- for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
- if (block_m != block_n)
- matrix.block(block_m,block_n).clear_rows(constrained_rows);
- }
- }
+ // clean up
+ matrix.compress (VectorOperation::insert);
+ solution.compress (VectorOperation::insert);
+ right_hand_side.compress (VectorOperation::insert);
}
- }
-
-
- void
- apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
- TrilinosWrappers::SparseMatrix &matrix,
- TrilinosWrappers::Vector &solution,
- TrilinosWrappers::Vector &right_hand_side,
- const bool eliminate_columns)
- {
- // simply redirect to the generic function
- // used for both trilinos matrix types
- internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
- right_hand_side, eliminate_columns);
- }
+ template <typename TrilinosMatrix, typename TrilinosBlockVector>
+ void
+ apply_block_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
+ TrilinosMatrix &matrix,
+ TrilinosBlockVector &solution,
+ TrilinosBlockVector &right_hand_side,
+ const bool eliminate_columns)
+ {
+ Assert (eliminate_columns == false, ExcNotImplemented());
+
+ Assert (matrix.n() == right_hand_side.size(),
+ ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+ Assert (matrix.n() == solution.size(),
+ ExcDimensionMismatch(matrix.n(), solution.size()));
+ Assert (matrix.n_block_rows() == matrix.n_block_cols(),
+ ExcNotQuadratic());
+
+ const unsigned int n_blocks = matrix.n_block_rows();
+
+ // We need to find the subdivision
+ // into blocks for the boundary values.
+ // To this end, generate a vector of
+ // maps with the respective indices.
+ std::vector<std::map<types::global_dof_index,TrilinosScalar> > block_boundary_values(n_blocks);
+ {
+ int block=0;
+ types::global_dof_index offset = 0;
+ for (std::map<types::global_dof_index,TrilinosScalar>::const_iterator
+ dof = boundary_values.begin();
+ dof != boundary_values.end();
+ ++dof)
+ {
+ if (dof->first >= matrix.block(block,0).m() + offset)
+ {
+ offset += matrix.block(block,0).m();
+ block++;
+ }
+ const types::global_dof_index index = dof->first - offset;
+ block_boundary_values[block].insert(
+ std::pair<types::global_dof_index, TrilinosScalar> (index,dof->second));
+ }
+ }
+ // Now call the non-block variants on
+ // the diagonal subblocks and the
+ // solution/rhs.
+ for (unsigned int block=0; block<n_blocks; ++block)
+ TrilinosWrappers::apply_boundary_values(block_boundary_values[block],
+ matrix.block(block,block),
+ solution.block(block),
+ right_hand_side.block(block),
+ eliminate_columns);
+
+ // Finally, we need to do something
+ // about the off-diagonal matrices. This
+ // is luckily not difficult. Just clear
+ // the whole row.
+ for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
+ {
+ const std::pair<types::global_dof_index, types::global_dof_index> local_range
+ = matrix.block(block_m,0).local_range();
- void
- apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
- TrilinosWrappers::SparseMatrix &matrix,
- TrilinosWrappers::MPI::Vector &solution,
- TrilinosWrappers::MPI::Vector &right_hand_side,
- const bool eliminate_columns)
- {
- // simply redirect to the generic function
- // used for both trilinos matrix types
- internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
- right_hand_side, eliminate_columns);
- }
-
-
-
- void
- apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
- TrilinosWrappers::BlockSparseMatrix &matrix,
- TrilinosWrappers::BlockVector &solution,
- TrilinosWrappers::BlockVector &right_hand_side,
- const bool eliminate_columns)
- {
- internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
- solution, right_hand_side,
- eliminate_columns);
- }
-
-
+ std::vector<types::global_dof_index> constrained_rows;
+ for (std::map<types::global_dof_index,TrilinosScalar>::const_iterator
+ dof = block_boundary_values[block_m].begin();
+ dof != block_boundary_values[block_m].end();
+ ++dof)
+ if ((dof->first >= local_range.first) &&
+ (dof->first < local_range.second))
+ constrained_rows.push_back (dof->first);
- void
- apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
- TrilinosWrappers::BlockSparseMatrix &matrix,
- TrilinosWrappers::MPI::BlockVector &solution,
- TrilinosWrappers::MPI::BlockVector &right_hand_side,
- const bool eliminate_columns)
- {
- internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
- solution, right_hand_side,
- eliminate_columns);
+ for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
+ if (block_m != block_n)
+ matrix.block(block_m,block_n).clear_rows(constrained_rows);
+ }
+ }
}
+}
+
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
+ TrilinosWrappers::SparseMatrix &matrix,
+ TrilinosWrappers::Vector &solution,
+ TrilinosWrappers::Vector &right_hand_side,
+ const bool eliminate_columns)
+{
+ // simply redirect to the generic function
+ // used for both trilinos matrix types
+ internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
+ right_hand_side, eliminate_columns);
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
+ TrilinosWrappers::SparseMatrix &matrix,
+ TrilinosWrappers::MPI::Vector &solution,
+ TrilinosWrappers::MPI::Vector &right_hand_side,
+ const bool eliminate_columns)
+{
+ // simply redirect to the generic function
+ // used for both trilinos matrix types
+ internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
+ right_hand_side, eliminate_columns);
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
+ TrilinosWrappers::BlockSparseMatrix &matrix,
+ TrilinosWrappers::BlockVector &solution,
+ TrilinosWrappers::BlockVector &right_hand_side,
+ const bool eliminate_columns)
+{
+ internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
+ solution, right_hand_side,
+ eliminate_columns);
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,TrilinosScalar> &boundary_values,
+ TrilinosWrappers::BlockSparseMatrix &matrix,
+ TrilinosWrappers::MPI::BlockVector &solution,
+ TrilinosWrappers::MPI::BlockVector &right_hand_side,
+ const bool eliminate_columns)
+{
+ internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
+ solution, right_hand_side,
+ eliminate_columns);
+}
#endif
projects / dealii.git / commitdiff