From a02a5f7733e37cc720454b0937c36d82d539c5dc Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Tue, 1 Mar 2016 11:48:16 -0500 Subject: [PATCH] wrap correctly in namespace fixes missing apply_boundary_values() implementation error --- source/numerics/matrix_tools_once.cc | 805 ++++++++++++++------------- 1 file changed, 404 insertions(+), 401 deletions(-) diff --git a/source/numerics/matrix_tools_once.cc b/source/numerics/matrix_tools_once.cc index d9bd2d1f18..d22fffbcf7 100644 --- a/source/numerics/matrix_tools_once.cc +++ b/source/numerics/matrix_tools_once.cc @@ -49,9 +49,6 @@ # include #endif -#include - - #include #include #include @@ -59,215 +56,220 @@ DEAL_II_NAMESPACE_OPEN +namespace MatrixTools +{ + #ifdef DEAL_II_WITH_PETSC -namespace internal -{ - namespace PETScWrappers + namespace internal { - template - void - apply_boundary_values (const std::map &boundary_values, - PETScMatrix &matrix, - PETScVector &solution, - PETScVector &right_hand_side, - const bool eliminate_columns) + namespace PETScWrappers { - (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 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 constrained_rows; - for (std::map::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 indices; - std::vector solution_values; - for (std::map::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 constrained_rows; - matrix.clear_rows (constrained_rows, 1.); - } + template + void + apply_boundary_values (const std::map &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 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 constrained_rows; + for (std::map::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 indices; + std::vector solution_values; + for (std::map::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 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 &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 &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 &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 &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 > block_boundary_values(n_blocks); + void + apply_boundary_values (const std::map &boundary_values, + PETScWrappers::MPI::SparseMatrix &matrix, + PETScWrappers::MPI::Vector &solution, + PETScWrappers::MPI::Vector &right_hand_side, + const bool eliminate_columns) { - int block = 0; - dealii::types::global_dof_index offset = 0; - for (std::map::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 (index,dof->second)); - } + // 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); } - // Now call the non-block variants on - // the diagonal subblocks and the - // solution/rhs. - for (unsigned int block=0; block local_range - = matrix.block(block_m,0).local_range(); - std::vector constrained_rows; + void + apply_boundary_values (const std::map &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 > block_boundary_values(n_blocks); + { + int block = 0; + dealii::types::global_dof_index offset = 0; for (std::map::const_iterator - dof = block_boundary_values[block_m].begin(); - dof != block_boundary_values[block_m].end(); + 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); - - for (unsigned int block_n=0; block_nfirst >= 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 (index,dof->second)); + } } -} + + // Now call the non-block variants on + // the diagonal subblocks and the + // solution/rhs. + for (unsigned int block=0; block local_range + = matrix.block(block_m,0).local_range(); + + std::vector constrained_rows; + for (std::map::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 &boundary #ifdef DEAL_II_WITH_TRILINOS -namespace internal -{ - namespace TrilinosWrappers + namespace internal { - template - void - apply_boundary_values (const std::map &boundary_values, - TrilinosMatrix &matrix, - TrilinosVector &solution, - TrilinosVector &right_hand_side, - const bool eliminate_columns) + namespace TrilinosWrappers { - 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) + template + void + apply_boundary_values (const std::map &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 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 constrained_rows; + for (std::map::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 indices; + std::vector solution_values; + for (std::map::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 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 + void + apply_block_boundary_values (const std::map &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 > block_boundary_values(n_blocks); { - const std::pair 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 constrained_rows; - for (std::map::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 indices; - std::vector solution_values; + int block=0; + types::global_dof_index offset = 0; for (std::map::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; ifirst >= 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 (index,dof->second)); + } } - else - { - // clear_rows() is a collective operation so we still have to call - // it: - std::vector 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 - void - apply_block_boundary_values (const std::map &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 > block_boundary_values(n_blocks); - { - int block=0; - types::global_dof_index offset = 0; - for (std::map::const_iterator - dof = boundary_values.begin(); - dof != boundary_values.end(); - ++dof) + // Now call the non-block variants on + // the diagonal subblocks and the + // solution/rhs. + for (unsigned int block=0; blockfirst >= 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 (index,dof->second)); + const std::pair local_range + = matrix.block(block_m,0).local_range(); + + std::vector constrained_rows; + for (std::map::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 local_range - = matrix.block(block_m,0).local_range(); - - std::vector constrained_rows; - for (std::map::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 &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 &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 &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 &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 &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 &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 &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); -} + void + apply_boundary_values (const std::map &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 +} DEAL_II_NAMESPACE_CLOSE -- 2.39.5