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