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