Reference documentation for deal.II version GIT 75f1417c0d 2023-02-03 16: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\}}\)
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
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 
16 #include <deal.II/base/function.h>
20 
23 #include <deal.II/dofs/dof_tools.h>
24 
25 #include <deal.II/fe/fe.h>
26 #include <deal.II/fe/fe_values.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 
58 namespace MatrixTools
59 {
60 #ifdef DEAL_II_WITH_PETSC
61  void
63  const std::map<types::global_dof_index, PetscScalar> &boundary_values,
65  PETScWrappers::VectorBase & solution,
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(),
267  ExcInternalError());
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
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(),
353  ExcNotQuadratic());
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:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1583
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNotQuadratic()
Expression fabs(const Expression &x)
@ matrix
Contents is actually a matrix.
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
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)
Definition: matrix_tools.cc:76
VectorType::value_type * begin(VectorType &V)
unsigned int global_dof_index
Definition: types.h:81