]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added apply_boundary_conditions() function to Trilinos Block matrices and vectors.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 2 Sep 2008 14:27:26 +0000 (14:27 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 2 Sep 2008 14:27:26 +0000 (14:27 +0000)
git-svn-id: https://svn.dealii.org/trunk@16707 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/source/numerics/matrices.all_dimensions.cc
deal.II/lac/include/lac/trilinos_block_sparse_matrix.h
deal.II/lac/include/lac/trilinos_block_vector.h
deal.II/lac/source/trilinos_block_sparse_matrix.cc
deal.II/lac/source/trilinos_sparse_matrix.cc

index 476afa305d2b2f388fd5ecd99d8db76357d17ebe..9b8db14cbab205dd7884ae614b92d2c900046fec 100644 (file)
@@ -65,6 +65,8 @@ namespace TrilinosWrappers
 {
   class SparseMatrix;
   class Vector;
+  class BlockSparseMatrix;
+  class BlockVector;
 }
 #endif
 
@@ -1137,10 +1139,6 @@ class MatrixTools : public MatrixCreator
                                      * is that it doesn't handle the case
                                      * where the matrix is distributed across
                                      * an MPI system.
-                                     *
-                                     * This function is used in
-                                     * @ref step_17 "step-17" and
-                                     * @ref step_18 "step-18".
                                      */
     static void
     apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
@@ -1148,6 +1146,18 @@ class MatrixTools : public MatrixCreator
                           TrilinosWrappers::Vector  &solution,
                           TrilinosWrappers::Vector  &right_hand_side,
                           const bool             eliminate_columns = true);
+    
+                                    /**
+                                     * This function does the same as
+                                     * the one above, except now working
+                                     * on block structures.
+                                     */
+    static void
+    apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+                          TrilinosWrappers::BlockSparseMatrix  &matrix,
+                          TrilinosWrappers::BlockVector  &solution,
+                          TrilinosWrappers::BlockVector  &right_hand_side,
+                          const bool                eliminate_columns = true);
 #endif
     
                                      /**
index 3034f8123651288e26e30297cdbde01880a1ae3a..9fcf748c6ce69c493cf4ffebe33c834044d5acbc 100644 (file)
@@ -29,6 +29,8 @@
 #ifdef DEAL_II_USE_TRILINOS
 #  include <lac/trilinos_sparse_matrix.h>
 #  include <lac/trilinos_vector.h>
+#  include <lac/trilinos_block_sparse_matrix.h>
+#  include <lac/trilinos_block_vector.h>
 #endif
 
 #include <algorithm>
@@ -700,7 +702,7 @@ namespace TrilinosWrappers
     Assert (matrix.n() == right_hand_side.size(),
             ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
     Assert (matrix.n() == solution.size(),
-            ExcDimensionMismatch(matrix.n(), solution.size()));
+            ExcDimensionMismatch(matrix.m(), solution.size()));
 
                                      // if no boundary values are to be applied
                                      // simply return
@@ -713,7 +715,6 @@ namespace TrilinosWrappers
             ExcInternalError());
     Assert (local_range == solution.local_range(),
             ExcInternalError());
-    
 
                                      // we have to read and write from this
                                      // matrix (in this order). this will only
@@ -786,7 +787,7 @@ namespace TrilinosWrappers
                                      // 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;
+      solution_values[i] *= matrix.diag_element(i);
 
     right_hand_side.set (indices, solution_values);
 
@@ -813,6 +814,81 @@ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
                                         right_hand_side, eliminate_columns);
 }
 
+void
+MatrixTools::
+apply_boundary_values (const std::map<unsigned int,double>  &boundary_values,
+                       TrilinosWrappers::BlockSparseMatrix  &matrix,
+                       TrilinosWrappers::BlockVector        &solution,
+                       TrilinosWrappers::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();
+
+  matrix.compress();
+
+                                  // 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<unsigned int,double> > block_boundary_values(n_blocks);
+  {
+    int offset = 0, block=0;
+    for (std::map<unsigned int,double>::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 unsigned int index = dof->first - offset;
+       block_boundary_values[block].insert(std::pair<unsigned int, double> (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<unsigned int, unsigned int> local_range
+       = matrix.block(block_m,0).local_range();
+      
+      std::vector<unsigned int> constrained_rows;
+      for (std::map<unsigned int,double>::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);
+    }
+}
+
 #endif
 
 
index cde0a7a12274de4f0e633c8cd5ffcf427cca2c78..586a478be5f0a193027897c7228493713ad4b992 100644 (file)
@@ -228,6 +228,13 @@ namespace TrilinosWrappers
                                         */
       void collect_sizes ();
 
