Reference documentation for deal.II version 9.3.3
matrix_tools_once.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2019 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
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>
28
30
33
38#include <deal.II/lac/vector.h>
39
41
42#ifdef DEAL_II_WITH_PETSC
47#endif
48
49#ifdef DEAL_II_WITH_TRILINOS
54#endif
55
56#include <algorithm>
57#include <cmath>
58
59
61
62namespace MatrixTools
63{
64#ifdef DEAL_II_WITH_PETSC
65 void
67 const std::map<types::global_dof_index, PetscScalar> &boundary_values,
70 PETScWrappers::VectorBase & right_hand_side,
71 const bool eliminate_columns)
72 {
73 (void)eliminate_columns;
74 Assert(eliminate_columns == false, ExcNotImplemented());
75
76 Assert(matrix.n() == right_hand_side.size(),
77 ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
78 Assert(matrix.n() == solution.size(),
79 ExcDimensionMismatch(matrix.n(), solution.size()));
80
81 // if no boundary values are to be applied, then
82 // jump straight to the compress() calls that we still have
83 // to perform because they are collective operations
84 if (boundary_values.size() > 0)
85 {
86 const std::pair<types::global_dof_index, types::global_dof_index>
87 local_range = matrix.local_range();
88 Assert(local_range == right_hand_side.local_range(),
90 Assert(local_range == solution.local_range(), ExcInternalError());
91
92 // determine the first nonzero diagonal
93 // entry from within the part of the
94 // matrix that we can see. if we can't
95 // find such an entry, take one
96 PetscScalar average_nonzero_diagonal_entry = 1;
97 for (types::global_dof_index i = local_range.first;
98 i < local_range.second;
99 ++i)
100 if (matrix.diag_element(i) != PetscScalar())
101 {
102 average_nonzero_diagonal_entry = std::abs(matrix.diag_element(i));
103 break;
104 }
105
106 // figure out which rows of the matrix we
107 // have to eliminate on this processor
108 std::vector<types::global_dof_index> constrained_rows;
109 for (const auto &boundary_value : boundary_values)
110 if ((boundary_value.first >= local_range.first) &&
111 (boundary_value.first < local_range.second))
112 constrained_rows.push_back(boundary_value.first);
113
114 // then eliminate these rows and set
115 // their diagonal entry to what we have
116 // determined above. note that for petsc
117 // matrices interleaving read with write
118 // operations is very expensive. thus, we
119 // here always replace the diagonal
120 // element, rather than first checking
121 // whether it is nonzero and in that case
122 // preserving it. this is different from
123 // the case of deal.II sparse matrices
124 // treated in the other functions.
125 matrix.clear_rows(constrained_rows, average_nonzero_diagonal_entry);
126
127 std::vector<types::global_dof_index> indices;
128 std::vector<PetscScalar> solution_values;
129 for (const auto &boundary_value : boundary_values)
130 if ((boundary_value.first >= local_range.first) &&
131 (boundary_value.first < local_range.second))
132 {
133 indices.push_back(boundary_value.first);
134 solution_values.push_back(boundary_value.second);
135 }
136 solution.set(indices, solution_values);
137
138 // now also set appropriate values for the rhs
139 for (auto &solution_value : solution_values)
140 solution_value *= average_nonzero_diagonal_entry;
141
142 right_hand_side.set(indices, solution_values);
143 }
144 else
145 {
146 // clear_rows() is a collective operation so we still have to call
147 // it:
148 std::vector<types::global_dof_index> constrained_rows;
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(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()));
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
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
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:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcNotImplemented()
std::size_t size() const
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcInternalError()
BlockType & block(const unsigned int i)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)