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