Reference documentation for deal.II version Git 4abc4a1666 2020-07-04 19:58:34 +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 - 2020 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.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 <int dim, int spacedim, typename SparsityPatternType>
564  void
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 DoFHandler<dim, spacedim>::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 <int dim, int spacedim, typename SparsityPatternType>
1022  void
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 DoFHandler<dim, spacedim>::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<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 <int dim, int spacedim>
1157  void
1159  const DoFHandler<dim, spacedim> & 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<dim, spacedim> &fe = dof_handler.get_fe();
1164  const unsigned int n_blocks = fe.n_blocks();
1165  const unsigned int n_levels =
1166  dof_handler.get_triangulation().n_global_levels();
1167 
1168  AssertDimension(dofs_per_block.size(), n_levels);
1169 
1170  for (unsigned int l = 0; l < n_levels; ++l)
1171  std::fill(dofs_per_block[l].begin(), dofs_per_block[l].end(), 0U);
1172  // If the empty vector was given as
1173  // default argument, set up this
1174  // vector as identity.
1175  if (target_block.size() == 0)
1176  {
1177  target_block.resize(n_blocks);
1178  for (unsigned int i = 0; i < n_blocks; ++i)
1179  target_block[i] = i;
1180  }
1181  Assert(target_block.size() == n_blocks,
1182  ExcDimensionMismatch(target_block.size(), n_blocks));
1183 
1184  const unsigned int max_block =
1185  *std::max_element(target_block.begin(), target_block.end());
1186  const unsigned int n_target_blocks = max_block + 1;
1187  (void)n_target_blocks;
1188 
1189  for (unsigned int l = 0; l < n_levels; ++l)
1190  AssertDimension(dofs_per_block[l].size(), n_target_blocks);
1191 
1192  // special case for only one
1193  // block. treat this first
1194  // since it does not require any
1195  // computations
1196  if (n_blocks == 1)
1197  {
1198  for (unsigned int l = 0; l < n_levels; ++l)
1199  dofs_per_block[l][0] = dof_handler.n_dofs(l);
1200  return;
1201  }
1202  // otherwise determine the number
1203  // of dofs in each block
1204  // separately. do so in parallel
1205  for (unsigned int l = 0; l < n_levels; ++l)
1206  {
1207  std::vector<std::vector<bool>> dofs_in_block(
1208  n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
1209  std::vector<BlockMask> block_select(n_blocks);
1210  Threads::TaskGroup<> tasks;
1211  for (unsigned int i = 0; i < n_blocks; ++i)
1212  {
1213  void (*fun_ptr)(const unsigned int level,
1214  const DoFHandler<dim, spacedim> &,
1215  const BlockMask &,
1216  std::vector<bool> &) =
1217  &DoFTools::extract_level_dofs<dim, spacedim>;
1218 
1219  std::vector<bool> tmp(n_blocks, false);
1220  tmp[i] = true;
1221  block_select[i] = tmp;
1222 
1223  tasks += Threads::new_task(
1224  fun_ptr, l, dof_handler, block_select[i], dofs_in_block[i]);
1225  }
1226  tasks.join_all();
1227 
1228  // next count what we got
1229  for (unsigned int block = 0; block < fe.n_blocks(); ++block)
1230  dofs_per_block[l][target_block[block]] +=
1231  std::count(dofs_in_block[block].begin(),
1232  dofs_in_block[block].end(),
1233  true);
1234  }
1235  }
1236 
1237 
1238 
1239  template <int dim, int spacedim>
1240  void
1242  const DoFHandler<dim, spacedim> &dof,
1243  const std::map<types::boundary_id, const Function<spacedim> *>
1244  & function_map,
1245  std::vector<std::set<types::global_dof_index>> &boundary_indices,
1246  const ComponentMask & component_mask)
1247  {
1248  Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
1249  ExcDimensionMismatch(boundary_indices.size(),
1250  dof.get_triangulation().n_global_levels()));
1251 
1252  std::set<types::boundary_id> boundary_ids;
1253  for (const auto &boundary_function : function_map)
1254  boundary_ids.insert(boundary_function.first);
1255 
1256  std::vector<IndexSet> boundary_indexset;
1257  make_boundary_list(dof, boundary_ids, boundary_indexset, component_mask);
1258  for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
1259  boundary_indices[i].insert(boundary_indexset[i].begin(),
1260  boundary_indexset[i].end());
1261  }
1262 
1263 
1264  template <int dim, int spacedim>
1265  void
1267  const std::map<types::boundary_id,
1268  const Function<spacedim> *> &function_map,
1269  std::vector<IndexSet> &boundary_indices,
1270  const ComponentMask & component_mask)
1271  {
1272  Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
1273  ExcDimensionMismatch(boundary_indices.size(),
1274  dof.get_triangulation().n_global_levels()));
1275 
1276  std::set<types::boundary_id> boundary_ids;
1277  for (const auto &boundary_function : function_map)
1278  boundary_ids.insert(boundary_function.first);
1279 
1280  make_boundary_list(dof, boundary_ids, boundary_indices, component_mask);
1281  }
1282 
1283 
1284 
1285  template <int dim, int spacedim>
1286  void
1288  const std::set<types::boundary_id> &boundary_ids,
1289  std::vector<IndexSet> & boundary_indices,
1290  const ComponentMask & component_mask)
1291  {
1292  boundary_indices.resize(dof.get_triangulation().n_global_levels());
1293 
1294  // if for whatever reason we were passed an empty set, return immediately
1295  if (boundary_ids.size() == 0)
1296  return;
1297 
1298  for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
1299  if (boundary_indices[i].size() == 0)
1300  boundary_indices[i] = IndexSet(dof.n_dofs(i));
1301 
1302  const unsigned int n_components = dof.get_fe_collection().n_components();
1303  const bool fe_is_system = (n_components != 1);
1304 
1305  std::vector<types::global_dof_index> local_dofs;
1306  local_dofs.reserve(dof.get_fe_collection().max_dofs_per_face());
1307  std::fill(local_dofs.begin(), local_dofs.end(), numbers::invalid_dof_index);
1308 
1309  std::vector<std::vector<types::global_dof_index>> dofs_by_level(
1310  dof.get_triangulation().n_levels());
1311 
1312  // First, deal with the simpler case when we have to identify all boundary
1313  // dofs
1314  if (component_mask.n_selected_components(n_components) == n_components)
1315  {
1316  for (const auto &cell : dof.cell_iterators())
1317  {
1318  if (dof.get_triangulation().locally_owned_subdomain() !=
1320  cell->level_subdomain_id() == numbers::artificial_subdomain_id)
1321  continue;
1322  const FiniteElement<dim> &fe = cell->get_fe();
1323  const unsigned int level = cell->level();
1324  local_dofs.resize(fe.dofs_per_face);
1325 
1326  for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
1327  if (cell->at_boundary(face_no) == true)
1328  {
1329  const typename DoFHandler<dim, spacedim>::face_iterator face =
1330  cell->face(face_no);
1331  const types::boundary_id bi = face->boundary_id();
1332  // Face is listed in boundary map
1333  if (boundary_ids.find(bi) != boundary_ids.end())
1334  {
1335  face->get_mg_dof_indices(level, local_dofs);
1336  dofs_by_level[level].insert(dofs_by_level[level].end(),
1337  local_dofs.begin(),
1338  local_dofs.end());
1339  }
1340  }
1341  }
1342  }
1343  else
1344  {
1345  Assert(component_mask.n_selected_components(n_components) > 0,
1346  ExcMessage(
1347  "It's probably worthwhile to select at least one component."));
1348 
1349  for (const auto &cell : dof.cell_iterators())
1350  if (dof.get_triangulation().locally_owned_subdomain() ==
1352  cell->level_subdomain_id() != numbers::artificial_subdomain_id)
1353  for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
1354  {
1355  if (cell->at_boundary(face_no) == false)
1356  continue;
1357 
1358  const FiniteElement<dim> &fe = cell->get_fe();
1359  const unsigned int level = cell->level();
1360 
1362  cell->face(face_no);
1363  const types::boundary_id boundary_component =
1364  face->boundary_id();
1365  if (boundary_ids.find(boundary_component) != boundary_ids.end())
1366  // we want to constrain this boundary
1367  {
1368  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
1369  ++i)
1370  {
1371  const ComponentMask &nonzero_component_array =
1372  cell->get_fe().get_nonzero_components(i);
1373  // if we want to constrain one of the nonzero
1374  // components, we have to constrain all of them
1375 
1376  bool selected = false;
1377  for (unsigned int c = 0; c < n_components; ++c)
1378  if (nonzero_component_array[c] == true &&
1379  component_mask[c] == true)
1380  {
1381  selected = true;
1382  break;
1383  }
1384  if (selected)
1385  for (unsigned int c = 0; c < n_components; ++c)
1386  Assert(
1387  nonzero_component_array[c] == false ||
1388  component_mask[c] == true,
1389  ExcMessage(
1390  "You are using a non-primitive FiniteElement "
1391  "and try to constrain just some of its components!"));
1392  }
1393 
1394  // get indices, physical location and boundary values of
1395  // dofs on this face
1396  local_dofs.resize(fe.dofs_per_face);
1397  face->get_mg_dof_indices(level, local_dofs);
1398  if (fe_is_system)
1399  {
1400  for (unsigned int i = 0; i < local_dofs.size(); ++i)
1401  {
1402  unsigned int component =
1404  if (fe.is_primitive())
1405  component =
1406  fe.face_system_to_component_index(i).first;
1407  else
1408  {
1409  // Just pick the first of the components
1410  // We already know that either all or none
1411  // of the components are selected
1412  const ComponentMask &nonzero_component_array =
1413  cell->get_fe().get_nonzero_components(i);
1414  for (unsigned int c = 0; c < n_components; ++c)
1415  if (nonzero_component_array[c] == true)
1416  {
1417  component = c;
1418  break;
1419  }
1420  }
1422  ExcInternalError());
1423  if (component_mask[component] == true)
1424  dofs_by_level[level].push_back(local_dofs[i]);
1425  }
1426  }
1427  else
1428  dofs_by_level[level].insert(dofs_by_level[level].end(),
1429  local_dofs.begin(),
1430  local_dofs.end());
1431  }
1432  }
1433  }
1434  for (unsigned int level = 0; level < dof.get_triangulation().n_levels();
1435  ++level)
1436  {
1437  std::sort(dofs_by_level[level].begin(), dofs_by_level[level].end());
1438  boundary_indices[level].add_indices(
1439  dofs_by_level[level].begin(),
1440  std::unique(dofs_by_level[level].begin(),
1441  dofs_by_level[level].end()));
1442  }
1443  }
1444 
1445 
1446  template <int dim, int spacedim>
1447  void
1449  const DoFHandler<dim, spacedim> & mg_dof_handler,
1450  std::vector<std::set<types::global_dof_index>> &non_interface_dofs)
1451  {
1452  Assert(non_interface_dofs.size() ==
1453  mg_dof_handler.get_triangulation().n_global_levels(),
1455  non_interface_dofs.size(),
1456  mg_dof_handler.get_triangulation().n_global_levels()));
1457 
1458  const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();
1459 
1460  const unsigned int dofs_per_cell = fe.dofs_per_cell;
1461  const unsigned int dofs_per_face = fe.dofs_per_face;
1462 
1463  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1464  std::vector<bool> cell_dofs(dofs_per_cell, false);
1465  std::vector<bool> cell_dofs_interface(dofs_per_cell, false);
1466 
1467  typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
1468  endc = mg_dof_handler.end();
1469 
1470 
1471  for (; cell != endc; ++cell)
1472  {
1473  if (mg_dof_handler.get_triangulation().locally_owned_subdomain() !=
1475  cell->level_subdomain_id() !=
1476  mg_dof_handler.get_triangulation().locally_owned_subdomain())
1477  continue;
1478 
1479  std::fill(cell_dofs.begin(), cell_dofs.end(), false);
1480  std::fill(cell_dofs_interface.begin(),
1481  cell_dofs_interface.end(),
1482  false);
1483 
1484  for (const unsigned int face_nr : GeometryInfo<dim>::face_indices())
1485  {
1486  const typename DoFHandler<dim, spacedim>::face_iterator face =
1487  cell->face(face_nr);
1488  if (!face->at_boundary() || cell->has_periodic_neighbor(face_nr))
1489  {
1490  // interior face
1491  const typename DoFHandler<dim>::cell_iterator neighbor =
1492  cell->neighbor_or_periodic_neighbor(face_nr);
1493 
1494  if ((neighbor->level() < cell->level()))
1495  {
1496  for (unsigned int j = 0; j < dofs_per_face; ++j)
1497  cell_dofs_interface[fe.face_to_cell_index(j, face_nr)] =
1498  true;
1499  }
1500  else
1501  {
1502  for (unsigned int j = 0; j < dofs_per_face; ++j)
1503  cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
1504  }
1505  }
1506  else
1507  {
1508  // boundary face
1509  for (unsigned int j = 0; j < dofs_per_face; ++j)
1510  cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
1511  }
1512  }
1513 
1514  const unsigned int level = cell->level();
1515  cell->get_mg_dof_indices(local_dof_indices);
1516 
1517  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1518  if (cell_dofs[i] && !cell_dofs_interface[i])
1519  non_interface_dofs[level].insert(local_dof_indices[i]);
1520  }
1521  }
1522 
1523 
1524  template <int dim, int spacedim>
1525  void
1527  std::vector<IndexSet> & interface_dofs)
1528  {
1529  Assert(interface_dofs.size() ==
1530  mg_dof_handler.get_triangulation().n_global_levels(),
1532  interface_dofs.size(),
1533  mg_dof_handler.get_triangulation().n_global_levels()));
1534 
1535  std::vector<std::vector<types::global_dof_index>> tmp_interface_dofs(
1536  interface_dofs.size());
1537 
1538  const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();
1539 
1540  const unsigned int dofs_per_cell = fe.dofs_per_cell;
1541  const unsigned int dofs_per_face = fe.dofs_per_face;
1542 
1543  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1544 
1545  std::vector<bool> cell_dofs(dofs_per_cell, false);
1546 
1547  typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
1548  endc = mg_dof_handler.end();
1549 
1550  for (; cell != endc; ++cell)
1551  {
1552  // Do not look at artificial level cells (in a serial computation we
1553  // need to ignore the level_subdomain_id() because it is never set).
1554  if (mg_dof_handler.get_triangulation().locally_owned_subdomain() !=
1556  cell->level_subdomain_id() == numbers::artificial_subdomain_id)
1557  continue;
1558 
1559  bool has_coarser_neighbor = false;
1560 
1561  std::fill(cell_dofs.begin(), cell_dofs.end(), false);
1562 
1563  for (const unsigned int face_nr : GeometryInfo<dim>::face_indices())
1564  {
1565  const typename DoFHandler<dim, spacedim>::face_iterator face =
1566  cell->face(face_nr);
1567  if (!face->at_boundary() || cell->has_periodic_neighbor(face_nr))
1568  {
1569  // interior face
1570  const typename DoFHandler<dim>::cell_iterator neighbor =
1571  cell->neighbor_or_periodic_neighbor(face_nr);
1572 
1573  // only process cell pairs if one or both of them are owned by
1574  // me (ignore if running in serial)
1575  if (mg_dof_handler.get_triangulation()
1576  .locally_owned_subdomain() !=
1578  neighbor->level_subdomain_id() ==
1580  continue;
1581 
1582  // Do refinement face from the coarse side
1583  if (neighbor->level() < cell->level())
1584  {
1585  for (unsigned int j = 0; j < dofs_per_face; ++j)
1586  cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
1587 
1588  has_coarser_neighbor = true;
1589  }
1590  }
1591  }
1592 
1593  if (has_coarser_neighbor == false)
1594  continue;
1595 
1596  const unsigned int level = cell->level();
1597  cell->get_mg_dof_indices(local_dof_indices);
1598 
1599  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1600  {
1601  if (cell_dofs[i])
1602  tmp_interface_dofs[level].push_back(local_dof_indices[i]);
1603  }
1604  }
1605 
1606  for (unsigned int l = 0;
1607  l < mg_dof_handler.get_triangulation().n_global_levels();
1608  ++l)
1609  {
1610  interface_dofs[l].clear();
1611  std::sort(tmp_interface_dofs[l].begin(), tmp_interface_dofs[l].end());
1612  interface_dofs[l].add_indices(tmp_interface_dofs[l].begin(),
1613  std::unique(tmp_interface_dofs[l].begin(),
1614  tmp_interface_dofs[l].end()));
1615  interface_dofs[l].compress();
1616  }
1617  }
1618 
1619 
1620 
1621  template <int dim, int spacedim>
1622  unsigned int
1624  {
1625  // Find minimum level for an active cell in
1626  // this locally owned subdomain
1627  // Note: with the way active cells are traversed,
1628  // the first locally owned cell we find will have
1629  // the lowest level in the particular subdomain.
1630  unsigned int min_level = tria.n_global_levels();
1632  cell = tria.begin_active(),
1633  endc = tria.end();
1634  for (; cell != endc; ++cell)
1635  if (cell->is_locally_owned())
1636  {
1637  min_level = cell->level();
1638  break;
1639  }
1640 
1641  unsigned int global_min = min_level;
1642  // If necessary, communicate to find minimum
1643  // level for an active cell over all subdomains
1645  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1646  &tria))
1647  global_min = Utilities::MPI::min(min_level, tr->get_communicator());
1648 
1649  AssertIndexRange(global_min, tria.n_global_levels());
1650 
1651  return global_min;
1652  }
1653 
1654 
1655 
1656  template <int dim, int spacedim>
1657  double
1659  {
1660  double workload_imbalance = 1.0;
1661 
1662  // It is only necessary to calculate the imbalance
1663  // on a distributed mesh. The imbalance is always
1664  // 1.0 for the serial case.
1666  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1667  &tria))
1668  {
1669  Assert(
1670  tr->is_multilevel_hierarchy_constructed(),
1671  ExcMessage(
1672  "We can only compute the workload imbalance if the multilevel hierarchy has been constructed!"));
1673 
1674  const unsigned int n_proc =
1675  Utilities::MPI::n_mpi_processes(tr->get_communicator());
1676  const unsigned int n_global_levels = tr->n_global_levels();
1677 
1678  // This value will represent the sum over the multigrid
1679  // levels of the maximum number of cells owned by any
1680  // one processesor on that level.
1681  types::global_dof_index work_estimate = 0;
1682 
1683  // Sum of all cells in the multigrid hierarchy
1684  types::global_dof_index total_cells_in_hierarchy = 0;
1685 
1686  for (int lvl = n_global_levels - 1; lvl >= 0; --lvl)
1687  {
1688  // Number of cells this processor owns on this level
1689  types::global_dof_index n_owned_cells_on_lvl = 0;
1690 
1691  for (const auto &cell : tr->cell_iterators_on_level(lvl))
1692  if (cell->is_locally_owned_on_level())
1693  ++n_owned_cells_on_lvl;
1694 
1695  work_estimate +=
1696  ::Utilities::MPI::max(n_owned_cells_on_lvl,
1697  tr->get_communicator());
1698 
1699  total_cells_in_hierarchy +=
1700  ::Utilities::MPI::sum(n_owned_cells_on_lvl,
1701  tr->get_communicator());
1702  }
1703 
1704  const double ideal_work =
1705  total_cells_in_hierarchy / static_cast<double>(n_proc);
1706  workload_imbalance = work_estimate / ideal_work;
1707  }
1708 
1709  return workload_imbalance;
1710  }
1711 } // namespace MGTools
1712 
1713 
1714 // explicit instantiations
1715 #include "mg_tools.inst"
1716 
const unsigned int first_hex_index
Definition: fe_base.h:256
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:565
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:191
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1560
const types::subdomain_id invalid_subdomain_id
Definition: types.h:285
void clear_user_flags()
Definition: tria.cc:11389
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:1158
cell_iterator begin(const unsigned int level=0) const
double workload_imbalance(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1658
Task< RT > new_task(const std::function< RT()> &function)
void load_user_flags(std::istream &in)
Definition: tria.cc:11449
cell_iterator end() const
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:1628
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:12228
static const char U
bool is_primitive() const
Definition: fe.h:3253
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)
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:1023
typename ActiveSelector::face_iterator face_iterator
Definition: dof_handler.h:389
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:12294
static ::ExceptionBase & ExcMessage(std::string arg1)
const unsigned int first_quad_index
Definition: fe_base.h:251
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1403
unsigned int element_multiplicity(const unsigned int index) const
Definition: fe.h:3085
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
unsigned int level
Definition: grid_out.cc:4355
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index) const
Definition: fe.h:3122
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:1242
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
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
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
const unsigned int dofs_per_cell
Definition: fe_base.h:280
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:3173
const types::subdomain_id artificial_subdomain_id
Definition: types.h:302
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:1241
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:3062
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:3205
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
virtual unsigned int n_global_levels() const
VectorType::value_type * begin(VectorType &V)
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:1448
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:273
const unsigned int first_line_index
Definition: fe_base.h:246
void extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition: mg_tools.cc:1526
void save_user_flags(std::ostream &out) const
Definition: tria.cc:11400
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:687
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:371
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:206
unsigned int n_base_elements() const
Definition: fe.h:3076
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:2410
T max(const T &t, const MPI_Comm &mpi_communicator)
unsigned int max_level_for_coarse_mesh(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1623
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
Definition: fe.cc:1078
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()