deal.II version GIT relicensing-2863-gb18dd3416c 2025-03-19 12:40: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 (void)n_dofs;
585
586 Assert(sparsity.n_rows() == n_dofs,
587 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
588 Assert(sparsity.n_cols() == n_dofs,
589 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
590
591 const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
592 std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
593 for (const auto &cell : dof.cell_iterators_on_level(level))
594 if (cell->is_locally_owned_on_level())
595 {
596 cell->get_mg_dof_indices(dofs_on_this_cell);
597 constraints.add_entries_local_to_global(dofs_on_this_cell,
598 sparsity,
599 keep_constrained_dofs);
600 }
601 }
602
603
604
605 template <int dim, int spacedim, typename number>
606 void
608 SparsityPatternBase &sparsity,
609 const unsigned int level,
610 const AffineConstraints<number> &constraints,
611 const bool keep_constrained_dofs)
612 {
613 const types::global_dof_index n_dofs = dof.n_dofs(level);
614 (void)n_dofs;
615
616 Assert(sparsity.n_rows() == n_dofs,
617 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
618 Assert(sparsity.n_cols() == n_dofs,
619 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
620
621 const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
622 std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
623 std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
624 for (const auto &cell : dof.cell_iterators_on_level(level))
625 {
626 if (!cell->is_locally_owned_on_level())
627 continue;
628
629 cell->get_mg_dof_indices(dofs_on_this_cell);
630 // make sparsity pattern for this cell
631 constraints.add_entries_local_to_global(dofs_on_this_cell,
632 sparsity,
633 keep_constrained_dofs);
634
635 // Loop over all interior neighbors
636 for (const unsigned int face : GeometryInfo<dim>::face_indices())
637 {
638 bool use_face = false;
639 if ((!cell->at_boundary(face)) &&
640 (static_cast<unsigned int>(cell->neighbor_level(face)) ==
641 level))
642 use_face = true;
643 else if (cell->has_periodic_neighbor(face) &&
644 (static_cast<unsigned int>(
645 cell->periodic_neighbor_level(face)) == level))
646 use_face = true;
647
648 if (use_face)
649 {
651 cell->neighbor_or_periodic_neighbor(face);
652 neighbor->get_mg_dof_indices(dofs_on_other_cell);
653 // only add one direction The other is taken care of by
654 // neighbor (except when the neighbor is not owned by the same
655 // processor)
656 constraints.add_entries_local_to_global(dofs_on_this_cell,
657 dofs_on_other_cell,
658 sparsity,
659 keep_constrained_dofs);
660
661 if (neighbor->is_locally_owned_on_level() == false)
662 {
663 constraints.add_entries_local_to_global(
664 dofs_on_other_cell, sparsity, keep_constrained_dofs);
665
666 constraints.add_entries_local_to_global(
667 dofs_on_other_cell,
668 dofs_on_this_cell,
669 sparsity,
670 keep_constrained_dofs);
671 }
672 }
673 }
674 }
675 }
676
677
678
679 template <int dim, int spacedim>
680 void
682 SparsityPatternBase &sparsity,
683 const unsigned int level)
684 {
685 Assert((level >= 1) && (level < dof.get_triangulation().n_global_levels()),
687
688 const types::global_dof_index fine_dofs = dof.n_dofs(level);
689 const types::global_dof_index coarse_dofs = dof.n_dofs(level - 1);
690 (void)fine_dofs;
691 (void)coarse_dofs;
692
693 // Matrix maps from fine level to coarse level
694
695 Assert(sparsity.n_rows() == coarse_dofs,
696 ExcDimensionMismatch(sparsity.n_rows(), coarse_dofs));
697 Assert(sparsity.n_cols() == fine_dofs,
698 ExcDimensionMismatch(sparsity.n_cols(), fine_dofs));
699
700 const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
701 std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
702 std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
703 for (const auto &cell : dof.cell_iterators_on_level(level))
704 {
705 if (!cell->is_locally_owned_on_level())
706 continue;
707
708 cell->get_mg_dof_indices(dofs_on_this_cell);
709 // Loop over all interior neighbors
710 for (const unsigned int face : GeometryInfo<dim>::face_indices())
711 {
712 // Neighbor is coarser
713 bool use_face = false;
714 if ((!cell->at_boundary(face)) &&
715 (static_cast<unsigned int>(cell->neighbor_level(face)) !=
716 level))
717 use_face = true;
718 else if (cell->has_periodic_neighbor(face) &&
719 (static_cast<unsigned int>(
720 cell->periodic_neighbor_level(face)) != level))
721 use_face = true;
722
723 if (use_face)
724 {
726 cell->neighbor_or_periodic_neighbor(face);
727 neighbor->get_mg_dof_indices(dofs_on_other_cell);
728
729 std::sort(dofs_on_this_cell.begin(), dofs_on_this_cell.end());
730 for (unsigned int i = 0; i < dofs_per_cell; ++i)
731 sparsity.add_row_entries(dofs_on_other_cell[i],
732 make_array_view(dofs_on_this_cell),
733 true);
734 }
735 }
736 }
737 }
738
739
740
741 template <int dim, int spacedim>
742 void
744 SparsityPatternBase &sparsity,
745 const unsigned int level,
746 const Table<2, DoFTools::Coupling> &int_mask,
747 const Table<2, DoFTools::Coupling> &flux_mask)
748 {
749 const FiniteElement<dim> &fe = dof.get_fe();
750 const types::global_dof_index n_dofs = dof.n_dofs(level);
751 const unsigned int n_comp = fe.n_components();
752 (void)n_dofs;
753 (void)n_comp;
754
755 Assert(sparsity.n_rows() == n_dofs,
756 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
757 Assert(sparsity.n_cols() == n_dofs,
758 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
759 Assert(int_mask.n_rows() == n_comp,
760 ExcDimensionMismatch(int_mask.n_rows(), n_comp));
761 Assert(int_mask.n_cols() == n_comp,
762 ExcDimensionMismatch(int_mask.n_cols(), n_comp));
763 Assert(flux_mask.n_rows() == n_comp,
764 ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
765 Assert(flux_mask.n_cols() == n_comp,
766 ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
767
768 const unsigned int total_dofs = fe.n_dofs_per_cell();
769 std::vector<types::global_dof_index> dofs_on_this_cell(total_dofs);
770 std::vector<types::global_dof_index> dofs_on_other_cell(total_dofs);
771 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
772 cell_entries;
773
774 Table<2, bool> support_on_face(total_dofs,
776
778 int_dof_mask =
780 flux_dof_mask =
782
783 for (unsigned int i = 0; i < total_dofs; ++i)
784 for (auto f : GeometryInfo<dim>::face_indices())
785 support_on_face(i, f) = fe.has_support_on_face(i, f);
786
787 std::vector<bool> face_touched(dim == 2 ?
790
791 for (const auto &cell : dof.cell_iterators_on_level(level))
792 {
793 if (!cell->is_locally_owned_on_level())
794 continue;
795
796 cell->get_mg_dof_indices(dofs_on_this_cell);
797 // make sparsity pattern for this cell
798 for (unsigned int i = 0; i < total_dofs; ++i)
799 for (unsigned int j = 0; j < total_dofs; ++j)
800 if (int_dof_mask[i][j] != DoFTools::none)
801 cell_entries.emplace_back(dofs_on_this_cell[i],
802 dofs_on_this_cell[j]);
803
804 // Loop over all interior neighbors
805 for (const unsigned int face : GeometryInfo<dim>::face_indices())
806 {
808 cell->face(face);
809 if (face_touched[cell_face->index()])
810 continue;
811
812 if (cell->at_boundary(face) && !cell->has_periodic_neighbor(face))
813 {
814 for (unsigned int i = 0; i < total_dofs; ++i)
815 {
816 const bool i_non_zero_i = support_on_face(i, face);
817 for (unsigned int j = 0; j < total_dofs; ++j)
818 {
819 const bool j_non_zero_i = support_on_face(j, face);
820
821 if (flux_dof_mask(i, j) == DoFTools::always)
822 cell_entries.emplace_back(dofs_on_this_cell[i],
823 dofs_on_this_cell[j]);
824 if (flux_dof_mask(i, j) == DoFTools::nonzero &&
825 i_non_zero_i && j_non_zero_i)
826 cell_entries.emplace_back(dofs_on_this_cell[i],
827 dofs_on_this_cell[j]);
828 }
829 }
830 }
831 else
832 {
834 cell->neighbor_or_periodic_neighbor(face);
835
836 if (neighbor->level() < cell->level())
837 continue;
838
839 unsigned int neighbor_face =
840 cell->has_periodic_neighbor(face) ?
841 cell->periodic_neighbor_of_periodic_neighbor(face) :
842 cell->neighbor_of_neighbor(face);
843
844 neighbor->get_mg_dof_indices(dofs_on_other_cell);
845 for (unsigned int i = 0; i < total_dofs; ++i)
846 {
847 const bool i_non_zero_i = support_on_face(i, face);
848 const bool i_non_zero_e = support_on_face(i, neighbor_face);
849 for (unsigned int j = 0; j < total_dofs; ++j)
850 {
851 const bool j_non_zero_i = support_on_face(j, face);
852 const bool j_non_zero_e =
853 support_on_face(j, neighbor_face);
854 if (flux_dof_mask(i, j) == DoFTools::always)
855 {
856 cell_entries.emplace_back(dofs_on_this_cell[i],
857 dofs_on_other_cell[j]);
858 cell_entries.emplace_back(dofs_on_other_cell[i],
859 dofs_on_this_cell[j]);
860 cell_entries.emplace_back(dofs_on_this_cell[i],
861 dofs_on_this_cell[j]);
862 cell_entries.emplace_back(dofs_on_other_cell[i],
863 dofs_on_other_cell[j]);
864 }
865 if (flux_dof_mask(i, j) == DoFTools::nonzero)
866 {
867 if (i_non_zero_i && j_non_zero_e)
868 cell_entries.emplace_back(dofs_on_this_cell[i],
869 dofs_on_other_cell[j]);
870 if (i_non_zero_e && j_non_zero_i)
871 cell_entries.emplace_back(dofs_on_other_cell[i],
872 dofs_on_this_cell[j]);
873 if (i_non_zero_i && j_non_zero_i)
874 cell_entries.emplace_back(dofs_on_this_cell[i],
875 dofs_on_this_cell[j]);
876 if (i_non_zero_e && j_non_zero_e)
877 cell_entries.emplace_back(dofs_on_other_cell[i],
878 dofs_on_other_cell[j]);
879 }
880
881 if (flux_dof_mask(j, i) == DoFTools::always)
882 {
883 cell_entries.emplace_back(dofs_on_this_cell[j],
884 dofs_on_other_cell[i]);
885 cell_entries.emplace_back(dofs_on_other_cell[j],
886 dofs_on_this_cell[i]);
887 cell_entries.emplace_back(dofs_on_this_cell[j],
888 dofs_on_this_cell[i]);
889 cell_entries.emplace_back(dofs_on_other_cell[j],
890 dofs_on_other_cell[i]);
891 }
892 if (flux_dof_mask(j, i) == DoFTools::nonzero)
893 {
894 if (j_non_zero_i && i_non_zero_e)
895 cell_entries.emplace_back(dofs_on_this_cell[j],
896 dofs_on_other_cell[i]);
897 if (j_non_zero_e && i_non_zero_i)
898 cell_entries.emplace_back(dofs_on_other_cell[j],
899 dofs_on_this_cell[i]);
900 if (j_non_zero_i && i_non_zero_i)
901 cell_entries.emplace_back(dofs_on_this_cell[j],
902 dofs_on_this_cell[i]);
903 if (j_non_zero_e && i_non_zero_e)
904 cell_entries.emplace_back(dofs_on_other_cell[j],
905 dofs_on_other_cell[i]);
906 }
907 }
908 }
909 face_touched[neighbor->face(neighbor_face)->index()] = true;
910 }
911 }
912 sparsity.add_entries(make_array_view(cell_entries));
913 cell_entries.clear();
914 }
915 }
916
917
918
919 template <int dim, int spacedim>
920 void
922 SparsityPatternBase &sparsity,
923 const unsigned int level,
924 const Table<2, DoFTools::Coupling> &flux_mask)
925 {
926 const FiniteElement<dim> &fe = dof.get_fe();
927 const unsigned int n_comp = fe.n_components();
928 (void)n_comp;
929
930 Assert((level >= 1) && (level < dof.get_triangulation().n_global_levels()),
932
933 const types::global_dof_index fine_dofs = dof.n_dofs(level);
934 const types::global_dof_index coarse_dofs = dof.n_dofs(level - 1);
935 (void)fine_dofs;
936 (void)coarse_dofs;
937
938 // Matrix maps from fine level to coarse level
939
940 Assert(sparsity.n_rows() == coarse_dofs,
941 ExcDimensionMismatch(sparsity.n_rows(), coarse_dofs));
942 Assert(sparsity.n_cols() == fine_dofs,
943 ExcDimensionMismatch(sparsity.n_cols(), fine_dofs));
944 Assert(flux_mask.n_rows() == n_comp,
945 ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
946 Assert(flux_mask.n_cols() == n_comp,
947 ExcDimensionMismatch(flux_mask.n_cols(), n_comp));
948
949 const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
950 std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
951 std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
952 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
953 cell_entries;
954
955 Table<2, bool> support_on_face(dofs_per_cell,
957
958 const Table<2, DoFTools::Coupling> flux_dof_mask =
960
961 for (unsigned int i = 0; i < dofs_per_cell; ++i)
962 for (auto f : GeometryInfo<dim>::face_indices())
963 support_on_face(i, f) = fe.has_support_on_face(i, f);
964
965 for (const auto &cell : dof.cell_iterators_on_level(level))
966 {
967 if (!cell->is_locally_owned_on_level())
968 continue;
969
970 cell->get_mg_dof_indices(dofs_on_this_cell);
971 // Loop over all interior neighbors
972 for (const unsigned int face : GeometryInfo<dim>::face_indices())
973 {
974 // Neighbor is coarser
975 bool use_face = false;
976 if ((!cell->at_boundary(face)) &&
977 (static_cast<unsigned int>(cell->neighbor_level(face)) !=
978 level))
979 use_face = true;
980 else if (cell->has_periodic_neighbor(face) &&
981 (static_cast<unsigned int>(
982 cell->periodic_neighbor_level(face)) != level))
983 use_face = true;
984
985 if (use_face)
986 {
988 cell->neighbor_or_periodic_neighbor(face);
989 neighbor->get_mg_dof_indices(dofs_on_other_cell);
990
991 for (unsigned int i = 0; i < dofs_per_cell; ++i)
992 {
993 for (unsigned int j = 0; j < dofs_per_cell; ++j)
994 {
995 if (flux_dof_mask(i, j) != DoFTools::none)
996 {
997 cell_entries.emplace_back(dofs_on_other_cell[i],
998 dofs_on_this_cell[j]);
999 cell_entries.emplace_back(dofs_on_other_cell[j],
1000 dofs_on_this_cell[i]);
1001 }
1002 }
1003 }
1004 }
1005 }
1006 sparsity.add_entries(make_array_view(cell_entries));
1007 cell_entries.clear();
1008 }
1009 }
1010
1011
1012
1013 template <int dim, int spacedim>
1014 void
1016 const MGConstrainedDoFs &mg_constrained_dofs,
1017 SparsityPatternBase &sparsity,
1018 const unsigned int level)
1019 {
1020 const types::global_dof_index n_dofs = dof.n_dofs(level);
1021 (void)n_dofs;
1022
1023 Assert(sparsity.n_rows() == n_dofs,
1024 ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
1025 Assert(sparsity.n_cols() == n_dofs,
1026 ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
1027
1028 const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
1029 std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
1030 std::vector<types::global_dof_index> cols;
1031 cols.reserve(dofs_per_cell);
1032
1033 for (const auto &cell : dof.cell_iterators_on_level(level))
1034 if (cell->is_locally_owned_on_level())
1035 {
1036 cell->get_mg_dof_indices(dofs_on_this_cell);
1037 std::sort(dofs_on_this_cell.begin(), dofs_on_this_cell.end());
1038 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1039 {
1040 for (unsigned int j = 0; j < dofs_per_cell; ++j)
1041 if (mg_constrained_dofs.is_interface_matrix_entry(
1042 level, dofs_on_this_cell[i], dofs_on_this_cell[j]))
1043 cols.push_back(dofs_on_this_cell[j]);
1044 sparsity.add_row_entries(dofs_on_this_cell[i],
1045 make_array_view(cols),
1046 true);
1047 cols.clear();
1048 }
1049 }
1050 }
1051
1052
1053
1054 template <int dim, int spacedim>
1055 void
1057 const DoFHandler<dim, spacedim> &dof_handler,
1058 std::vector<std::vector<types::global_dof_index>> &result,
1059 bool only_once,
1060 std::vector<unsigned int> target_component)
1061 {
1062 const FiniteElement<dim> &fe = dof_handler.get_fe();
1063 const unsigned int n_components = fe.n_components();
1064 const unsigned int nlevels =
1065 dof_handler.get_triangulation().n_global_levels();
1066
1067 Assert(result.size() == nlevels,
1068 ExcDimensionMismatch(result.size(), nlevels));
1069
1070 if (target_component.empty())
1071 {
1072 target_component.resize(n_components);
1073 for (unsigned int i = 0; i < n_components; ++i)
1074 target_component[i] = i;
1075 }
1076
1077 Assert(target_component.size() == n_components,
1078 ExcDimensionMismatch(target_component.size(), n_components));
1079
1080 for (unsigned int l = 0; l < nlevels; ++l)
1081 {
1082 result[l].resize(n_components);
1083 std::fill(result[l].begin(), result[l].end(), 0U);
1084
1085 // special case for only one
1086 // component. treat this first
1087 // since it does not require any
1088 // computations
1089 if (n_components == 1)
1090 {
1091 result[l][0] = dof_handler.n_dofs(l);
1092 }
1093 else
1094 {
1095 // otherwise determine the number
1096 // of dofs in each component
1097 // separately. do so in parallel
1098 std::vector<std::vector<bool>> dofs_in_component(
1099 n_components, std::vector<bool>(dof_handler.n_dofs(l), false));
1100 std::vector<ComponentMask> component_select(n_components);
1102 for (unsigned int i = 0; i < n_components; ++i)
1103 {
1104 void (*fun_ptr)(const unsigned int level,
1106 const ComponentMask &,
1107 std::vector<bool> &) =
1108 &DoFTools::extract_level_dofs<dim, spacedim>;
1109
1110 std::vector<bool> tmp(n_components, false);
1111 tmp[i] = true;
1112 component_select[i] = ComponentMask(tmp);
1113
1114 tasks += Threads::new_task(fun_ptr,
1115 l,
1116 dof_handler,
1117 component_select[i],
1118 dofs_in_component[i]);
1119 }
1120 tasks.join_all();
1121
1122 // next count what we got
1123 unsigned int component = 0;
1124 for (unsigned int b = 0; b < fe.n_base_elements(); ++b)
1125 {
1126 const FiniteElement<dim> &base = fe.base_element(b);
1127 // Dimension of base element
1128 unsigned int d = base.n_components();
1129
1130 for (unsigned int m = 0; m < fe.element_multiplicity(b); ++m)
1131 {
1132 for (unsigned int dd = 0; dd < d; ++dd)
1133 {
1134 if (base.is_primitive() || (!only_once || dd == 0))
1135 result[l][target_component[component]] +=
1136 std::count(dofs_in_component[component].begin(),
1137 dofs_in_component[component].end(),
1138 true);
1139 ++component;
1140 }
1141 }
1142 }
1143 // finally sanity check
1144 Assert(!dof_handler.get_fe().is_primitive() ||
1145 std::accumulate(result[l].begin(),
1146 result[l].end(),
1148 dof_handler.n_dofs(l),
1150 }
1151 }
1152 }
1153
1154
1155
1156 template <int dim, int spacedim>
1157 void
1159 const DoFHandler<dim, spacedim> &dof_handler,
1160 std::vector<std::vector<types::global_dof_index>> &dofs_per_block,
1161 std::vector<unsigned int> target_block)
1162 {
1163 const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
1164 const unsigned int n_blocks = fe.n_blocks();
1165 const unsigned int n_levels =
1166 dof_handler.get_triangulation().n_global_levels();
1167
1168 AssertDimension(dofs_per_block.size(), n_levels);
1169
1170 for (unsigned int l = 0; l < n_levels; ++l)
1171 std::fill(dofs_per_block[l].begin(), dofs_per_block[l].end(), 0U);
1172 // If the empty vector was given as
1173 // default argument, set up this
1174 // vector as identity.
1175 if (target_block.empty())
1176 {
1177 target_block.resize(n_blocks);
1178 for (unsigned int i = 0; i < n_blocks; ++i)
1179 target_block[i] = i;
1180 }
1181 Assert(target_block.size() == n_blocks,
1182 ExcDimensionMismatch(target_block.size(), n_blocks));
1183
1184 const unsigned int max_block =
1185 *std::max_element(target_block.begin(), target_block.end());
1186 const unsigned int n_target_blocks = max_block + 1;
1187 (void)n_target_blocks;
1188
1189 for (unsigned int l = 0; l < n_levels; ++l)
1190 AssertDimension(dofs_per_block[l].size(), n_target_blocks);
1191
1192 // special case for only one
1193 // block. treat this first
1194 // since it does not require any
1195 // computations
1196 if (n_blocks == 1)
1197 {
1198 for (unsigned int l = 0; l < n_levels; ++l)
1199 dofs_per_block[l][0] = dof_handler.n_dofs(l);
1200 return;
1201 }
1202 // otherwise determine the number
1203 // of dofs in each block
1204 // separately. do so in parallel
1205 for (unsigned int l = 0; l < n_levels; ++l)
1206 {
1207 std::vector<std::vector<bool>> dofs_in_block(
1208 n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
1209 std::vector<BlockMask> block_select(n_blocks);
1211 for (unsigned int i = 0; i < n_blocks; ++i)
1212 {
1213 void (*fun_ptr)(const unsigned int level,
1215 const BlockMask &,
1216 std::vector<bool> &) =
1217 &DoFTools::extract_level_dofs<dim, spacedim>;
1218
1219 std::vector<bool> tmp(n_blocks, false);
1220 tmp[i] = true;
1221 block_select[i] = BlockMask(tmp);
1222
1223 tasks += Threads::new_task(
1224 fun_ptr, l, dof_handler, block_select[i], dofs_in_block[i]);
1225 }
1226 tasks.join_all();
1227
1228 // next count what we got
1229 for (unsigned int block = 0; block < fe.n_blocks(); ++block)
1230 dofs_per_block[l][target_block[block]] +=
1231 std::count(dofs_in_block[block].begin(),
1232 dofs_in_block[block].end(),
1233 true);
1234 }
1235 }
1236
1237
1238
1239 template <int dim, int spacedim>
1240 void
1242 const DoFHandler<dim, spacedim> &dof,
1243 const std::map<types::boundary_id, const Function<spacedim> *>
1244 &function_map,
1245 std::vector<std::set<types::global_dof_index>> &boundary_indices,
1246 const ComponentMask &component_mask)
1247 {
1248 Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
1249 ExcDimensionMismatch(boundary_indices.size(),
1251
1252 std::set<types::boundary_id> boundary_ids;
1253 for (const auto &boundary_function : function_map)
1254 boundary_ids.insert(boundary_function.first);
1255
1256 std::vector<IndexSet> boundary_indexset;
1257 make_boundary_list(dof, boundary_ids, boundary_indexset, component_mask);
1258 for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
1259 boundary_indices[i].insert(boundary_indexset[i].begin(),
1260 boundary_indexset[i].end());
1261 }
1262
1263
1264 template <int dim, int spacedim>
1265 void
1267 const std::map<types::boundary_id,
1268 const Function<spacedim> *> &function_map,
1269 std::vector<IndexSet> &boundary_indices,
1270 const ComponentMask &component_mask)
1271 {
1272 Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
1273 ExcDimensionMismatch(boundary_indices.size(),
1275
1276 std::set<types::boundary_id> boundary_ids;
1277 for (const auto &boundary_function : function_map)
1278 boundary_ids.insert(boundary_function.first);
1279
1280 make_boundary_list(dof, boundary_ids, boundary_indices, component_mask);
1281 }
1282
1283
1284
1285 template <int dim, int spacedim>
1286 void
1288 const std::set<types::boundary_id> &boundary_ids,
1289 std::vector<IndexSet> &boundary_indices,
1290 const ComponentMask &component_mask)
1291 {
1292 boundary_indices.resize(dof.get_triangulation().n_global_levels());
1293
1294 // if for whatever reason we were passed an empty set, return immediately
1295 if (boundary_ids.empty())
1296 return;
1297
1298 for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
1299 if (boundary_indices[i].size() == 0)
1300 boundary_indices[i] = IndexSet(dof.n_dofs(i));
1301
1302 const unsigned int n_components = dof.get_fe_collection().n_components();
1303 const bool fe_is_system = (n_components != 1);
1304
1305 std::vector<types::global_dof_index> local_dofs;
1306 local_dofs.reserve(dof.get_fe_collection().max_dofs_per_face());
1307 std::fill(local_dofs.begin(), local_dofs.end(), numbers::invalid_dof_index);
1308
1309 std::vector<std::vector<types::global_dof_index>> dofs_by_level(
1310 dof.get_triangulation().n_levels());
1311
1312 // First, deal with the simpler case when we have to identify all boundary
1313 // dofs
1314 if (component_mask.n_selected_components(n_components) == n_components)
1315 {
1316 for (const auto &cell : dof.cell_iterators())
1317 {
1318 if (cell->is_artificial_on_level())
1319 continue;
1320 const FiniteElement<dim> &fe = cell->get_fe();
1321 const unsigned int level = cell->level();
1322
1323 for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
1324 if (cell->at_boundary(face_no) == true)
1325 {
1326 const typename DoFHandler<dim, spacedim>::face_iterator face =
1327 cell->face(face_no);
1328 const types::boundary_id bi = face->boundary_id();
1329 // Face is listed in boundary map
1330 if (boundary_ids.find(bi) != boundary_ids.end())
1331 {
1332 local_dofs.resize(fe.n_dofs_per_face(face_no));
1333 face->get_mg_dof_indices(level, local_dofs);
1334 dofs_by_level[level].insert(dofs_by_level[level].end(),
1335 local_dofs.begin(),
1336 local_dofs.end());
1337 }
1338 }
1339 }
1340 }
1341 else
1342 {
1343 Assert(component_mask.n_selected_components(n_components) > 0,
1344 ExcMessage(
1345 "It's probably worthwhile to select at least one component."));
1346
1347 for (const auto &cell : dof.cell_iterators())
1348 if (!cell->is_artificial_on_level())
1349 for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
1350 {
1351 if (cell->at_boundary(face_no) == false)
1352 continue;
1353
1354 const FiniteElement<dim> &fe = cell->get_fe();
1355 const unsigned int level = cell->level();
1356
1358 cell->face(face_no);
1359 const types::boundary_id boundary_component =
1360 face->boundary_id();
1361 if (boundary_ids.find(boundary_component) != boundary_ids.end())
1362 // we want to constrain this boundary
1363 {
1364 for (unsigned int i = 0;
1365 i < cell->get_fe().n_dofs_per_cell();
1366 ++i)
1367 {
1368 const ComponentMask &nonzero_component_array =
1369 cell->get_fe().get_nonzero_components(i);
1370 // if we want to constrain one of the nonzero
1371 // components, we have to constrain all of them
1372
1373 bool selected = false;
1374 for (unsigned int c = 0; c < n_components; ++c)
1375 if (nonzero_component_array[c] == true &&
1376 component_mask[c] == true)
1377 {
1378 selected = true;
1379 break;
1380 }
1381 if (selected)
1382 for (unsigned int c = 0; c < n_components; ++c)
1383 Assert(
1384 nonzero_component_array[c] == false ||
1385 component_mask[c] == true,
1386 ExcMessage(
1387 "You are using a non-primitive FiniteElement "
1388 "and try to constrain just some of its components!"));
1389 }
1390
1391 // get indices, physical location and boundary values of
1392 // dofs on this face
1393 local_dofs.resize(fe.n_dofs_per_face(face_no));
1394 face->get_mg_dof_indices(level, local_dofs);
1395 if (fe_is_system)
1396 {
1397 for (unsigned int i = 0; i < local_dofs.size(); ++i)
1398 {
1399 unsigned int component =
1401 if (fe.is_primitive())
1402 component =
1403 fe.face_system_to_component_index(i, face_no)
1404 .first;
1405 else
1406 {
1407 // Just pick the first of the components
1408 // We already know that either all or none
1409 // of the components are selected
1410 const ComponentMask &nonzero_component_array =
1411 cell->get_fe().get_nonzero_components(i);
1412 for (unsigned int c = 0; c < n_components; ++c)
1413 if (nonzero_component_array[c] == true)
1414 {
1415 component = c;
1416 break;
1417 }
1418 }
1421 if (component_mask[component] == true)
1422 dofs_by_level[level].push_back(local_dofs[i]);
1423 }
1424 }
1425 else
1426 dofs_by_level[level].insert(dofs_by_level[level].end(),
1427 local_dofs.begin(),
1428 local_dofs.end());
1429 }
1430 }
1431 }
1432 for (unsigned int level = 0; level < dof.get_triangulation().n_levels();
1433 ++level)
1434 {
1435 std::sort(dofs_by_level[level].begin(), dofs_by_level[level].end());
1436 boundary_indices[level].add_indices(dofs_by_level[level].begin(),
1437 dofs_by_level[level].end());
1438 }
1439 }
1440
1441
1442
1443 template <int dim, int spacedim>
1444 void
1446 std::vector<IndexSet> &interface_dofs)
1447 {
1448 Assert(interface_dofs.size() ==
1449 mg_dof_handler.get_triangulation().n_global_levels(),
1451 interface_dofs.size(),
1452 mg_dof_handler.get_triangulation().n_global_levels()));
1453
1454 std::vector<std::vector<types::global_dof_index>> tmp_interface_dofs(
1455 interface_dofs.size());
1456
1457 const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();
1458
1459 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1460
1461 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1462
1463 std::vector<bool> cell_dofs(dofs_per_cell, false);
1464
1465 for (const auto &cell : mg_dof_handler.cell_iterators())
1466 {
1467 // Do not look at artificial level cells (in a serial computation we
1468 // need to ignore the level_subdomain_id() because it is never set).
1469 if (cell->is_artificial_on_level())
1470 continue;
1471
1472 bool has_coarser_neighbor = false;
1473
1474 std::fill(cell_dofs.begin(), cell_dofs.end(), false);
1475
1476 for (const unsigned int face_nr : GeometryInfo<dim>::face_indices())
1477 {
1478 const typename DoFHandler<dim, spacedim>::face_iterator face =
1479 cell->face(face_nr);
1480 if (!face->at_boundary() || cell->has_periodic_neighbor(face_nr))
1481 {
1482 // interior face
1483 const typename DoFHandler<dim>::cell_iterator neighbor =
1484 cell->neighbor_or_periodic_neighbor(face_nr);
1485
1486 // only process cell pairs if one or both of them are owned by
1487 // me (ignore if running in serial)
1488 if (neighbor->is_artificial_on_level())
1489 continue;
1490
1491 // Do refinement face from the coarse side
1492 if (neighbor->level() < cell->level())
1493 {
1494 for (unsigned int j = 0; j < fe.n_dofs_per_face(face_nr);
1495 ++j)
1496 cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
1497
1498 has_coarser_neighbor = true;
1499 }
1500 }
1501 }
1502
1503 if (has_coarser_neighbor == false)
1504 continue;
1505
1506 const unsigned int level = cell->level();
1507 cell->get_mg_dof_indices(local_dof_indices);
1508
1509 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1510 {
1511 if (cell_dofs[i])
1512 tmp_interface_dofs[level].push_back(local_dof_indices[i]);
1513 }
1514 }
1515
1516 for (unsigned int l = 0;
1517 l < mg_dof_handler.get_triangulation().n_global_levels();
1518 ++l)
1519 {
1520 interface_dofs[l].clear();
1521 std::sort(tmp_interface_dofs[l].begin(), tmp_interface_dofs[l].end());
1522 interface_dofs[l].add_indices(tmp_interface_dofs[l].begin(),
1523 tmp_interface_dofs[l].end());
1524 interface_dofs[l].compress();
1525 }
1526 }
1527
1528
1529
1530 template <int dim, int spacedim>
1531 unsigned int
1533 {
1534 // Find minimum level for an active cell in
1535 // this locally owned subdomain
1536 // Note: with the way active cells are traversed,
1537 // the first locally owned cell we find will have
1538 // the lowest level in the particular subdomain.
1539 unsigned int min_level = tria.n_global_levels();
1540 for (const auto &cell : tria.active_cell_iterators())
1541 if (cell->is_locally_owned())
1542 {
1543 min_level = cell->level();
1544 break;
1545 }
1546
1547 unsigned int global_min = min_level;
1548 // If necessary, communicate to find minimum
1549 // level for an active cell over all subdomains
1551 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1552 &tria))
1553 global_min = Utilities::MPI::min(min_level, tr->get_mpi_communicator());
1554
1555 AssertIndexRange(global_min, tria.n_global_levels());
1556
1557 return global_min;
1558 }
1559
1560
1561
1562 namespace internal
1563 {
1564 double
1566 const std::vector<types::global_dof_index> &n_cells_on_levels,
1567 const MPI_Comm comm)
1568 {
1569 std::vector<types::global_dof_index> n_cells_on_levels_max(
1570 n_cells_on_levels.size());
1571 std::vector<types::global_dof_index> n_cells_on_levels_sum(
1572 n_cells_on_levels.size());
1573
1574 Utilities::MPI::max(n_cells_on_levels, comm, n_cells_on_levels_max);
1575 Utilities::MPI::sum(n_cells_on_levels, comm, n_cells_on_levels_sum);
1576
1577 const unsigned int n_proc = Utilities::MPI::n_mpi_processes(comm);
1578
1579 const double ideal_work = std::accumulate(n_cells_on_levels_sum.begin(),
1580 n_cells_on_levels_sum.end(),
1581 0) /
1582 static_cast<double>(n_proc);
1583 const double workload_imbalance =
1584 std::accumulate(n_cells_on_levels_max.begin(),
1585 n_cells_on_levels_max.end(),
1586 0) /
1587 ideal_work;
1588
1589 return workload_imbalance;
1590 }
1591
1592
1593
1594 double
1596 const std::vector<
1597 std::pair<types::global_dof_index, types::global_dof_index>> &cells,
1598 const MPI_Comm comm)
1599 {
1600 std::vector<types::global_dof_index> cells_local(cells.size());
1601 std::vector<types::global_dof_index> cells_remote(cells.size());
1602
1603 for (unsigned int i = 0; i < cells.size(); ++i)
1604 {
1605 cells_local[i] = cells[i].first;
1606 cells_remote[i] = cells[i].second;
1607 }
1608
1609 std::vector<types::global_dof_index> cells_local_sum(cells_local.size());
1610 Utilities::MPI::sum(cells_local, comm, cells_local_sum);
1611
1612 std::vector<types::global_dof_index> cells_remote_sum(
1613 cells_remote.size());
1614 Utilities::MPI::sum(cells_remote, comm, cells_remote_sum);
1615
1616 const auto n_cells_local =
1617 std::accumulate(cells_local_sum.begin(), cells_local_sum.end(), 0);
1618 const auto n_cells_remote =
1619 std::accumulate(cells_remote_sum.begin(), cells_remote_sum.end(), 0);
1620
1621 return static_cast<double>(n_cells_local) /
1622 (n_cells_local + n_cells_remote);
1623 }
1624 } // namespace internal
1625
1626
1627
1628 template <int dim, int spacedim>
1629 std::vector<types::global_dof_index>
1631 {
1633 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1634 &tria))
1635 Assert(
1636 tr->is_multilevel_hierarchy_constructed(),
1637 ExcMessage(
1638 "We can only compute the workload imbalance if the multilevel hierarchy has been constructed!"));
1639
1640 const unsigned int n_global_levels = tria.n_global_levels();
1641
1642 std::vector<types::global_dof_index> n_cells_on_levels(n_global_levels);
1643
1644 for (unsigned int lvl = 0; lvl < n_global_levels; ++lvl)
1645 for (const auto &cell : tria.cell_iterators_on_level(lvl))
1646 if (cell->is_locally_owned_on_level())
1647 ++n_cells_on_levels[lvl];
1648
1649 return n_cells_on_levels;
1650 }
1651
1652
1653
1654 template <int dim, int spacedim>
1655 std::vector<types::global_dof_index>
1657 const std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
1658 &trias)
1659 {
1660 const unsigned int n_global_levels = trias.size();
1661
1662 std::vector<types::global_dof_index> n_cells_on_levels(n_global_levels);
1663
1664 for (unsigned int lvl = 0; lvl < n_global_levels; ++lvl)
1665 for (const auto &cell : trias[lvl]->active_cell_iterators())
1666 if (cell->is_locally_owned())
1667 ++n_cells_on_levels[lvl];
1668
1669 return n_cells_on_levels;
1670 }
1671
1672
1673
1674 template <int dim, int spacedim>
1675 double
1681
1682
1683
1684 template <int dim, int spacedim>
1685 double
1687 const std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
1688 &trias)
1689 {
1691 trias.back()->get_mpi_communicator());
1692 }
1693
1694
1695
1696 template <int dim, int spacedim>
1697 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1699 {
1700 const unsigned int n_global_levels = tria.n_global_levels();
1701
1702 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1703 cells(n_global_levels);
1704
1705 const MPI_Comm communicator = tria.get_mpi_communicator();
1706
1707 const unsigned int my_rank = Utilities::MPI::this_mpi_process(communicator);
1708
1709 for (unsigned int lvl = 0; lvl < n_global_levels - 1; ++lvl)
1710 for (const auto &cell : tria.cell_iterators_on_level(lvl))
1711 if (cell->is_locally_owned_on_level() && cell->has_children())
1712 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1713 ++i)
1714 {
1715 const auto level_subdomain_id =
1716 cell->child(i)->level_subdomain_id();
1717 if (level_subdomain_id == my_rank)
1718 ++cells[lvl + 1].first;
1719 else if (level_subdomain_id != numbers::invalid_unsigned_int)
1720 ++cells[lvl + 1].second;
1721 else
1723 }
1724
1725 return cells;
1726 }
1727
1728
1729
1730 template <int dim, int spacedim>
1731 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1733 const std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
1734 &trias)
1735 {
1736 const unsigned int n_global_levels = trias.size();
1737
1738 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1739 cells(n_global_levels);
1740
1741 const MPI_Comm communicator = trias.back()->get_mpi_communicator();
1742
1743 const unsigned int my_rank = Utilities::MPI::this_mpi_process(communicator);
1744
1745 for (unsigned int lvl = 0; lvl < n_global_levels - 1; ++lvl)
1746 {
1747 const auto &tria_coarse = *trias[lvl];
1748 const auto &tria_fine = *trias[lvl + 1];
1749
1750 const unsigned int n_coarse_cells = tria_fine.n_global_coarse_cells();
1751 const unsigned int n_global_levels = tria_fine.n_global_levels();
1752
1753 const ::internal::CellIDTranslator<dim> cell_id_translator(
1754 n_coarse_cells, n_global_levels);
1755
1756 IndexSet is_fine_owned(cell_id_translator.size());
1757 IndexSet is_fine_required(cell_id_translator.size());
1758
1759 for (const auto &cell : tria_fine.active_cell_iterators())
1760 if (!cell->is_artificial() && cell->is_locally_owned())
1761 is_fine_owned.add_index(cell_id_translator.translate(cell));
1762
1763 for (const auto &cell : tria_coarse.active_cell_iterators())
1764 if (!cell->is_artificial() && cell->is_locally_owned())
1765 {
1766 if (cell->level() + 1u == tria_fine.n_global_levels())
1767 continue;
1768
1769 for (unsigned int i = 0;
1770 i < GeometryInfo<dim>::max_children_per_cell;
1771 ++i)
1772 is_fine_required.add_index(
1773 cell_id_translator.translate(cell, i));
1774 }
1775
1776 const std::vector<unsigned int> is_fine_required_ranks =
1778 is_fine_required,
1779 communicator);
1780
1781 for (unsigned i = 0; i < is_fine_required.n_elements(); ++i)
1782 if (is_fine_required_ranks[i] == my_rank)
1783 ++cells[lvl + 1].first;
1784 else if (is_fine_required_ranks[i] != numbers::invalid_unsigned_int)
1785 ++cells[lvl + 1].second;
1786 }
1787
1788 return cells;
1789 }
1790
1791
1792
1793 template <int dim, int spacedim>
1794 double
1800
1801
1802
1803 template <int dim, int spacedim>
1804 double
1806 const std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
1807 &trias)
1808 {
1811 trias.back()->get_mpi_communicator());
1812 }
1813
1814} // namespace MGTools
1815
1816
1817// explicit instantiations
1818#include "multigrid/mg_tools.inst"
1819
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:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#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:929
std::size_t size
Definition mpi.cc:745
const MPI_Comm comm
Definition mpi.cc:924
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:1565
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:1595
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:607
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:1056
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:1698
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:1241
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:1015
std::vector< types::global_dof_index > local_workload(const Triangulation< dim, spacedim > &tria)
Definition mg_tools.cc:1630
void make_flux_sparsity_pattern_edge(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity, const unsigned int level)
Definition mg_tools.cc:681
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:1445
double workload_imbalance(const Triangulation< dim, spacedim > &tria)
Definition mg_tools.cc:1676
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:1158
double vertical_communication_efficiency(const Triangulation< dim, spacedim > &tria)
Definition mg_tools.cc:1795
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:1532
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:99
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:1831
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:114
constexpr types::global_dof_index invalid_dof_index
Definition types.h:263
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()