Reference documentation for deal.II version 9.5.0
\(\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_once.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2023 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
20
24
25#include <deal.II/fe/fe.h>
27
29
34#include <deal.II/lac/vector.h>
35
37
38#ifdef DEAL_II_WITH_PETSC
43#endif
44
45#ifdef DEAL_II_WITH_TRILINOS
50#endif
51
52#include <algorithm>
53#include <cmath>
54
55
57
58namespace MatrixTools
59{
60#ifdef DEAL_II_WITH_PETSC
61 void
63 const std::map<types::global_dof_index, PetscScalar> &boundary_values,
66 PETScWrappers::VectorBase & right_hand_side,
67 const bool eliminate_columns)
68 {
69 Assert(matrix.n() == right_hand_side.size(),
70 ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
71 Assert(matrix.n() == solution.size(),
72 ExcDimensionMismatch(matrix.n(), solution.size()));
73
74 // if no boundary values are to be applied, then
75 // jump straight to the compress() calls that we still have
76 // to perform because they are collective operations
77 if (boundary_values.size() > 0)
78 {
79 const std::pair<types::global_dof_index, types::global_dof_index>
80 local_range = matrix.local_range();
81 Assert(local_range == right_hand_side.local_range(),
83 Assert(local_range == solution.local_range(), ExcInternalError());
84
85 // determine the first nonzero diagonal
86 // entry from within the part of the
87 // matrix that we can see. if we can't
88 // find such an entry, take one
89 PetscScalar average_nonzero_diagonal_entry = 1;
90 for (types::global_dof_index i = local_range.first;
91 i < local_range.second;
92 ++i)
93 if (matrix.diag_element(i) != PetscScalar())
94 {
95 average_nonzero_diagonal_entry = std::abs(matrix.diag_element(i));
96 break;
97 }
98
99 // figure out which rows of the matrix we
100 // have to eliminate on this processor
101 std::vector<types::global_dof_index> constrained_rows;
102 for (const auto &boundary_value : boundary_values)
103 if ((boundary_value.first >= local_range.first) &&
104 (boundary_value.first < local_range.second))
105 constrained_rows.push_back(boundary_value.first);
106
107 // then eliminate these rows and set
108 // their diagonal entry to what we have
109 // determined above. note that for petsc
110 // matrices interleaving read with write
111 // operations is very expensive. thus, we
112 // here always replace the diagonal
113 // element, rather than first checking
114 // whether it is nonzero and in that case
115 // preserving it. this is different from
116 // the case of deal.II sparse matrices
117 // treated in the other functions.
118 if (eliminate_columns)
119 matrix.clear_rows_columns(constrained_rows,
120 average_nonzero_diagonal_entry);
121 else
122 matrix.clear_rows(constrained_rows, average_nonzero_diagonal_entry);
123
124 std::vector<types::global_dof_index> indices;
125 std::vector<PetscScalar> solution_values;
126 for (const auto &boundary_value : boundary_values)
127 if ((boundary_value.first >= local_range.first) &&
128 (boundary_value.first < local_range.second))
129 {
130 indices.push_back(boundary_value.first);
131 solution_values.push_back(boundary_value.second);
132 }
133 solution.set(indices, solution_values);
134
135 // now also set appropriate values for the rhs
136 for (auto &solution_value : solution_values)
137 solution_value *= average_nonzero_diagonal_entry;
138
139 right_hand_side.set(indices, solution_values);
140 }
141 else
142 {
143 // clear_rows() is a collective operation so we still have to call
144 // it:
145 std::vector<types::global_dof_index> constrained_rows;
146 if (eliminate_columns)
147 matrix.clear_rows_columns(constrained_rows, 1.);
148 else
149 matrix.clear_rows(constrained_rows, 1.);
150 }
151
152 // clean up
154 right_hand_side.compress(VectorOperation::insert);
155 }
156
157
158 void
160 const std::map<types::global_dof_index, PetscScalar> &boundary_values,
163 PETScWrappers::MPI::BlockVector & right_hand_side,
164 const bool eliminate_columns)
165 {
166 Assert(eliminate_columns == false, ExcNotImplemented());
167 Assert(matrix.n() == right_hand_side.size(),
168 ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
169 Assert(matrix.n() == solution.size(),
170 ExcDimensionMismatch(matrix.n(), solution.size()));
171 Assert(matrix.n_block_rows() == matrix.n_block_cols(), ExcNotQuadratic());
172
173 const unsigned int n_blocks = matrix.n_block_rows();
174
175 // We need to find the subdivision
176 // into blocks for the boundary values.
177 // To this end, generate a vector of
178 // maps with the respective indices.
179 std::vector<std::map<::types::global_dof_index, PetscScalar>>
180 block_boundary_values(n_blocks);
181 {
182 int block = 0;
183 ::types::global_dof_index offset = 0;
184 for (const auto &boundary_value : boundary_values)
185 {
186 if (boundary_value.first >= matrix.block(block, 0).m() + offset)
187 {
188 offset += matrix.block(block, 0).m();
189 block++;
190 }
191 const types::global_dof_index index = boundary_value.first - offset;
192 block_boundary_values[block].insert(
193 std::pair<types::global_dof_index, PetscScalar>(
194 index, boundary_value.second));
195 }
196 }
197
198 // Now call the non-block variants on
199 // the diagonal subblocks and the
200 // solution/rhs.
201 for (unsigned int block = 0; block < n_blocks; ++block)
202 apply_boundary_values(block_boundary_values[block],
203 matrix.block(block, block),
204 solution.block(block),
205 right_hand_side.block(block),
206 eliminate_columns);
207
208 // Finally, we need to do something
209 // about the off-diagonal matrices. This
210 // is luckily not difficult. Just clear
211 // the whole row.
212 for (unsigned int block_m = 0; block_m < n_blocks; ++block_m)
213 {
214 const std::pair<types::global_dof_index, types::global_dof_index>
215 local_range = matrix.block(block_m, 0).local_range();
216
217 std::vector<types::global_dof_index> constrained_rows;
218 for (std::map<types::global_dof_index, PetscScalar>::const_iterator
219 dof = block_boundary_values[block_m].begin();
220 dof != block_boundary_values[block_m].end();
221 ++dof)
222 if ((dof->first >= local_range.first) &&
223 (dof->first < local_range.second))
224 constrained_rows.push_back(dof->first);
225
226 for (unsigned int block_n = 0; block_n < n_blocks; ++block_n)
227 if (block_m != block_n)
228 matrix.block(block_m, block_n).clear_rows(constrained_rows);
229 }
230 }
231
232#endif
233
234
235
236#ifdef DEAL_II_WITH_TRILINOS
237
238 namespace internal
239 {
241 {
242 template <typename TrilinosMatrix, typename TrilinosVector>
243 void
245 TrilinosScalar> &boundary_values,
246 TrilinosMatrix & matrix,
247 TrilinosVector & solution,
248 TrilinosVector & right_hand_side,
249 const bool eliminate_columns)
250 {
251 Assert(eliminate_columns == false, ExcNotImplemented());
252 (void)eliminate_columns;
253
254 Assert(matrix.n() == right_hand_side.size(),
255 ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
256 Assert(matrix.n() == solution.size(),
257 ExcDimensionMismatch(matrix.m(), solution.size()));
258
259 // if no boundary values are to be applied, then
260 // jump straight to the compress() calls that we still have
261 // to perform because they are collective operations
262 if (boundary_values.size() > 0)
263 {
264 const std::pair<types::global_dof_index, types::global_dof_index>
265 local_range = matrix.local_range();
266 Assert(local_range == right_hand_side.local_range(),
268 Assert(local_range == solution.local_range(), ExcInternalError());
269
270 // determine the first nonzero diagonal
271 // entry from within the part of the
272 // matrix that we can see. if we can't
273 // find such an entry, take one
274 TrilinosScalar average_nonzero_diagonal_entry = 1;
275 for (types::global_dof_index i = local_range.first;
276 i < local_range.second;
277 ++i)
278 if (matrix.diag_element(i) != 0)
279 {
280 average_nonzero_diagonal_entry =
281 std::fabs(matrix.diag_element(i));
282 break;
283 }
284
285 // figure out which rows of the matrix we
286 // have to eliminate on this processor
287 std::vector<types::global_dof_index> constrained_rows;
288 for (const auto &boundary_value : boundary_values)
289 if ((boundary_value.first >= local_range.first) &&
290 (boundary_value.first < local_range.second))
291 constrained_rows.push_back(boundary_value.first);
292
293 // then eliminate these rows and
294 // set their diagonal entry to
295 // what we have determined
296 // above. if the value already is
297 // nonzero, it will be preserved,
298 // in accordance with the basic
299 // matrix classes in deal.II.
300 matrix.clear_rows(constrained_rows, average_nonzero_diagonal_entry);
301
302 std::vector<types::global_dof_index> indices;
303 std::vector<TrilinosScalar> solution_values;
304 for (const auto &boundary_value : boundary_values)
305 if ((boundary_value.first >= local_range.first) &&
306 (boundary_value.first < local_range.second))
307 {
308 indices.push_back(boundary_value.first);
309 solution_values.push_back(boundary_value.second);
310 }
311 solution.set(indices, solution_values);
312
313 // now also set appropriate
314 // values for the rhs
315 for (unsigned int i = 0; i < solution_values.size(); ++i)
316 solution_values[i] *= matrix.diag_element(indices[i]);
317
318 right_hand_side.set(indices, solution_values);
319 }
320 else
321 {
322 // clear_rows() is a collective operation so we still have to call
323 // it:
324 std::vector<types::global_dof_index> constrained_rows;
325 matrix.clear_rows(constrained_rows, 1.);
326 }
327
328 // clean up
329 matrix.compress(VectorOperation::insert);
330 solution.compress(VectorOperation::insert);
331 right_hand_side.compress(VectorOperation::insert);
332 }
333
334
335
336 template <typename TrilinosMatrix, typename TrilinosBlockVector>
337 void
339 const std::map<types::global_dof_index, TrilinosScalar>
340 & boundary_values,
341 TrilinosMatrix & matrix,
342 TrilinosBlockVector &solution,
343 TrilinosBlockVector &right_hand_side,
344 const bool eliminate_columns)
345 {
346 Assert(eliminate_columns == false, ExcNotImplemented());
347
348 Assert(matrix.n() == right_hand_side.size(),
349 ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
350 Assert(matrix.n() == solution.size(),
351 ExcDimensionMismatch(matrix.n(), solution.size()));
352 Assert(matrix.n_block_rows() == matrix.n_block_cols(),
354
355 const unsigned int n_blocks = matrix.n_block_rows();
356
357 // We need to find the subdivision
358 // into blocks for the boundary values.
359 // To this end, generate a vector of
360 // maps with the respective indices.
361 std::vector<std::map<types::global_dof_index, TrilinosScalar>>
362 block_boundary_values(n_blocks);
363 {
364 int block = 0;
365 types::global_dof_index offset = 0;
366 for (const auto &boundary_value : boundary_values)
367 {
368 if (boundary_value.first >= matrix.block(block, 0).m() + offset)
369 {
370 offset += matrix.block(block, 0).m();
371 block++;
372 }
373 const types::global_dof_index index =
374 boundary_value.first - offset;
375 block_boundary_values[block].insert(
376 std::pair<types::global_dof_index, TrilinosScalar>(
377 index, boundary_value.second));
378 }
379 }
380
381 // Now call the non-block variants on
382 // the diagonal subblocks and the
383 // solution/rhs.
384 for (unsigned int block = 0; block < n_blocks; ++block)
385 TrilinosWrappers::apply_boundary_values(block_boundary_values[block],
386 matrix.block(block, block),
387 solution.block(block),
388 right_hand_side.block(block),
389 eliminate_columns);
390
391 // Finally, we need to do something
392 // about the off-diagonal matrices. This
393 // is luckily not difficult. Just clear
394 // the whole row.
395 for (unsigned int block_m = 0; block_m < n_blocks; ++block_m)
396 {
397 const std::pair<types::global_dof_index, types::global_dof_index>
398 local_range = matrix.block(block_m, 0).local_range();
399
400 std::vector<types::global_dof_index> constrained_rows;
401 for (std::map<types::global_dof_index,
402 TrilinosScalar>::const_iterator dof =
403 block_boundary_values[block_m].begin();
404 dof != block_boundary_values[block_m].end();
405 ++dof)
406 if ((dof->first >= local_range.first) &&
407 (dof->first < local_range.second))
408 constrained_rows.push_back(dof->first);
409
410 for (unsigned int block_n = 0; block_n < n_blocks; ++block_n)
411 if (block_m != block_n)
412 matrix.block(block_m, block_n).clear_rows(constrained_rows);
413 }
414 }
415 } // namespace TrilinosWrappers
416 } // namespace internal
417
418
419
420 void
422 const std::map<types::global_dof_index, TrilinosScalar> &boundary_values,
425 TrilinosWrappers::MPI::Vector & right_hand_side,
426 const bool eliminate_columns)
427 {
428 // simply redirect to the generic function
429 // used for both trilinos matrix types
431 boundary_values, matrix, solution, right_hand_side, eliminate_columns);
432 }
433
434
435
436 void
438 const std::map<types::global_dof_index, TrilinosScalar> &boundary_values,
441 TrilinosWrappers::MPI::BlockVector & right_hand_side,
442 const bool eliminate_columns)
443 {
445 boundary_values, matrix, solution, right_hand_side, eliminate_columns);
446 }
447
448#endif
449
450} // namespace MatrixTools
451
std::size_t size() const
BlockType & block(const unsigned int i)
std::pair< size_type, size_type > local_range() const
void compress(const VectorOperation::values operation)
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
void apply_block_boundary_values(const std::map< types::global_dof_index, TrilinosScalar > &boundary_values, TrilinosMatrix &matrix, TrilinosBlockVector &solution, TrilinosBlockVector &right_hand_side, const bool eliminate_columns)
void apply_boundary_values(const std::map< types::global_dof_index, TrilinosScalar > &boundary_values, TrilinosMatrix &matrix, TrilinosVector &solution, TrilinosVector &right_hand_side, 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)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
double TrilinosScalar
Definition types.h:175