deal.II version GIT relicensing-3535-g02565b63bb 2025-06-19 13:50: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_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(),
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;
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 Assert(matrix.n() == right_hand_side.size(),
252 ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
253 Assert(matrix.n() == solution.size(),
254 ExcDimensionMismatch(matrix.m(), solution.size()));
255
256 // if no boundary values are to be applied, then
257 // jump straight to the compress() calls that we still have
258 // to perform because they are collective operations
259 if (boundary_values.size() > 0)
260 {
261 const std::pair<types::global_dof_index, types::global_dof_index>
262 local_range = matrix.local_range();
263 Assert(local_range == right_hand_side.local_range(),
265 Assert(local_range == solution.local_range(), ExcInternalError());
266
267 // determine the first nonzero diagonal
268 // entry from within the part of the
269 // matrix that we can see. if we can't
270 // find such an entry, take one
271 TrilinosScalar average_nonzero_diagonal_entry = 1;
273 i < local_range.second;
274 ++i)
275 if (matrix.diag_element(i) != 0)
276 {
277 average_nonzero_diagonal_entry =
278 std::fabs(matrix.diag_element(i));
279 break;
280 }
281
282 // figure out which rows of the matrix we
283 // have to eliminate on this processor
284 std::vector<types::global_dof_index> constrained_rows;
285 for (const auto &boundary_value : boundary_values)
286 if ((boundary_value.first >= local_range.first) &&
287 (boundary_value.first < local_range.second))
288 constrained_rows.push_back(boundary_value.first);
289
290 // then eliminate these rows and
291 // set their diagonal entry to
292 // what we have determined
293 // above. if the value already is
294 // nonzero, it will be preserved,
295 // in accordance with the basic
296 // matrix classes in deal.II.
297 matrix.clear_rows(constrained_rows, average_nonzero_diagonal_entry);
298
299 std::vector<types::global_dof_index> indices;
300 std::vector<TrilinosScalar> solution_values;
301 for (const auto &boundary_value : boundary_values)
302 if ((boundary_value.first >= local_range.first) &&
303 (boundary_value.first < local_range.second))
304 {
305 indices.push_back(boundary_value.first);
306 solution_values.push_back(boundary_value.second);
307 }
308 solution.set(indices, solution_values);
309
310 // now also set appropriate
311 // values for the rhs
312 for (unsigned int i = 0; i < solution_values.size(); ++i)
313 solution_values[i] *= matrix.diag_element(indices[i]);
314
315 right_hand_side.set(indices, solution_values);
316 }
317 else
318 {
319 // clear_rows() is a collective operation so we still have to call
320 // it:
321 std::vector<types::global_dof_index> constrained_rows;
322 matrix.clear_rows(constrained_rows, 1.);
323 }
324
325 // clean up
326 matrix.compress(VectorOperation::insert);
327 solution.compress(VectorOperation::insert);
328 right_hand_side.compress(VectorOperation::insert);
329 }
330
331
332
333 template <typename TrilinosMatrix, typename TrilinosBlockVector>
334 void
336 const std::map<types::global_dof_index, TrilinosScalar>
337 &boundary_values,
338 TrilinosMatrix &matrix,
339 TrilinosBlockVector &solution,
340 TrilinosBlockVector &right_hand_side,
341 const bool eliminate_columns)
342 {
343 Assert(eliminate_columns == false, ExcNotImplemented());
344
345 Assert(matrix.n() == right_hand_side.size(),
346 ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
347 Assert(matrix.n() == solution.size(),
348 ExcDimensionMismatch(matrix.n(), solution.size()));
349 Assert(matrix.n_block_rows() == matrix.n_block_cols(),
351
352 const unsigned int n_blocks = matrix.n_block_rows();
353
354 // We need to find the subdivision
355 // into blocks for the boundary values.
356 // To this end, generate a vector of
357 // maps with the respective indices.
358 std::vector<std::map<types::global_dof_index, TrilinosScalar>>
359 block_boundary_values(n_blocks);
360 {
361 int block = 0;
362 types::global_dof_index offset = 0;
363 for (const auto &boundary_value : boundary_values)
364 {
365 if (boundary_value.first >= matrix.block(block, 0).m() + offset)
366 {
367 offset += matrix.block(block, 0).m();
368 ++block;
369 }
370 const types::global_dof_index index =
371 boundary_value.first - offset;
372 block_boundary_values[block].insert(
373 std::pair<types::global_dof_index, TrilinosScalar>(
374 index, boundary_value.second));
375 }
376 }
377
378 // Now call the non-block variants on
379 // the diagonal subblocks and the
380 // solution/rhs.
381 for (unsigned int block = 0; block < n_blocks; ++block)
382 TrilinosWrappers::apply_boundary_values(block_boundary_values[block],
383 matrix.block(block, block),
384 solution.block(block),
385 right_hand_side.block(block),
386 eliminate_columns);
387
388 // Finally, we need to do something
389 // about the off-diagonal matrices. This
390 // is luckily not difficult. Just clear
391 // the whole row.
392 for (unsigned int block_m = 0; block_m < n_blocks; ++block_m)
393 {
394 const std::pair<types::global_dof_index, types::global_dof_index>
395 local_range = matrix.block(block_m, 0).local_range();
396
397 std::vector<types::global_dof_index> constrained_rows;
398 for (std::map<types::global_dof_index,
399 TrilinosScalar>::const_iterator dof =
400 block_boundary_values[block_m].begin();
401 dof != block_boundary_values[block_m].end();
402 ++dof)
403 if ((dof->first >= local_range.first) &&
404 (dof->first < local_range.second))
405 constrained_rows.push_back(dof->first);
406
407 for (unsigned int block_n = 0; block_n < n_blocks; ++block_n)
408 if (block_m != block_n)
409 matrix.block(block_m, block_n).clear_rows(constrained_rows);
410 }
411 }
412 } // namespace TrilinosWrappers
413 } // namespace internal
414
415
416
417 void
419 const std::map<types::global_dof_index, TrilinosScalar> &boundary_values,
422 TrilinosWrappers::MPI::Vector &right_hand_side,
423 const bool eliminate_columns)
424 {
425 // simply redirect to the generic function
426 // used for both trilinos matrix types
428 boundary_values, matrix, solution, right_hand_side, eliminate_columns);
429 }
430
431
432
433 void
435 const std::map<types::global_dof_index, TrilinosScalar> &boundary_values,
438 TrilinosWrappers::MPI::BlockVector &right_hand_side,
439 const bool eliminate_columns)
440 {
442 boundary_values, matrix, solution, right_hand_side, eliminate_columns);
443 }
444
445#endif
446
447} // namespace MatrixTools
448
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:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
std::pair< types::global_dof_index, types::global_dof_index > local_range
Definition mpi.cc:830
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:198