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