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