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