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