Reference documentation for deal.II version Git 2276b94619 2020-02-28 10:11:55 +0100
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mg_tools.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2018 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 #include <deal.II/base/logstream.h>
17 #include <deal.II/base/mg_level_object.h>
18 #include <deal.II/base/thread_management.h>
19 
20 #include <deal.II/dofs/dof_accessor.h>
21 #include <deal.II/dofs/dof_handler.h>
22 #include <deal.II/dofs/dof_tools.h>
23 
24 #include <deal.II/fe/component_mask.h>
25 #include <deal.II/fe/fe.h>
26 
27 #include <deal.II/grid/tria.h>
28 #include <deal.II/grid/tria_iterator.h>
29 
30 #include <deal.II/lac/block_sparse_matrix.h>
31 #include <deal.II/lac/block_sparsity_pattern.h>
32 #include <deal.II/lac/block_vector.h>
33 #include <deal.II/lac/dynamic_sparsity_pattern.h>
34 #include <deal.II/lac/sparse_matrix.h>
35 #include <deal.II/lac/sparsity_pattern.h>
36 
37 #include <deal.II/multigrid/mg_base.h>
38 #include <deal.II/multigrid/mg_constrained_dofs.h>
39 #include <deal.II/multigrid/mg_tools.h>
40 
41 #include <deal.II/numerics/matrix_tools.h>
42 
43 #include <algorithm>
44 #include <numeric>
45 #include <vector>
46 
47 DEAL_II_NAMESPACE_OPEN
48 
49 
50 namespace MGTools
51 {
52  // specializations for 1D
53  template <>
54  void
56  const unsigned int,
57  std::vector<unsigned int> &,
58  const DoFTools::Coupling)
59  {
60  Assert(false, ExcNotImplemented());
61  }
62 
63 
64 
65  template <>
66  void
68  const unsigned int,
69  std::vector<unsigned int> &,
72  {
73  Assert(false, ExcNotImplemented());
74  }
75 
76 
77 
78  template <>
79  void
81  const unsigned int,
82  std::vector<unsigned int> &,
83  const DoFTools::Coupling)
84  {
85  Assert(false, ExcNotImplemented());
86  }
87 
88 
89  template <>
90  void
92  const unsigned int,
93  std::vector<unsigned int> &,
96  {
97  Assert(false, ExcNotImplemented());
98  }
99 
100 
101 
102  // Template for 2D and 3D. For 1D see specialization above
103  template <int dim, int spacedim>
104  void
106  const unsigned int level,
107  std::vector<unsigned int> & row_lengths,
108  const DoFTools::Coupling flux_coupling)
109  {
110  Assert(row_lengths.size() == dofs.n_dofs(),
111  ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs()));
112 
113  // Function starts here by
114  // resetting the counters.
115  std::fill(row_lengths.begin(), row_lengths.end(), 0);
116  // We need the user flags, so we
117  // save them for later restoration
118  std::vector<bool> old_flags;
119  // We need a non-constant
120  // triangulation for the user
121  // flags. Since we restore them in
122  // the end, this cast is safe.
123  Triangulation<dim, spacedim> &user_flags_triangulation =
124  const_cast<Triangulation<dim, spacedim> &>(dofs.get_triangulation());
125  user_flags_triangulation.save_user_flags(old_flags);
126  user_flags_triangulation.clear_user_flags();
127 
128  std::vector<types::global_dof_index> cell_indices;
129  std::vector<types::global_dof_index> neighbor_indices;
130 
131  // We loop over cells and go from
132  // cells to lower dimensional
133  // objects. This is the only way to
134  // cope with the fact, that an
135  // unknown number of cells may
136  // share an object of dimension
137  // smaller than dim-1.
138  for (const auto &cell : dofs.cell_iterators_on_level(level))
139  {
140  const FiniteElement<dim> &fe = cell->get_fe();
141  cell_indices.resize(fe.dofs_per_cell);
142  cell->get_mg_dof_indices(cell_indices);
143  unsigned int i = 0;
144  // First, dofs on
145  // vertices. We assume that
146  // each vertex dof couples
147  // with all dofs on
148  // adjacent grid cells.
149 
150  // Adding all dofs of the cells
151  // will add dofs of the faces
152  // of the cell adjacent to the
153  // vertex twice. Therefore, we
154  // subtract these here and add
155  // them in a loop over the
156  // faces below.
157 
158  // in 1d, faces and vertices
159  // are identical. Nevertheless,
160  // this will only work if
161  // dofs_per_face is zero and
162  // dofs_per_vertex is
163  // arbitrary, not the other way
164  // round.
165  // TODO: This assumes that even in hp context, the dofs per face
166  // coincide!
167  unsigned int increment = fe.dofs_per_cell - dim * fe.dofs_per_face;
168  while (i < fe.first_line_index)
169  row_lengths[cell_indices[i++]] += increment;
170  // From now on, if an object is
171  // a cell, its dofs only couple
172  // inside the cell. Since the
173  // faces are handled below, we
174  // have to subtract ALL faces
175  // in this case.
176 
177  // In all other cases we
178  // subtract adjacent faces to be
179  // added in the loop below.
180  increment = (dim > 1) ?
181  fe.dofs_per_cell - (dim - 1) * fe.dofs_per_face :
182  fe.dofs_per_cell -
184  while (i < fe.first_quad_index)
185  row_lengths[cell_indices[i++]] += increment;
186 
187  // Now quads in 2D and 3D
188  increment = (dim > 2) ?
189  fe.dofs_per_cell - (dim - 2) * fe.dofs_per_face :
190  fe.dofs_per_cell -
192  while (i < fe.first_hex_index)
193  row_lengths[cell_indices[i++]] += increment;
194  // Finally, cells in 3D
195  increment = fe.dofs_per_cell -
197  while (i < fe.dofs_per_cell)
198  row_lengths[cell_indices[i++]] += increment;
199 
200  // At this point, we have
201  // counted all dofs
202  // contributiong from cells
203  // coupled topologically to the
204  // adjacent cells, but we
205  // subtracted some faces.
206 
207  // Now, let's go by the faces
208  // and add the missing
209  // contribution as well as the
210  // flux contributions.
211  for (const unsigned int iface : GeometryInfo<dim>::face_indices())
212  {
213  bool level_boundary = cell->at_boundary(iface);
215  if (!level_boundary)
216  {
217  neighbor = cell->neighbor(iface);
218  if (static_cast<unsigned int>(neighbor->level()) != level)
219  level_boundary = true;
220  }
221 
222  if (level_boundary)
223  {
224  for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
225  ++local_dof)
226  row_lengths[cell_indices[local_dof]] += fe.dofs_per_face;
227  continue;
228  }
229 
230  const FiniteElement<dim> &nfe = neighbor->get_fe();
232  cell->face(iface);
233 
234  // Flux couplings are
235  // computed from both sides
236  // for simplicity.
237 
238  // The dofs on the common face
239  // will be handled below,
240  // therefore, we subtract them
241  // here.
242  if (flux_coupling != DoFTools::none)
243  {
244  const unsigned int dof_increment =
245  nfe.dofs_per_cell - nfe.dofs_per_face;
246  for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
247  ++local_dof)
248  row_lengths[cell_indices[local_dof]] += dof_increment;
249  }
250 
251  // Do this only once per
252  // face.
253  if (face->user_flag_set())
254  continue;
255  face->set_user_flag();
256  // At this point, we assume
257  // that each cell added its
258  // dofs minus the face to
259  // the couplings of the
260  // face dofs. Since we
261  // subtracted two faces, we
262  // have to re-add one.
263 
264  // If one side of the face
265  // is refined, all the fine
266  // face dofs couple with
267  // the coarse one.
268  neighbor_indices.resize(nfe.dofs_per_cell);
269  neighbor->get_mg_dof_indices(neighbor_indices);
270  for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
271  ++local_dof)
272  row_lengths[cell_indices[local_dof]] += nfe.dofs_per_face;
273  for (unsigned int local_dof = 0; local_dof < nfe.dofs_per_cell;
274  ++local_dof)
275  row_lengths[neighbor_indices[local_dof]] += fe.dofs_per_face;
276  }
277  }
278  user_flags_triangulation.load_user_flags(old_flags);
279  }
280 
281 
282  // This is the template for 2D and 3D. See version for 1D above
283  template <int dim, int spacedim>
284  void
286  const unsigned int level,
287  std::vector<unsigned int> & row_lengths,
288  const Table<2, DoFTools::Coupling> &couplings,
289  const Table<2, DoFTools::Coupling> &flux_couplings)
290  {
291  Assert(row_lengths.size() == dofs.n_dofs(),
292  ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs()));
293 
294  // Function starts here by
295  // resetting the counters.
296  std::fill(row_lengths.begin(), row_lengths.end(), 0);
297  // We need the user flags, so we
298  // save them for later restoration
299  std::vector<bool> old_flags;
300  // We need a non-constant
301  // triangulation for the user
302  // flags. Since we restore them in
303  // the end, this cast is safe.
304  Triangulation<dim, spacedim> &user_flags_triangulation =
305  const_cast<Triangulation<dim, spacedim> &>(dofs.get_triangulation());
306  user_flags_triangulation.save_user_flags(old_flags);
307  user_flags_triangulation.clear_user_flags();
308 
309  std::vector<types::global_dof_index> cell_indices;
310  std::vector<types::global_dof_index> neighbor_indices;
311 
312  // We have to translate the
313  // couplings from components to
314  // blocks, so this works for
315  // nonprimitive elements as well.
316  std::vector<Table<2, DoFTools::Coupling>> couple_cell;
317  std::vector<Table<2, DoFTools::Coupling>> couple_face;
318  DoFTools::convert_couplings_to_blocks(dofs, couplings, couple_cell);
319  DoFTools::convert_couplings_to_blocks(dofs, flux_couplings, couple_face);
320 
321  // We loop over cells and go from
322  // cells to lower dimensional
323  // objects. This is the only way to
324  // cope with the fact, that an
325  // unknown number of cells may
326  // share an object of dimension
327  // smaller than dim-1.
328  for (const auto &cell : dofs.cell_iterators_on_level(level))
329  {
330  const FiniteElement<dim> &fe = cell->get_fe();
331  const unsigned int fe_index = cell->active_fe_index();
332 
333  Assert(couplings.n_rows() == fe.n_components(),
334  ExcDimensionMismatch(couplings.n_rows(), fe.n_components()));
335  Assert(couplings.n_cols() == fe.n_components(),
336  ExcDimensionMismatch(couplings.n_cols(), fe.n_components()));
337  Assert(flux_couplings.n_rows() == fe.n_components(),
338  ExcDimensionMismatch(flux_couplings.n_rows(),
339  fe.n_components()));
340  Assert(flux_couplings.n_cols() == fe.n_components(),
341  ExcDimensionMismatch(flux_couplings.n_cols(),
342  fe.n_components()));
343 
344  cell_indices.resize(fe.dofs_per_cell);
345  cell->get_mg_dof_indices(cell_indices);
346  unsigned int i = 0;
347  // First, dofs on
348  // vertices. We assume that
349  // each vertex dof couples
350  // with all dofs on
351  // adjacent grid cells.
352 
353  // Adding all dofs of the cells
354  // will add dofs of the faces
355  // of the cell adjacent to the
356  // vertex twice. Therefore, we
357  // subtract these here and add
358  // them in a loop over the
359  // faces below.
360 
361  // in 1d, faces and vertices
362  // are identical. Nevertheless,
363  // this will only work if
364  // dofs_per_face is zero and
365  // dofs_per_vertex is
366  // arbitrary, not the other way
367  // round.
368  unsigned int increment;
369  while (i < fe.first_line_index)
370  {
371  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
372  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
373  ++mult)
374  if (couple_cell[fe_index](fe.system_to_block_index(i).first,
375  fe.first_block_of_base(base) +
376  mult) != DoFTools::none)
377  {
378  increment = fe.base_element(base).dofs_per_cell -
379  dim * fe.base_element(base).dofs_per_face;
380  row_lengths[cell_indices[i]] += increment;
381  }
382  ++i;
383  }
384  // From now on, if an object is
385  // a cell, its dofs only couple
386  // inside the cell. Since the
387  // faces are handled below, we
388  // have to subtract ALL faces
389  // in this case.
390 
391  // In all other cases we
392  // subtract adjacent faces to be
393  // added in the loop below.
394  while (i < fe.first_quad_index)
395  {
396  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
397  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
398  ++mult)
399  if (couple_cell[fe_index](fe.system_to_block_index(i).first,
400  fe.first_block_of_base(base) +
401  mult) != DoFTools::none)
402  {
403  increment =
404  fe.base_element(base).dofs_per_cell -
405  ((dim > 1) ? (dim - 1) :
407  fe.base_element(base).dofs_per_face;
408  row_lengths[cell_indices[i]] += increment;
409  }
410  ++i;
411  }
412 
413  // Now quads in 2D and 3D
414  while (i < fe.first_hex_index)
415  {
416  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
417  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
418  ++mult)
419  if (couple_cell[fe_index](fe.system_to_block_index(i).first,
420  fe.first_block_of_base(base) +
421  mult) != DoFTools::none)
422  {
423  increment =
424  fe.base_element(base).dofs_per_cell -
425  ((dim > 2) ? (dim - 2) :
427  fe.base_element(base).dofs_per_face;
428  row_lengths[cell_indices[i]] += increment;
429  }
430  ++i;
431  }
432 
433  // Finally, cells in 3D
434  while (i < fe.dofs_per_cell)
435  {
436  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
437  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
438  ++mult)
439  if (couple_cell[fe_index](fe.system_to_block_index(i).first,
440  fe.first_block_of_base(base) +
441  mult) != DoFTools::none)
442  {
443  increment = fe.base_element(base).dofs_per_cell -
445  fe.base_element(base).dofs_per_face;
446  row_lengths[cell_indices[i]] += increment;
447  }
448  ++i;
449  }
450 
451  // At this point, we have
452  // counted all dofs
453  // contributiong from cells
454  // coupled topologically to the
455  // adjacent cells, but we
456  // subtracted some faces.
457 
458  // Now, let's go by the faces
459  // and add the missing
460  // contribution as well as the
461  // flux contributions.
462  for (const unsigned int iface : GeometryInfo<dim>::face_indices())
463  {
464  bool level_boundary = cell->at_boundary(iface);
466  if (!level_boundary)
467  {
468  neighbor = cell->neighbor(iface);
469  if (static_cast<unsigned int>(neighbor->level()) != level)
470  level_boundary = true;
471  }
472 
473  if (level_boundary)
474  {
475  for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
476  ++local_dof)
477  row_lengths[cell_indices[local_dof]] += fe.dofs_per_face;
478  continue;
479  }
480 
481  const FiniteElement<dim> &nfe = neighbor->get_fe();
483  cell->face(iface);
484 
485  // Flux couplings are
486  // computed from both sides
487  // for simplicity.
488 
489  // The dofs on the common face
490  // will be handled below,
491  // therefore, we subtract them
492  // here.
493  for (unsigned int base = 0; base < nfe.n_base_elements(); ++base)
494  for (unsigned int mult = 0; mult < nfe.element_multiplicity(base);
495  ++mult)
496  for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
497  ++local_dof)
498  if (couple_face[fe_index](
499  fe.system_to_block_index(local_dof).first,
500  nfe.first_block_of_base(base) + mult) != DoFTools::none)
501  {
502  const unsigned int dof_increment =
503  nfe.base_element(base).dofs_per_cell -
504  nfe.base_element(base).dofs_per_face;
505  row_lengths[cell_indices[local_dof]] += dof_increment;
506  }
507 
508  // Do this only once per
509  // face and not on the
510  // hanging faces.
511  if (face->user_flag_set())
512  continue;
513  face->set_user_flag();
514  // At this point, we assume
515  // that each cell added its
516  // dofs minus the face to
517  // the couplings of the
518  // face dofs. Since we
519  // subtracted two faces, we
520  // have to re-add one.
521 
522  // If one side of the face
523  // is refined, all the fine
524  // face dofs couple with
525  // the coarse one.
526 
527  // Wolfgang, do they couple
528  // with each other by
529  // constraints?
530 
531  // This will not work with
532  // different couplings on
533  // different cells.
534  neighbor_indices.resize(nfe.dofs_per_cell);
535  neighbor->get_mg_dof_indices(neighbor_indices);
536  for (unsigned int base = 0; base < nfe.n_base_elements(); ++base)
537  for (unsigned int mult = 0; mult < nfe.element_multiplicity(base);
538  ++mult)
539  for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
540  ++local_dof)
541  if (couple_cell[fe_index](
542  fe.system_to_component_index(local_dof).first,
543  nfe.first_block_of_base(base) + mult) != DoFTools::none)
544  row_lengths[cell_indices[local_dof]] +=
545  nfe.base_element(base).dofs_per_face;
546  for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
547  for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
548  ++mult)
549  for (unsigned int local_dof = 0; local_dof < nfe.dofs_per_cell;
550  ++local_dof)
551  if (couple_cell[fe_index](
552  nfe.system_to_component_index(local_dof).first,
553  fe.first_block_of_base(base) + mult) != DoFTools::none)
554  row_lengths[neighbor_indices[local_dof]] +=
555  fe.base_element(base).dofs_per_face;
556  }
557  }
558  user_flags_triangulation.load_user_flags(old_flags);
559  }
560 
561 
562 
563  template <typename DoFHandlerType, typename SparsityPatternType>
564  void
565  make_sparsity_pattern(const DoFHandlerType &dof,
566  SparsityPatternType & sparsity,
567  const unsigned int level)
568  {
569  const types::global_dof_index n_dofs = dof.n_dofs(level);
570  (void)n_dofs;
571 
572  Assert(sparsity.n_rows() == n_dofs,
573  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
574  Assert(sparsity.n_cols() == n_dofs,
575  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
576 
577  const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
578  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
579  typename DoFHandlerType::cell_iterator cell = dof.begin(level),
580  endc = dof.end(level);
581  for (; cell != endc; ++cell)
582  if (dof.get_triangulation().locally_owned_subdomain() ==
584  cell->level_subdomain_id() ==
585  dof.get_triangulation().locally_owned_subdomain())
586  {
587  cell->get_mg_dof_indices(dofs_on_this_cell);
588  // make sparsity pattern for this cell
589  for (unsigned int i = 0; i < dofs_per_cell; ++i)
590  for (unsigned int j = 0; j < dofs_per_cell; ++j)
591  sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);
592  }
593  }
594 
595 
596 
597  template <int dim, typename SparsityPatternType, int spacedim>
598  void
600  SparsityPatternType & sparsity,
601  const unsigned int level)
602  {
603  const types::global_dof_index n_dofs = dof.n_dofs(level);
604  (void)n_dofs;
605 
606  Assert(sparsity.n_rows() == n_dofs,
607  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
608  Assert(sparsity.n_cols() == n_dofs,
609  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
610 
611  const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
612  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
613  std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
614  typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(level),
615  endc = dof.end(level);
616  for (; cell != endc; ++cell)
617  {
618  if (!cell->is_locally_owned_on_level())
619  continue;
620 
621  cell->get_mg_dof_indices(dofs_on_this_cell);
622  // make sparsity pattern for this cell
623  for (unsigned int i = 0; i < dofs_per_cell; ++i)
624  for (unsigned int j = 0; j < dofs_per_cell; ++j)
625  sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);
626 
627  // Loop over all interior neighbors
628  for (const unsigned int face : GeometryInfo<dim>::face_indices())
629  {
630  bool use_face = false;
631  if ((!cell->at_boundary(face)) &&
632  (static_cast<unsigned int>(cell->neighbor_level(face)) ==
633  level))
634  use_face = true;
635  else if (cell->has_periodic_neighbor(face) &&
636  (static_cast<unsigned int>(
637  cell->periodic_neighbor_level(face)) == level))
638  use_face = true;
639 
640  if (use_face)
641  {
642  typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
643  cell->neighbor_or_periodic_neighbor(face);
644  neighbor->get_mg_dof_indices(dofs_on_other_cell);
645  // only add one direction The other is taken care of by
646  // neighbor (except when the neighbor is not owned by the same
647  // processor)
648  for (unsigned int i = 0; i < dofs_per_cell; ++i)
649  {
650  for (unsigned int j = 0; j < dofs_per_cell; ++j)
651  {
652  sparsity.add(dofs_on_this_cell[i],
653  dofs_on_other_cell[j]);
654  }
655  }
656  if (neighbor->is_locally_owned_on_level() == false)
657  for (unsigned int i = 0; i < dofs_per_cell; ++i)
658  for (unsigned int j = 0; j < dofs_per_cell; ++j)
659  {
660  sparsity.add(dofs_on_other_cell[i],
661  dofs_on_other_cell[j]);
662  sparsity.add(dofs_on_other_cell[i],
663  dofs_on_this_cell[j]);
664  }
665  }
666  }
667  }
668  }
669 
670 
671 
672  template <int dim, typename SparsityPatternType, int spacedim>
673  void
675  SparsityPatternType & sparsity,
676  const unsigned int level)
677  {
678  Assert((level >= 1) && (level < dof.get_triangulation().n_global_levels()),
679  ExcIndexRange(level, 1, dof.get_triangulation().n_global_levels()));
680 
681  const types::global_dof_index fine_dofs = dof.n_dofs(level);
682  const types::global_dof_index coarse_dofs = dof.n_dofs(level - 1);
683  (void)fine_dofs;
684  (void)coarse_dofs;
685 
686  // Matrix maps from fine level to coarse level
687 
688  Assert(sparsity.n_rows() == coarse_dofs,
689  ExcDimensionMismatch(sparsity.n_rows(), coarse_dofs));
690  Assert(sparsity.n_cols() == fine_dofs,
691  ExcDimensionMismatch(sparsity.n_cols(), fine_dofs));
692 
693  const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
694  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
695  std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
696  typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(level),
697  endc = dof.end(level);
698  for (; cell != endc; ++cell)
699  {
700  if (!cell->is_locally_owned_on_level())
701  continue;
702 
703  cell->get_mg_dof_indices(dofs_on_this_cell);
704  // Loop over all interior neighbors
705  for (const unsigned int face : GeometryInfo<dim>::face_indices())
706  {
707  // Neighbor is coarser
708  bool use_face = false;
709  if ((!cell->at_boundary(face)) &&
710  (static_cast<unsigned int>(cell->neighbor_level(face)) !=
711  level))
712  use_face = true;
713  else if (cell->has_periodic_neighbor(face) &&
714  (static_cast<unsigned int>(
715  cell->periodic_neighbor_level(face)) != level))
716  use_face = true;
717 
718  if (use_face)
719  {
720  typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
721  cell->neighbor_or_periodic_neighbor(face);
722  neighbor->get_mg_dof_indices(dofs_on_other_cell);
723 
724  for (unsigned int i = 0; i < dofs_per_cell; ++i)
725  {
726  for (unsigned int j = 0; j < dofs_per_cell; ++j)
727  {
728  sparsity.add(dofs_on_other_cell[i],
729  dofs_on_this_cell[j]);
730  sparsity.add(dofs_on_other_cell[j],
731  dofs_on_this_cell[i]);
732  }
733  }
734  }
735  }
736  }
737  }
738 
739 
740 
741  template <int dim, typename SparsityPatternType, int spacedim>
742  void
744  SparsityPatternType & sparsity,
745  const unsigned int level,
746  const Table<2, DoFTools::Coupling> &int_mask,
747  const Table<2, DoFTools::Coupling> &flux_mask)
748  {
749  const FiniteElement<dim> & fe = dof.get_fe();
750  const types::global_dof_index n_dofs = dof.n_dofs(level);
751  const unsigned int n_comp = fe.n_components();
752  (void)n_dofs;
753  (void)n_comp;
754 
755  Assert(sparsity.n_rows() == n_dofs,
756  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
757  Assert(sparsity.n_cols() == n_dofs,
758  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
759  Assert(int_mask.n_rows() == n_comp,
760  ExcDimensionMismatch(int_mask.n_rows(), n_comp));
761  Assert(int_mask.n_cols() == n_comp,
762  ExcDimensionMismatch(int_mask.n_cols(), n_comp));
763  Assert(flux_mask.n_rows() == n_comp,
764  ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
765  Assert(flux_mask.n_cols() == n_comp,
766  ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
767 
768  const unsigned int total_dofs = fe.dofs_per_cell;
769  std::vector<types::global_dof_index> dofs_on_this_cell(total_dofs);
770  std::vector<types::global_dof_index> dofs_on_other_cell(total_dofs);
771  Table<2, bool> support_on_face(total_dofs,
773 
774  typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(level),
775  endc = dof.end(level);
776 
778  int_dof_mask =
780  flux_dof_mask =
782 
783  for (unsigned int i = 0; i < total_dofs; ++i)
784  for (auto f : GeometryInfo<dim>::face_indices())
785  support_on_face(i, f) = fe.has_support_on_face(i, f);
786 
787  // Clear user flags because we will
788  // need them. But first we save
789  // them and make sure that we
790  // restore them later such that at
791  // the end of this function the
792  // Triangulation will be in the
793  // same state as it was at the
794  // beginning of this function.
795  std::vector<bool> user_flags;
796  dof.get_triangulation().save_user_flags(user_flags);
798  .clear_user_flags();
799 
800  for (; cell != endc; ++cell)
801  {
802  if (!cell->is_locally_owned_on_level())
803  continue;
804 
805  cell->get_mg_dof_indices(dofs_on_this_cell);
806  // make sparsity pattern for this cell
807  for (unsigned int i = 0; i < total_dofs; ++i)
808  for (unsigned int j = 0; j < total_dofs; ++j)
809  if (int_dof_mask[i][j] != DoFTools::none)
810  sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);
811 
812  // Loop over all interior neighbors
813  for (const unsigned int face : GeometryInfo<dim>::face_indices())
814  {
815  typename DoFHandler<dim, spacedim>::face_iterator cell_face =
816  cell->face(face);
817  if (cell_face->user_flag_set())
818  continue;
819 
820  if (cell->at_boundary(face) && !cell->has_periodic_neighbor(face))
821  {
822  for (unsigned int i = 0; i < total_dofs; ++i)
823  {
824  const bool i_non_zero_i = support_on_face(i, face);
825  for (unsigned int j = 0; j < total_dofs; ++j)
826  {
827  const bool j_non_zero_i = support_on_face(j, face);
828 
829  if (flux_dof_mask(i, j) == DoFTools::always)
830  sparsity.add(dofs_on_this_cell[i],
831  dofs_on_this_cell[j]);
832  if (flux_dof_mask(i, j) == DoFTools::nonzero &&
833  i_non_zero_i && j_non_zero_i)
834  sparsity.add(dofs_on_this_cell[i],
835  dofs_on_this_cell[j]);
836  }
837  }
838  }
839  else
840  {
841  typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
842  cell->neighbor_or_periodic_neighbor(face);
843 
844  if (neighbor->level() < cell->level())
845  continue;
846 
847  unsigned int neighbor_face =
848  cell->has_periodic_neighbor(face) ?
849  cell->periodic_neighbor_of_periodic_neighbor(face) :
850  cell->neighbor_of_neighbor(face);
851 
852  neighbor->get_mg_dof_indices(dofs_on_other_cell);
853  for (unsigned int i = 0; i < total_dofs; ++i)
854  {
855  const bool i_non_zero_i = support_on_face(i, face);
856  const bool i_non_zero_e = support_on_face(i, neighbor_face);
857  for (unsigned int j = 0; j < total_dofs; ++j)
858  {
859  const bool j_non_zero_i = support_on_face(j, face);
860  const bool j_non_zero_e =
861  support_on_face(j, neighbor_face);
862  if (flux_dof_mask(i, j) == DoFTools::always)
863  {
864  sparsity.add(dofs_on_this_cell[i],
865  dofs_on_other_cell[j]);
866  sparsity.add(dofs_on_other_cell[i],
867  dofs_on_this_cell[j]);
868  sparsity.add(dofs_on_this_cell[i],
869  dofs_on_this_cell[j]);
870  sparsity.add(dofs_on_other_cell[i],
871  dofs_on_other_cell[j]);
872  }
873  if (flux_dof_mask(i, j) == DoFTools::nonzero)
874  {
875  if (i_non_zero_i && j_non_zero_e)
876  sparsity.add(dofs_on_this_cell[i],
877  dofs_on_other_cell[j]);
878  if (i_non_zero_e && j_non_zero_i)
879  sparsity.add(dofs_on_other_cell[i],
880  dofs_on_this_cell[j]);
881  if (i_non_zero_i && j_non_zero_i)
882  sparsity.add(dofs_on_this_cell[i],
883  dofs_on_this_cell[j]);
884  if (i_non_zero_e && j_non_zero_e)
885  sparsity.add(dofs_on_other_cell[i],
886  dofs_on_other_cell[j]);
887  }
888 
889  if (flux_dof_mask(j, i) == DoFTools::always)
890  {
891  sparsity.add(dofs_on_this_cell[j],
892  dofs_on_other_cell[i]);
893  sparsity.add(dofs_on_other_cell[j],
894  dofs_on_this_cell[i]);
895  sparsity.add(dofs_on_this_cell[j],
896  dofs_on_this_cell[i]);
897  sparsity.add(dofs_on_other_cell[j],
898  dofs_on_other_cell[i]);
899  }
900  if (flux_dof_mask(j, i) == DoFTools::nonzero)
901  {
902  if (j_non_zero_i && i_non_zero_e)
903  sparsity.add(dofs_on_this_cell[j],
904  dofs_on_other_cell[i]);
905  if (j_non_zero_e && i_non_zero_i)
906  sparsity.add(dofs_on_other_cell[j],
907  dofs_on_this_cell[i]);
908  if (j_non_zero_i && i_non_zero_i)
909  sparsity.add(dofs_on_this_cell[j],
910  dofs_on_this_cell[i]);
911  if (j_non_zero_e && i_non_zero_e)
912  sparsity.add(dofs_on_other_cell[j],
913  dofs_on_other_cell[i]);
914  }
915  }
916  }
917  neighbor->face(neighbor_face)->set_user_flag();
918  }
919  }
920  }
921 
922  // finally restore the user flags
924  .load_user_flags(user_flags);
925  }
926 
927 
928 
929  template <int dim, typename SparsityPatternType, int spacedim>
930  void
932  SparsityPatternType & sparsity,
933  const unsigned int level,
934  const Table<2, DoFTools::Coupling> &flux_mask)
935  {
936  const FiniteElement<dim> &fe = dof.get_fe();
937  const unsigned int n_comp = fe.n_components();
938  (void)n_comp;
939 
940  Assert((level >= 1) && (level < dof.get_triangulation().n_global_levels()),
941  ExcIndexRange(level, 1, dof.get_triangulation().n_global_levels()));
942 
943  const types::global_dof_index fine_dofs = dof.n_dofs(level);
944  const types::global_dof_index coarse_dofs = dof.n_dofs(level - 1);
945  (void)fine_dofs;
946  (void)coarse_dofs;
947 
948  // Matrix maps from fine level to coarse level
949 
950  Assert(sparsity.n_rows() == coarse_dofs,
951  ExcDimensionMismatch(sparsity.n_rows(), coarse_dofs));
952  Assert(sparsity.n_cols() == fine_dofs,
953  ExcDimensionMismatch(sparsity.n_cols(), fine_dofs));
954  Assert(flux_mask.n_rows() == n_comp,
955  ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
956  Assert(flux_mask.n_cols() == n_comp,
957  ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
958 
959  const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
960  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
961  std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
962  Table<2, bool> support_on_face(dofs_per_cell,
964 
965  typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(level),
966  endc = dof.end(level);
967 
968  const Table<2, DoFTools::Coupling> flux_dof_mask =
970 
971  for (unsigned int i = 0; i < dofs_per_cell; ++i)
972  for (auto f : GeometryInfo<dim>::face_indices())
973  support_on_face(i, f) = fe.has_support_on_face(i, f);
974 
975  for (; cell != endc; ++cell)
976  {
977  if (!cell->is_locally_owned_on_level())
978  continue;
979 
980  cell->get_mg_dof_indices(dofs_on_this_cell);
981  // Loop over all interior neighbors
982  for (const unsigned int face : GeometryInfo<dim>::face_indices())
983  {
984  // Neighbor is coarser
985  bool use_face = false;
986  if ((!cell->at_boundary(face)) &&
987  (static_cast<unsigned int>(cell->neighbor_level(face)) !=
988  level))
989  use_face = true;
990  else if (cell->has_periodic_neighbor(face) &&
991  (static_cast<unsigned int>(
992  cell->periodic_neighbor_level(face)) != level))
993  use_face = true;
994 
995  if (use_face)
996  {
997  typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
998  cell->neighbor_or_periodic_neighbor(face);
999  neighbor->get_mg_dof_indices(dofs_on_other_cell);
1000 
1001  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1002  {
1003  for (unsigned int j = 0; j < dofs_per_cell; ++j)
1004  {
1005  if (flux_dof_mask(i, j) != DoFTools::none)
1006  {
1007  sparsity.add(dofs_on_other_cell[i],
1008  dofs_on_this_cell[j]);
1009  sparsity.add(dofs_on_other_cell[j],
1010  dofs_on_this_cell[i]);
1011  }
1012  }
1013  }
1014  }
1015  }
1016  }
1017  }
1018 
1019 
1020 
1021  template <typename DoFHandlerType, typename SparsityPatternType>
1022  void
1023  make_interface_sparsity_pattern(const DoFHandlerType & dof,
1024  const MGConstrainedDoFs &mg_constrained_dofs,
1025  SparsityPatternType & sparsity,
1026  const unsigned int level)
1027  {
1028  const types::global_dof_index n_dofs = dof.n_dofs(level);
1029  (void)n_dofs;
1030 
1031  Assert(sparsity.n_rows() == n_dofs,
1032  ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1033  Assert(sparsity.n_cols() == n_dofs,
1034  ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1035 
1036  const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
1037  std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
1038  typename DoFHandlerType::cell_iterator cell = dof.begin(level),
1039  endc = dof.end(level);
1040  for (; cell != endc; ++cell)
1041  if (cell->is_locally_owned_on_level())
1042  {
1043  cell->get_mg_dof_indices(dofs_on_this_cell);
1044  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1045  for (unsigned int j = 0; j < dofs_per_cell; ++j)
1046  if (mg_constrained_dofs.is_interface_matrix_entry(
1047  level, dofs_on_this_cell[i], dofs_on_this_cell[j]))
1048  sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);
1049  }
1050  }
1051 
1052 
1053 
1054  template <int dim, int spacedim>
1055  void
1057  const DoFHandler<dim, spacedim> & dof_handler,
1058  std::vector<std::vector<types::global_dof_index>> &result,
1059  bool only_once,
1060  std::vector<unsigned int> target_component)
1061  {
1062  const FiniteElement<dim> &fe = dof_handler.get_fe();
1063  const unsigned int n_components = fe.n_components();
1064  const unsigned int nlevels =
1065  dof_handler.get_triangulation().n_global_levels();
1066 
1067  Assert(result.size() == nlevels,
1068  ExcDimensionMismatch(result.size(), nlevels));
1069 
1070  if (target_component.size() == 0)
1071  {
1072  target_component.resize(n_components);
1073  for (unsigned int i = 0; i < n_components; ++i)
1074  target_component[i] = i;
1075  }
1076 
1077  Assert(target_component.size() == n_components,
1078  ExcDimensionMismatch(target_component.size(), n_components));
1079 
1080  for (unsigned int l = 0; l < nlevels; ++l)
1081  {
1082  result[l].resize(n_components);
1083  std::fill(result[l].begin(), result[l].end(), 0U);
1084 
1085  // special case for only one
1086  // component. treat this first
1087  // since it does not require any
1088  // computations
1089  if (n_components == 1)
1090  {
1091  result[l][0] = dof_handler.n_dofs(l);
1092  }
1093  else
1094  {
1095  // otherwise determine the number
1096  // of dofs in each component
1097  // separately. do so in parallel
1098  std::vector<std::vector<bool>> dofs_in_component(
1099  n_components, std::vector<bool>(dof_handler.n_dofs(l), false));
1100  std::vector<ComponentMask> component_select(n_components);
1101  Threads::TaskGroup<> tasks;
1102  for (unsigned int i = 0; i < n_components; ++i)
1103  {
1104  void (*fun_ptr)(const unsigned int level,
1105  const DoFHandler<dim, spacedim> &,
1106  const ComponentMask &,
1107  std::vector<bool> &) =
1108  &DoFTools::extract_level_dofs<DoFHandler<dim, spacedim>>;
1109 
1110  std::vector<bool> tmp(n_components, false);
1111  tmp[i] = true;
1112  component_select[i] = ComponentMask(tmp);
1113 
1114  tasks += Threads::new_task(fun_ptr,
1115  l,
1116  dof_handler,
1117  component_select[i],
1118  dofs_in_component[i]);
1119  }
1120  tasks.join_all();
1121 
1122  // next count what we got
1123  unsigned int component = 0;
1124  for (unsigned int b = 0; b < fe.n_base_elements(); ++b)
1125  {
1126  const FiniteElement<dim> &base = fe.base_element(b);
1127  // Dimension of base element
1128  unsigned int d = base.n_components();
1129 
1130  for (unsigned int m = 0; m < fe.element_multiplicity(b); ++m)
1131  {
1132  for (unsigned int dd = 0; dd < d; ++dd)
1133  {
1134  if (base.is_primitive() || (!only_once || dd == 0))
1135  result[l][target_component[component]] +=
1136  std::count(dofs_in_component[component].begin(),
1137  dofs_in_component[component].end(),
1138  true);
1139  ++component;
1140  }
1141  }
1142  }
1143  // finally sanity check
1144  Assert(!dof_handler.get_fe().is_primitive() ||
1145  std::accumulate(result[l].begin(),
1146  result[l].end(),
1148  dof_handler.n_dofs(l),
1149  ExcInternalError());
1150  }
1151  }
1152  }
1153 
1154 
1155 
1156  template <typename DoFHandlerType>
1157  void
1159  const DoFHandlerType & dof_handler,
1160  std::vector<std::vector<types::global_dof_index>> &dofs_per_block,
1161  std::vector<unsigned int> target_block)
1162  {
1163  const FiniteElement<DoFHandlerType::dimension,
1164  DoFHandlerType::space_dimension> &fe =
1165  dof_handler.get_fe();
1166  const unsigned int n_blocks = fe.n_blocks();
1167  const unsigned int n_levels =
1168  dof_handler.get_triangulation().n_global_levels();
1169 
1170  AssertDimension(dofs_per_block.size(), n_levels);
1171 
1172  for (unsigned int l = 0; l < n_levels; ++l)
1173  std::fill(dofs_per_block[l].begin(), dofs_per_block[l].end(), 0U);
1174  // If the empty vector was given as
1175  // default argument, set up this
1176  // vector as identity.
1177  if (target_block.size() == 0)
1178  {
1179  target_block.resize(n_blocks);
1180  for (unsigned int i = 0; i < n_blocks; ++i)
1181  target_block[i] = i;
1182  }
1183  Assert(target_block.size() == n_blocks,
1184  ExcDimensionMismatch(target_block.size(), n_blocks));
1185 
1186  const unsigned int max_block =
1187  *std::max_element(target_block.begin(), target_block.end());
1188  const unsigned int n_target_blocks = max_block + 1;
1189  (void)n_target_blocks;
1190 
1191  for (unsigned int l = 0; l < n_levels; ++l)
1192  AssertDimension(dofs_per_block[l].size(), n_target_blocks);
1193 
1194  // special case for only one
1195  // block. treat this first
1196  // since it does not require any
1197  // computations
1198  if (n_blocks == 1)
1199  {
1200  for (unsigned int l = 0; l < n_levels; ++l)
1201  dofs_per_block[l][0] = dof_handler.n_dofs(l);
1202  return;
1203  }
1204  // otherwise determine the number
1205  // of dofs in each block
1206  // separately. do so in parallel
1207  for (unsigned int l = 0; l < n_levels; ++l)
1208  {
1209  std::vector<std::vector<bool>> dofs_in_block(
1210  n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
1211  std::vector<BlockMask> block_select(n_blocks);
1212  Threads::TaskGroup<> tasks;
1213  for (unsigned int i = 0; i < n_blocks; ++i)
1214  {
1215  void (*fun_ptr)(const unsigned int level,
1216  const DoFHandlerType &,
1217  const BlockMask &,
1218  std::vector<bool> &) =
1219  &DoFTools::extract_level_dofs<DoFHandlerType>;
1220 
1221  std::vector<bool> tmp(n_blocks, false);
1222  tmp[i] = true;
1223  block_select[i] = tmp;
1224 
1225  tasks += Threads::new_task(
1226  fun_ptr, l, dof_handler, block_select[i], dofs_in_block[i]);
1227  }
1228  tasks.join_all();
1229 
1230  // next count what we got
1231  for (unsigned int block = 0; block < fe.n_blocks(); ++block)
1232  dofs_per_block[l][target_block[block]] +=
1233  std::count(dofs_in_block[block].begin(),
1234  dofs_in_block[block].end(),
1235  true);
1236  }
1237  }
1238 
1239 
1240 
1241  template <int dim, int spacedim>
1242  void
1244  const DoFHandler<dim, spacedim> &dof,
1245  const std::map<types::boundary_id, const Function<spacedim> *>
1246  & function_map,
1247  std::vector<std::set<types::global_dof_index>> &boundary_indices,
1248  const ComponentMask & component_mask)
1249  {
1250  Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
1251  ExcDimensionMismatch(boundary_indices.size(),
1252  dof.get_triangulation().n_global_levels()));
1253 
1254  std::set<types::boundary_id> boundary_ids;
1255  for (const auto &boundary_function : function_map)
1256  boundary_ids.insert(boundary_function.first);
1257 
1258  std::vector<IndexSet> boundary_indexset;
1259  make_boundary_list(dof, boundary_ids, boundary_indexset, component_mask);
1260  for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
1261  boundary_indices[i].insert(boundary_indexset[i].begin(),
1262  boundary_indexset[i].end());
1263  }
1264 
1265 
1266  template <int dim, int spacedim>
1267  void
1269  const std::map<types::boundary_id,
1270  const Function<spacedim> *> &function_map,
1271  std::vector<IndexSet> &boundary_indices,
1272  const ComponentMask & component_mask)
1273  {
1274  Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
1275  ExcDimensionMismatch(boundary_indices.size(),
1276  dof.get_triangulation().n_global_levels()));
1277 
1278  std::set<types::boundary_id> boundary_ids;
1279  for (const auto &boundary_function : function_map)
1280  boundary_ids.insert(boundary_function.first);
1281 
1282  make_boundary_list(dof, boundary_ids, boundary_indices, component_mask);
1283  }
1284 
1285 
1286 
1287  template <int dim, int spacedim>
1288  void
1290  const std::set<types::boundary_id> &boundary_ids,
1291  std::vector<IndexSet> & boundary_indices,
1292  const ComponentMask & component_mask)
1293  {
1294  boundary_indices.resize(dof.get_triangulation().n_global_levels());
1295 
1296  // if for whatever reason we were passed an empty set, return immediately
1297  if (boundary_ids.size() == 0)
1298  return;
1299 
1300  for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
1301  if (boundary_indices[i].size() == 0)
1302  boundary_indices[i] = IndexSet(dof.n_dofs(i));
1303 
1304  const unsigned int n_components = DoFTools::n_components(dof);
1305  const bool fe_is_system = (n_components != 1);
1306 
1307  std::vector<types::global_dof_index> local_dofs;
1308  local_dofs.reserve(DoFTools::max_dofs_per_face(dof));
1309  std::fill(local_dofs.begin(), local_dofs.end(), numbers::invalid_dof_index);
1310 
1311  std::vector<std::vector<types::global_dof_index>> dofs_by_level(
1312  dof.get_triangulation().n_levels());
1313 
1314  // First, deal with the simpler case when we have to identify all boundary
1315  // dofs
1316  if (component_mask.n_selected_components(n_components) == n_components)
1317  {
1318  for (const auto &cell : dof.cell_iterators())
1319  {
1320  if (dof.get_triangulation().locally_owned_subdomain() !=
1322  cell->level_subdomain_id() == numbers::artificial_subdomain_id)
1323  continue;
1324  const FiniteElement<dim> &fe = cell->get_fe();
1325  const unsigned int level = cell->level();
1326  local_dofs.resize(fe.dofs_per_face);
1327 
1328  for (unsigned int face_no = 0;
1329  face_no < GeometryInfo<dim>::faces_per_cell;
1330  ++face_no)
1331  if (cell->at_boundary(face_no) == true)
1332  {
1333  const typename DoFHandler<dim, spacedim>::face_iterator face =
1334  cell->face(face_no);
1335  const types::boundary_id bi = face->boundary_id();
1336  // Face is listed in boundary map
1337  if (boundary_ids.find(bi) != boundary_ids.end())
1338  {
1339  face->get_mg_dof_indices(level, local_dofs);
1340  dofs_by_level[level].insert(dofs_by_level[level].end(),
1341  local_dofs.begin(),
1342  local_dofs.end());
1343  }
1344  }
1345  }
1346  }
1347  else
1348  {
1349  Assert(component_mask.n_selected_components(n_components) > 0,
1350  ExcMessage(
1351  "It's probably worthwhile to select at least one component."));
1352 
1353  for (const auto &cell : dof.cell_iterators())
1354  if (dof.get_triangulation().locally_owned_subdomain() ==
1356  cell->level_subdomain_id() != numbers::artificial_subdomain_id)
1357  for (unsigned int face_no = 0;
1358  face_no < GeometryInfo<dim>::faces_per_cell;
1359  ++face_no)
1360  {
1361  if (cell->at_boundary(face_no) == false)
1362  continue;
1363 
1364  const FiniteElement<dim> &fe = cell->get_fe();
1365  const unsigned int level = cell->level();
1366 
1368  cell->face(face_no);
1369  const types::boundary_id boundary_component =
1370  face->boundary_id();
1371  if (boundary_ids.find(boundary_component) != boundary_ids.end())
1372  // we want to constrain this boundary
1373  {
1374  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
1375  ++i)
1376  {
1377  const ComponentMask &nonzero_component_array =
1378  cell->get_fe().get_nonzero_components(i);
1379  // if we want to constrain one of the nonzero
1380  // components, we have to constrain all of them
1381 
1382  bool selected = false;
1383  for (unsigned int c = 0; c < n_components; ++c)
1384  if (nonzero_component_array[c] == true &&
1385  component_mask[c] == true)
1386  {
1387  selected = true;
1388  break;
1389  }
1390  if (selected)
1391  for (unsigned int c = 0; c < n_components; ++c)
1392  Assert(
1393  nonzero_component_array[c] == false ||
1394  component_mask[c] == true,
1395  ExcMessage(
1396  "You are using a non-primitive FiniteElement "
1397  "and try to constrain just some of its components!"));
1398  }
1399 
1400  // get indices, physical location and boundary values of
1401  // dofs on this face
1402  local_dofs.resize(fe.dofs_per_face);
1403  face->get_mg_dof_indices(level, local_dofs);
1404  if (fe_is_system)
1405  {
1406  for (unsigned int i = 0; i < local_dofs.size(); ++i)
1407  {
1408  unsigned int component =
1410  if (fe.is_primitive())
1411  component =
1412  fe.face_system_to_component_index(i).first;
1413  else
1414  {
1415  // Just pick the first of the components
1416  // We already know that either all or none
1417  // of the components are selected
1418  const ComponentMask &nonzero_component_array =
1419  cell->get_fe().get_nonzero_components(i);
1420  for (unsigned int c = 0; c < n_components; ++c)
1421  if (nonzero_component_array[c] == true)
1422  {
1423  component = c;
1424  break;
1425  }
1426  }
1428  ExcInternalError());
1429  if (component_mask[component] == true)
1430  dofs_by_level[level].push_back(local_dofs[i]);
1431  }
1432  }
1433  else
1434  dofs_by_level[level].insert(dofs_by_level[level].end(),
1435  local_dofs.begin(),
1436  local_dofs.end());
1437  }
1438  }
1439  }
1440  for (unsigned int level = 0; level < dof.get_triangulation().n_levels();
1441  ++level)
1442  {
1443  std::sort(dofs_by_level[level].begin(), dofs_by_level[level].end());
1444  boundary_indices[level].add_indices(
1445  dofs_by_level[level].begin(),
1446  std::unique(dofs_by_level[level].begin(),
1447  dofs_by_level[level].end()));
1448  }
1449  }
1450 
1451 
1452  template <int dim, int spacedim>
1453  void
1455  const DoFHandler<dim, spacedim> & mg_dof_handler,
1456  std::vector<std::set<types::global_dof_index>> &non_interface_dofs)
1457  {
1458  Assert(non_interface_dofs.size() ==
1459  mg_dof_handler.get_triangulation().n_global_levels(),
1461  non_interface_dofs.size(),
1462  mg_dof_handler.get_triangulation().n_global_levels()));
1463 
1464  const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();
1465 
1466  const unsigned int dofs_per_cell = fe.dofs_per_cell;
1467  const unsigned int dofs_per_face = fe.dofs_per_face;
1468 
1469  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1470  std::vector<bool> cell_dofs(dofs_per_cell, false);
1471  std::vector<bool> cell_dofs_interface(dofs_per_cell, false);
1472 
1473  typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
1474  endc = mg_dof_handler.end();
1475 
1476 
1477  for (; cell != endc; ++cell)
1478  {
1479  if (mg_dof_handler.get_triangulation().locally_owned_subdomain() !=
1481  cell->level_subdomain_id() !=
1482  mg_dof_handler.get_triangulation().locally_owned_subdomain())
1483  continue;
1484 
1485  std::fill(cell_dofs.begin(), cell_dofs.end(), false);
1486  std::fill(cell_dofs_interface.begin(),
1487  cell_dofs_interface.end(),
1488  false);
1489 
1490  for (unsigned int face_nr = 0;
1491  face_nr < GeometryInfo<dim>::faces_per_cell;
1492  ++face_nr)
1493  {
1494  const typename DoFHandler<dim, spacedim>::face_iterator face =
1495  cell->face(face_nr);
1496  if (!face->at_boundary() || cell->has_periodic_neighbor(face_nr))
1497  {
1498  // interior face
1499  const typename DoFHandler<dim>::cell_iterator neighbor =
1500  cell->neighbor_or_periodic_neighbor(face_nr);
1501 
1502  if ((neighbor->level() < cell->level()))
1503  {
1504  for (unsigned int j = 0; j < dofs_per_face; ++j)
1505  cell_dofs_interface[fe.face_to_cell_index(j, face_nr)] =
1506  true;
1507  }
1508  else
1509  {
1510  for (unsigned int j = 0; j < dofs_per_face; ++j)
1511  cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
1512  }
1513  }
1514  else
1515  {
1516  // boundary face
1517  for (unsigned int j = 0; j < dofs_per_face; ++j)
1518  cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
1519  }
1520  }
1521 
1522  const unsigned int level = cell->level();
1523  cell->get_mg_dof_indices(local_dof_indices);
1524 
1525  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1526  if (cell_dofs[i] && !cell_dofs_interface[i])
1527  non_interface_dofs[level].insert(local_dof_indices[i]);
1528  }
1529  }
1530 
1531 
1532  template <int dim, int spacedim>
1533  void
1535  std::vector<IndexSet> & interface_dofs)
1536  {
1537  Assert(interface_dofs.size() ==
1538  mg_dof_handler.get_triangulation().n_global_levels(),
1540  interface_dofs.size(),
1541  mg_dof_handler.get_triangulation().n_global_levels()));
1542 
1543  std::vector<std::vector<types::global_dof_index>> tmp_interface_dofs(
1544  interface_dofs.size());
1545 
1546  const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();
1547 
1548  const unsigned int dofs_per_cell = fe.dofs_per_cell;
1549  const unsigned int dofs_per_face = fe.dofs_per_face;
1550 
1551  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1552 
1553  std::vector<bool> cell_dofs(dofs_per_cell, false);
1554 
1555  typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
1556  endc = mg_dof_handler.end();
1557 
1558  for (; cell != endc; ++cell)
1559  {
1560  // Do not look at artificial level cells (in a serial computation we
1561  // need to ignore the level_subdomain_id() because it is never set).
1562  if (mg_dof_handler.get_triangulation().locally_owned_subdomain() !=
1564  cell->level_subdomain_id() == numbers::artificial_subdomain_id)
1565  continue;
1566 
1567  bool has_coarser_neighbor = false;
1568 
1569  std::fill(cell_dofs.begin(), cell_dofs.end(), false);
1570 
1571  for (unsigned int face_nr = 0;
1572  face_nr < GeometryInfo<dim>::faces_per_cell;
1573  ++face_nr)
1574  {
1575  const typename DoFHandler<dim, spacedim>::face_iterator face =
1576  cell->face(face_nr);
1577  if (!face->at_boundary() || cell->has_periodic_neighbor(face_nr))
1578  {
1579  // interior face
1580  const typename DoFHandler<dim>::cell_iterator neighbor =
1581  cell->neighbor_or_periodic_neighbor(face_nr);
1582 
1583  // only process cell pairs if one or both of them are owned by
1584  // me (ignore if running in serial)
1585  if (mg_dof_handler.get_triangulation()
1586  .locally_owned_subdomain() !=
1588  neighbor->level_subdomain_id() ==
1590  continue;
1591 
1592  // Do refinement face from the coarse side
1593  if (neighbor->level() < cell->level())
1594  {
1595  for (unsigned int j = 0; j < dofs_per_face; ++j)
1596  cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
1597 
1598  has_coarser_neighbor = true;
1599  }
1600  }
1601  }
1602 
1603  if (has_coarser_neighbor == false)
1604  continue;
1605 
1606  const unsigned int level = cell->level();
1607  cell->get_mg_dof_indices(local_dof_indices);
1608 
1609  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1610  {
1611  if (cell_dofs[i])
1612  tmp_interface_dofs[level].push_back(local_dof_indices[i]);
1613  }
1614  }
1615 
1616  for (unsigned int l = 0;
1617  l < mg_dof_handler.get_triangulation().n_global_levels();
1618  ++l)
1619  {
1620  interface_dofs[l].clear();
1621  std::sort(tmp_interface_dofs[l].begin(), tmp_interface_dofs[l].end());
1622  interface_dofs[l].add_indices(tmp_interface_dofs[l].begin(),
1623  std::unique(tmp_interface_dofs[l].begin(),
1624  tmp_interface_dofs[l].end()));
1625  interface_dofs[l].compress();
1626  }
1627  }
1628 
1629 
1630 
1631  template <int dim, int spacedim>
1632  unsigned int
1634  {
1635  // Find minimum level for an active cell in
1636  // this locally owned subdomain
1637  // Note: with the way active cells are traversed,
1638  // the first locally owned cell we find will have
1639  // the lowest level in the particular subdomain.
1640  unsigned int min_level = tria.n_global_levels();
1642  cell = tria.begin_active(),
1643  endc = tria.end();
1644  for (; cell != endc; ++cell)
1645  if (cell->is_locally_owned())
1646  {
1647  min_level = cell->level();
1648  break;
1649  }
1650 
1651  unsigned int global_min = min_level;
1652  // If necessary, communicate to find minimum
1653  // level for an active cell over all subdomains
1655  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1656  &tria))
1657  global_min = Utilities::MPI::min(min_level, tr->get_communicator());
1658 
1659  AssertIndexRange(global_min, tria.n_global_levels());
1660 
1661  return global_min;
1662  }
1663 
1664 
1665 
1666  template <int dim, int spacedim>
1667  double
1669  {
1670  double workload_imbalance = 1.0;
1671 
1672  // It is only necessary to calculate the imbalance
1673  // on a distributed mesh. The imbalance is always
1674  // 1.0 for the serial case.
1676  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1677  &tria))
1678  {
1679  const unsigned int n_proc =
1680  Utilities::MPI::n_mpi_processes(tr->get_communicator());
1681  const unsigned int n_global_levels = tr->n_global_levels();
1682 
1683  // This value will represent the sum over the multigrid
1684  // levels of the maximum number of cells owned by any
1685  // one processesor on that level.
1686  types::global_dof_index work_estimate = 0;
1687 
1688  // Sum of all cells in the multigrid hierarchy
1689  types::global_dof_index total_cells_in_hierarchy = 0;
1690 
1691  for (int lvl = n_global_levels - 1; lvl >= 0; --lvl)
1692  {
1693  // Number of cells this processor owns on this level
1694  types::global_dof_index n_owned_cells_on_lvl = 0;
1695 
1696  for (const auto &cell : tr->cell_iterators_on_level(lvl))
1697  if (cell->is_locally_owned_on_level())
1698  ++n_owned_cells_on_lvl;
1699 
1700  work_estimate +=
1701  ::Utilities::MPI::max(n_owned_cells_on_lvl,
1702  tr->get_communicator());
1703 
1704  total_cells_in_hierarchy +=
1705  ::Utilities::MPI::sum(n_owned_cells_on_lvl,
1706  tr->get_communicator());
1707  }
1708 
1709  const double ideal_work =
1710  total_cells_in_hierarchy / static_cast<double>(n_proc);
1711  workload_imbalance = work_estimate / ideal_work;
1712  }
1713 
1714  return workload_imbalance;
1715  }
1716 } // namespace MGTools
1717 
1718 
1719 // explicit instantiations
1720 #include "mg_tools.inst"
1721 
1722 DEAL_II_NAMESPACE_CLOSE
const unsigned int first_hex_index
Definition: fe_base.h:258
const Triangulation< dim, spacedim > & get_triangulation() const
void count_dofs_per_component(const DoFHandler< dim, spacedim > &mg_dof, std::vector< std::vector< types::global_dof_index >> &result, const bool only_once=false, std::vector< unsigned int > target_component={})
Definition: mg_tools.cc:1056
static const unsigned int invalid_unsigned_int
Definition: types.h:187
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1571
const types::subdomain_id invalid_subdomain_id
Definition: types.h:279
void clear_user_flags()
Definition: tria.cc:11213
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:931
double workload_imbalance(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1668
Task< RT > new_task(const std::function< RT()> &function)
void load_user_flags(std::istream &in)
Definition: tria.cc:11273
cell_iterator end() const
Definition: dof_handler.cc:960
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:599
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:565
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:1158
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:12052
bool is_primitive() const
Definition: fe.h:3301
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
Table< 2, Coupling > dof_couplings_from_component_couplings(const FiniteElement< dim, spacedim > &fe, const Table< 2, Coupling > &component_couplings)
typename ActiveSelector::face_iterator face_iterator
Definition: dof_handler.h:375
bool is_interface_matrix_entry(const unsigned int level, const types::global_dof_index i, const types::global_dof_index j) const
cell_iterator end() const
Definition: tria.cc:12118
void make_interface_sparsity_pattern(const DoFHandlerType &dof_handler, const MGConstrainedDoFs &mg_constrained_dofs, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:1023
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int max_dofs_per_face(const DoFHandler< dim, spacedim > &dh)
const unsigned int first_quad_index
Definition: fe_base.h:253
void convert_couplings_to_blocks(const hp::DoFHandler< dim, spacedim > &dof_handler, const Table< 2, Coupling > &table_by_component, std::vector< Table< 2, Coupling >> &tables_by_block)
Definition: dof_tools.cc:2481
#define Assert(cond, exc)
Definition: exceptions.h:1411
unsigned int element_multiplicity(const unsigned int index) const
Definition: fe.h:3133
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index) const
Definition: fe.h:3170
types::global_dof_index n_dofs() const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
Definition: fe.cc:1276
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
Definition: fe.cc:537
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:100
const unsigned int dofs_per_cell
Definition: fe_base.h:282
types::global_dof_index first_block_of_base(const unsigned int b) const
Definition: fe.h:3221
unsigned int global_dof_index
Definition: types.h:89
const types::subdomain_id artificial_subdomain_id
Definition: types.h:296
void make_boundary_list(const DoFHandler< dim, spacedim > &mg_dof, const std::map< types::boundary_id, const Function< spacedim > *> &function_map, std::vector< std::set< types::global_dof_index >> &boundary_indices, const ComponentMask &component_mask=ComponentMask())
Definition: mg_tools.cc:1243
void compute_row_length_vector(const DoFHandler< dim, spacedim > &dofs, const unsigned int level, std::vector< unsigned int > &row_lengths, const DoFTools::Coupling flux_couplings=DoFTools::none)
Definition: mg_tools.cc:105
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3110
unsigned int n_components() const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
Definition: fe.h:3253
virtual unsigned int n_global_levels() const
T min(const T &t, const MPI_Comm &mpi_communicator)
void extract_non_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< std::set< types::global_dof_index >> &non_interface_dofs)
Definition: mg_tools.cc:1454
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const unsigned int dofs_per_face
Definition: fe_base.h:275
const unsigned int first_line_index
Definition: fe_base.h:248
void extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition: mg_tools.cc:1534
void save_user_flags(std::ostream &out) const
Definition: tria.cc:11224
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:699
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:357
void extract_level_dofs(const unsigned int level, const DoFHandlerType &dof, const ComponentMask &component_mask, std::vector< bool > &selected_dofs)
Definition: dof_tools.cc:566
void make_flux_sparsity_pattern_edge(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:674
const types::global_dof_index invalid_dof_index
Definition: types.h:202
unsigned int n_base_elements() const
Definition: fe.h:3124
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
unsigned int boundary_id
Definition: types.h:125
IteratorRange< cell_iterator > cell_iterators() const
unsigned int max_level_for_coarse_mesh(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1633
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
Definition: fe.cc:1112
static ::ExceptionBase & ExcInternalError()