]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Patch by Michal Wichrowski: fix the same bug in a second place.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 15 Jun 2014 15:46:13 +0000 (15:46 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 15 Jun 2014 15:46:13 +0000 (15:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@33047 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/numerics/matrix_tools.cc

index dbb7d854f6a586148b72fa7625b7b7969f62a29f..b836e73cd8e6f4e3a96419bc3adc268f711df1d2 100644 (file)
@@ -2553,77 +2553,85 @@ namespace MatrixTools
         Assert (matrix.n() == solution.size(),
                 ExcDimensionMismatch(matrix.m(), solution.size()));
 
-        // if no boundary values are to be applied
-        // simply return
-        if (boundary_values.size() == 0)
-          return;
-
-        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());
-
-        // we have to read and write from this
-        // matrix (in this order). this will only
-        // work if we compress the matrix first,
-        // done here
-        matrix.compress ();
-
-        // 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,double>::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,double>::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);
-
+        // 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());
+
+           // we have to read and write from this
+           // matrix (in this order). this will only
+           // work if we compress the matrix first,
+           // done here
+           matrix.compress ();
+
+           // 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,double>::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,double>::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 ();
         solution.compress (VectorOperation::insert);

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.