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