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