Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mg_transfer_component.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/function.h>
18 #include <deal.II/base/logstream.h>
19 
20 #include <deal.II/dofs/dof_accessor.h>
21 #include <deal.II/dofs/dof_tools.h>
22 
23 #include <deal.II/fe/fe.h>
24 
25 #include <deal.II/grid/tria.h>
26 #include <deal.II/grid/tria_iterator.h>
27 
28 #include <deal.II/lac/block_indices.h>
29 #include <deal.II/lac/block_sparse_matrix.h>
30 #include <deal.II/lac/block_vector.h>
31 #include <deal.II/lac/sparse_matrix.h>
32 #include <deal.II/lac/vector.h>
33 
34 #include <deal.II/multigrid/mg_tools.h>
35 #include <deal.II/multigrid/mg_transfer_component.h>
36 #include <deal.II/multigrid/mg_transfer_component.templates.h>
37 
38 #include <algorithm>
39 #include <iostream>
40 #include <numeric>
41 
42 DEAL_II_NAMESPACE_OPEN
43 
44 
45 namespace
46 {
81  template <int dim, typename number, int spacedim>
82  void
83  reinit_vector_by_components(
84  const ::DoFHandler<dim, spacedim> & mg_dof,
86  const std::vector<bool> & sel,
87  const std::vector<unsigned int> & target_comp,
88  std::vector<std::vector<types::global_dof_index>> &ndofs)
89  {
90  std::vector<bool> selected = sel;
91  std::vector<unsigned int> target_component = target_comp;
92  const unsigned int ncomp = mg_dof.get_fe(0).n_components();
93 
94  // If the selected and
95  // target_component have size 0,
96  // they must be replaced by default
97  // values.
98  //
99  // Since we already made copies
100  // directly after this function was
101  // called, we use the arguments
102  // directly.
103  if (target_component.size() == 0)
104  {
105  target_component.resize(ncomp);
106  for (unsigned int i = 0; i < ncomp; ++i)
107  target_component[i] = i;
108  }
109 
110  // If selected is an empty vector,
111  // all components are selected.
112  if (selected.size() == 0)
113  {
114  selected.resize(target_component.size());
115  std::fill_n(selected.begin(), ncomp, false);
116  for (const unsigned int component : target_component)
117  selected[component] = true;
118  }
119 
120  Assert(selected.size() == target_component.size(),
121  ExcDimensionMismatch(selected.size(), target_component.size()));
122 
123  // Compute the number of blocks needed
124  const unsigned int n_selected =
125  std::accumulate(selected.begin(), selected.end(), 0u);
126 
127  if (ndofs.size() == 0)
128  {
129  std::vector<std::vector<types::global_dof_index>> new_dofs(
130  mg_dof.get_triangulation().n_levels(),
131  std::vector<types::global_dof_index>(target_component.size()));
132  std::swap(ndofs, new_dofs);
133  MGTools::count_dofs_per_block(mg_dof, ndofs, target_component);
134  }
135 
136  for (unsigned int level = v.min_level(); level <= v.max_level(); ++level)
137  {
138  v[level].reinit(n_selected, 0);
139  unsigned int k = 0;
140  for (unsigned int i = 0;
141  i < selected.size() && (k < v[level].n_blocks());
142  ++i)
143  {
144  if (selected[i])
145  {
146  v[level].block(k++).reinit(ndofs[level][i]);
147  }
148  v[level].collect_sizes();
149  }
150  }
151  }
152 
153 
180  template <int dim, typename number, int spacedim>
181  void
182  reinit_vector_by_components(
183  const ::DoFHandler<dim, spacedim> & mg_dof,
185  const ComponentMask & component_mask,
186  const std::vector<unsigned int> & target_component,
187  std::vector<std::vector<types::global_dof_index>> &ndofs)
188  {
189  Assert(component_mask.represents_n_components(target_component.size()),
190  ExcMessage("The component mask does not have the correct size."));
191 
192  unsigned int selected_block = 0;
193  for (unsigned int i = 0; i < target_component.size(); ++i)
194  if (component_mask[i])
195  selected_block = target_component[i];
196 
197  if (ndofs.size() == 0)
198  {
199  std::vector<std::vector<types::global_dof_index>> new_dofs(
200  mg_dof.get_triangulation().n_levels(),
201  std::vector<types::global_dof_index>(target_component.size()));
202  std::swap(ndofs, new_dofs);
203  MGTools::count_dofs_per_block(mg_dof, ndofs, target_component);
204  }
205 
206  for (unsigned int level = v.min_level(); level <= v.max_level(); ++level)
207  {
208  v[level].reinit(ndofs[level][selected_block]);
209  }
210  }
211 } // namespace
212 
213 
214 template <typename number>
215 template <int dim, class InVector, int spacedim>
216 void
218  const DoFHandler<dim, spacedim> &mg_dof_handler,
220  const InVector & src) const
221 {
222  dst = 0;
223 
224  Assert(sizes.size() == mg_dof_handler.get_triangulation().n_levels(),
225  ExcMatricesNotBuilt());
226 
227  reinit_vector_by_components(
228  mg_dof_handler, dst, mg_component_mask, mg_target_component, sizes);
229 
230  // traverse the grid top-down
231  // (i.e. starting with the most
232  // refined grid). this way, we can
233  // always get that part of one
234  // level of the output vector which
235  // corresponds to a region which is
236  // more refined, by restriction of
237  // the respective vector on the
238  // next finer level, which we then
239  // already have built.
240 
241  for (unsigned int level = mg_dof_handler.get_triangulation().n_levels();
242  level != 0;)
243  {
244  --level;
245 
246  using IT = std::vector<
247  std::pair<types::global_dof_index, unsigned int>>::const_iterator;
248  for (IT i = copy_to_and_from_indices[level].begin();
249  i != copy_to_and_from_indices[level].end();
250  ++i)
251  dst[level](i->second) = src(i->first);
252  }
253 }
254 
255 
256 template <int dim, int spacedim>
257 void
259  const DoFHandler<dim, spacedim> &mg_dof)
260 {
261  // Fill target component with
262  // standard values (identity) if it
263  // is empty
264  if (target_component.size() == 0)
265  {
266  target_component.resize(mg_dof.get_fe(0).n_components());
267  for (unsigned int i = 0; i < target_component.size(); ++i)
268  target_component[i] = i;
269  }
270  else
271  {
272  // otherwise, check it for consistency
273  Assert(target_component.size() == mg_dof.get_fe(0).n_components(),
275  mg_dof.get_fe(0).n_components()));
276 
277  for (unsigned int i = 0; i < target_component.size(); ++i)
278  {
279  Assert(i < target_component.size(),
280  ExcIndexRange(i, 0, target_component.size()));
281  }
282  }
283  // Do the same for the multilevel
284  // components. These may be
285  // different.
286  if (mg_target_component.size() == 0)
287  {
288  mg_target_component.resize(mg_dof.get_fe(0).n_components());
289  for (unsigned int i = 0; i < mg_target_component.size(); ++i)
291  }
292  else
293  {
294  Assert(mg_target_component.size() == mg_dof.get_fe(0).n_components(),
296  mg_dof.get_fe(0).n_components()));
297 
298  for (unsigned int i = 0; i < mg_target_component.size(); ++i)
299  {
300  Assert(i < mg_target_component.size(),
301  ExcIndexRange(i, 0, mg_target_component.size()));
302  }
303  }
304 
305  const FiniteElement<dim> &fe = mg_dof.get_fe();
306 
307  // Effective number of components
308  // is the maximum entry in
309  // mg_target_component. This
310  // assumes that the values in that
311  // vector don't have holes.
312  const unsigned int n_components =
313  *std::max_element(mg_target_component.begin(), mg_target_component.end()) +
314  1;
315  const unsigned int dofs_per_cell = fe.dofs_per_cell;
316  const unsigned int n_levels = mg_dof.get_triangulation().n_levels();
317 
319  ExcMessage("Component mask has wrong size."));
320 
321  // Compute the lengths of all blocks
322  sizes.resize(n_levels);
323  for (unsigned int l = 0; l < n_levels; ++l)
324  sizes[l].resize(n_components);
325 
327 
328  // Fill some index vectors
329  // for later use.
331  // Compute start indices from sizes
332  for (auto &level_components : mg_component_start)
333  {
335  for (types::global_dof_index &level_component_start : level_components)
336  {
337  const types::global_dof_index t = level_component_start;
338  level_component_start = k;
339  k += t;
340  }
341  }
342 
343  component_start.resize(
344  *std::max_element(target_component.begin(), target_component.end()) + 1);
346 
348  for (types::global_dof_index &first_index : component_start)
349  {
350  const types::global_dof_index t = first_index;
351  first_index = k;
352  k += t;
353  }
354 
355  // Build index vectors for
356  // copy_to_mg and
357  // copy_from_mg. These vectors must
358  // be prebuilt, since the
359  // get_dof_indices functions are
360  // too slow
361 
362  copy_to_and_from_indices.resize(n_levels);
363 
364  // Building the prolongation matrices starts here!
365 
366  // reset the size of the array of
367  // matrices. call resize(0) first,
368  // in order to delete all elements
369  // and clear their memory. then
370  // repopulate these arrays
371  //
372  // note that on resize(0), the
373  // shared_ptr class takes care of
374  // deleting the object it points to
375  // by itself
376  prolongation_matrices.resize(0);
377  prolongation_sparsities.resize(0);
378  prolongation_matrices.reserve(n_levels - 1);
379  prolongation_sparsities.reserve(n_levels - 1);
380 
381  for (unsigned int i = 0; i < n_levels - 1; ++i)
382  {
383  prolongation_sparsities.emplace_back(new BlockSparsityPattern);
385  }
386 
387  // two fields which will store the
388  // indices of the multigrid dofs
389  // for a cell and one of its children
390  std::vector<types::global_dof_index> dof_indices_parent(dofs_per_cell);
391  std::vector<types::global_dof_index> dof_indices_child(dofs_per_cell);
392 
393  // for each level: first build the
394  // sparsity pattern of the matrices
395  // and then build the matrices
396  // themselves. note that we only
397  // need to take care of cells on
398  // the coarser level which have
399  // children
400  for (unsigned int level = 0; level < n_levels - 1; ++level)
401  {
402  // reset the dimension of the
403  // structure. note that for
404  // the number of entries per
405  // row, the number of parent
406  // dofs coupling to a child dof
407  // is necessary. this, is the
408  // number of degrees of freedom
409  // per cell
410  prolongation_sparsities[level]->reinit(n_components, n_components);
411  for (unsigned int i = 0; i < n_components; ++i)
412  for (unsigned int j = 0; j < n_components; ++j)
413  if (i == j)
414  prolongation_sparsities[level]->block(i, j).reinit(
415  sizes[level + 1][i], sizes[level][j], dofs_per_cell + 1);
416  else
417  prolongation_sparsities[level]->block(i, j).reinit(
418  sizes[level + 1][i], sizes[level][j], 0);
419 
420  prolongation_sparsities[level]->collect_sizes();
421 
422  for (typename DoFHandler<dim, spacedim>::cell_iterator cell =
423  mg_dof.begin(level);
424  cell != mg_dof.end(level);
425  ++cell)
426  if (cell->has_children())
427  {
428  cell->get_mg_dof_indices(dof_indices_parent);
429 
430  for (unsigned int child = 0; child < cell->n_children(); ++child)
431  {
432  // set an alias to the
433  // prolongation matrix for
434  // this child
435  const FullMatrix<double> &prolongation =
436  mg_dof.get_fe().get_prolongation_matrix(
437  child, cell->refinement_case());
438 
439  cell->child(child)->get_mg_dof_indices(dof_indices_child);
440 
441  // now tag the entries in the
442  // matrix which will be used
443  // for this pair of parent/child
444  for (unsigned int i = 0; i < dofs_per_cell; ++i)
445  for (unsigned int j = 0; j < dofs_per_cell; ++j)
446  if (prolongation(i, j) != 0)
447  {
448  const unsigned int icomp =
449  fe.system_to_component_index(i).first;
450  const unsigned int jcomp =
451  fe.system_to_component_index(j).first;
452  if ((icomp == jcomp) && mg_component_mask[icomp])
453  prolongation_sparsities[level]->add(
454  dof_indices_child[i], dof_indices_parent[j]);
455  }
456  }
457  }
458  prolongation_sparsities[level]->compress();
459 
460  prolongation_matrices[level]->reinit(*prolongation_sparsities[level]);
461  // now actually build the matrices
462  for (typename DoFHandler<dim, spacedim>::cell_iterator cell =
463  mg_dof.begin(level);
464  cell != mg_dof.end(level);
465  ++cell)
466  if (cell->has_children())
467  {
468  cell->get_mg_dof_indices(dof_indices_parent);
469 
470  for (unsigned int child = 0; child < cell->n_children(); ++child)
471  {
472  // set an alias to the
473  // prolongation matrix for
474  // this child
475  const FullMatrix<double> &prolongation =
476  mg_dof.get_fe().get_prolongation_matrix(
477  child, cell->refinement_case());
478 
479  cell->child(child)->get_mg_dof_indices(dof_indices_child);
480 
481  // now set the entries in the
482  // matrix
483  for (unsigned int i = 0; i < dofs_per_cell; ++i)
484  for (unsigned int j = 0; j < dofs_per_cell; ++j)
485  if (prolongation(i, j) != 0)
486  {
487  const unsigned int icomp =
488  fe.system_to_component_index(i).first;
489  const unsigned int jcomp =
490  fe.system_to_component_index(j).first;
491  if ((icomp == jcomp) && mg_component_mask[icomp])
492  prolongation_matrices[level]->set(
493  dof_indices_child[i],
494  dof_indices_parent[j],
495  prolongation(i, j));
496  }
497  }
498  }
499  }
500  // impose boundary conditions
501  // but only in the column of
502  // the prolongation matrix
503  // TODO: this way is not very efficient
504 
505  if (boundary_indices.size() != 0)
506  {
507  std::vector<std::vector<types::global_dof_index>> dofs_per_component(
508  mg_dof.get_triangulation().n_levels(),
509  std::vector<types::global_dof_index>(n_components));
510 
512  dofs_per_component,
514  for (unsigned int level = 0; level < n_levels - 1; ++level)
515  {
516  if (boundary_indices[level].size() == 0)
517  continue;
518 
519  for (unsigned int iblock = 0; iblock < n_components; ++iblock)
520  for (unsigned int jblock = 0; jblock < n_components; ++jblock)
521  if (iblock == jblock)
522  {
523  const types::global_dof_index n_dofs =
524  prolongation_matrices[level]->block(iblock, jblock).m();
525  for (types::global_dof_index i = 0; i < n_dofs; ++i)
526  {
528  anfang = prolongation_matrices[level]
529  ->block(iblock, jblock)
530  .begin(i),
531  ende = prolongation_matrices[level]
532  ->block(iblock, jblock)
533  .end(i);
534  for (; anfang != ende; ++anfang)
535  {
536  const types::global_dof_index column_number =
537  anfang->column();
538 
539  // convert global indices into local ones
540  const BlockIndices block_indices_coarse(
541  dofs_per_component[level]);
542  const types::global_dof_index global_j =
543  block_indices_coarse.local_to_global(iblock,
544  column_number);
545 
546  std::set<types::global_dof_index>::const_iterator
547  found_dof = boundary_indices[level].find(global_j);
548 
549  const bool is_boundary_index =
550  (found_dof != boundary_indices[level].end());
551 
552  if (is_boundary_index)
553  {
554  prolongation_matrices[level]
555  ->block(iblock, jblock)
556  .set(i, column_number, 0);
557  }
558  }
559  }
560  }
561  }
562  }
563 }
564 
565 
566 template <typename number>
567 template <int dim, int spacedim>
568 void
570  const DoFHandler<dim, spacedim> & dof,
571  const DoFHandler<dim, spacedim> & mg_dof,
572  unsigned int select,
573  unsigned int mg_select,
574  const std::vector<unsigned int> & t_component,
575  const std::vector<unsigned int> & mg_t_component,
576  const std::vector<std::set<types::global_dof_index>> &bdry_indices)
577 {
578  const FiniteElement<dim> &fe = mg_dof.get_fe();
579  unsigned int ncomp = mg_dof.get_fe(0).n_components();
580 
581  target_component = t_component;
582  mg_target_component = mg_t_component;
583  boundary_indices = bdry_indices;
584 
585  selected_component = select;
586  mg_selected_component = mg_select;
587 
588  {
589  std::vector<bool> tmp(ncomp, false);
590  for (unsigned int c = 0; c < ncomp; ++c)
591  if (t_component[c] == selected_component)
592  tmp[c] = true;
593  component_mask = ComponentMask(tmp);
594  }
595 
596  {
597  std::vector<bool> tmp(ncomp, false);
598  for (unsigned int c = 0; c < ncomp; ++c)
599  if (mg_t_component[c] == mg_selected_component)
600  tmp[c] = true;
601  mg_component_mask = ComponentMask(tmp);
602  }
603 
604  // If components are renumbered,
605  // find the first original
606  // component corresponding to the
607  // target component.
608  for (unsigned int i = 0; i < target_component.size(); ++i)
609  {
610  if (target_component[i] == select)
611  {
612  selected_component = i;
613  break;
614  }
615  }
616 
617  for (unsigned int i = 0; i < mg_target_component.size(); ++i)
618  {
619  if (mg_target_component[i] == mg_select)
620  {
621  mg_selected_component = i;
622  break;
623  }
624  }
625 
627 
628  interface_dofs.resize(mg_dof.get_triangulation().n_levels());
629  for (unsigned int l = 0; l < mg_dof.get_triangulation().n_levels(); ++l)
630  {
631  interface_dofs[l].clear();
632  interface_dofs[l].set_size(mg_dof.n_dofs(l));
633  }
634  MGTools::extract_inner_interface_dofs(mg_dof, interface_dofs);
635 
636  // use a temporary vector to create the
637  // relation between global and level dofs
638  std::vector<types::global_dof_index> temp_copy_indices;
639  std::vector<types::global_dof_index> global_dof_indices(fe.dofs_per_cell);
640  std::vector<types::global_dof_index> level_dof_indices(fe.dofs_per_cell);
641  for (int level = dof.get_triangulation().n_levels() - 1; level >= 0; --level)
642  {
643  copy_to_and_from_indices[level].clear();
645  mg_dof.begin_active(level);
646  const typename DoFHandler<dim, spacedim>::active_cell_iterator level_end =
647  mg_dof.end_active(level);
648 
649  temp_copy_indices.resize(0);
650  temp_copy_indices.resize(mg_dof.n_dofs(level),
652 
653  // Compute coarse level right hand side
654  // by restricting from fine level.
655  for (; level_cell != level_end; ++level_cell)
656  {
657  // get the dof numbers of
658  // this cell for the global
659  // and the level-wise
660  // numbering
661  level_cell->get_dof_indices(global_dof_indices);
662  level_cell->get_mg_dof_indices(level_dof_indices);
663 
664  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
665  {
666  const unsigned int component =
667  fe.system_to_component_index(i).first;
668  if (component_mask[component] &&
669  !interface_dofs[level].is_element(level_dof_indices[i]))
670  {
671  const types::global_dof_index level_start =
672  mg_component_start[level][mg_target_component[component]];
673  const types::global_dof_index global_start =
674  component_start[target_component[component]];
675  temp_copy_indices[level_dof_indices[i] - level_start] =
676  global_dof_indices[i] - global_start;
677  }
678  }
679  }
680 
681  // write indices from vector into the map from
682  // global to level dofs
683  const types::global_dof_index n_active_dofs =
684  std::count_if(temp_copy_indices.begin(),
685  temp_copy_indices.end(),
686  std::bind(std::not_equal_to<types::global_dof_index>(),
687  std::placeholders::_1,
689  copy_to_and_from_indices[level].resize(n_active_dofs);
690  types::global_dof_index counter = 0;
691  for (types::global_dof_index i = 0; i < temp_copy_indices.size(); ++i)
692  if (temp_copy_indices[i] != numbers::invalid_dof_index)
693  copy_to_and_from_indices[level][counter++] =
694  std::pair<types::global_dof_index, unsigned int>(
695  temp_copy_indices[i], i);
696  Assert(counter == n_active_dofs, ExcInternalError());
697  }
698 }
699 
700 
701 
702 // explicit instantiations
703 #include "mg_transfer_component.inst"
704 
705 
706 DEAL_II_NAMESPACE_CLOSE
const Triangulation< dim, spacedim > & get_triangulation() const
std::vector< std::shared_ptr< BlockSparseMatrix< double > > > prolongation_matrices
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:930
std::vector< std::set< types::global_dof_index > > boundary_indices
cell_iterator end() const
Definition: dof_handler.cc:959
std::vector< unsigned int > target_component
bool represents_n_components(const unsigned int n) const
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
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:943
static ::ExceptionBase & ExcMessage(std::string arg1)
void build_matrices(const DoFHandler< dim, spacedim > &dof, const DoFHandler< dim, spacedim > &mg_dof, unsigned int selected, unsigned int mg_selected, const std::vector< unsigned int > &target_component=std::vector< unsigned int >(), const std::vector< unsigned int > &mg_target_component=std::vector< unsigned int >(), const std::vector< std::set< types::global_dof_index >> &boundary_indices=std::vector< std::set< types::global_dof_index >>())
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< types::global_dof_index > component_start
active_cell_iterator end_active(const unsigned int level) const
Definition: dof_handler.cc:979
std::vector< std::vector< std::pair< types::global_dof_index, unsigned int > > > copy_to_and_from_indices
types::global_dof_index n_dofs() const
size_type local_to_global(const unsigned int block, const size_type index) const
const unsigned int dofs_per_cell
Definition: fe_base.h:282
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
unsigned int global_dof_index
Definition: types.h:89
void do_copy_to_mg(const DoFHandler< dim, spacedim > &mg_dof, MGLevelObject< Vector< number >> &dst, const InVector &src) const
std::vector< unsigned int > mg_target_component
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3110
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:1953
unsigned int n_components() const
void extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition: mg_tools.cc:1536
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:352
std::vector< std::vector< types::global_dof_index > > mg_component_start
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:324
const types::global_dof_index invalid_dof_index
Definition: types.h:188
void build_matrices(const DoFHandler< dim, spacedim > &dof, const DoFHandler< dim, spacedim > &mg_dof)
std::vector< std::vector< types::global_dof_index > > sizes
static ::ExceptionBase & ExcInternalError()