Reference documentation for deal.II version 8.4.1
matrix_tools_once.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2015 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/function.h>
17 #include <deal.II/base/quadrature.h>
18 #include <deal.II/base/work_stream.h>
19 #include <deal.II/base/geometry_info.h>
20 #include <deal.II/base/quadrature.h>
21 #include <deal.II/dofs/dof_handler.h>
22 #include <deal.II/dofs/dof_accessor.h>
23 #include <deal.II/dofs/dof_tools.h>
24 #include <deal.II/fe/fe.h>
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping_q1.h>
27 #include <deal.II/grid/tria_iterator.h>
28 #include <deal.II/hp/fe_values.h>
29 #include <deal.II/hp/mapping_collection.h>
30 #include <deal.II/numerics/matrix_tools.h>
31 #include <deal.II/lac/vector.h>
32 #include <deal.II/lac/block_vector.h>
33 #include <deal.II/lac/full_matrix.h>
34 #include <deal.II/lac/sparse_matrix.h>
35 #include <deal.II/lac/block_sparse_matrix.h>
36 
37 #ifdef DEAL_II_WITH_PETSC
38 # include <deal.II/lac/petsc_parallel_sparse_matrix.h>
39 # include <deal.II/lac/petsc_sparse_matrix.h>
40 # include <deal.II/lac/petsc_parallel_vector.h>
41 # include <deal.II/lac/petsc_vector.h>
42 # include <deal.II/lac/petsc_parallel_block_sparse_matrix.h>
43 #endif
44 
45 #ifdef DEAL_II_WITH_TRILINOS
46 # include <deal.II/lac/trilinos_sparse_matrix.h>
47 # include <deal.II/lac/trilinos_vector.h>
48 # include <deal.II/lac/trilinos_block_sparse_matrix.h>
49 # include <deal.II/lac/trilinos_block_vector.h>
50 #endif
51 
52 #include <algorithm>
53 #include <set>
54 #include <cmath>
55 
56 
57 DEAL_II_NAMESPACE_OPEN
58 
59 namespace MatrixTools
60 {
61 
62 #ifdef DEAL_II_WITH_PETSC
63 
64  namespace internal
65  {
66  namespace PETScWrappers
67  {
68  template <typename PETScMatrix, typename PETScVector>
69  void
70  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
71  PETScMatrix &matrix,
72  PETScVector &solution,
73  PETScVector &right_hand_side,
74  const bool eliminate_columns)
75  {
76  (void)eliminate_columns;
77  Assert (eliminate_columns == false, ExcNotImplemented());
78 
79  Assert (matrix.n() == right_hand_side.size(),
80  ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
81  Assert (matrix.n() == solution.size(),
82  ExcDimensionMismatch(matrix.n(), solution.size()));
83 
84  // if no boundary values are to be applied, then
85  // jump straight to the compress() calls that we still have
86  // to perform because they are collective operations
87  if (boundary_values.size() > 0)
88  {
89  const std::pair<types::global_dof_index, types::global_dof_index> local_range
90  = matrix.local_range();
91  Assert (local_range == right_hand_side.local_range(),
92  ExcInternalError());
93  Assert (local_range == solution.local_range(),
94  ExcInternalError());
95 
96  // determine the first nonzero diagonal
97  // entry from within the part of the
98  // matrix that we can see. if we can't
99  // find such an entry, take one
100  PetscScalar average_nonzero_diagonal_entry = 1;
101  for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
102  if (matrix.diag_element(i) != PetscScalar ())
103  {
104  average_nonzero_diagonal_entry = std::abs(matrix.diag_element(i));
105  break;
106  }
107 
108  // figure out which rows of the matrix we
109  // have to eliminate on this processor
110  std::vector<types::global_dof_index> constrained_rows;
111  for (std::map<types::global_dof_index,double>::const_iterator
112  dof = boundary_values.begin();
113  dof != boundary_values.end();
114  ++dof)
115  if ((dof->first >= local_range.first) &&
116  (dof->first < local_range.second))
117  constrained_rows.push_back (dof->first);
118 
119  // then eliminate these rows and set
120  // their diagonal entry to what we have
121  // determined above. note that for petsc
122  // matrices interleaving read with write
123  // operations is very expensive. thus, we
124  // here always replace the diagonal
125  // element, rather than first checking
126  // whether it is nonzero and in that case
127  // preserving it. this is different from
128  // the case of deal.II sparse matrices
129  // treated in the other functions.
130  matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
131 
132  std::vector<types::global_dof_index> indices;
133  std::vector<PetscScalar> solution_values;
134  for (std::map<types::global_dof_index,double>::const_iterator
135  dof = boundary_values.begin();
136  dof != boundary_values.end();
137  ++dof)
138  if ((dof->first >= local_range.first) &&
139  (dof->first < local_range.second))
140  {
141  indices.push_back (dof->first);
142  solution_values.push_back (dof->second);
143  }
144  solution.set (indices, solution_values);
145 
146  // now also set appropriate values for
147  // the rhs
148  for (unsigned int i=0; i<solution_values.size(); ++i)
149  solution_values[i] *= average_nonzero_diagonal_entry;
150 
151  right_hand_side.set (indices, solution_values);
152  }
153  else
154  {
155  // clear_rows() is a collective operation so we still have to call
156  // it:
157  std::vector<types::global_dof_index> constrained_rows;
158  matrix.clear_rows (constrained_rows, 1.);
159  }
160 
161  // clean up
162  solution.compress (VectorOperation::insert);
163  right_hand_side.compress (VectorOperation::insert);
164  }
165  }
166  }
167 
168 
169 
170 
171 
172  void
173  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
175  PETScWrappers::Vector &solution,
176  PETScWrappers::Vector &right_hand_side,
177  const bool eliminate_columns)
178  {
179  // simply redirect to the generic function
180  // used for both petsc matrix types
181  internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
182  right_hand_side, eliminate_columns);
183  }
184 
185 
186 
187  void
188  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
190  PETScWrappers::MPI::Vector &solution,
191  PETScWrappers::MPI::Vector &right_hand_side,
192  const bool eliminate_columns)
193  {
194  // simply redirect to the generic function
195  // used for both petsc matrix types
196  internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
197  right_hand_side, eliminate_columns);
198  }
199 
200 
201  void
202  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
205  PETScWrappers::MPI::BlockVector &right_hand_side,
206  const bool eliminate_columns)
207  {
208  Assert (matrix.n() == right_hand_side.size(),
209  ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
210  Assert (matrix.n() == solution.size(),
211  ExcDimensionMismatch(matrix.n(), solution.size()));
212  Assert (matrix.n_block_rows() == matrix.n_block_cols(),
213  ExcNotQuadratic());
214 
215  const unsigned int n_blocks = matrix.n_block_rows();
216 
217  // We need to find the subdivision
218  // into blocks for the boundary values.
219  // To this end, generate a vector of
220  // maps with the respective indices.
221  std::vector<std::map<::types::global_dof_index,double> > block_boundary_values(n_blocks);
222  {
223  int block = 0;
224  ::types::global_dof_index offset = 0;
225  for (std::map<types::global_dof_index,double>::const_iterator
226  dof = boundary_values.begin();
227  dof != boundary_values.end();
228  ++dof)
229  {
230  if (dof->first >= matrix.block(block,0).m() + offset)
231  {
232  offset += matrix.block(block,0).m();
233  block++;
234  }
235  const types::global_dof_index index = dof->first - offset;
236  block_boundary_values[block].insert(std::pair<types::global_dof_index, double> (index,dof->second));
237  }
238  }
239 
240  // Now call the non-block variants on
241  // the diagonal subblocks and the
242  // solution/rhs.
243  for (unsigned int block=0; block<n_blocks; ++block)
244  internal::PETScWrappers::apply_boundary_values(block_boundary_values[block],
245  matrix.block(block,block),
246  solution.block(block),
247  right_hand_side.block(block),
248  eliminate_columns);
249 
250  // Finally, we need to do something
251  // about the off-diagonal matrices. This
252  // is luckily not difficult. Just clear
253  // the whole row.
254  for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
255  {
256  const std::pair<types::global_dof_index, types::global_dof_index> local_range
257  = matrix.block(block_m,0).local_range();
258 
259  std::vector<types::global_dof_index> constrained_rows;
260  for (std::map<types::global_dof_index,double>::const_iterator
261  dof = block_boundary_values[block_m].begin();
262  dof != block_boundary_values[block_m].end();
263  ++dof)
264  if ((dof->first >= local_range.first) &&
265  (dof->first < local_range.second))
266  constrained_rows.push_back (dof->first);
267 
268  for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
269  if (block_m != block_n)
270  matrix.block(block_m,block_n).clear_rows(constrained_rows);
271  }
272  }
273 
274 #endif
275 
276 
277 
278 #ifdef DEAL_II_WITH_TRILINOS
279 
280  namespace internal
281  {
282  namespace TrilinosWrappers
283  {
284  template <typename TrilinosMatrix, typename TrilinosVector>
285  void
286  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
287  TrilinosMatrix &matrix,
288  TrilinosVector &solution,
289  TrilinosVector &right_hand_side,
290  const bool eliminate_columns)
291  {
292  Assert (eliminate_columns == false, ExcNotImplemented());
293  (void)eliminate_columns;
294 
295  Assert (matrix.n() == right_hand_side.size(),
296  ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
297  Assert (matrix.n() == solution.size(),
298  ExcDimensionMismatch(matrix.m(), solution.size()));
299 
300  // if no boundary values are to be applied, then
301  // jump straight to the compress() calls that we still have
302  // to perform because they are collective operations
303  if (boundary_values.size() > 0)
304  {
305  const std::pair<types::global_dof_index, types::global_dof_index> local_range
306  = matrix.local_range();
307  Assert (local_range == right_hand_side.local_range(),
308  ExcInternalError());
309  Assert (local_range == solution.local_range(),
310  ExcInternalError());
311 
312  // determine the first nonzero diagonal
313  // entry from within the part of the
314  // matrix that we can see. if we can't
315  // find such an entry, take one
316  TrilinosScalar average_nonzero_diagonal_entry = 1;
317  for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
318  if (matrix.diag_element(i) != 0)
319  {
320  average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
321  break;
322  }
323 
324  // figure out which rows of the matrix we
325  // have to eliminate on this processor
326  std::vector<types::global_dof_index> constrained_rows;
327  for (std::map<types::global_dof_index,double>::const_iterator
328  dof = boundary_values.begin();
329  dof != boundary_values.end();
330  ++dof)
331  if ((dof->first >= local_range.first) &&
332  (dof->first < local_range.second))
333  constrained_rows.push_back (dof->first);
334 
335  // then eliminate these rows and
336  // set their diagonal entry to
337  // what we have determined
338  // above. if the value already is
339  // nonzero, it will be preserved,
340  // in accordance with the basic
341  // matrix classes in deal.II.
342  matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
343 
344  std::vector<types::global_dof_index> indices;
345  std::vector<TrilinosScalar> solution_values;
346  for (std::map<types::global_dof_index,double>::const_iterator
347  dof = boundary_values.begin();
348  dof != boundary_values.end();
349  ++dof)
350  if ((dof->first >= local_range.first) &&
351  (dof->first < local_range.second))
352  {
353  indices.push_back (dof->first);
354  solution_values.push_back (dof->second);
355  }
356  solution.set (indices, solution_values);
357 
358  // now also set appropriate
359  // values for the rhs
360  for (unsigned int i=0; i<solution_values.size(); ++i)
361  solution_values[i] *= matrix.diag_element(indices[i]);
362 
363  right_hand_side.set (indices, solution_values);
364  }
365  else
366  {
367  // clear_rows() is a collective operation so we still have to call
368  // it:
369  std::vector<types::global_dof_index> constrained_rows;
370  matrix.clear_rows (constrained_rows, 1.);
371  }
372 
373  // clean up
374  matrix.compress (VectorOperation::insert);
375  solution.compress (VectorOperation::insert);
376  right_hand_side.compress (VectorOperation::insert);
377  }
378 
379 
380 
381  template <typename TrilinosMatrix, typename TrilinosBlockVector>
382  void
383  apply_block_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
384  TrilinosMatrix &matrix,
385  TrilinosBlockVector &solution,
386  TrilinosBlockVector &right_hand_side,
387  const bool eliminate_columns)
388  {
389  Assert (eliminate_columns == false, ExcNotImplemented());
390 
391  Assert (matrix.n() == right_hand_side.size(),
392  ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
393  Assert (matrix.n() == solution.size(),
394  ExcDimensionMismatch(matrix.n(), solution.size()));
395  Assert (matrix.n_block_rows() == matrix.n_block_cols(),
396  ExcNotQuadratic());
397 
398  const unsigned int n_blocks = matrix.n_block_rows();
399 
400  // We need to find the subdivision
401  // into blocks for the boundary values.
402  // To this end, generate a vector of
403  // maps with the respective indices.
404  std::vector<std::map<types::global_dof_index,double> > block_boundary_values(n_blocks);
405  {
406  int block=0;
407  types::global_dof_index offset = 0;
408  for (std::map<types::global_dof_index,double>::const_iterator
409  dof = boundary_values.begin();
410  dof != boundary_values.end();
411  ++dof)
412  {
413  if (dof->first >= matrix.block(block,0).m() + offset)
414  {
415  offset += matrix.block(block,0).m();
416  block++;
417  }
418  const types::global_dof_index index = dof->first - offset;
419  block_boundary_values[block].insert(
420  std::pair<types::global_dof_index, double> (index,dof->second));
421  }
422  }
423 
424  // Now call the non-block variants on
425  // the diagonal subblocks and the
426  // solution/rhs.
427  for (unsigned int block=0; block<n_blocks; ++block)
428  TrilinosWrappers::apply_boundary_values(block_boundary_values[block],
429  matrix.block(block,block),
430  solution.block(block),
431  right_hand_side.block(block),
432  eliminate_columns);
433 
434  // Finally, we need to do something
435  // about the off-diagonal matrices. This
436  // is luckily not difficult. Just clear
437  // the whole row.
438  for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
439  {
440  const std::pair<types::global_dof_index, types::global_dof_index> local_range
441  = matrix.block(block_m,0).local_range();
442 
443  std::vector<types::global_dof_index> constrained_rows;
444  for (std::map<types::global_dof_index,double>::const_iterator
445  dof = block_boundary_values[block_m].begin();
446  dof != block_boundary_values[block_m].end();
447  ++dof)
448  if ((dof->first >= local_range.first) &&
449  (dof->first < local_range.second))
450  constrained_rows.push_back (dof->first);
451 
452  for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
453  if (block_m != block_n)
454  matrix.block(block_m,block_n).clear_rows(constrained_rows);
455  }
456  }
457  }
458  }
459 
460 
461 
462 
463  void
464  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
466  TrilinosWrappers::Vector &solution,
467  TrilinosWrappers::Vector &right_hand_side,
468  const bool eliminate_columns)
469  {
470  // simply redirect to the generic function
471  // used for both trilinos matrix types
472  internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
473  right_hand_side, eliminate_columns);
474  }
475 
476 
477 
478  void
479  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
482  TrilinosWrappers::MPI::Vector &right_hand_side,
483  const bool eliminate_columns)
484  {
485  // simply redirect to the generic function
486  // used for both trilinos matrix types
487  internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
488  right_hand_side, eliminate_columns);
489  }
490 
491 
492 
493  void
494  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
497  TrilinosWrappers::BlockVector &right_hand_side,
498  const bool eliminate_columns)
499  {
500  internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
501  solution, right_hand_side,
502  eliminate_columns);
503  }
504 
505 
506 
507  void
508  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
511  TrilinosWrappers::MPI::BlockVector &right_hand_side,
512  const bool eliminate_columns)
513  {
514  internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
515  solution, right_hand_side,
516  eliminate_columns);
517  }
518 
519 #endif
520 
521 }
522 
523 DEAL_II_NAMESPACE_CLOSE
std::size_t size() const
size_type n() const
size_type m() const
unsigned int n_block_cols() const
unsigned int global_dof_index
Definition: types.h:88
#define Assert(cond, exc)
Definition: exceptions.h:294
BlockType & block(const unsigned int row, const unsigned int column)
BlockType & block(const unsigned int i)
void apply_boundary_values(const std::map< types::global_dof_index, double > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Definition: matrix_tools.cc:77
unsigned int n_block_rows() const