+                                       /**
+                                        * Return the number of nonzero
+                                        * elements of this
+                                        * matrix. 
+                                        */
+      unsigned int n_nonzero_elements () const;
+
                                        /**
                                         * Matrix-vector multiplication:
                                         * let $dst = M*src$ with $M$
index 61547947b0a82e5f075ef87a40dc8eddc63d89c7..d7b3c2d0e2a3e86a1181d0d0032efd0fba6359f2 100644 (file)
@@ -171,7 +171,8 @@ namespace TrilinosWrappers
                                           * since they may be routed to
                                           * the wrong block.
                                           */
-        void reinit (const BlockVector &V);
+        void reinit (const BlockVector &V,
+                    const bool fast=false);
 
                                          /**
                                           * Change the number of blocks to
@@ -328,14 +329,15 @@ namespace TrilinosWrappers
 
     inline
     void
-    BlockVector::reinit (const BlockVector& v)
+    BlockVector::reinit (const BlockVector& v,
+                        const bool fast)
     {
       this->block_indices = v.get_block_indices();
       if (this->components.size() != this->n_blocks())
         this->components.resize(this->n_blocks());
   
       for (unsigned int i=0;i<this->n_blocks();++i)
-        block(i).reinit(v.block(i));
+        block(i).reinit(v.block(i), fast);
     }
 
 
index 04f43fe0f5c7709b7042a04be9b8686c239fc9a9..90731a9cc46956f782bd39071ce8e98af40cb1a6 100644 (file)
@@ -156,16 +156,31 @@ namespace TrilinosWrappers
        this->block(r,c).compress();
   }
 
+
+
   void
   BlockSparseMatrix::collect_sizes ()
   {
     compress();
     BaseClass::collect_sizes ();
   }
-  
+
+
+
+  unsigned int
+  BlockSparseMatrix::n_nonzero_elements () const
+  {
+    unsigned int n_nonzero = 0;
+    for (unsigned int rows = 0; rows<this->n_block_rows(); ++rows)
+      for (unsigned int cols = 0; cols<this->n_block_cols(); ++cols)
+       n_nonzero += this->block(rows,cols).n_nonzero_elements();
+
+    return n_nonzero;
+  }
+
 }
 
-  
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index f83268fe4e9ee83d556ef883720eaffb9441d07c..b48139bc9f3b2cb0b5f2289fda462dc438d48905 100755 (executable)
@@ -492,10 +492,10 @@ namespace TrilinosWrappers
        int diag_index = (int)(diag_find - col_indices);
 
        for (int j=0; j<num_entries; ++j)
-         if (diag_index != j)
+         if (diag_index != j || new_diag_value == 0)
            values[j] = 0.;
 
-       if (diag_find && std::fabs(values[diag_index]) == 0.)
+       if (diag_find && std::fabs(values[diag_index]) == 0 && new_diag_value!=0)
          values[diag_index] = new_diag_value;
       }
   }
@@ -730,9 +730,9 @@ namespace TrilinosWrappers
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
     
-    Assert (col_map.SameAs(dst.map),
+    Assert (col_map.SameAs(src.map),
            ExcMessage ("Column map of matrix does not fit with vector map!"));
-    Assert (row_map.SameAs(src.map),
+    Assert (row_map.SameAs(dst.map),
            ExcMessage ("Row map of matrix does not fit with vector map!"));
 
     if (!matrix->Filled())
@@ -750,9 +750,9 @@ namespace TrilinosWrappers
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
 
-    Assert (row_map.SameAs(dst.map),
+    Assert (row_map.SameAs(src.map),
            ExcMessage ("Row map of matrix does not fit with vector map!"));
-    Assert (col_map.SameAs(src.map),
+    Assert (col_map.SameAs(dst.map),
            ExcMessage ("Column map of matrix does not fit with vector map!"));
 
     if (!matrix->Filled())

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.