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