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