Reference documentation for deal.II version 9.3.3
\(\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
19
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
50namespace MGTools
51{
52 // specializations for 1D
53 template <>
54 void
56 const unsigned int,
57 std::vector<unsigned int> &,
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> &,
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 =
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 =
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);
610 endc = dof.end(level);
611 for (; cell != endc; ++cell)
614 cell->level_subdomain_id() ==
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);
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 {
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()),
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);
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 {
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
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 {
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 {
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()),
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
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 {
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);
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);
1131 for (unsigned int i = 0; i < n_components; ++i)
1132 {
1133 void (*fun_ptr)(const unsigned int level,
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),
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);
1240 for (unsigned int i = 0; i < n_blocks; ++i)
1241 {
1242 void (*fun_ptr)(const unsigned int level,
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(),
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(),
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 {
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())
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 }
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()
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 add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
cell_iterator end() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
cell_iterator begin(const unsigned int level=0) const
void clear()
unsigned int get_first_line_index() const
unsigned int n_dofs_per_cell() const
unsigned int get_first_quad_index(const unsigned int quad_no=0) const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_blocks() const
unsigned int n_components() const
ReferenceCell reference_cell() const
unsigned int get_first_hex_index() const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
bool is_primitive() const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
unsigned int element_multiplicity(const unsigned int index) const
unsigned int n_base_elements() const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
types::global_dof_index first_block_of_base(const unsigned int b) const
bool is_interface_matrix_entry(const unsigned int level, const types::global_dof_index i, const types::global_dof_index j) const
virtual types::subdomain_id locally_owned_subdomain() const
void save_user_flags(std::ostream &out) const
unsigned int n_levels() const
cell_iterator end() const
void clear_user_flags()
virtual unsigned int n_global_levels() const
void load_user_flags(std::istream &in)
active_cell_iterator begin_active(const unsigned int level=0) const
unsigned int max_dofs_per_face() const
unsigned int n_components() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
unsigned int level
Definition: grid_out.cc:4590
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
Task< RT > new_task(const std::function< RT()> &function)
void convert_couplings_to_blocks(const DoFHandler< dim, spacedim > &dof_handler, const Table< 2, Coupling > &table_by_component, std::vector< Table< 2, Coupling > > &tables_by_block)
Definition: dof_tools.cc:2462
Table< 2, Coupling > dof_couplings_from_component_couplings(const FiniteElement< dim, spacedim > &fe, const Table< 2, Coupling > &component_couplings)
static const char U
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
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
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:628
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
void extract_inner_interface_dofs(const DoFHandler< dim, spacedim > &mg_dof_handler, std::vector< IndexSet > &interface_dofs)
Definition: mg_tools.cc:1480
double workload_imbalance(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1612
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
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
void make_flux_sparsity_pattern_edge(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:703
unsigned int max_level_for_coarse_mesh(const Triangulation< dim, spacedim > &tria)
Definition: mg_tools.cc:1577
std::enable_if< IsBlockVector< VectorType >::value, unsignedint >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T min(const T &t, const MPI_Comm &mpi_communicator)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
const types::subdomain_id invalid_subdomain_id
Definition: types.h:276
static const unsigned int invalid_unsigned_int
Definition: types.h:196
const types::global_dof_index invalid_dof_index
Definition: types.h:211