Reference documentation for deal.II version 8.4.1
matrix_tools.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2015 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/function.h>
17 #include <deal.II/base/quadrature.h>
18 #include <deal.II/base/work_stream.h>
19 #include <deal.II/base/geometry_info.h>
20 #include <deal.II/base/quadrature.h>
21 #include <deal.II/dofs/dof_handler.h>
22 #include <deal.II/dofs/dof_accessor.h>
23 #include <deal.II/dofs/dof_tools.h>
24 #include <deal.II/fe/fe.h>
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping_q1.h>
27 #include <deal.II/grid/tria_iterator.h>
28 #include <deal.II/hp/fe_values.h>
29 #include <deal.II/hp/mapping_collection.h>
30 #include <deal.II/numerics/matrix_tools.h>
31 #include <deal.II/lac/vector.h>
32 #include <deal.II/lac/block_vector.h>
33 #include <deal.II/lac/full_matrix.h>
34 #include <deal.II/lac/sparse_matrix.h>
35 #include <deal.II/lac/block_sparse_matrix.h>
36 
37 #ifdef DEAL_II_WITH_PETSC
38 # include <deal.II/lac/petsc_parallel_sparse_matrix.h>
39 # include <deal.II/lac/petsc_sparse_matrix.h>
40 # include <deal.II/lac/petsc_parallel_vector.h>
41 # include <deal.II/lac/petsc_vector.h>
42 # include <deal.II/lac/petsc_parallel_block_sparse_matrix.h>
43 #endif
44 
45 #ifdef DEAL_II_WITH_TRILINOS
46 # include <deal.II/lac/trilinos_sparse_matrix.h>
47 # include <deal.II/lac/trilinos_vector.h>
48 # include <deal.II/lac/trilinos_block_sparse_matrix.h>
49 # include <deal.II/lac/trilinos_block_vector.h>
50 #endif
51 
52 #include <algorithm>
53 
54 
55 #include <algorithm>
56 #include <set>
57 #include <cmath>
58 
59 
60 DEAL_II_NAMESPACE_OPEN
61 
62 namespace MatrixTools
63 {
64  namespace
65  {
66  template <typename Iterator>
67  bool column_less_than(const typename Iterator::value_type p,
68  const unsigned int column)
69  {
70  return (p.column() < column);
71  }
72  }
73 
74 //TODO:[WB] I don't think that the optimized storage of diagonals is needed (GK)
75  template <typename number>
76  void
77  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
78  SparseMatrix<number> &matrix,
79  Vector<number> &solution,
80  Vector<number> &right_hand_side,
81  const bool eliminate_columns)
82  {
83  Assert (matrix.n() == right_hand_side.size(),
84  ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
85  Assert (matrix.n() == solution.size(),
86  ExcDimensionMismatch(matrix.n(), solution.size()));
87  Assert (matrix.n() == matrix.m(),
88  ExcDimensionMismatch(matrix.n(), matrix.m()));
89 
90  // if no boundary values are to be applied
91  // simply return
92  if (boundary_values.size() == 0)
93  return;
94 
95 
96  const types::global_dof_index n_dofs = matrix.m();
97 
98  // if a diagonal entry is zero
99  // later, then we use another
100  // number instead. take it to be
101  // the first nonzero diagonal
102  // element of the matrix, or 1 if
103  // there is no such thing
104  number first_nonzero_diagonal_entry = 1;
105  for (unsigned int i=0; i<n_dofs; ++i)
106  if (matrix.diag_element(i) != 0)
107  {
108  first_nonzero_diagonal_entry = matrix.diag_element(i);
109  break;
110  }
111 
112 
113  std::map<types::global_dof_index,double>::const_iterator dof = boundary_values.begin(),
114  endd = boundary_values.end();
115  for (; dof != endd; ++dof)
116  {
117  Assert (dof->first < n_dofs, ExcInternalError());
118 
119  const types::global_dof_index dof_number = dof->first;
120  // for each boundary dof:
121 
122  // set entries of this line to zero except for the diagonal
123  // entry
124  for (typename SparseMatrix<number>::iterator
125  p = matrix.begin(dof_number);
126  p != matrix.end(dof_number); ++p)
127  if (p->column() != dof_number)
128  p->value() = 0.;
129 
130  // set right hand side to
131  // wanted value: if main diagonal
132  // entry nonzero, don't touch it
133  // and scale rhs accordingly. If
134  // zero, take the first main
135  // diagonal entry we can find, or
136  // one if no nonzero main diagonal
137  // element exists. Normally, however,
138  // the main diagonal entry should
139  // not be zero.
140  //
141  // store the new rhs entry to make
142  // the gauss step more efficient
143  number new_rhs;
144  if (matrix.diag_element(dof_number) != 0.0)
145  {
146  new_rhs = dof->second * matrix.diag_element(dof_number);
147  right_hand_side(dof_number) = new_rhs;
148  }
149  else
150  {
151  matrix.set (dof_number, dof_number,
152  first_nonzero_diagonal_entry);
153  new_rhs = dof->second * first_nonzero_diagonal_entry;
154  right_hand_side(dof_number) = new_rhs;
155  }
156 
157 
158  // if the user wants to have
159  // the symmetry of the matrix
160  // preserved, and if the
161  // sparsity pattern is
162  // symmetric, then do a Gauss
163  // elimination step with the
164  // present row
165  if (eliminate_columns)
166  {
167  // store the only nonzero entry
168  // of this line for the Gauss
169  // elimination step
170  const number diagonal_entry = matrix.diag_element(dof_number);
171 
172  // we have to loop over all rows of the matrix which have
173  // a nonzero entry in the column which we work in
174  // presently. if the sparsity pattern is symmetric, then
175  // we can get the positions of these rows cheaply by
176  // looking at the nonzero column numbers of the present
177  // row. we need not look at the first entry of each row,
178  // since that is the diagonal element and thus the present
179  // row
180  for (typename SparseMatrix<number>::iterator
181  q = matrix.begin(dof_number)+1;
182  q != matrix.end(dof_number); ++q)
183  {
184  const types::global_dof_index row = q->column();
185 
186  // find the position of
187  // element
188  // (row,dof_number)
189  bool (*comp)(const typename SparseMatrix<number>::iterator::value_type p,
190  const unsigned int column)
191  = &column_less_than<typename SparseMatrix<number>::iterator>;
192  const typename SparseMatrix<number>::iterator
193  p = Utilities::lower_bound(matrix.begin(row)+1,
194  matrix.end(row),
195  dof_number,
196  comp);
197 
198  // check whether this line has an entry in the
199  // regarding column (check for ==dof_number and !=
200  // next_row, since if row==dof_number-1, *p is a
201  // past-the-end pointer but points to dof_number
202  // anyway...)
203  //
204  // there should be such an entry! we know this because
205  // we have assumed that the sparsity pattern is
206  // symmetric and we only walk over those rows for
207  // which the current row has a column entry
208  Assert ((p != matrix.end(row))
209  &&
210  (p->column() == dof_number),
211  ExcMessage("This function is trying to access an element of the "
212  "matrix that doesn't seem to exist. Are you using a "
213  "nonsymmetric sparsity pattern? If so, you are not "
214  "allowed to set the eliminate_column argument of this "
215  "function, see the documentation."));
216 
217  // correct right hand side
218  right_hand_side(row) -= p->value() /
219  diagonal_entry * new_rhs;
220 
221  // set matrix entry to zero
222  p->value() = 0.;
223  }
224  }
225 
226  // preset solution vector
227  solution(dof_number) = dof->second;
228  }
229  }
230 
231 
232 
233  template <typename number>
234  void
235  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
237  BlockVector<number> &solution,
238  BlockVector<number> &right_hand_side,
239  const bool eliminate_columns)
240  {
241  const unsigned int blocks = matrix.n_block_rows();
242 
243  Assert (matrix.n() == right_hand_side.size(),
244  ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
245  Assert (matrix.n() == solution.size(),
246  ExcDimensionMismatch(matrix.n(), solution.size()));
247  Assert (matrix.n_block_rows() == matrix.n_block_cols(),
248  ExcNotQuadratic());
251  ExcNotQuadratic());
253  solution.get_block_indices (),
254  ExcBlocksDontMatch ());
256  right_hand_side.get_block_indices (),
257  ExcBlocksDontMatch ());
258 
259  // if no boundary values are to be applied
260  // simply return
261  if (boundary_values.size() == 0)
262  return;
263 
264 
265  const types::global_dof_index n_dofs = matrix.m();
266 
267  // if a diagonal entry is zero
268  // later, then we use another
269  // number instead. take it to be
270  // the first nonzero diagonal
271  // element of the matrix, or 1 if
272  // there is no such thing
273  number first_nonzero_diagonal_entry = 0;
274  for (unsigned int diag_block=0; diag_block<blocks; ++diag_block)
275  {
276  for (unsigned int i=0; i<matrix.block(diag_block,diag_block).n(); ++i)
277  if (matrix.block(diag_block,diag_block).diag_element(i) != 0)
278  {
279  first_nonzero_diagonal_entry
280  = matrix.block(diag_block,diag_block).diag_element(i);
281  break;
282  }
283  // check whether we have found
284  // something in the present
285  // block
286  if (first_nonzero_diagonal_entry != 0)
287  break;
288  }
289  // nothing found on all diagonal
290  // blocks? if so, use 1.0 instead
291  if (first_nonzero_diagonal_entry == 0)
292  first_nonzero_diagonal_entry = 1;
293 
294 
295  std::map<types::global_dof_index,double>::const_iterator dof = boundary_values.begin(),
296  endd = boundary_values.end();
297  const BlockSparsityPattern &
298  sparsity_pattern = matrix.get_sparsity_pattern();
299 
300  // pointer to the mapping between
301  // global and block indices. since
302  // the row and column mappings are
303  // equal, store a pointer on only
304  // one of them
305  const BlockIndices &
306  index_mapping = sparsity_pattern.get_column_indices();
307 
308  // now loop over all boundary dofs
309  for (; dof != endd; ++dof)
310  {
311  Assert (dof->first < n_dofs, ExcInternalError());
312  (void)n_dofs;
313 
314  // get global index and index
315  // in the block in which this
316  // dof is located
317  const types::global_dof_index dof_number = dof->first;
318  const std::pair<unsigned int,types::global_dof_index>
319  block_index = index_mapping.global_to_local (dof_number);
320 
321  // for each boundary dof:
322 
323  // set entries of this line
324  // to zero except for the diagonal
325  // entry. Note that the diagonal
326  // entry is always the first one
327  // in a row for square matrices
328  for (unsigned int block_col=0; block_col<blocks; ++block_col)
329  for (typename SparseMatrix<number>::iterator
330  p = (block_col == block_index.first ?
331  matrix.block(block_index.first,block_col).begin(block_index.second) + 1 :
332  matrix.block(block_index.first,block_col).begin(block_index.second));
333  p != matrix.block(block_index.first,block_col).end(block_index.second);
334  ++p)
335  p->value() = 0;
336 
337  // set right hand side to
338  // wanted value: if main diagonal
339  // entry nonzero, don't touch it
340  // and scale rhs accordingly. If
341  // zero, take the first main
342  // diagonal entry we can find, or
343  // one if no nonzero main diagonal
344  // element exists. Normally, however,
345  // the main diagonal entry should
346  // not be zero.
347  //
348  // store the new rhs entry to make
349  // the gauss step more efficient
350  number new_rhs;
351  if (matrix.block(block_index.first, block_index.first)
352  .diag_element(block_index.second) != 0.0)
353  new_rhs = dof->second *
354  matrix.block(block_index.first, block_index.first)
355  .diag_element(block_index.second);
356  else
357  {
358  matrix.block(block_index.first, block_index.first)
359  .diag_element(block_index.second)
360  = first_nonzero_diagonal_entry;
361  new_rhs = dof->second * first_nonzero_diagonal_entry;
362  }
363  right_hand_side.block(block_index.first)(block_index.second)
364  = new_rhs;
365 
366 
367  // if the user wants to have
368  // the symmetry of the matrix
369  // preserved, and if the
370  // sparsity pattern is
371  // symmetric, then do a Gauss
372  // elimination step with the
373  // present row. this is a
374  // little more complicated for
375  // block matrices.
376  if (eliminate_columns)
377  {
378  // store the only nonzero entry
379  // of this line for the Gauss
380  // elimination step
381  const number diagonal_entry
382  = matrix.block(block_index.first,block_index.first)
383  .diag_element(block_index.second);
384 
385  // we have to loop over all
386  // rows of the matrix which
387  // have a nonzero entry in
388  // the column which we work
389  // in presently. if the
390  // sparsity pattern is
391  // symmetric, then we can
392  // get the positions of
393  // these rows cheaply by
394  // looking at the nonzero
395  // column numbers of the
396  // present row.
397  //
398  // note that if we check
399  // whether row @p{row} in
400  // block (r,c) is non-zero,
401  // then we have to check
402  // for the existence of
403  // column @p{row} in block
404  // (c,r), i.e. of the
405  // transpose block
406  for (unsigned int block_row=0; block_row<blocks; ++block_row)
407  {
408  // get pointers to the sparsity patterns of this block and of
409  // the transpose one
410  const SparsityPattern &this_sparsity
411  = sparsity_pattern.block (block_row, block_index.first);
412 
413  SparseMatrix<number> &this_matrix
414  = matrix.block(block_row, block_index.first);
415  SparseMatrix<number> &transpose_matrix
416  = matrix.block(block_index.first, block_row);
417 
418  // traverse the row of the transpose block to find the
419  // interesting rows in the present block. don't use the
420  // diagonal element of the diagonal block
421  for (typename SparseMatrix<number>::iterator
422  q = (block_index.first == block_row ?
423  transpose_matrix.begin(block_index.second)+1 :
424  transpose_matrix.begin(block_index.second));
425  q != transpose_matrix.end(block_index.second);
426  ++q)
427  {
428  // get the number of the column in this row in which a
429  // nonzero entry is. this is also the row of the transpose
430  // block which has an entry in the interesting row
431  const types::global_dof_index row = q->column();
432 
433  // find the position of element (row,dof_number) in this
434  // block (not in the transpose one). note that we have to
435  // take care of special cases with square sub-matrices
436  bool (*comp)(typename SparseMatrix<number>::iterator::value_type p,
437  const unsigned int column)
438  = &column_less_than<typename SparseMatrix<number>::iterator>;
439 
440  typename SparseMatrix<number>::iterator p = this_matrix.end();
441 
442  if (this_sparsity.n_rows() == this_sparsity.n_cols())
443  {
444  if (this_matrix.begin(row)->column()
445  ==
446  block_index.second)
447  p = this_matrix.begin(row);
448  else
449  p = Utilities::lower_bound(this_matrix.begin(row)+1,
450  this_matrix.end(row),
451  block_index.second,
452  comp);
453  }
454  else
455  p = Utilities::lower_bound(this_matrix.begin(row),
456  this_matrix.end(row),
457  block_index.second,
458  comp);
459 
460  // check whether this line has an entry in the
461  // regarding column (check for ==dof_number and !=
462  // next_row, since if row==dof_number-1, *p is a
463  // past-the-end pointer but points to dof_number
464  // anyway...)
465  //
466  // there should be such an entry! we know this because
467  // we have assumed that the sparsity pattern is
468  // symmetric and we only walk over those rows for
469  // which the current row has a column entry
470  Assert ((p->column() == block_index.second) &&
471  (p != this_matrix.end(row)),
472  ExcInternalError());
473 
474  // correct right hand side
475  right_hand_side.block(block_row)(row)
476  -= p->value() /
477  diagonal_entry * new_rhs;
478 
479  // set matrix entry to zero
480  p->value() = 0.;
481  }
482  }
483  }
484 
485  // preset solution vector
486  solution.block(block_index.first)(block_index.second) = dof->second;
487  }
488  }
489 
490 
491 
492 
493 
494 
495  void
496  local_apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
497  const std::vector<types::global_dof_index> &local_dof_indices,
498  FullMatrix<double> &local_matrix,
499  Vector<double> &local_rhs,
500  const bool eliminate_columns)
501  {
502  Assert (local_dof_indices.size() == local_matrix.m(),
503  ExcDimensionMismatch(local_dof_indices.size(),
504  local_matrix.m()));
505  Assert (local_dof_indices.size() == local_matrix.n(),
506  ExcDimensionMismatch(local_dof_indices.size(),
507  local_matrix.n()));
508  Assert (local_dof_indices.size() == local_rhs.size(),
509  ExcDimensionMismatch(local_dof_indices.size(),
510  local_rhs.size()));
511 
512  // if there is nothing to do, then exit
513  // right away
514  if (boundary_values.size() == 0)
515  return;
516 
517  // otherwise traverse all the dofs used in
518  // the local matrices and vectors and see
519  // what's there to do
520 
521  // if we need to treat an entry, then we
522  // set the diagonal entry to its absolute
523  // value. if it is zero, we used to set it
524  // to one, which is a really terrible
525  // choice that can lead to hours of
526  // searching for bugs in programs (I
527  // experienced this :-( ) if the matrix
528  // entries are otherwise very large. this
529  // is so since iterative solvers would
530  // simply not correct boundary nodes for
531  // their correct values since the residual
532  // contributions of their rows of the
533  // linear system is almost zero if the
534  // diagonal entry is one. thus, set it to
535  // the average absolute value of the
536  // nonzero diagonal elements.
537  //
538  // we only compute this value lazily the
539  // first time we need it.
540  double average_diagonal = 0;
541  const unsigned int n_local_dofs = local_dof_indices.size();
542  for (unsigned int i=0; i<n_local_dofs; ++i)
543  {
544  const std::map<types::global_dof_index, double>::const_iterator
545  boundary_value = boundary_values.find (local_dof_indices[i]);
546  if (boundary_value != boundary_values.end())
547  {
548  // remove this row, except for the
549  // diagonal element
550  for (unsigned int j=0; j<n_local_dofs; ++j)
551  if (i != j)
552  local_matrix(i,j) = 0;
553 
554  // replace diagonal entry by its
555  // absolute value to make sure that
556  // everything remains positive, or
557  // by the average diagonal value if
558  // zero
559  if (local_matrix(i,i) == 0.)
560  {
561  // if average diagonal hasn't
562  // yet been computed, do so now
563  if (average_diagonal == 0.)
564  {
565  unsigned int nonzero_diagonals = 0;
566  for (unsigned int k=0; k<n_local_dofs; ++k)
567  if (local_matrix(k,k) != 0.)
568  {
569  average_diagonal += std::fabs(local_matrix(k,k));
570  ++nonzero_diagonals;
571  }
572  if (nonzero_diagonals != 0)
573  average_diagonal /= nonzero_diagonals;
574  else
575  average_diagonal = 0;
576  }
577 
578  // only if all diagonal entries
579  // are zero, then resort to the
580  // last measure: choose one
581  if (average_diagonal == 0.)
582  average_diagonal = 1.;
583 
584  local_matrix(i,i) = average_diagonal;
585  }
586  else
587  local_matrix(i,i) = std::fabs(local_matrix(i,i));
588 
589  // and replace rhs entry by correct
590  // value
591  local_rhs(i) = local_matrix(i,i) * boundary_value->second;
592 
593  // finally do the elimination step
594  // if requested
595  if (eliminate_columns == true)
596  {
597  for (unsigned int row=0; row<n_local_dofs; ++row)
598  if (row != i)
599  {
600  local_rhs(row) -= local_matrix(row,i) * boundary_value->second;
601  local_matrix(row,i) = 0;
602  }
603  }
604  }
605  }
606  }
607 }
608 
609 
610 
611 // explicit instantiations
612 #include "matrix_tools.inst"
613 
614 
615 DEAL_II_NAMESPACE_CLOSE
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:604
const BlockIndices & get_row_indices() const
std::size_t size() const
size_type m() const
::ExceptionBase & ExcMessage(std::string arg1)
void set(const size_type i, const size_type j, const number value)
size_type n() const
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:54
SparsityPatternType & block(const size_type row, const size_type column)
const_iterator begin() const
void local_apply_boundary_values(const std::map< types::global_dof_index, double > &boundary_values, const std::vector< types::global_dof_index > &local_dof_indices, FullMatrix< double > &local_matrix, Vector< double > &local_rhs, const bool eliminate_columns)
size_type n() const
number diag_element(const size_type i) const
unsigned int global_dof_index
Definition: types.h:88
#define Assert(cond, exc)
Definition: exceptions.h:294
const BlockIndices & get_block_indices() const
size_type n_cols() const
std::size_t size() const
size_type m() const
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
const Accessor< number, Constness > & value_type
const_iterator end() const
const BlockSparsityPattern & get_sparsity_pattern() const
const BlockIndices & get_column_indices() const
BlockType & block(const unsigned int row, const unsigned int column)
size_type n_rows() const
BlockType & block(const unsigned int i)
void apply_boundary_values(const std::map< types::global_dof_index, double > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Definition: matrix_tools.cc:77