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