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