]> https://gitweb.dealii.org/ - dealii.git/commitdiff
fix indentation 2273/head
authorTimo Heister <timo.heister@gmail.com>
Mon, 29 Feb 2016 16:23:42 +0000 (11:23 -0500)
committerTimo Heister <timo.heister@gmail.com>
Mon, 29 Feb 2016 16:23:42 +0000 (11:23 -0500)
include/deal.II/numerics/matrix_creator.templates.h
source/numerics/matrix_tools_once.cc

index 4541403de0889966d9989af777e32071283584b2..72c90ddc3640e7675165e76b20b9cc8fee2b0508 100644 (file)
@@ -1196,9 +1196,9 @@ namespace MatrixCreator
 
   }
 
-  
-  
-  
+
+
+
 
 
 
@@ -1260,7 +1260,7 @@ namespace MatrixCreator
                                      std_cxx11::cref(component_mapping))),
                     static_cast<std_cxx11::function<void (MatrixCreator::internal::AssemblerBoundary
                                                           ::CopyData<DoFHandler<dim,spacedim>,number> const &)> > (std_cxx11::bind(
-                                                                                                                    &internal::copy_boundary_mass_matrix_1<dim,spacedim,number>,
+                                                                &internal::copy_boundary_mass_matrix_1<dim,spacedim,number>,
                                                                 std_cxx11::_1,
                                                                 std_cxx11::cref(boundary_functions),
                                                                 std_cxx11::cref(dof_to_boundary_mapping),
index 4ccee230c86c751b6d83601a2e0c380e066ba362..cf2d90a50305fdd620c32782374da7537cea88f4 100644 (file)
@@ -61,213 +61,213 @@ DEAL_II_NAMESPACE_OPEN
 
 #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
 
@@ -275,244 +275,244 @@ DEAL_II_NAMESPACE_OPEN
 
 #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
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.