Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mg_transfer_block.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 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 
17 #include <deal.II/base/logstream.h>
18 
19 #include <deal.II/dofs/dof_accessor.h>
20 #include <deal.II/dofs/dof_tools.h>
21 
22 #include <deal.II/fe/fe.h>
23 
24 #include <deal.II/grid/tria.h>
25 #include <deal.II/grid/tria_iterator.h>
26 
27 #include <deal.II/lac/block_sparse_matrix.h>
28 #include <deal.II/lac/block_vector.h>
29 #include <deal.II/lac/sparse_matrix.h>
30 #include <deal.II/lac/vector.h>
31 
32 #include <deal.II/multigrid/mg_tools.h>
33 #include <deal.II/multigrid/mg_transfer_block.h>
34 #include <deal.II/multigrid/mg_transfer_block.templates.h>
35 
36 #include <algorithm>
37 #include <iostream>
38 #include <numeric>
39 #include <utility>
40 
41 DEAL_II_NAMESPACE_OPEN
42 
43 namespace
44 {
52  template <int dim, typename number, int spacedim>
53  void
54  reinit_vector_by_blocks(
55  const ::DoFHandler<dim, spacedim> & mg_dof,
57  const std::vector<bool> & sel,
58  std::vector<std::vector<types::global_dof_index>> &ndofs)
59  {
60  std::vector<bool> selected = sel;
61  // Compute the number of blocks needed
62  const unsigned int n_selected =
63  std::accumulate(selected.begin(), selected.end(), 0u);
64 
65  if (ndofs.size() == 0)
66  {
67  std::vector<std::vector<types::global_dof_index>> new_dofs(
68  mg_dof.get_triangulation().n_levels(),
69  std::vector<types::global_dof_index>(selected.size()));
70  std::swap(ndofs, new_dofs);
71  MGTools::count_dofs_per_block(mg_dof, ndofs);
72  }
73 
74  for (unsigned int level = v.min_level(); level <= v.max_level(); ++level)
75  {
76  v[level].reinit(n_selected, 0);
77  unsigned int k = 0;
78  for (unsigned int i = 0;
79  i < selected.size() && (k < v[level].n_blocks());
80  ++i)
81  {
82  if (selected[i])
83  {
84  v[level].block(k++).reinit(ndofs[level][i]);
85  }
86  v[level].collect_sizes();
87  }
88  }
89  }
90 
91 
98  template <int dim, typename number, int spacedim>
99  void
100  reinit_vector_by_blocks(
101  const ::DoFHandler<dim, spacedim> & mg_dof,
103  const unsigned int selected_block,
104  std::vector<std::vector<types::global_dof_index>> &ndofs)
105  {
106  const unsigned int n_blocks = mg_dof.get_fe().n_blocks();
107  Assert(selected_block < n_blocks,
108  ExcIndexRange(selected_block, 0, n_blocks));
109 
110  std::vector<bool> selected(n_blocks, false);
111  selected[selected_block] = true;
112 
113  if (ndofs.size() == 0)
114  {
115  std::vector<std::vector<types::global_dof_index>> new_dofs(
116  mg_dof.get_triangulation().n_levels(),
117  std::vector<types::global_dof_index>(selected.size()));
118  std::swap(ndofs, new_dofs);
119  MGTools::count_dofs_per_block(mg_dof, ndofs);
120  }
121 
122  for (unsigned int level = v.min_level(); level <= v.max_level(); ++level)
123  {
124  v[level].reinit(ndofs[level][selected_block]);
125  }
126  }
127 } // namespace
128 
129 
130 template <typename number>
131 template <int dim, typename number2, int spacedim>
132 void
134  const DoFHandler<dim, spacedim> &mg_dof_handler,
136  const BlockVector<number2> & src) const
137 {
138  reinit_vector_by_blocks(mg_dof_handler, dst, selected_block, sizes);
139  // For MGTransferBlockSelect, the
140  // multilevel block is always the
141  // first, since only one block is
142  // selected.
143  for (unsigned int level = mg_dof_handler.get_triangulation().n_levels();
144  level != 0;)
145  {
146  --level;
147  for (IT i = copy_indices[selected_block][level].begin();
148  i != copy_indices[selected_block][level].end();
149  ++i)
150  dst[level](i->second) = src.block(selected_block)(i->first);
151  }
152 }
153 
154 
155 
156 template <typename number>
157 template <int dim, typename number2, int spacedim>
158 void
160  const DoFHandler<dim, spacedim> &mg_dof_handler,
162  const Vector<number2> & src) const
163 {
164  reinit_vector_by_blocks(mg_dof_handler, dst, selected_block, sizes);
165  // For MGTransferBlockSelect, the
166  // multilevel block is always the
167  // first, since only one block is selected.
168  for (unsigned int level = mg_dof_handler.get_triangulation().n_levels();
169  level != 0;)
170  {
171  --level;
172  for (IT i = copy_indices[selected_block][level].begin();
173  i != copy_indices[selected_block][level].end();
174  ++i)
175  dst[level](i->second) = src(i->first);
176  }
177 }
178 
179 
180 
181 template <typename number>
182 template <int dim, typename number2, int spacedim>
183 void
185  const DoFHandler<dim, spacedim> & mg_dof_handler,
187  const BlockVector<number2> & src) const
188 {
189  reinit_vector_by_blocks(mg_dof_handler, dst, selected, sizes);
190  for (unsigned int level = mg_dof_handler.get_triangulation().n_levels();
191  level != 0;)
192  {
193  --level;
194  for (unsigned int block = 0; block < selected.size(); ++block)
195  if (selected[block])
196  for (IT i = copy_indices[block][level].begin();
197  i != copy_indices[block][level].end();
198  ++i)
199  dst[level].block(mg_block[block])(i->second) =
200  src.block(block)(i->first);
201  }
202 }
203 
204 
205 
206 template <int dim, int spacedim>
207 void
209  const DoFHandler<dim, spacedim> &mg_dof)
210 {
211  const FiniteElement<dim> &fe = mg_dof.get_fe();
212  const unsigned int n_blocks = fe.n_blocks();
213  const unsigned int dofs_per_cell = fe.dofs_per_cell;
214  const unsigned int n_levels = mg_dof.get_triangulation().n_levels();
215 
216  Assert(selected.size() == n_blocks,
217  ExcDimensionMismatch(selected.size(), n_blocks));
218 
219  // Compute the mapping between real
220  // blocks and blocks used for
221  // multigrid computations.
222  mg_block.resize(n_blocks);
223  n_mg_blocks = 0;
224  for (unsigned int i = 0; i < n_blocks; ++i)
225  if (selected[i])
226  mg_block[i] = n_mg_blocks++;
227  else
229 
230  // Compute the lengths of all blocks
231  sizes.clear();
232  sizes.resize(n_levels, std::vector<types::global_dof_index>(fe.n_blocks()));
234 
235  // Fill some index vectors
236  // for later use.
238  // Compute start indices from sizes
239  for (auto &level_blocks : mg_block_start)
240  {
242  for (types::global_dof_index &level_block_start : level_blocks)
243  {
244  const types::global_dof_index t = level_block_start;
245  level_block_start = k;
246  k += t;
247  }
248  }
249 
250  block_start.resize(n_blocks);
252  static_cast<const DoFHandler<dim, spacedim> &>(mg_dof), block_start);
253 
255  for (types::global_dof_index &first_index : block_start)
256  {
257  const types::global_dof_index t = first_index;
258  first_index = k;
259  k += t;
260  }
261  // Build index vectors for
262  // copy_to_mg and
263  // copy_from_mg. These vectors must
264  // be prebuilt, since the
265  // get_dof_indices functions are
266  // too slow
267  copy_indices.resize(n_blocks);
268  for (unsigned int block = 0; block < n_blocks; ++block)
269  if (selected[block])
270  copy_indices[block].resize(n_levels);
271 
272  // Building the prolongation matrices starts here!
273 
274  // reset the size of the array of
275  // matrices. call resize(0) first,
276  // in order to delete all elements
277  // and clear their memory. then
278  // repopulate these arrays
279  //
280  // note that on resize(0), the
281  // shared_ptr class takes care of
282  // deleting the object it points to
283  // by itself
284  prolongation_matrices.resize(0);
285  prolongation_sparsities.resize(0);
286  prolongation_matrices.reserve(n_levels - 1);
287  prolongation_sparsities.reserve(n_levels - 1);
288 
289  for (unsigned int i = 0; i < n_levels - 1; ++i)
290  {
291  prolongation_sparsities.emplace_back(new BlockSparsityPattern);
293  }
294 
295  // two fields which will store the
296  // indices of the multigrid dofs
297  // for a cell and one of its children
298  std::vector<types::global_dof_index> dof_indices_parent(dofs_per_cell);
299  std::vector<types::global_dof_index> dof_indices_child(dofs_per_cell);
300 
301  // for each level: first build the
302  // sparsity pattern of the matrices
303  // and then build the matrices
304  // themselves. note that we only
305  // need to take care of cells on
306  // the coarser level which have
307  // children
308 
309  for (unsigned int level = 0; level < n_levels - 1; ++level)
310  {
311  // reset the dimension of the
312  // structure. note that for
313  // the number of entries per
314  // row, the number of parent
315  // dofs coupling to a child dof
316  // is necessary. this, is the
317  // number of degrees of freedom
318  // per cell
319  prolongation_sparsities[level]->reinit(n_blocks, n_blocks);
320  for (unsigned int i = 0; i < n_blocks; ++i)
321  for (unsigned int j = 0; j < n_blocks; ++j)
322  if (i == j)
323  prolongation_sparsities[level]->block(i, j).reinit(
324  sizes[level + 1][i], sizes[level][j], dofs_per_cell + 1);
325  else
326  prolongation_sparsities[level]->block(i, j).reinit(
327  sizes[level + 1][i], sizes[level][j], 0);
328 
329  prolongation_sparsities[level]->collect_sizes();
330 
331  for (typename DoFHandler<dim, spacedim>::cell_iterator cell =
332  mg_dof.begin(level);
333  cell != mg_dof.end(level);
334  ++cell)
335  if (cell->has_children())
336  {
337  cell->get_mg_dof_indices(dof_indices_parent);
338 
339  Assert(cell->n_children() ==
342  for (unsigned int child = 0; child < cell->n_children(); ++child)
343  {
344  // set an alias to the
345  // prolongation matrix for
346  // this child
347  const FullMatrix<double> &prolongation =
348  mg_dof.get_fe().get_prolongation_matrix(
349  child, cell->refinement_case());
350 
351  cell->child(child)->get_mg_dof_indices(dof_indices_child);
352 
353  // now tag the entries in the
354  // matrix which will be used
355  // for this pair of parent/child
356  for (unsigned int i = 0; i < dofs_per_cell; ++i)
357  for (unsigned int j = 0; j < dofs_per_cell; ++j)
358  if (prolongation(i, j) != 0)
359  {
360  const unsigned int icomp =
361  fe.system_to_block_index(i).first;
362  const unsigned int jcomp =
363  fe.system_to_block_index(j).first;
364  if ((icomp == jcomp) && selected[icomp])
365  prolongation_sparsities[level]->add(
366  dof_indices_child[i], dof_indices_parent[j]);
367  }
368  }
369  }
370  prolongation_sparsities[level]->compress();
371 
372  prolongation_matrices[level]->reinit(*prolongation_sparsities[level]);
373  // now actually build the matrices
374  for (typename DoFHandler<dim, spacedim>::cell_iterator cell =
375  mg_dof.begin(level);
376  cell != mg_dof.end(level);
377  ++cell)
378  if (cell->has_children())
379  {
380  cell->get_mg_dof_indices(dof_indices_parent);
381 
382  Assert(cell->n_children() ==
385  for (unsigned int child = 0; child < cell->n_children(); ++child)
386  {
387  // set an alias to the
388  // prolongation matrix for
389  // this child
390  const FullMatrix<double> &prolongation =
391  mg_dof.get_fe().get_prolongation_matrix(
392  child, cell->refinement_case());
393 
394  cell->child(child)->get_mg_dof_indices(dof_indices_child);
395 
396  // now set the entries in the
397  // matrix
398  for (unsigned int i = 0; i < dofs_per_cell; ++i)
399  for (unsigned int j = 0; j < dofs_per_cell; ++j)
400  if (prolongation(i, j) != 0)
401  {
402  const unsigned int icomp =
403  fe.system_to_block_index(i).first;
404  const unsigned int jcomp =
405  fe.system_to_block_index(j).first;
406  if ((icomp == jcomp) && selected[icomp])
407  prolongation_matrices[level]->set(
408  dof_indices_child[i],
409  dof_indices_parent[j],
410  prolongation(i, j));
411  }
412  }
413  }
414  }
415  // impose boundary conditions
416  // but only in the column of
417  // the prolongation matrix
418  if (mg_constrained_dofs != nullptr &&
420  {
421  std::vector<types::global_dof_index> constrain_indices;
422  std::vector<std::vector<bool>> constraints_per_block(n_blocks);
423  for (int level = n_levels - 2; level >= 0; --level)
424  {
426  0)
427  continue;
428 
429  // need to delete all the columns in the
430  // matrix that are on the boundary. to achieve
431  // this, create an array as long as there are
432  // matrix columns, and find which columns we
433  // need to filter away.
434  constrain_indices.resize(0);
435  constrain_indices.resize(prolongation_matrices[level]->n(), 0);
439  for (; dof != endd; ++dof)
440  constrain_indices[*dof] = 1;
441 
442  unsigned int index = 0;
443  for (unsigned int block = 0; block < n_blocks; ++block)
444  {
445  const types::global_dof_index n_dofs =
446  prolongation_matrices[level]->block(block, block).m();
447  constraints_per_block[block].resize(0);
448  constraints_per_block[block].resize(n_dofs, false);
449  for (types::global_dof_index i = 0; i < n_dofs; ++i, ++index)
450  constraints_per_block[block][i] =
451  (constrain_indices[index] == 1);
452 
453  for (types::global_dof_index i = 0; i < n_dofs; ++i)
454  {
456  start_row = prolongation_matrices[level]
457  ->block(block, block)
458  .begin(i),
459  end_row =
460  prolongation_matrices[level]->block(block, block).end(i);
461  for (; start_row != end_row; ++start_row)
462  {
463  if (constraints_per_block[block][start_row->column()])
464  start_row->value() = 0.;
465  }
466  }
467  }
468  }
469  }
470 }
471 
472 
473 
474 template <typename number>
475 template <int dim, int spacedim>
476 void
478  const DoFHandler<dim, spacedim> &dof,
479  const DoFHandler<dim, spacedim> &mg_dof,
480  unsigned int select)
481 {
482  const FiniteElement<dim> &fe = mg_dof.get_fe();
483  unsigned int n_blocks = mg_dof.get_fe().n_blocks();
484 
485  selected_block = select;
486  selected.resize(n_blocks, false);
487  selected[select] = true;
488 
490 
491  std::vector<types::global_dof_index> temp_copy_indices;
492  std::vector<types::global_dof_index> global_dof_indices(fe.dofs_per_cell);
493  std::vector<types::global_dof_index> level_dof_indices(fe.dofs_per_cell);
494 
495  for (int level = dof.get_triangulation().n_levels() - 1; level >= 0; --level)
496  {
498  mg_dof.begin_active(level);
499  const typename DoFHandler<dim, spacedim>::active_cell_iterator level_end =
500  mg_dof.end_active(level);
501 
502  temp_copy_indices.resize(0);
503  temp_copy_indices.resize(sizes[level][selected_block],
505 
506  // Compute coarse level right hand side
507  // by restricting from fine level.
508  for (; level_cell != level_end; ++level_cell)
509  {
510  // get the dof numbers of
511  // this cell for the global
512  // and the level-wise
513  // numbering
514  level_cell->get_dof_indices(global_dof_indices);
515  level_cell->get_mg_dof_indices(level_dof_indices);
516 
517  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
518  {
519  const unsigned int block = fe.system_to_block_index(i).first;
520  if (selected[block])
521  {
522  if (mg_constrained_dofs != nullptr)
523  {
525  level, level_dof_indices[i]))
526  temp_copy_indices[level_dof_indices[i] -
527  mg_block_start[level][block]] =
528  global_dof_indices[i] - block_start[block];
529  }
530  else
531  temp_copy_indices[level_dof_indices[i] -
532  mg_block_start[level][block]] =
533  global_dof_indices[i] - block_start[block];
534  }
535  }
536  }
537 
538  // now all the active dofs got a valid entry,
539  // the other ones have an invalid entry. Count
540  // the invalid entries and then resize the
541  // copy_indices object. Then, insert the pairs
542  // of global index and level index into
543  // copy_indices.
544  const types::global_dof_index n_active_dofs =
545  std::count_if(temp_copy_indices.begin(),
546  temp_copy_indices.end(),
547  [](const types::global_dof_index index) {
548  return index != numbers::invalid_dof_index;
549  });
550  copy_indices[selected_block][level].resize(n_active_dofs);
551  types::global_dof_index counter = 0;
552  for (types::global_dof_index i = 0; i < temp_copy_indices.size(); ++i)
553  if (temp_copy_indices[i] != numbers::invalid_dof_index)
554  copy_indices[selected_block][level][counter++] =
555  std::pair<types::global_dof_index, unsigned int>(
556  temp_copy_indices[i], i);
557  Assert(counter == n_active_dofs, ExcInternalError());
558  }
559 }
560 
561 
562 
563 template <typename number>
564 template <int dim, int spacedim>
565 void
567  const DoFHandler<dim, spacedim> &mg_dof,
568  const std::vector<bool> & sel)
569 {
570  const FiniteElement<dim> &fe = mg_dof.get_fe();
571  unsigned int n_blocks = mg_dof.get_fe().n_blocks();
572 
573  if (sel.size() != 0)
574  {
575  Assert(sel.size() == n_blocks,
576  ExcDimensionMismatch(sel.size(), n_blocks));
577  selected = sel;
578  }
579  if (selected.size() == 0)
580  selected = std::vector<bool>(n_blocks, true);
581 
583 
584  std::vector<std::vector<types::global_dof_index>> temp_copy_indices(n_blocks);
585  std::vector<types::global_dof_index> global_dof_indices(fe.dofs_per_cell);
586  std::vector<types::global_dof_index> level_dof_indices(fe.dofs_per_cell);
587  for (int level = dof.get_triangulation().n_levels() - 1; level >= 0; --level)
588  {
590  mg_dof.begin_active(level);
591  const typename DoFHandler<dim, spacedim>::active_cell_iterator level_end =
592  mg_dof.end_active(level);
593 
594  for (unsigned int block = 0; block < n_blocks; ++block)
595  if (selected[block])
596  {
597  temp_copy_indices[block].resize(0);
598  temp_copy_indices[block].resize(sizes[level][block],
600  }
601 
602  // Compute coarse level right hand side
603  // by restricting from fine level.
604  for (; level_cell != level_end; ++level_cell)
605  {
606  // get the dof numbers of
607  // this cell for the global
608  // and the level-wise
609  // numbering
610  level_cell->get_dof_indices(global_dof_indices);
611  level_cell->get_mg_dof_indices(level_dof_indices);
612 
613  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
614  {
615  const unsigned int block = fe.system_to_block_index(i).first;
616  if (selected[block])
617  temp_copy_indices[block][level_dof_indices[i] -
618  mg_block_start[level][block]] =
619  global_dof_indices[i] - block_start[block];
620  }
621  }
622 
623  for (unsigned int block = 0; block < n_blocks; ++block)
624  if (selected[block])
625  {
626  const types::global_dof_index n_active_dofs =
627  std::count_if(temp_copy_indices[block].begin(),
628  temp_copy_indices[block].end(),
629  [](const types::global_dof_index index) {
630  return index != numbers::invalid_dof_index;
631  });
632  copy_indices[block][level].resize(n_active_dofs);
633  types::global_dof_index counter = 0;
634  for (types::global_dof_index i = 0;
635  i < temp_copy_indices[block].size();
636  ++i)
637  if (temp_copy_indices[block][i] != numbers::invalid_dof_index)
638  copy_indices[block][level][counter++] =
639  std::pair<types::global_dof_index, unsigned int>(
640  temp_copy_indices[block][i], i);
641  Assert(counter == n_active_dofs, ExcInternalError());
642  }
643  }
644 }
645 
646 
647 
648 // explicit instantiations
649 #include "mg_transfer_block.inst"
650 
651 
652 DEAL_II_NAMESPACE_CLOSE
void build_matrices(const DoFHandler< dim, spacedim > &dof, const DoFHandler< dim, spacedim > &mg_dof)
const Triangulation< dim, spacedim > & get_triangulation() const
static const unsigned int invalid_unsigned_int
Definition: types.h:187
ElementIterator end() const
Definition: index_set.h:1546
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:930
cell_iterator end() const
Definition: dof_handler.cc:959
void count_dofs_per_block(const DoFHandlerType &dof_handler, std::vector< std::vector< types::global_dof_index >> &dofs_per_block, std::vector< unsigned int > target_block={})
Definition: mg_tools.cc:1170
std::vector< types::global_dof_index > block_start
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
SmartPointer< const MGConstrainedDoFs, MGTransferBlockBase > mg_constrained_dofs
unsigned int n_blocks() const
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:943
const IndexSet & get_boundary_indices(const unsigned int level) const
std::vector< std::vector< types::global_dof_index > > sizes
std::vector< std::shared_ptr< BlockSparseMatrix< double > > > prolongation_matrices
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
active_cell_iterator end_active(const unsigned int level) const
Definition: dof_handler.cc:979
const unsigned int dofs_per_cell
Definition: fe_base.h:282
std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > copy_indices
std::vector< unsigned int > mg_block
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
unsigned int global_dof_index
Definition: types.h:89
bool at_refinement_edge(const unsigned int level, const types::global_dof_index index) const
void build_matrices(const DoFHandler< dim, spacedim > &dof, const DoFHandler< dim, spacedim > &mg_dof, unsigned int selected)
std::vector< std::vector< types::global_dof_index > > mg_block_start
void count_dofs_per_block(const DoFHandlerType &dof, std::vector< types::global_dof_index > &dofs_per_block, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
Definition: dof_tools.cc:1894
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
Definition: fe.h:3258
void copy_to_mg(const DoFHandler< dim, spacedim > &mg_dof, MGLevelObject< Vector< number >> &dst, const Vector< number2 > &src) const
void build_matrices(const DoFHandler< dim, spacedim > &dof, const DoFHandler< dim, spacedim > &mg_dof, const std::vector< bool > &selected)
static ::ExceptionBase & ExcNotImplemented()
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:356
void copy_to_mg(const DoFHandler< dim, spacedim > &mg_dof, MGLevelObject< BlockVector< number >> &dst, const BlockVector< number2 > &src) const
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:328
bool have_boundary_indices() const
const types::global_dof_index invalid_dof_index
Definition: types.h:202
ElementIterator begin() const
Definition: index_set.h:1483
BlockType & block(const unsigned int i)
size_type n_elements() const
Definition: index_set.h:1799
static ::ExceptionBase & ExcInternalError()