Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
matrix_tools_once.cc
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>
17 #include <deal.II/base/geometry_info.h>
18 #include <deal.II/base/quadrature.h>
19 #include <deal.II/base/work_stream.h>
20 
21 #include <deal.II/dofs/dof_accessor.h>
22 #include <deal.II/dofs/dof_handler.h>
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 #include <deal.II/fe/mapping_q1.h>
28 
29 #include <deal.II/grid/tria_iterator.h>
30 
31 #include <deal.II/hp/fe_values.h>
32 #include <deal.II/hp/mapping_collection.h>
33 
34 #include <deal.II/lac/block_sparse_matrix.h>
35 #include <deal.II/lac/block_vector.h>
36 #include <deal.II/lac/full_matrix.h>
37 #include <deal.II/lac/sparse_matrix.h>
38 #include <deal.II/lac/vector.h>
39 
40 #include <deal.II/numerics/matrix_tools.h>
41 
42 #ifdef DEAL_II_WITH_PETSC
43 # include <deal.II/lac/petsc_block_sparse_matrix.h>
44 # include <deal.II/lac/petsc_block_vector.h>
45 # include <deal.II/lac/petsc_matrix_base.h>
46 # include <deal.II/lac/petsc_vector_base.h>
47 #endif
48 
49 #ifdef DEAL_II_WITH_TRILINOS
50 # include <deal.II/lac/trilinos_block_sparse_matrix.h>
51 # include <deal.II/lac/trilinos_parallel_block_vector.h>
52 # include <deal.II/lac/trilinos_sparse_matrix.h>
53 # include <deal.II/lac/trilinos_vector.h>
54 #endif
55 
56 #include <algorithm>
57 #include <cmath>
58 
59 
60 DEAL_II_NAMESPACE_OPEN
61 
62 namespace MatrixTools
63 {
64 #ifdef DEAL_II_WITH_PETSC
65  void
67  const std::map<types::global_dof_index, PetscScalar> &boundary_values,
69  PETScWrappers::VectorBase & solution,
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()));
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  {
239  namespace TrilinosWrappers
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(),
266  ExcInternalError());
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
337  apply_block_boundary_values(
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(),
352  ExcNotQuadratic());
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
429  internal::TrilinosWrappers::apply_boundary_values(
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  {
443  internal::TrilinosWrappers::apply_block_boundary_values(
444  boundary_values, matrix, solution, right_hand_side, eliminate_columns);
445  }
446 
447 #endif
448 
449 } // namespace MatrixTools
450 
451 DEAL_II_NAMESPACE_CLOSE
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:81
Contents is actually a matrix.
std::pair< size_type, size_type > local_range() const
std::size_t size() const
void compress(const VectorOperation::values operation)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
unsigned int global_dof_index
Definition: types.h:89
static ::ExceptionBase & ExcNotImplemented()
BlockType & block(const unsigned int i)
static ::ExceptionBase & ExcInternalError()
void set(const std::vector< size_type > &indices, const std::vector< PetscScalar > &values)