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