]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new functions for multigrid boundaries
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 3 Jul 2006 07:57:23 +0000 (07:57 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 3 Jul 2006 07:57:23 +0000 (07:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@13324 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc [new file with mode: 0644]

diff --git a/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc b/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc
new file mode 100644 (file)
index 0000000..d11b7bd
--- /dev/null
@@ -0,0 +1,445 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+#include <base/logstream.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/block_vector.h>
+#include <multigrid/mg_tools.h>
+#include <multigrid/mg_base.h>
+#include <multigrid/mg_level_object.h>
+
+#include <vector>
+#include <algorithm>
+#include <numeric>
+
+
+template <typename number>
+void
+MGTools::apply_boundary_values (
+  const std::vector<unsigned int> &boundary_values,
+  SparseMatrix<number>& matrix,
+  const bool preserve_symmetry)
+{
+                                  // if no boundary values are to be applied
+                                  // simply return
+  if (boundary_values.size() == 0)
+    return;
+
+
+  const unsigned int n_dofs = matrix.m();
+
+                                  // if a diagonal entry is zero
+                                  // later, then we use another
+                                  // number instead. take it to be
+                                  // the first nonzero diagonal
+                                  // element of the matrix, or 1 if
+                                  // there is no such thing
+  number first_nonzero_diagonal_entry = 1;
+  for (unsigned int i=0; i<n_dofs; ++i)
+    if (matrix.diag_element(i) != 0)
+      {
+       first_nonzero_diagonal_entry = matrix.diag_element(i);
+       break;
+      }
+
+  
+  std::vector<unsigned int>::const_iterator dof  = boundary_values.begin(),
+                                           endd = boundary_values.end();
+  const SparsityPattern    &sparsity    = matrix.get_sparsity_pattern();
+  const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
+  const unsigned int *sparsity_colnums  = sparsity.get_column_numbers();
+  for (; dof != endd; ++dof)
+    {
+      Assert (*dof < n_dofs, ExcInternalError());
+      
+      const unsigned int dof_number = *dof;
+                                      // for each boundary dof:
+      
+                                      // set entries of this line
+                                      // to zero except for the diagonal
+                                      // entry. Note that the diagonal
+                                      // entry is always the first one
+                                      // for square matrices, i.e.
+                                      // we shall not set
+                                      // matrix.global_entry(
+                                      //     sparsity_rowstart[dof.first])
+      const unsigned int last = sparsity_rowstart[dof_number+1];
+      for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+       matrix.global_entry(j) = 0.;
+
+
+                                      // set right hand side to
+                                      // wanted value: if main diagonal
+                                      // entry nonzero, don't touch it
+                                      // and scale rhs accordingly. If
+                                      // zero, take the first main
+                                      // diagonal entry we can find, or
+                                      // one if no nonzero main diagonal
+                                      // element exists. Normally, however,
+                                      // the main diagonal entry should
+                                      // not be zero.
+                                      //
+                                      // store the new rhs entry to make
+                                      // the gauss step more efficient
+      matrix.set (dof_number, dof_number,
+                 first_nonzero_diagonal_entry);
+
+
+                                      // if the user wants to have
+                                      // the symmetry of the matrix
+                                      // preserved, and if the
+                                      // sparsity pattern is
+                                      // symmetric, then do a Gauss
+                                      // elimination step with the
+                                      // present row
+      if (preserve_symmetry)
+       {
+                                          // we have to loop over all
+                                          // rows of the matrix which
+                                          // have a nonzero entry in
+                                          // the column which we work
+                                          // in presently. if the
+                                          // sparsity pattern is
+                                          // symmetric, then we can
+                                          // get the positions of
+                                          // these rows cheaply by
+                                          // looking at the nonzero
+                                          // column numbers of the
+                                          // present row. we need not
+                                          // look at the first entry,
+                                          // since that is the
+                                          // diagonal element and
+                                          // thus the present row
+         for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+           {
+             const unsigned int row = sparsity_colnums[j];
+
+                                              // find the position of
+                                              // element
+                                              // (row,dof_number)
+             const unsigned int *
+               p = std::lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
+                                    &sparsity_colnums[sparsity_rowstart[row+1]],
+                                    dof_number);
+
+                                              // check whether this line has
+                                              // an entry in the regarding column
+                                              // (check for ==dof_number and
+                                              // != next_row, since if
+                                              // row==dof_number-1, *p is a
+                                              // past-the-end pointer but points
+                                              // to dof_number anyway...)
+                                              //
+                                              // there should be such an entry!
+             Assert ((*p == dof_number) &&
+                     (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
+                     ExcInternalError());
+
+             const unsigned int global_entry
+               = (p - &sparsity_colnums[sparsity_rowstart[0]]);
+             
+                                              // correct right hand side
+                                              // set matrix entry to zero
+             matrix.global_entry(global_entry) = 0.;
+           }
+       }
+    }
+}
+
+
+
+template <typename number>
+void
+MGTools::apply_boundary_values (
+  const std::vector<unsigned int>& boundary_values,
+  BlockSparseMatrix<number>& matrix,
+  const bool preserve_symmetry)
+{
+  const unsigned int blocks = matrix.n_block_rows();
+  
+  Assert (matrix.n_block_rows() == matrix.n_block_cols(),
+         ExcNotQuadratic());
+  Assert (matrix.get_sparsity_pattern().get_row_indices() == 
+         matrix.get_sparsity_pattern().get_column_indices(),
+         ExcNotQuadratic());
+
+  for (unsigned int i=0; i<blocks; ++i)
+    Assert (matrix.block(i,i).get_sparsity_pattern().optimize_diagonal(),
+           SparsityPattern::ExcDiagonalNotOptimized());
+
+  
+                                  // if no boundary values are to be applied
+                                  // simply return
+  if (boundary_values.size() == 0)
+    return;
+
+
+  const unsigned int n_dofs = matrix.m();
+
+                                  // if a diagonal entry is zero
+                                  // later, then we use another
+                                  // number instead. take it to be
+                                  // the first nonzero diagonal
+                                  // element of the matrix, or 1 if
+                                  // there is no such thing
+  number first_nonzero_diagonal_entry = 0;
+  for (unsigned int diag_block=0; diag_block<blocks; ++diag_block)
+    {
+      for (unsigned int i=0; i<matrix.block(diag_block,diag_block).n(); ++i)
+       if (matrix.block(diag_block,diag_block).diag_element(i) != 0)
+         {
+           first_nonzero_diagonal_entry 
+             = matrix.block(diag_block,diag_block).diag_element(i);
+           break;
+         }
+                                      // check whether we have found
+                                      // something in the present
+                                      // block
+      if (first_nonzero_diagonal_entry != 0)
+       break;
+    }
+                                  // nothing found on all diagonal
+                                  // blocks? if so, use 1.0 instead
+  if (first_nonzero_diagonal_entry == 0)
+    first_nonzero_diagonal_entry = 1;
+  
+
+  std::vector<unsigned int>::const_iterator dof  = boundary_values.begin(),
+                                           endd = boundary_values.end();
+  const BlockSparsityPattern &
+    sparsity_pattern = matrix.get_sparsity_pattern();
+  
+                                  // pointer to the mapping between
+                                  // global and block indices. since
+                                  // the row and column mappings are
+                                  // equal, store a pointer on only
+                                  // one of them
+  const BlockIndices &
+    index_mapping = sparsity_pattern.get_column_indices();
+  
+                                  // now loop over all boundary dofs
+  for (; dof != endd; ++dof)
+    {
+      Assert (*dof < n_dofs, ExcInternalError());
+
+                                      // get global index and index
+                                      // in the block in which this
+                                      // dof is located
+      const unsigned int dof_number = *dof;
+      const std::pair<unsigned int,unsigned int>
+       block_index = index_mapping.global_to_local (dof_number);
+
+                                      // for each boundary dof:
+      
+                                      // set entries of this line
+                                      // to zero except for the diagonal
+                                      // entry. Note that the diagonal
+                                      // entry is always the first one
+                                      // for square matrices, i.e.
+                                      // we shall not set
+                                      // matrix.global_entry(
+                                      //     sparsity_rowstart[dof.first])
+                                      // of the diagonal block
+      for (unsigned int block_col=0; block_col<blocks; ++block_col)
+       {
+         const SparsityPattern &
+           local_sparsity = sparsity_pattern.block(block_index.first,
+                                                   block_col);
+
+                                          // find first and last
+                                          // entry in the present row
+                                          // of the present
+                                          // block. exclude the main
+                                          // diagonal element, which
+                                          // is the diagonal element
+                                          // of a diagonal block,
+                                          // which must be a square
+                                          // matrix so the diagonal
+                                          // element is the first of
+                                          // this row.
+         const unsigned int 
+           last  = local_sparsity.get_rowstart_indices()[block_index.second+1],
+           first = (block_col == block_index.first ?
+                    local_sparsity.get_rowstart_indices()[block_index.second]+1 :
+                    local_sparsity.get_rowstart_indices()[block_index.second]);
+         
+         for (unsigned int j=first; j<last; ++j)
+           matrix.block(block_index.first,block_col).global_entry(j) = 0.;
+       }
+      
+      matrix.block(block_index.first, block_index.first)
+       .diag_element(block_index.second)
+       = first_nonzero_diagonal_entry;
+
+                                      // if the user wants to have
+                                      // the symmetry of the matrix
+                                      // preserved, and if the
+                                      // sparsity pattern is
+                                      // symmetric, then do a Gauss
+                                      // elimination step with the
+                                      // present row. this is a
+                                      // little more complicated for
+                                      // block matrices.
+      if (preserve_symmetry)
+       {
+                                          // we have to loop over all
+                                          // rows of the matrix which
+                                          // have a nonzero entry in
+                                          // the column which we work
+                                          // in presently. if the
+                                          // sparsity pattern is
+                                          // symmetric, then we can
+                                          // get the positions of
+                                          // these rows cheaply by
+                                          // looking at the nonzero
+                                          // column numbers of the
+                                          // present row.
+                                          //
+                                          // note that if we check
+                                          // whether row @p{row} in
+                                          // block (r,c) is non-zero,
+                                          // then we have to check
+                                          // for the existence of
+                                          // column @p{row} in block
+                                          // (c,r), i.e. of the
+                                          // transpose block
+         for (unsigned int block_row=0; block_row<blocks; ++block_row)
+           {
+                                              // get pointers to the
+                                              // sparsity patterns of
+                                              // this block and of
+                                              // the transpose one
+             const SparsityPattern &this_sparsity
+               = sparsity_pattern.block (block_row, block_index.first);
+             const SparsityPattern &transpose_sparsity
+               = sparsity_pattern.block (block_index.first, block_row);
+             
+                                              // traverse the row of
+                                              // the transpose block
+                                              // to find the
+                                              // interesting rows in
+                                              // the present block.
+                                              // don't use the
+                                              // diagonal element of
+                                              // the diagonal block
+             const unsigned int
+               first = (block_index.first == block_row ?
+                        transpose_sparsity.get_rowstart_indices()[block_index.second]+1 :
+                        transpose_sparsity.get_rowstart_indices()[block_index.second]),
+               last  = transpose_sparsity.get_rowstart_indices()[block_index.second+1];
+             
+             for (unsigned int j=first; j<last; ++j)
+               {
+                                                  // get the number
+                                                  // of the column in
+                                                  // this row in
+                                                  // which a nonzero
+                                                  // entry is. this
+                                                  // is also the row
+                                                  // of the transpose
+                                                  // block which has
+                                                  // an entry in the
+                                                  // interesting row
+                 const unsigned int row = transpose_sparsity.get_column_numbers()[j];
+
+                                                  // find the
+                                                  // position of
+                                                  // element
+                                                  // (row,dof_number)
+                                                  // in this block
+                                                  // (not in the
+                                                  // transpose
+                                                  // one). note that
+                                                  // we have to take
+                                                  // care of special
+                                                  // cases with
+                                                  // square
+                                                  // sub-matrices
+                 const unsigned int *p = 0;
+                 if (this_sparsity.n_rows() == this_sparsity.n_cols())
+                   {
+                     if (this_sparsity.get_column_numbers()
+                         [this_sparsity.get_rowstart_indices()[row]]
+                         ==
+                         block_index.second)
+                       p = &this_sparsity.get_column_numbers()
+                           [this_sparsity.get_rowstart_indices()[row]];
+                     else
+                       p = std::lower_bound(&this_sparsity.get_column_numbers()
+                                            [this_sparsity.get_rowstart_indices()[row]+1],
+                                            &this_sparsity.get_column_numbers()
+                                            [this_sparsity.get_rowstart_indices()[row+1]],
+                                            block_index.second);
+                   }
+                 else
+                   p = std::lower_bound(&this_sparsity.get_column_numbers()
+                                        [this_sparsity.get_rowstart_indices()[row]],
+                                        &this_sparsity.get_column_numbers()
+                                        [this_sparsity.get_rowstart_indices()[row+1]],
+                                        block_index.second);
+
+                                                  // check whether this line has
+                                                  // an entry in the regarding column
+                                                  // (check for ==dof_number and
+                                                  // != next_row, since if
+                                                  // row==dof_number-1, *p is a
+                                                  // past-the-end pointer but points
+                                                  // to dof_number anyway...)
+                                                  //
+                                                  // there should be
+                                                  // such an entry!
+                                                  // note, however,
+                                                  // that this
+                                                  // assertion will
+                                                  // fail sometimes
+                                                  // if the sparsity
+                                                  // pattern is not
+                                                  // symmetric!
+                 Assert ((*p == block_index.second) &&
+                         (p != &this_sparsity.get_column_numbers()
+                          [this_sparsity.get_rowstart_indices()[row+1]]),
+                         ExcInternalError());
+                 
+                 const unsigned int global_entry
+                   = (p
+                      -
+                      &this_sparsity.get_column_numbers()
+                      [this_sparsity.get_rowstart_indices()[0]]);
+
+                                                  // set matrix entry to zero
+                 matrix.block(block_row,block_index.first).global_entry(global_entry) = 0.;
+               }
+           }
+       }
+    }
+}
+
+
+template void MGTools::apply_boundary_values (
+  const std::vector<unsigned int>& boundary_values,
+  SparseMatrix<float>& matrix,
+  const bool preserve_symmetry);
+template void MGTools::apply_boundary_values (
+  const std::vector<unsigned int>& boundary_values,
+  SparseMatrix<double>& matrix,
+  const bool preserve_symmetry);
+template void MGTools::apply_boundary_values (
+  const std::vector<unsigned int>& boundary_values,
+  BlockSparseMatrix<float>& matrix,
+  const bool preserve_symmetry);
+template void MGTools::apply_boundary_values (
+  const std::vector<unsigned int>& boundary_values,
+  BlockSparseMatrix<double>& matrix,
+  const bool preserve_symmetry);

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.