deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
dof_tools_constraints.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2013 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#include <deal.II/base/table.h>
14#include <deal.II/base/template_constraints.h>
15#include <deal.II/base/utilities.h>
16#include <deal.II/base/work_stream.h>
17
18#include <deal.II/dofs/dof_accessor.h>
19#include <deal.II/dofs/dof_handler.h>
20#include <deal.II/dofs/dof_tools.h>
21
22#include <deal.II/fe/fe.h>
23#include <deal.II/fe/fe_tools.h>
24#include <deal.II/fe/fe_values.h>
25
26#include <deal.II/grid/filtered_iterator.h>
27#include <deal.II/grid/grid_tools.h>
28#include <deal.II/grid/intergrid_map.h>
29#include <deal.II/grid/tria.h>
30#include <deal.II/grid/tria_iterator.h>
31
32#include <deal.II/hp/fe_collection.h>
33#include <deal.II/hp/fe_values.h>
34
35#include <deal.II/lac/affine_constraints.h>
36#include <deal.II/lac/vector.h>
37
38#ifdef DEAL_II_WITH_MPI
39# include <deal.II/lac/la_parallel_vector.h>
40#endif
41
42#include <algorithm>
43#include <array>
44#include <complex>
45#include <memory>
46#include <numeric>
47
48DEAL_II_NAMESPACE_OPEN
49
50
51
52namespace DoFTools
53{
54 namespace internal
55 {
56 namespace
57 {
58 bool
59 check_primary_dof_list(
60 const FullMatrix<double> &face_interpolation_matrix,
61 const std::vector<types::global_dof_index> &primary_dof_list)
62 {
63 const unsigned int N = primary_dof_list.size();
64
65 FullMatrix<double> tmp(N, N);
66 for (unsigned int i = 0; i < N; ++i)
67 for (unsigned int j = 0; j < N; ++j)
68 tmp(i, j) = face_interpolation_matrix(primary_dof_list[i], j);
69
70 // then use the algorithm from FullMatrix::gauss_jordan on this matrix
71 // to find out whether it is singular. the algorithm there does pivoting
72 // and at the end swaps rows back into their proper order -- we omit
73 // this step here, since we don't care about the inverse matrix, all we
74 // care about is whether the matrix is regular or singular
75
76 // first get an estimate of the size of the elements of this matrix, for
77 // later checks whether the pivot element is large enough, or whether we
78 // have to fear that the matrix is not regular
79 double diagonal_sum = 0;
80 for (unsigned int i = 0; i < N; ++i)
81 diagonal_sum += std::fabs(tmp(i, i));
82 const double typical_diagonal_element = diagonal_sum / N;
83
84 // initialize the array that holds the permutations that we find during
85 // pivot search
86 std::vector<unsigned int> p(N);
87 for (unsigned int i = 0; i < N; ++i)
88 p[i] = i;
89
90 for (unsigned int j = 0; j < N; ++j)
91 {
92 // pivot search: search that part of the line on and right of the
93 // diagonal for the largest element
94 double max = std::fabs(tmp(j, j));
95 unsigned int r = j;
96 for (unsigned int i = j + 1; i < N; ++i)
97 {
98 if (std::fabs(tmp(i, j)) > max)
99 {
100 max = std::fabs(tmp(i, j));
101 r = i;
102 }
103 }
104 // check whether the pivot is too small. if that is the case, then
105 // the matrix is singular and we shouldn't use this set of primary
106 // dofs
107 if (max < 1.e-12 * typical_diagonal_element)
108 return false;
109
110 // row interchange
111 if (r > j)
112 {
113 for (unsigned int k = 0; k < N; ++k)
114 std::swap(tmp(j, k), tmp(r, k));
115
116 std::swap(p[j], p[r]);
117 }
118
119 // transformation
120 const double hr = 1. / tmp(j, j);
121 tmp(j, j) = hr;
122 for (unsigned int k = 0; k < N; ++k)
123 {
124 if (k == j)
125 continue;
126 for (unsigned int i = 0; i < N; ++i)
127 {
128 if (i == j)
129 continue;
130 tmp(i, k) -= tmp(i, j) * tmp(j, k) * hr;
131 }
132 }
133 for (unsigned int i = 0; i < N; ++i)
134 {
135 tmp(i, j) *= hr;
136 tmp(j, i) *= -hr;
137 }
138 tmp(j, j) = hr;
139 }
140
141 // everything went fine, so we can accept this set of primary dofs (at
142 // least as far as they have already been collected)
143 return true;
144 }
145
146
147
169 template <int dim, int spacedim>
170 void
171 select_primary_dofs_for_face_restriction(
172 const FiniteElement<dim, spacedim> &fe1,
173 const FiniteElement<dim, spacedim> &fe2,
174 const FullMatrix<double> &face_interpolation_matrix,
175 std::vector<bool> &primary_dof_mask)
176 {
177 // TODO: the implementation makes the assumption that all faces have the
178 // same number of dofs
179 AssertDimension(fe1.n_unique_faces(), 1);
180 AssertDimension(fe2.n_unique_faces(), 1);
181 const unsigned int face_no = 0;
182 (void)face_no;
183
184 Assert(fe1.n_dofs_per_face(face_no) >= fe2.n_dofs_per_face(face_no),
185 ExcInternalError());
186 AssertDimension(primary_dof_mask.size(), fe1.n_dofs_per_face(face_no));
187
188 Assert(fe2.n_dofs_per_vertex() <= fe1.n_dofs_per_vertex(),
189 ExcInternalError());
190 Assert(fe2.n_dofs_per_line() <= fe1.n_dofs_per_line(),
191 ExcInternalError());
192 Assert((dim < 3) ||
193 (fe2.n_dofs_per_quad(face_no) <= fe1.n_dofs_per_quad(face_no)),
194 ExcInternalError());
195
196 // the idea here is to designate as many DoFs in fe1 per object (vertex,
197 // line, quad) as primary as there are such dofs in fe2 (indices are
198 // int, because we want to avoid the 'unsigned int < 0 is always false
199 // warning for the cases at the bottom in 1d and 2d)
200 //
201 // as mentioned in the paper, it is not always easy to find a set of
202 // primary dofs that produces an invertible matrix. to this end, we
203 // check in each step whether the matrix is still invertible and simply
204 // discard this dof if the matrix is not invertible anymore.
205 //
206 // the cases where we did have trouble in the past were with adding more
207 // quad dofs when Q3 and Q4 elements meet at a refined face in 3d (see
208 // the hp/crash_12 test that tests that we can do exactly this, and
209 // failed before we had code to compensate for this case). the other
210 // case are system elements: if we have say a Q1Q2 vs a Q2Q3 element,
211 // then we can't just take all primary dofs on a line from a single base
212 // element, since the shape functions of that base element are
213 // independent of that of the other one. this latter case shows up when
214 // running hp/hp_constraints_q_system_06
215
216 std::vector<types::global_dof_index> primary_dof_list;
217 unsigned int index = 0;
218 for (int v = 0;
219 v < static_cast<signed int>(GeometryInfo<dim>::vertices_per_face);
220 ++v)
221 {
222 unsigned int dofs_added = 0;
223 unsigned int i = 0;
224 while (dofs_added < fe2.n_dofs_per_vertex())
225 {
226 // make sure that we were able to find a set of primary dofs and
227 // that the code down below didn't just reject all our efforts
228 Assert(i < fe1.n_dofs_per_vertex(), ExcInternalError());
229
230 // tentatively push this vertex dof
231 primary_dof_list.push_back(index + i);
232
233 // then see what happens. if it succeeds, fine
234 if (check_primary_dof_list(face_interpolation_matrix,
235 primary_dof_list) == true)
236 ++dofs_added;
237 else
238 // well, it didn't. simply pop that dof from the list again
239 // and try with the next dof
240 primary_dof_list.pop_back();
241
242 // forward counter by one
243 ++i;
244 }
245 index += fe1.n_dofs_per_vertex();
246 }
247
248 for (int l = 0;
249 l < static_cast<signed int>(GeometryInfo<dim>::lines_per_face);
250 ++l)
251 {
252 // same algorithm as above
253 unsigned int dofs_added = 0;
254 unsigned int i = 0;
255 while (dofs_added < fe2.n_dofs_per_line())
256 {
257 Assert(i < fe1.n_dofs_per_line(), ExcInternalError());
258
259 primary_dof_list.push_back(index + i);
260 if (check_primary_dof_list(face_interpolation_matrix,
261 primary_dof_list) == true)
262 ++dofs_added;
263 else
264 primary_dof_list.pop_back();
265
266 ++i;
267 }
268 index += fe1.n_dofs_per_line();
269 }
270
271 for (int q = 0;
272 q < static_cast<signed int>(GeometryInfo<dim>::quads_per_face);
273 ++q)
274 {
275 // same algorithm as above
276 unsigned int dofs_added = 0;
277 unsigned int i = 0;
278 while (dofs_added < fe2.n_dofs_per_quad(q))
279 {
280 Assert(i < fe1.n_dofs_per_quad(q), ExcInternalError());
281
282 primary_dof_list.push_back(index + i);
283 if (check_primary_dof_list(face_interpolation_matrix,
284 primary_dof_list) == true)
285 ++dofs_added;
286 else
287 primary_dof_list.pop_back();
288
289 ++i;
290 }
291 index += fe1.n_dofs_per_quad(q);
292 }
293
294 AssertDimension(index, fe1.n_dofs_per_face(face_no));
295 AssertDimension(primary_dof_list.size(), fe2.n_dofs_per_face(face_no));
296
297 // finally copy the list into the mask
298 std::fill(primary_dof_mask.begin(), primary_dof_mask.end(), false);
299 for (const auto dof : primary_dof_list)
300 primary_dof_mask[dof] = true;
301 }
302
303
304
309 template <int dim, int spacedim>
310 void
311 ensure_existence_of_primary_dof_mask(
312 const FiniteElement<dim, spacedim> &fe1,
313 const FiniteElement<dim, spacedim> &fe2,
314 const FullMatrix<double> &face_interpolation_matrix,
315 std::unique_ptr<std::vector<bool>> &primary_dof_mask)
316 {
317 // TODO: the implementation makes the assumption that all faces have the
318 // same number of dofs
319 AssertDimension(fe1.n_unique_faces(), 1);
320 AssertDimension(fe2.n_unique_faces(), 1);
321 const unsigned int face_no = 0;
322
323 if (primary_dof_mask == nullptr)
324 {
325 primary_dof_mask =
326 std::make_unique<std::vector<bool>>(fe1.n_dofs_per_face(face_no));
327 select_primary_dofs_for_face_restriction(fe1,
328 fe2,
329 face_interpolation_matrix,
330 *primary_dof_mask);
331 }
332 }
333
334
335
341 template <int dim, int spacedim>
342 void
343 ensure_existence_of_face_matrix(
344 const FiniteElement<dim, spacedim> &fe1,
345 const FiniteElement<dim, spacedim> &fe2,
346 std::unique_ptr<FullMatrix<double>> &matrix)
347 {
348 // TODO: the implementation makes the assumption that all faces have the
349 // same number of dofs
350 AssertDimension(fe1.n_unique_faces(), 1);
351 AssertDimension(fe2.n_unique_faces(), 1);
352 const unsigned int face_no = 0;
353
354 if (matrix == nullptr)
355 {
356 matrix = std::make_unique<FullMatrix<double>>(
357 fe2.n_dofs_per_face(face_no), fe1.n_dofs_per_face(face_no));
358 fe1.get_face_interpolation_matrix(fe2, *matrix, face_no);
359 }
360 }
361
362
363
367 template <int dim, int spacedim>
368 void
369 ensure_existence_of_subface_matrix(
370 const FiniteElement<dim, spacedim> &fe1,
371 const FiniteElement<dim, spacedim> &fe2,
372 const unsigned int subface,
373 std::unique_ptr<FullMatrix<double>> &matrix)
374 {
375 // TODO: the implementation makes the assumption that all faces have the
376 // same number of dofs
377 AssertDimension(fe1.n_unique_faces(), 1);
378 AssertDimension(fe2.n_unique_faces(), 1);
379 const unsigned int face_no = 0;
380
381 if (matrix == nullptr)
382 {
383 matrix = std::make_unique<FullMatrix<double>>(
384 fe2.n_dofs_per_face(face_no), fe1.n_dofs_per_face(face_no));
385 fe1.get_subface_interpolation_matrix(fe2,
386 subface,
387 *matrix,
388 face_no);
389 }
390 }
391
392
393
399 void
400 ensure_existence_of_split_face_matrix(
401 const FullMatrix<double> &face_interpolation_matrix,
402 const std::vector<bool> &primary_dof_mask,
403 std::unique_ptr<std::pair<FullMatrix<double>, FullMatrix<double>>>
404 &split_matrix)
405 {
406 AssertDimension(primary_dof_mask.size(), face_interpolation_matrix.m());
407 Assert(std::count(primary_dof_mask.begin(),
408 primary_dof_mask.end(),
409 true) ==
410 static_cast<signed int>(face_interpolation_matrix.n()),
411 ExcInternalError());
412
413 if (split_matrix == nullptr)
414 {
415 split_matrix = std::make_unique<
416 std::pair<FullMatrix<double>, FullMatrix<double>>>();
417
418 const unsigned int n_primary_dofs = face_interpolation_matrix.n();
419 const unsigned int n_dofs = face_interpolation_matrix.m();
420
421 Assert(n_primary_dofs <= n_dofs, ExcInternalError());
422
423 // copy and invert the primary component, copy the dependent
424 // component
425 split_matrix->first.reinit(n_primary_dofs, n_primary_dofs);
426 split_matrix->second.reinit(n_dofs - n_primary_dofs,
427 n_primary_dofs);
428
429 unsigned int nth_primary_dof = 0, nth_dependent_dof = 0;
430
431 for (unsigned int i = 0; i < n_dofs; ++i)
432 if (primary_dof_mask[i] == true)
433 {
434 for (unsigned int j = 0; j < n_primary_dofs; ++j)
435 split_matrix->first(nth_primary_dof, j) =
436 face_interpolation_matrix(i, j);
437 ++nth_primary_dof;
438 }
439 else
440 {
441 for (unsigned int j = 0; j < n_primary_dofs; ++j)
442 split_matrix->second(nth_dependent_dof, j) =
443 face_interpolation_matrix(i, j);
444 ++nth_dependent_dof;
445 }
446
447 AssertDimension(nth_primary_dof, n_primary_dofs);
448 AssertDimension(nth_dependent_dof, n_dofs - n_primary_dofs);
449
450 // TODO[WB]: We should make sure very small entries are removed
451 // after inversion
452 split_matrix->first.gauss_jordan();
453 }
454 }
455
456
462 template <int dim, int spacedim>
463 unsigned int
464 n_finite_elements(const DoFHandler<dim, spacedim> &dof_handler)
465 {
466 if (dof_handler.has_hp_capabilities() == true)
467 return dof_handler.get_fe_collection().size();
468 else
469 return 1;
470 }
471
472
473
484 template <typename number1, typename number2>
485 void
486 filter_constraints(
487 const std::vector<types::global_dof_index> &primary_dofs,
488 const std::vector<types::global_dof_index> &dependent_dofs,
489 const FullMatrix<number1> &face_constraints,
490 AffineConstraints<number2> &constraints)
491 {
492 Assert(face_constraints.n() == primary_dofs.size(),
493 ExcDimensionMismatch(primary_dofs.size(), face_constraints.n()));
494 Assert(face_constraints.m() == dependent_dofs.size(),
495 ExcDimensionMismatch(dependent_dofs.size(),
496 face_constraints.m()));
497
498 const unsigned int n_primary_dofs = primary_dofs.size();
499 const unsigned int n_dependent_dofs = dependent_dofs.size();
500
501 // check for a couple conditions that happened in parallel distributed
502 // mode
503 for (unsigned int row = 0; row != n_dependent_dofs; ++row)
504 Assert(dependent_dofs[row] != numbers::invalid_dof_index,
505 ExcInternalError());
506 for (unsigned int col = 0; col != n_primary_dofs; ++col)
507 Assert(primary_dofs[col] != numbers::invalid_dof_index,
508 ExcInternalError());
509
510 // Build constraints in a vector of pairs that can be
511 // arbitrarily large, but that holds up to 25 elements without
512 // external memory allocation. This is good enough for hanging
513 // node constraints of Q4 elements in 3d, so covers most
514 // common cases. Sort the primary dofs to add a sorted list to the
515 // affine constraints, which increases performance there.
516 using size_type = typename AffineConstraints<number2>::size_type;
517 boost::container::small_vector<std::pair<size_type, size_type>, 25>
518 sorted_primary_dofs;
519 sorted_primary_dofs.reserve(n_primary_dofs);
520 for (unsigned int i = 0; i < n_primary_dofs; ++i)
521 sorted_primary_dofs.emplace_back(primary_dofs[i], i);
522 std::sort(sorted_primary_dofs.begin(), sorted_primary_dofs.end());
523
524 boost::container::small_vector<std::pair<size_type, number2>, 25>
525 entries;
526 entries.reserve(n_primary_dofs);
527 for (unsigned int row = 0; row != n_dependent_dofs; ++row)
528 if (constraints.is_constrained(dependent_dofs[row]) == false)
529 {
530 // Check if we have an identity constraint, i.e.,
531 // something of the form
532 // U(dependent_dof[row])==U(primary_dof[row]),
533 // where
534 // dependent_dof[row] == primary_dof[row].
535 // This can happen in the hp context where we have previously
536 // unified DoF indices, for example, the middle node on the
537 // face of a Q4 element will have gotten the same index
538 // as the middle node of the Q2 element on the neighbor
539 // cell. But because the other Q4 nodes will still have to be
540 // constrained, so the middle node shows up again here.
541 //
542 // If we find such a constraint, then it is trivially
543 // satisfied, and we can move on to the next dependent
544 // DoF (row). The only thing we should make sure is that the
545 // row of the matrix really just contains this one entry.
546 {
547 bool is_trivial_constraint = false;
548
549 for (unsigned int i = 0; i < n_primary_dofs; ++i)
550 if (face_constraints(row, i) == 1.0)
551 if (dependent_dofs[row] == primary_dofs[i])
552 {
553 is_trivial_constraint = true;
554
555 for (unsigned int ii = 0; ii < n_primary_dofs; ++ii)
556 if (ii != i)
557 Assert(face_constraints(row, ii) == 0.0,
558 ExcInternalError());
559
560 break;
561 }
562
563 if (is_trivial_constraint == true)
564 continue;
565 }
566
567 // then enter those constraints that are larger than
568 // 1e-14; since numbers are normalized for the subface
569 // interpolation matrices, we do not need to normalize here.
570 // everything else probably originated from
571 // inexact inversion of matrices and similar effects. having
572 // those constraints in here will only lead to problems because
573 // it makes sparsity patterns fuller than necessary without
574 // producing any significant effect. do this in two steps, first
575 // filling a vector and then adding to the constraints in order
576 // to reduce the number of memory allocations.
577 entries.clear();
578 for (const auto &[dof_index, unsorted_index] :
579 sorted_primary_dofs)
580 if (std::fabs(face_constraints(row, unsorted_index)) >= 1e-14)
581 entries.emplace_back(dof_index,
582 face_constraints(row, unsorted_index));
583 constraints.add_constraint(dependent_dofs[row],
584 entries,
585 /* inhomogeneity= */ 0.);
586 }
587 }
588
589 } // namespace
590
591
592
593 template <typename number, int spacedim>
594 void
595 make_hp_hanging_node_constraints(
596 const DoFHandler<1, spacedim> & /*dof_handler*/,
597 AffineConstraints<number> & /*constraints*/)
598 {
599 // nothing to do for dof handlers in 1d
600 }
601
602
603
604 template <typename number, int spacedim>
605 void
606 make_hanging_node_constraints_nedelec(
607 const ::DoFHandler<1, spacedim> & /*dof_handler*/,
608 AffineConstraints<number> & /*constraints*/,
609 std::integral_constant<int, 1>)
610 {
611 // nothing to do for dof handlers in 1d
612 }
613
614
615
616 template <typename number, int spacedim>
617 void
618 make_oldstyle_hanging_node_constraints(
619 const DoFHandler<1, spacedim> & /*dof_handler*/,
620 AffineConstraints<number> & /*constraints*/,
621 std::integral_constant<int, 1>)
622 {
623 // nothing to do for dof handlers in 1d
624 }
625
626
627
628 template <int dim_, int spacedim, typename number>
629 void
630 make_oldstyle_hanging_node_constraints(
631 const DoFHandler<dim_, spacedim> &dof_handler,
632 AffineConstraints<number> &constraints,
633 std::integral_constant<int, 2>)
634 {
635 const unsigned int dim = 2;
636
637 std::vector<types::global_dof_index> dofs_on_mother;
638 std::vector<types::global_dof_index> dofs_on_children;
639
640 // Build constraints in a vector of pairs that can be
641 // arbitrarily large, but that holds up to 25 elements without
642 // external memory allocation. This is good enough for hanging
643 // node constraints of Q4 elements in 3d, so covers most
644 // common cases.
645 boost::container::small_vector<
646 std::pair<typename AffineConstraints<number>::size_type, number>,
647 25>
648 constraint_entries;
649
650 // loop over all lines; only on lines there can be constraints. We do so
651 // by looping over all active cells and checking whether any of the faces
652 // are refined which can only be from the neighboring cell because this
653 // one is active. In that case, the face is subject to constraints
654 //
655 // note that even though we may visit a face twice if the neighboring
656 // cells are equally refined, we can only visit each face with hanging
657 // nodes once
658 for (const auto &cell : dof_handler.active_cell_iterators())
659 {
660 // artificial cells can at best neighbor ghost cells, but we're not
661 // interested in these interfaces
662 if (cell->is_artificial())
663 continue;
664
665 for (const unsigned int face : cell->face_indices())
666 if (cell->face(face)->has_children())
667 {
668 // in any case, faces can have at most two active FE indices,
669 // but here the face can have only one (namely the same as that
670 // from the cell we're sitting on), and each of the children can
671 // have only one as well. check this
672 Assert(cell->face(face)->n_active_fe_indices() == 1,
673 ExcInternalError());
674 Assert(cell->face(face)->fe_index_is_active(
675 cell->active_fe_index()) == true,
676 ExcInternalError());
677 for (unsigned int c = 0; c < cell->face(face)->n_children();
678 ++c)
679 if (!cell->neighbor_child_on_subface(face, c)
680 ->is_artificial())
681 Assert(cell->face(face)->child(c)->n_active_fe_indices() ==
682 1,
683 ExcInternalError());
684
685 // right now, all that is implemented is the case that both
686 // sides use the same FE
687 for (unsigned int c = 0; c < cell->face(face)->n_children();
688 ++c)
689 if (!cell->neighbor_child_on_subface(face, c)
690 ->is_artificial())
691 Assert(cell->face(face)->child(c)->fe_index_is_active(
692 cell->active_fe_index()) == true,
693 ExcNotImplemented());
694
695 // ok, start up the work
696 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
697 const types::fe_index fe_index = cell->active_fe_index();
698
699 const unsigned int n_dofs_on_mother =
700 2 * fe.n_dofs_per_vertex() +
701 fe.n_dofs_per_line(),
702 n_dofs_on_children =
703 fe.n_dofs_per_vertex() +
704 2 * fe.n_dofs_per_line();
705
706 dofs_on_mother.resize(n_dofs_on_mother);
707 // we might not use all of those in case of artificial cells, so
708 // do not resize(), but reserve() and use push_back later.
709 dofs_on_children.clear();
710 dofs_on_children.reserve(n_dofs_on_children);
711
712 Assert(n_dofs_on_mother == fe.constraints().n(),
713 ExcDimensionMismatch(n_dofs_on_mother,
714 fe.constraints().n()));
715 Assert(n_dofs_on_children == fe.constraints().m(),
716 ExcDimensionMismatch(n_dofs_on_children,
717 fe.constraints().m()));
718
719 const typename DoFHandler<dim_, spacedim>::line_iterator
720 this_face = cell->face(face);
721
722 // fill the dofs indices. Use same enumeration scheme as in
723 // @p{FiniteElement::constraints()}
724 unsigned int next_index = 0;
725 for (unsigned int vertex = 0; vertex < 2; ++vertex)
726 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
727 ++dof)
728 dofs_on_mother[next_index++] =
729 this_face->vertex_dof_index(vertex, dof, fe_index);
730 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof)
731 dofs_on_mother[next_index++] =
732 this_face->dof_index(dof, fe_index);
733 AssertDimension(next_index, dofs_on_mother.size());
734
735 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); ++dof)
736 dofs_on_children.push_back(
737 this_face->child(0)->vertex_dof_index(1, dof, fe_index));
738 for (unsigned int child = 0; child < 2; ++child)
739 {
740 // skip artificial cells
741 if (cell->neighbor_child_on_subface(face, child)
742 ->is_artificial())
743 continue;
744 for (unsigned int dof = 0; dof != fe.n_dofs_per_line();
745 ++dof)
746 dofs_on_children.push_back(
747 this_face->child(child)->dof_index(dof, fe_index));
748 }
749 // note: can get fewer DoFs when we have artificial cells
750 Assert(dofs_on_children.size() <= n_dofs_on_children,
751 ExcInternalError());
752
753 // for each row in the AffineConstraints object for this line:
754 for (unsigned int row = 0; row != dofs_on_children.size();
755 ++row)
756 {
757 constraint_entries.clear();
758 constraint_entries.reserve(dofs_on_mother.size());
759 for (unsigned int i = 0; i != dofs_on_mother.size(); ++i)
760 constraint_entries.emplace_back(dofs_on_mother[i],
761 fe.constraints()(row, i));
762
763 constraints.add_constraint(dofs_on_children[row],
764 constraint_entries,
765 0.);
766 }
767 }
768 else
769 {
770 // this face has no children, but it could still be that it is
771 // shared by two cells that use a different FE index. check a
772 // couple of things, but ignore the case that the neighbor is an
773 // artificial cell
774 if (!cell->at_boundary(face) &&
775 !cell->neighbor(face)->is_artificial())
776 {
777 Assert(cell->face(face)->n_active_fe_indices() == 1,
778 ExcNotImplemented());
779 Assert(cell->face(face)->fe_index_is_active(
780 cell->active_fe_index()) == true,
781 ExcInternalError());
782 }
783 }
784 }
785 }
786
787
788
789 template <int dim_, int spacedim, typename number>
790 void
791 make_oldstyle_hanging_node_constraints(
792 const DoFHandler<dim_, spacedim> &dof_handler,
793 AffineConstraints<number> &constraints,
794 std::integral_constant<int, 3>)
795 {
796 const unsigned int dim = 3;
797
798 std::vector<types::global_dof_index> dofs_on_mother;
799 std::vector<types::global_dof_index> dofs_on_children;
800
801 // Build constraints in a vector of pairs that can be
802 // arbitrarily large, but that holds up to 25 elements without
803 // external memory allocation. This is good enough for hanging
804 // node constraints of Q4 elements in 3d, so covers most
805 // common cases.
806 boost::container::small_vector<
807 std::pair<typename AffineConstraints<number>::size_type, number>,
808 25>
809 constraint_entries;
810
811 // loop over all quads; only on quads there can be constraints. We do so
812 // by looping over all active cells and checking whether any of the faces
813 // are refined which can only be from the neighboring cell because this
814 // one is active. In that case, the face is subject to constraints
815 //
816 // note that even though we may visit a face twice if the neighboring
817 // cells are equally refined, we can only visit each face with hanging
818 // nodes once
819 for (const auto &cell : dof_handler.active_cell_iterators())
820 {
821 // artificial cells can at best neighbor ghost cells, but we're not
822 // interested in these interfaces
823 if (cell->is_artificial())
824 continue;
825
826 for (const unsigned int face : cell->face_indices())
827 if (cell->face(face)->has_children())
828 {
829 // first of all, make sure that we treat a case which is
830 // possible, i.e. either no dofs on the face at all or no
831 // anisotropic refinement
832 if (cell->get_fe().n_dofs_per_face(face) == 0)
833 continue;
834
835 Assert(cell->face(face)->refinement_case() ==
836 RefinementCase<dim - 1>::isotropic_refinement,
837 ExcNotImplemented());
838
839 // in any case, faces can have at most two active FE indices,
840 // but here the face can have only one (namely the same as that
841 // from the cell we're sitting on), and each of the children can
842 // have only one as well. check this
843 AssertDimension(cell->face(face)->n_active_fe_indices(), 1);
844 Assert(cell->face(face)->fe_index_is_active(
845 cell->active_fe_index()) == true,
846 ExcInternalError());
847 for (unsigned int c = 0; c < cell->face(face)->n_children();
848 ++c)
849 if (!cell->neighbor_child_on_subface(face, c)
850 ->is_artificial())
851 AssertDimension(
852 cell->face(face)->child(c)->n_active_fe_indices(), 1);
853
854 // right now, all that is implemented is the case that both
855 // sides use the same fe, and not only that but also that all
856 // lines bounding this face and the children have the same FE
857 for (unsigned int c = 0; c < cell->face(face)->n_children();
858 ++c)
859 if (!cell->neighbor_child_on_subface(face, c)
860 ->is_artificial())
861 {
862 Assert(cell->face(face)->child(c)->fe_index_is_active(
863 cell->active_fe_index()) == true,
864 ExcNotImplemented());
865 for (unsigned int e = 0; e < 4; ++e)
866 {
867 Assert(cell->face(face)
868 ->child(c)
869 ->line(e)
870 ->n_active_fe_indices() == 1,
871 ExcNotImplemented());
872 Assert(cell->face(face)
873 ->child(c)
874 ->line(e)
875 ->fe_index_is_active(
876 cell->active_fe_index()) == true,
877 ExcNotImplemented());
878 }
879 }
880 for (unsigned int e = 0; e < 4; ++e)
881 {
882 Assert(cell->face(face)->line(e)->n_active_fe_indices() ==
883 1,
884 ExcNotImplemented());
885 Assert(cell->face(face)->line(e)->fe_index_is_active(
886 cell->active_fe_index()) == true,
887 ExcNotImplemented());
888 }
889
890 // ok, start up the work
891 const FiniteElement<dim> &fe = cell->get_fe();
892 const types::fe_index fe_index = cell->active_fe_index();
893
894 const unsigned int n_dofs_on_mother = fe.n_dofs_per_face(face);
895 const unsigned int n_dofs_on_children =
896 (5 * fe.n_dofs_per_vertex() + 12 * fe.n_dofs_per_line() +
897 4 * fe.n_dofs_per_quad(face));
898
899 // TODO[TL]: think about this and the following in case of
900 // anisotropic refinement
901
902 dofs_on_mother.resize(n_dofs_on_mother);
903 // we might not use all of those in case of artificial cells, so
904 // do not resize(), but reserve() and use push_back later.
905 dofs_on_children.clear();
906 dofs_on_children.reserve(n_dofs_on_children);
907
908 Assert(n_dofs_on_mother == fe.constraints().n(),
909 ExcDimensionMismatch(n_dofs_on_mother,
910 fe.constraints().n()));
911 Assert(n_dofs_on_children == fe.constraints().m(),
912 ExcDimensionMismatch(n_dofs_on_children,
913 fe.constraints().m()));
914
915 const typename DoFHandler<dim_, spacedim>::face_iterator
916 this_face = cell->face(face);
917
918 // fill the dofs indices. Use same enumeration scheme as in
919 // @p{FiniteElement::constraints()}
920 unsigned int next_index = 0;
921 for (unsigned int vertex = 0; vertex < 4; ++vertex)
922 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
923 ++dof)
924 dofs_on_mother[next_index++] =
925 this_face->vertex_dof_index(vertex, dof, fe_index);
926 for (unsigned int line = 0; line < 4; ++line)
927 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof)
928 dofs_on_mother[next_index++] =
929 this_face->line(line)->dof_index(dof, fe_index);
930 for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face);
931 ++dof)
932 dofs_on_mother[next_index++] =
933 this_face->dof_index(dof, fe_index);
934 AssertDimension(next_index, dofs_on_mother.size());
935
936 // TODO: assert some consistency assumptions
937
938 // TODO[TL]: think about this in case of anisotropic refinement
939
940 Assert(dof_handler.get_triangulation()
941 .get_anisotropic_refinement_flag() ||
942 ((this_face->child(0)->vertex_index(3) ==
943 this_face->child(1)->vertex_index(2)) &&
944 (this_face->child(0)->vertex_index(3) ==
945 this_face->child(2)->vertex_index(1)) &&
946 (this_face->child(0)->vertex_index(3) ==
947 this_face->child(3)->vertex_index(0))),
948 ExcInternalError());
949
950 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); ++dof)
951 dofs_on_children.push_back(
952 this_face->child(0)->vertex_dof_index(3, dof));
953
954 // dof numbers on the centers of the lines bounding this face
955 for (unsigned int line = 0; line < 4; ++line)
956 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
957 ++dof)
958 dofs_on_children.push_back(
959 this_face->line(line)->child(0)->vertex_dof_index(
960 1, dof, fe_index));
961
962 // next the dofs on the lines interior to the face; the order of
963 // these lines is laid down in the FiniteElement class
964 // documentation
965 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
966 dofs_on_children.push_back(
967 this_face->child(0)->line(1)->dof_index(dof, fe_index));
968 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
969 dofs_on_children.push_back(
970 this_face->child(2)->line(1)->dof_index(dof, fe_index));
971 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
972 dofs_on_children.push_back(
973 this_face->child(0)->line(3)->dof_index(dof, fe_index));
974 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
975 dofs_on_children.push_back(
976 this_face->child(1)->line(3)->dof_index(dof, fe_index));
977
978 // dofs on the bordering lines
979 for (unsigned int line = 0; line < 4; ++line)
980 for (unsigned int child = 0; child < 2; ++child)
981 {
982 for (unsigned int dof = 0; dof != fe.n_dofs_per_line();
983 ++dof)
984 dofs_on_children.push_back(
985 this_face->line(line)->child(child)->dof_index(
986 dof, fe_index));
987 }
988
989 // finally, for the dofs interior to the four child faces
990 for (unsigned int child = 0; child < 4; ++child)
991 {
992 // skip artificial cells
993 if (cell->neighbor_child_on_subface(face, child)
994 ->is_artificial())
995 continue;
996 for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face);
997 ++dof)
998 dofs_on_children.push_back(
999 this_face->child(child)->dof_index(dof, fe_index));
1000 }
1001
1002 // note: can get fewer DoFs when we have artificial cells:
1003 Assert(dofs_on_children.size() <= n_dofs_on_children,
1004 ExcInternalError());
1005
1006 // For each row in the AffineConstraints object for
1007 // this line, add the constraint. Ignore rows that
1008 // have already been added (e.g., in 3d degrees of
1009 // freedom on edges with hanging nodes will be visited
1010 // more than once).
1011 for (unsigned int row = 0; row != dofs_on_children.size();
1012 ++row)
1013 if (constraints.is_constrained(dofs_on_children[row]) ==
1014 false)
1015 {
1016 constraint_entries.clear();
1017 constraint_entries.reserve(dofs_on_mother.size());
1018 for (unsigned int i = 0; i != dofs_on_mother.size(); ++i)
1019 constraint_entries.emplace_back(dofs_on_mother[i],
1020 fe.constraints()(row,
1021 i));
1022
1023 constraints.add_constraint(dofs_on_children[row],
1024 constraint_entries,
1025 0.);
1026 }
1027 }
1028 else
1029 {
1030 // this face has no children, but it could still be that it is
1031 // shared by two cells that use a different FE index. check a
1032 // couple of things, but ignore the case that the neighbor is an
1033 // artificial cell
1034 if (!cell->at_boundary(face) &&
1035 !cell->neighbor(face)->is_artificial())
1036 {
1037 Assert(cell->face(face)->n_active_fe_indices() == 1,
1038 ExcNotImplemented());
1039 Assert(cell->face(face)->fe_index_is_active(
1040 cell->active_fe_index()) == true,
1041 ExcInternalError());
1042 }
1043 }
1044 }
1045 }
1046
1047
1048
1049 template <int dim_, int spacedim, typename number>
1050 void
1051 make_hanging_node_constraints_nedelec(
1052 const DoFHandler<dim_, spacedim> &dof_handler,
1053 AffineConstraints<number> &constraints,
1054 std::integral_constant<int, 2>)
1055 {
1056 // Parts of this function are very similar to
1057 // make_oldstyle_hanging_node_constraints.
1058 // Therefore, only the parts that differ from the
1059 // make_oldstyle_hanging_node_constraints are commented on.
1060
1061 const unsigned int dim = 2;
1062
1063 std::vector<types::global_dof_index> face_dof_indices;
1064 std::map<types::global_dof_index, std::set<types::global_dof_index>>
1065 depends_on;
1066
1067 // loop over all lines
1068 for (const auto &cell : dof_handler.active_cell_iterators())
1069 {
1070 // skip artificial cells
1071 if (cell->is_artificial())
1072 continue;
1073
1074 // loop over all faces:
1075 for (const unsigned int f : cell->face_indices())
1076 {
1077 // check if the neighbor is refined; if so, we need to
1078 // treat the constraints on this interface
1079 if (!cell->face(f)->has_children())
1080 continue;
1081
1082 Assert(cell->face(f)->n_active_fe_indices() == 1,
1083 ExcInternalError());
1084 Assert(cell->face(f)->fe_index_is_active(
1085 cell->active_fe_index()) == true,
1086 ExcInternalError());
1087
1088 if constexpr (running_in_debug_mode())
1089 {
1090 for (unsigned int c = 0; c < cell->face(f)->n_children(); ++c)
1091 {
1092 if (cell->neighbor_child_on_subface(f, c)
1093 ->is_artificial())
1094 continue;
1095
1096 Assert(cell->face(f)->child(c)->n_active_fe_indices() ==
1097 1,
1098 ExcInternalError());
1099
1100 Assert(cell->face(f)->child(c)->fe_index_is_active(
1101 cell->active_fe_index()) == true,
1102 ExcNotImplemented());
1103 }
1104 } // DEBUG
1105
1106 // Ok, start up the work:
1107 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
1108
1109 const unsigned int n_dofs = fe.n_dofs_per_line();
1110 face_dof_indices.resize(n_dofs);
1111
1112 cell->face(f)->get_dof_indices(face_dof_indices);
1113 const std::vector<types::global_dof_index> dof_on_mother_face =
1114 face_dof_indices;
1115
1116 cell->face(f)->child(0)->get_dof_indices(face_dof_indices);
1117 const std::vector<types::global_dof_index> dof_on_child_face_0 =
1118 face_dof_indices;
1119
1120 cell->face(f)->child(1)->get_dof_indices(face_dof_indices);
1121 const std::vector<types::global_dof_index> dof_on_child_face_1 =
1122 face_dof_indices;
1123
1124 // As the Nedelec elements are oriented, we need to take care of
1125 // the orientation of the lines.
1126 // Remark: "false" indicates the line is not flipped.
1127 // "true" indicates the line is flipped.
1128
1129 // get the orientation of the faces
1130 const bool direction_mother = (cell->face(f)->vertex_index(0) >
1131 cell->face(f)->vertex_index(1)) ?
1132 false :
1133 true;
1134 const bool direction_child_0 =
1135 (cell->face(f)->child(0)->vertex_index(0) >
1136 cell->face(f)->child(0)->vertex_index(1)) ?
1137 false :
1138 true;
1139 const bool direction_child_1 =
1140 (cell->face(f)->child(1)->vertex_index(0) >
1141 cell->face(f)->child(1)->vertex_index(1)) ?
1142 false :
1143 true;
1144
1145 for (unsigned int row = 0; row < n_dofs; ++row)
1146 {
1147 constraints.add_line(dof_on_child_face_0[row]);
1148 constraints.add_line(dof_on_child_face_1[row]);
1149 }
1150
1151 for (unsigned int row = 0; row < n_dofs; ++row)
1152 {
1153 for (unsigned int dof_i_on_mother = 0;
1154 dof_i_on_mother < n_dofs;
1155 ++dof_i_on_mother)
1156 {
1157 // We need to keep in mind that, if we use a FE_System
1158 // with multiple FE_NedelecSZ blocks inside, we need
1159 // to consider, that n_dofs depends on the number
1160 // of FE_NedelecSZ blocks used.
1161 unsigned int shift_0 =
1162 (direction_mother == direction_child_0) ? 0 : n_dofs;
1163 constraints.add_entry(dof_on_child_face_0[row],
1164 dof_on_mother_face[dof_i_on_mother],
1165 fe.constraints()(row + shift_0,
1166 dof_i_on_mother));
1167
1168 unsigned int shift_1 =
1169 (direction_mother == direction_child_1) ? 0 : n_dofs;
1170 constraints.add_entry(dof_on_child_face_1[row],
1171 dof_on_mother_face[dof_i_on_mother],
1172 fe.constraints()(row + shift_1,
1173 dof_i_on_mother));
1174 }
1175 }
1176 }
1177 }
1178 }
1179
1180
1181 template <int dim_, int spacedim, typename number>
1182 void
1183 make_hanging_node_constraints_nedelec(
1184 const DoFHandler<dim_, spacedim> &dof_handler,
1185 AffineConstraints<number> &constraints,
1186 std::integral_constant<int, 3>)
1187 {
1188 // Parts of this function are very similar to
1189 // make_oldstyle_hanging_node_constraints.
1190 // Therefore, only the parts that differ from the
1191 // make_oldstyle_hanging_node_constraints are commented on.
1192
1193 const unsigned int dim = 3;
1194
1195 // In order to find all hanging edges in a reasonable
1196 // computing time, even when five or more cells share
1197 // one edge, we pre-compute the edge-to-cell map.
1198 // Unfortunately, we need to get rid of a 'const'...
1199 Triangulation<dim, spacedim> &tria =
1200 const_cast<Triangulation<dim, spacedim> &>(
1201 dof_handler.get_triangulation());
1202 tria.compute_line_to_adjacent_cells_map();
1203
1204 std::vector<types::global_dof_index> dofs_on_mother;
1205 std::vector<types::global_dof_index> dofs_on_children;
1206
1207 // loop over all quads
1208 for (const auto &cell : dof_handler.active_cell_iterators())
1209 {
1210 // skip artificial cells
1211 if (cell->is_artificial())
1212 continue;
1213
1214 // loop over all faces
1215 for (const unsigned int face : cell->face_indices())
1216 {
1217 // skip cells without children
1218 if (cell->face(face)->has_children() == false)
1219 continue;
1220
1221 if (cell->get_fe().n_dofs_per_face(face) == 0)
1222 continue;
1223
1224 Assert(cell->face(face)->refinement_case() ==
1225 RefinementCase<dim - 1>::isotropic_refinement,
1226 ExcNotImplemented());
1227
1228 AssertDimension(cell->face(face)->n_active_fe_indices(), 1);
1229
1230 Assert(cell->face(face)->fe_index_is_active(
1231 cell->active_fe_index()) == true,
1232 ExcInternalError());
1233
1234 if constexpr (running_in_debug_mode())
1235 {
1236 for (unsigned int c = 0; c < cell->face(face)->n_children();
1237 ++c)
1238 {
1239 if (cell->neighbor_child_on_subface(face, c)
1240 ->is_artificial())
1241 continue;
1242
1243 AssertDimension(
1244 cell->face(face)->child(c)->n_active_fe_indices(), 1);
1245
1246 Assert(cell->face(face)->child(c)->fe_index_is_active(
1247 cell->active_fe_index()) == true,
1248 ExcNotImplemented());
1249
1250 for (unsigned int e = 0;
1251 e < GeometryInfo<dim>::vertices_per_face;
1252 ++e)
1253 {
1254 Assert(cell->face(face)
1255 ->child(c)
1256 ->line(e)
1257 ->n_active_fe_indices() == 1,
1258 ExcNotImplemented());
1259
1260 Assert(cell->face(face)
1261 ->child(c)
1262 ->line(e)
1263 ->fe_index_is_active(
1264 cell->active_fe_index()) == true,
1265 ExcNotImplemented());
1266 }
1267 }
1268
1269 for (unsigned int e = 0;
1270 e < GeometryInfo<dim>::vertices_per_face;
1271 ++e)
1272 {
1273 Assert(cell->face(face)->line(e)->n_active_fe_indices() ==
1274 1,
1275 ExcNotImplemented());
1276
1277 Assert(cell->face(face)->line(e)->fe_index_is_active(
1278 cell->active_fe_index()) == true,
1279 ExcNotImplemented());
1280 }
1281 } // DEBUG
1282
1283 // Ok, start up the work
1284 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
1285 const unsigned int fe_index = cell->active_fe_index();
1286
1287 // get the polynomial degree
1288 unsigned int degree(fe.degree);
1289
1290 // get the number of DoFs on mother and children;
1291 // number of DoFs on the mother
1292 const unsigned int n_dofs_on_mother = fe.n_dofs_per_face(face);
1293 dofs_on_mother.resize(n_dofs_on_mother);
1294
1295 const unsigned int n_lines_on_mother =
1296 GeometryInfo<dim>::lines_per_face;
1297
1298 // number of internal lines of the children;
1299 // for more details see description of the
1300 // GeometryInfo<dim> class
1301 // .................
1302 // . | .
1303 // . c2 1 c3 .
1304 // . | .
1305 // .---2---+---3---.
1306 // . | .
1307 // . c0 0 c1 .
1308 // . | .
1309 // .................
1310 const unsigned int n_internal_lines_on_children = 4;
1311
1312 // number of external lines of the children
1313 // +---6--------7--+
1314 // | . |
1315 // 1 c2 . c3 3
1316 // | . |
1317 // |...............|
1318 // | . |
1319 // 0 c0 . c1 2
1320 // | . |
1321 // +---4---+---5---+
1322 const unsigned int n_external_lines_on_children = 8;
1323
1324 const unsigned int n_lines_on_children =
1325 n_internal_lines_on_children + n_external_lines_on_children;
1326
1327 // we only consider the isotropic case here
1328 const unsigned int n_children_per_face =
1329 GeometryInfo<dim>::max_children_per_face;
1330 const unsigned int n_children_per_line =
1331 GeometryInfo<dim - 1>::max_children_per_face;
1332
1333 // number of DoFs on the children
1334 // Remark: Nedelec elements have no DoFs on the vertices,
1335 // therefore we skip the vertices
1336 const unsigned int n_dofs_on_children =
1337 (n_lines_on_children * fe.n_dofs_per_line() +
1338 n_children_per_face * fe.n_dofs_per_quad(face));
1339
1340 dofs_on_children.clear();
1341 dofs_on_children.reserve(n_dofs_on_children);
1342
1343 AssertDimension(n_dofs_on_mother, fe.constraints().n());
1344 AssertDimension(n_dofs_on_children, fe.constraints().m());
1345
1346 // get the current face
1347 const typename DoFHandler<dim, dim>::face_iterator this_face =
1348 cell->face(face);
1349
1350 // fill the DoFs on the mother:
1351 unsigned int next_index = 0;
1352
1353 // DoFs on vertices:
1354 // Nedelec elements have no DoFs on the vertices
1355
1356 // DoFs on lines:
1357 for (unsigned int line = 0;
1358 line < GeometryInfo<dim>::lines_per_face;
1359 ++line)
1360 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof)
1361 dofs_on_mother[next_index++] =
1362 this_face->line(line)->dof_index(dof, fe_index);
1363
1364 // DoFs on the face:
1365 for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face); ++dof)
1366 dofs_on_mother[next_index++] =
1367 this_face->dof_index(dof, fe_index);
1368
1369 // check that we have added all DoFs
1370 AssertDimension(next_index, dofs_on_mother.size());
1371
1372 // the implementation does not support anisotropic refinement
1373 Assert(!dof_handler.get_triangulation()
1374 .get_anisotropic_refinement_flag(),
1375 ExcInternalError());
1376
1377 // fill the DoF on the children:
1378 // DoFs on vertices:
1379 // Nedelec elements have no DoFs on the vertices
1380
1381 // DoFs on lines:
1382 // the DoFs on the interior lines to the children; the order
1383 // of these lines is shown above (see
1384 // n_internal_lines_on_children)
1385 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
1386 dofs_on_children.push_back(
1387 this_face->child(0)->line(1)->dof_index(dof, fe_index));
1388
1389 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
1390 dofs_on_children.push_back(
1391 this_face->child(2)->line(1)->dof_index(dof, fe_index));
1392
1393 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
1394 dofs_on_children.push_back(
1395 this_face->child(0)->line(3)->dof_index(dof, fe_index));
1396
1397 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
1398 dofs_on_children.push_back(
1399 this_face->child(1)->line(3)->dof_index(dof, fe_index));
1400
1401 // DoFs on the bordering lines:
1402 // DoFs on the exterior lines to the children; the order of
1403 // these lines is shown above (see n_external_lines_on_children)
1404 for (unsigned int line = 0;
1405 line < GeometryInfo<dim>::lines_per_face;
1406 ++line)
1407 for (unsigned int child = 0; child < n_children_per_line;
1408 ++child)
1409 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
1410 dofs_on_children.push_back(
1411 this_face->line(line)->child(child)->dof_index(dof,
1412 fe_index));
1413
1414 // DoFs on the faces of the four children:
1415 for (unsigned int child = 0; child < n_children_per_face; ++child)
1416 {
1417 // skip artificial cells
1418 if (cell->neighbor_child_on_subface(face, child)
1419 ->is_artificial())
1420 continue;
1421
1422 for (unsigned int dof = 0; dof < fe.n_dofs_per_quad(face);
1423 ++dof)
1424 dofs_on_children.push_back(
1425 this_face->child(child)->dof_index(dof, fe_index));
1426 } // rof: child
1427
1428 // consistency check:
1429 // note: we can get fewer DoFs when we have artificial cells
1430 Assert(dofs_on_children.size() <= n_dofs_on_children,
1431 ExcInternalError());
1432
1433 // As the Nedelec elements are oriented, we need to take care of
1434 // the orientation of the lines.
1435 // Remark: "false" indicates the line is not flipped.
1436 // "true" indicates the line is flipped.
1437
1438 // Orientation - Lines:
1439 // get the orientation from the edges from the mother cell
1440 std::vector<bool> direction_mother(
1441 GeometryInfo<dim>::lines_per_face, false);
1442 for (unsigned int line = 0;
1443 line < GeometryInfo<dim>::lines_per_face;
1444 ++line)
1445 if (this_face->line(line)->vertex_index(0) >
1446 this_face->line(line)->vertex_index(1))
1447 direction_mother[line] = true;
1448
1449 // get the orientation from the intern edges of the children
1450 std::vector<bool> direction_child_intern(
1451 n_internal_lines_on_children, false);
1452
1453 // get the global vertex index of vertex in the center;
1454 // we need this vertex index, to compute the direction
1455 // of the internal edges
1456 unsigned int center = this_face->child(0)->vertex_index(3);
1457
1458 // compute the direction of the internal edges
1459 for (unsigned int line = 0; line < n_internal_lines_on_children;
1460 ++line)
1461 if (line % 2 == 0)
1462 {
1463 direction_child_intern[line] =
1464 this_face->line(line)->child(0)->vertex_index(1) <
1465 center ?
1466 false :
1467 true;
1468 }
1469 else
1470 {
1471 direction_child_intern[line] =
1472 this_face->line(line)->child(0)->vertex_index(1) >
1473 center ?
1474 false :
1475 true;
1476 }
1477
1478 // compute the direction of the outer edges
1479 std::vector<bool> direction_child(n_external_lines_on_children,
1480 false);
1481 for (unsigned int line = 0;
1482 line < GeometryInfo<dim>::lines_per_face;
1483 ++line)
1484 {
1485 if (this_face->line(line)->child(0)->vertex_index(0) >
1486 this_face->line(line)->child(0)->vertex_index(1))
1487 direction_child[2 * line] = true;
1488 if (this_face->line(line)->child(1)->vertex_index(0) >
1489 this_face->line(line)->child(1)->vertex_index(1))
1490 direction_child[2 * line + 1] = true;
1491 }
1492
1493
1494 // Orientation - Faces:
1495 bool mother_flip_x = false;
1496 bool mother_flip_y = false;
1497 bool mother_flip_xy = false;
1498 std::vector<bool> child_flip_x(n_children_per_face, false);
1499 std::vector<bool> child_flip_y(n_children_per_face, false);
1500 std::vector<bool> child_flip_xy(n_children_per_face, false);
1501 const unsigned int
1502 vertices_adjacent_on_face[GeometryInfo<dim>::vertices_per_face]
1503 [2] = {{1, 2}, {0, 3}, {3, 0}, {2, 1}};
1504
1505 {
1506 // Mother
1507 // get the position of the vertex with the highest number
1508 unsigned int current_glob = cell->face(face)->vertex_index(0);
1509 unsigned int current_max = 0;
1510 for (unsigned int v = 1;
1511 v < GeometryInfo<dim>::vertices_per_face;
1512 ++v)
1513 if (current_glob < this_face->vertex_index(v))
1514 {
1515 current_max = v;
1516 current_glob = this_face->vertex_index(v);
1517 }
1518
1519 // if the vertex with the highest DoF index is in the lower row
1520 // of the face, the face is flipped in y direction
1521 if (current_max < 2)
1522 mother_flip_y = true;
1523
1524 // if the vertex with the highest DoF index is on the left side
1525 // of the face is flipped in x direction
1526 if (current_max % 2 == 0)
1527 mother_flip_x = true;
1528
1529 // get the minor direction of the face of the mother
1530 if (this_face->vertex_index(
1531 vertices_adjacent_on_face[current_max][0]) <
1532 this_face->vertex_index(
1533 vertices_adjacent_on_face[current_max][1]))
1534 mother_flip_xy = true;
1535 }
1536
1537 // Children:
1538 // get the orientation of the faces of the children
1539 for (unsigned int child = 0; child < n_children_per_face; ++child)
1540 {
1541 unsigned int current_max = 0;
1542 unsigned int current_glob =
1543 this_face->child(child)->vertex_index(0);
1544
1545 for (unsigned int v = 1;
1546 v < GeometryInfo<dim>::vertices_per_face;
1547 ++v)
1548 if (current_glob < this_face->child(child)->vertex_index(v))
1549 {
1550 current_max = v;
1551 current_glob = this_face->child(child)->vertex_index(v);
1552 }
1553
1554 if (current_max < 2)
1555 child_flip_y[child] = true;
1556
1557 if (current_max % 2 == 0)
1558 child_flip_x[child] = true;
1559
1560 if (this_face->child(child)->vertex_index(
1561 vertices_adjacent_on_face[current_max][0]) <
1562 this_face->child(child)->vertex_index(
1563 vertices_adjacent_on_face[current_max][1]))
1564 child_flip_xy[child] = true;
1565
1566 child_flip_xy[child] = mother_flip_xy;
1567 }
1568
1569 // copy the constraint matrix, since we need to modify that matrix
1570 std::vector<std::vector<double>> constraints_matrix(
1571 n_lines_on_children * fe.n_dofs_per_line() +
1572 n_children_per_face * fe.n_dofs_per_quad(),
1573 std::vector<double>(dofs_on_mother.size(), 0));
1574
1575 {
1576 // copy the constraint matrix
1577 // internal lines
1578 for (unsigned int line = 0; line < n_internal_lines_on_children;
1579 ++line)
1580 {
1581 unsigned int row_start = line * fe.n_dofs_per_line();
1582 unsigned int line_mother = line / 2;
1583 unsigned int row_mother =
1584 (line_mother * 2) * fe.n_dofs_per_line();
1585 for (unsigned int row = 0; row < fe.n_dofs_per_line();
1586 ++row)
1587 for (unsigned int i = 0;
1588 i < n_lines_on_mother * fe.n_dofs_per_line();
1589 ++i)
1590 constraints_matrix[row + row_start][i] =
1591 fe.constraints()(row + row_mother, i);
1592 }
1593
1594 for (unsigned int line = 0; line < n_internal_lines_on_children;
1595 ++line)
1596 {
1597 unsigned int row_start = line * fe.n_dofs_per_line();
1598 unsigned int line_mother = line / 2;
1599 unsigned int row_mother =
1600 (line_mother * 2) * fe.n_dofs_per_line();
1601 for (unsigned int row = 0; row < fe.n_dofs_per_line();
1602 ++row)
1603 for (unsigned int i =
1604 n_lines_on_mother * fe.n_dofs_per_line();
1605 i < dofs_on_mother.size();
1606 ++i)
1607 constraints_matrix[row + row_start][i] =
1608 fe.constraints()(row + row_mother, i);
1609 }
1610
1611 // external lines
1612 unsigned int row_offset =
1613 n_internal_lines_on_children * fe.n_dofs_per_line();
1614 for (unsigned int line = 0; line < n_external_lines_on_children;
1615 line++)
1616 {
1617 unsigned int row_start = line * fe.n_dofs_per_line();
1618 unsigned int line_mother = line / 2;
1619 unsigned int row_mother =
1620 (line_mother * 2) * fe.n_dofs_per_line();
1621 for (unsigned int row = row_offset;
1622 row < row_offset + fe.n_dofs_per_line();
1623 ++row)
1624 for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
1625 constraints_matrix[row + row_start][i] =
1626 fe.constraints()(row + row_mother, i);
1627 }
1628
1629 // copy the weights for the faces
1630 row_offset = n_lines_on_children * fe.n_dofs_per_line();
1631 for (unsigned int face = 0; face < n_children_per_face; ++face)
1632 {
1633 unsigned int row_start = face * fe.n_dofs_per_quad();
1634 for (unsigned int row = row_offset;
1635 row < row_offset + fe.n_dofs_per_quad();
1636 row++)
1637 for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
1638 constraints_matrix[row + row_start][i] =
1639 fe.constraints()(row, i);
1640 }
1641 }
1642
1643 // Modify the matrix
1644 // Edge - Edge:
1645 // Interior edges: the interior edges have support on the
1646 // corresponding edges and faces loop over all 4 intern edges
1647 for (unsigned int i = 0;
1648 i < n_internal_lines_on_children * fe.n_dofs_per_line();
1649 ++i)
1650 {
1651 unsigned int line_i = i / fe.n_dofs_per_line();
1652 unsigned int tmp_i = i % degree;
1653
1654 // loop over the edges of the mother cell
1655 for (unsigned int j = 0;
1656 j < n_lines_on_mother * fe.n_dofs_per_line();
1657 ++j)
1658 {
1659 unsigned int line_j = j / fe.n_dofs_per_line();
1660 unsigned int tmp_j = j % degree;
1661
1662 if ((line_i < 2 && line_j < 2) ||
1663 (line_i >= 2 && line_j >= 2))
1664 {
1665 if (direction_child_intern[line_i] !=
1666 direction_mother[line_j])
1667 {
1668 if ((tmp_i + tmp_j) % 2 == 1)
1669 { // anti-symmetric
1670 constraints_matrix[i][j] *= -1.0;
1671 }
1672 }
1673 }
1674 else
1675 {
1676 if (direction_mother[line_i])
1677 {
1678 if ((tmp_i + tmp_j) % 2 == 1)
1679 { // anti-symmetric
1680 constraints_matrix[i][j] *= -1.0;
1681 }
1682 }
1683 }
1684 }
1685 }
1686
1687 // Exterior edges:
1688 for (unsigned int i =
1689 n_internal_lines_on_children * fe.n_dofs_per_line();
1690 i < n_lines_on_children * fe.n_dofs_per_line();
1691 ++i)
1692 {
1693 unsigned int line_i = (i / fe.n_dofs_per_line()) - 4;
1694 unsigned int tmp_i = i % degree;
1695
1696 // loop over the edges of the mother cell
1697 for (unsigned int j = 0;
1698 j < n_lines_on_mother * fe.n_dofs_per_line();
1699 ++j)
1700 {
1701 unsigned int line_j = j / fe.n_dofs_per_line();
1702 unsigned int tmp_j = j % degree;
1703
1704 if (direction_child[line_i] != direction_mother[line_j])
1705 {
1706 if ((tmp_i + tmp_j) % 2 == 1)
1707 { // anti-symmetric
1708 constraints_matrix[i][j] *= -1.0;
1709 }
1710 }
1711 }
1712 }
1713
1714 // Note:
1715 // We need to keep in mind that, if we use a FE_System
1716 // with multiple FE_NedelecSZ blocks inside, we need
1717 // to consider, that fe.n_dofs_per_line() depends on the number
1718 // of FE_NedelecSZ blocks used.
1719 const unsigned int n_blocks = fe.n_dofs_per_line() / degree;
1720
1721 // Edge - Face
1722 // Interior edges: for x-direction
1723 for (unsigned int i = 0; i < 2 * fe.n_dofs_per_line(); ++i)
1724 {
1725 unsigned int line_i = i / fe.n_dofs_per_line();
1726 unsigned int tmp_i = i % degree;
1727
1728 unsigned int start_j =
1729 n_lines_on_mother * fe.n_dofs_per_line();
1730
1731 for (unsigned int block = 0; block < n_blocks; ++block)
1732 {
1733 // Type 1:
1734 for (unsigned int jy = 0; jy < degree - 1; ++jy)
1735 for (unsigned int jx = 0; jx < degree - 1; ++jx)
1736 {
1737 unsigned int j = start_j + jx + (jy * (degree - 1));
1738 if (direction_child_intern[line_i] != mother_flip_y)
1739 {
1740 if ((jy + tmp_i) % 2 == 0)
1741 { // anti-symmetric case
1742 constraints_matrix[i][j] *= -1.0;
1743 }
1744 }
1745 }
1746
1747 start_j += (degree - 1) * (degree - 1);
1748
1749 // Type 2:
1750 for (unsigned int jy = 0; jy < degree - 1; ++jy)
1751 for (unsigned int jx = 0; jx < degree - 1; ++jx)
1752 {
1753 unsigned int j = start_j + jx + (jy * (degree - 1));
1754
1755 if (direction_child_intern[line_i] != mother_flip_y)
1756 {
1757 if ((jy + tmp_i) % 2 == 0)
1758 { // anti-symmetric case
1759 constraints_matrix[i][j] *= -1.0;
1760 }
1761 }
1762 }
1763 start_j += (degree - 1) * (degree - 1);
1764
1765 // Type 3.1:
1766 // nothing to do
1767 start_j += degree - 1;
1768
1769 // Type 3.2:
1770 // nothing to do
1771 start_j += degree - 1;
1772 }
1773 }
1774
1775 // Interior edges: for y-direction
1776 for (unsigned int i = 2 * fe.n_dofs_per_line();
1777 i < 4 * fe.n_dofs_per_line();
1778 i++)
1779 {
1780 unsigned int line_i = i / fe.n_dofs_per_line();
1781 unsigned int tmp_i = i % degree;
1782
1783 unsigned int start_j =
1784 n_lines_on_mother * fe.n_dofs_per_line();
1785
1786 for (unsigned int block = 0; block < n_blocks; block++)
1787 {
1788 // Type 1:
1789 for (unsigned int jy = 0; jy < degree - 1; ++jy)
1790 for (unsigned int jx = 0; jx < degree - 1; ++jx)
1791 {
1792 unsigned int j = start_j + jx + (jy * (degree - 1));
1793 if (direction_child_intern[line_i] != mother_flip_x)
1794 {
1795 if ((jx + tmp_i) % 2 == 0)
1796 { // anti-symmetric case
1797 constraints_matrix[i][j] *= -1.0;
1798 }
1799 }
1800 }
1801
1802 start_j += (degree - 1) * (degree - 1);
1803
1804 // Type 2:
1805 for (unsigned int jy = 0; jy < degree - 1; ++jy)
1806 for (unsigned int jx = 0; jx < degree - 1; ++jx)
1807 {
1808 unsigned int j = start_j + jx + (jy * (degree - 1));
1809 if (direction_child_intern[line_i] != mother_flip_x)
1810 {
1811 if ((jx + tmp_i) % 2 == 0)
1812 { // anti-symmetric case
1813 constraints_matrix[i][j] *= -1.0;
1814 }
1815 }
1816 }
1817 start_j += (degree - 1) * (degree - 1);
1818
1819 // Type 3.1:
1820 // nothing to do
1821 start_j += degree - 1;
1822
1823 // Type 3.2:
1824 // nothing to do
1825 start_j += degree - 1;
1826 }
1827 }
1828
1829 // Face - Face
1830 unsigned int degree_square = (degree - 1) * (degree - 1);
1831 {
1832 // Face
1833 unsigned int i = n_lines_on_children * fe.n_dofs_per_line();
1834 for (unsigned int child_face = 0;
1835 child_face < n_children_per_face;
1836 ++child_face)
1837 for (unsigned int block = 0; block < n_blocks; ++block)
1838 {
1839 unsigned int block_size = fe.n_dofs_per_quad() / n_blocks;
1840
1841 // check if the counting of the DoFs is correct:
1842 Assert((block == 0 &&
1843 i != n_lines_on_children * fe.n_dofs_per_line() +
1844 child_face * fe.n_dofs_per_quad()) ==
1845 false,
1846 ExcInternalError());
1847
1848 // Type 1:
1849 for (unsigned int iy = 0; iy < degree - 1; ++iy)
1850 for (unsigned int ix = 0; ix < degree - 1; ++ix)
1851 {
1852 // Type 1 on mother:
1853 unsigned int j =
1854 n_lines_on_mother * fe.n_dofs_per_line() +
1855 block * block_size;
1856 for (unsigned int jy = 0; jy < degree - 1; ++jy)
1857 for (unsigned int jx = 0; jx < degree - 1; ++jx)
1858 {
1859 if (child_flip_x[child_face] !=
1860 mother_flip_x) // x - direction (x-flip)
1861 {
1862 if ((ix + jx) % 2 == 1)
1863 { // anti-symmetric in x
1864 constraints_matrix[i][j] *= -1.0;
1865 }
1866 }
1867
1868 if (child_flip_y[child_face] !=
1869 mother_flip_y) // y - direction (y-flip)
1870 {
1871 if ((iy + jy) % 2 == 1)
1872 { // anti-symmetric in y
1873 constraints_matrix[i][j] *= -1.0;
1874 }
1875 }
1876
1877 j++;
1878 }
1879 i++;
1880 }
1881
1882 // Type 2:
1883 for (unsigned int iy = 0; iy < degree - 1; ++iy)
1884 for (unsigned int ix = 0; ix < degree - 1; ++ix)
1885 {
1886 // Type 2 on mother:
1887 unsigned int j =
1888 n_lines_on_mother * fe.n_dofs_per_line() +
1889 degree_square + block * block_size;
1890 for (unsigned int jy = 0; jy < degree - 1; ++jy)
1891 for (unsigned int jx = 0; jx < degree - 1; ++jx)
1892 {
1893 if (child_flip_x[child_face] !=
1894 mother_flip_x) // x - direction (x-flip)
1895 {
1896 if ((ix + jx) % 2 == 1)
1897 { // anti-symmetric in x
1898 constraints_matrix[i][j] *= -1.0;
1899 }
1900 }
1901
1902 if (child_flip_y[child_face] !=
1903 mother_flip_y) // y - direction (y-flip)
1904 {
1905 if ((iy + jy) % 2 == 1)
1906 { // anti-symmetric in y
1907 constraints_matrix[i][j] *= -1.0;
1908 }
1909 }
1910
1911 j++;
1912 }
1913
1914 i++;
1915 }
1916
1917
1918 // Type 3 (y):
1919 for (unsigned int iy = 0; iy < degree - 1; ++iy)
1920 {
1921 // Type 2 on mother:
1922 unsigned int j =
1923 n_lines_on_mother * fe.n_dofs_per_line() +
1924 degree_square + block * block_size;
1925 for (unsigned int jy = 0; jy < degree - 1; ++jy)
1926 for (unsigned int jx = 0; jx < degree - 1; ++jx)
1927 {
1928 if (child_flip_x[child_face] !=
1929 mother_flip_x) // x - direction (x-flip)
1930 {
1931 if ((jx) % 2 == 0)
1932 { // anti-symmetric in x
1933 constraints_matrix[i][j] *= -1.0;
1934 }
1935 }
1936
1937 if (child_flip_y[child_face] !=
1938 mother_flip_y) // y - direction (y-flip)
1939 {
1940 if ((iy + jy) % 2 == 1)
1941 { // anti-symmetric in y
1942 constraints_matrix[i][j] *= -1.0;
1943 }
1944 }
1945
1946 j++;
1947 }
1948
1949 // Type 3 on mother:
1950 j = n_lines_on_mother * fe.n_dofs_per_line() +
1951 2 * degree_square + block * block_size;
1952 for (unsigned int jy = 0; jy < degree - 1; ++jy)
1953 {
1954 if (child_flip_y[child_face] !=
1955 mother_flip_y) // y - direction (y-flip)
1956 {
1957 if ((iy + jy) % 2 == 1)
1958 { // anti-symmetric in y
1959 constraints_matrix[i][j] *= -1.0;
1960 }
1961 }
1962
1963 j++;
1964 }
1965 i++;
1966 }
1967
1968 // Type 3 (x):
1969 for (unsigned int ix = 0; ix < degree - 1; ++ix)
1970 {
1971 // Type 2 on mother:
1972 unsigned int j =
1973 n_lines_on_mother * fe.n_dofs_per_line() +
1974 degree_square + block * block_size;
1975 for (unsigned int jy = 0; jy < degree - 1; ++jy)
1976 for (unsigned int jx = 0; jx < degree - 1; ++jx)
1977 {
1978 if (child_flip_x[child_face] !=
1979 mother_flip_x) // x - direction (x-flip)
1980 {
1981 if ((ix + jx) % 2 == 1)
1982 { // anti-symmetric in x
1983 constraints_matrix[i][j] *= -1.0;
1984 }
1985 }
1986
1987 if (child_flip_y[child_face] !=
1988 mother_flip_y) // y - direction (y-flip)
1989 {
1990 if ((jy) % 2 == 0)
1991 { // anti-symmetric in y
1992 constraints_matrix[i][j] *= -1.0;
1993 }
1994 }
1995
1996 j++;
1997 } // rof: Dof j
1998
1999 // Type 3 on mother:
2000 j = n_lines_on_mother * fe.n_dofs_per_line() +
2001 2 * degree_square + (degree - 1) +
2002 block * block_size;
2003 for (unsigned int jx = 0; jx < degree - 1; ++jx)
2004 {
2005 if (child_flip_x[child_face] !=
2006 mother_flip_x) // x - direction (x-flip)
2007 {
2008 if ((ix + jx) % 2 == 1)
2009 { // anti-symmetric in x
2010 constraints_matrix[i][j] *= -1.0;
2011 }
2012 }
2013
2014 j++;
2015 }
2016 i++;
2017 }
2018 }
2019 }
2020
2021 // Next, after we have adapted the signs in the constraint matrix,
2022 // based on the directions of the edges, we need to modify the
2023 // constraint matrix based on the orientation of the faces (i.e.
2024 // if x and y direction are exchanged on the face)
2025
2026 // interior edges:
2027 for (unsigned int i = 0;
2028 i < n_internal_lines_on_children * fe.n_dofs_per_line();
2029 ++i)
2030 {
2031 // check if x and y are permuted on the parent's face
2032 if (mother_flip_xy)
2033 {
2034 // copy the constraints:
2035 std::vector<double> constraints_matrix_old(
2036 dofs_on_mother.size(), 0);
2037 for (unsigned int j = 0; j < dofs_on_mother.size(); ++j)
2038 {
2039 constraints_matrix_old[j] = constraints_matrix[i][j];
2040 }
2041
2042 unsigned int j_start =
2043 n_lines_on_mother * fe.n_dofs_per_line();
2044 for (unsigned block = 0; block < n_blocks; block++)
2045 {
2046 // Type 1
2047 for (unsigned int jy = 0; jy < degree - 1; ++jy)
2048 for (unsigned int jx = 0; jx < degree - 1; ++jx)
2049 {
2050 unsigned int j_old =
2051 j_start + jx + (jy * (degree - 1));
2052 unsigned int j_new =
2053 j_start + jy + (jx * (degree - 1));
2054 constraints_matrix[i][j_new] =
2055 constraints_matrix_old[j_old];
2056 }
2057 j_start += degree_square;
2058
2059 // Type 2
2060 for (unsigned int jy = 0; jy < degree - 1; ++jy)
2061 for (unsigned int jx = 0; jx < degree - 1; ++jx)
2062 {
2063 unsigned int j_old =
2064 j_start + jx + (jy * (degree - 1));
2065 unsigned int j_new =
2066 j_start + jy + (jx * (degree - 1));
2067 constraints_matrix[i][j_new] =
2068 -constraints_matrix_old[j_old];
2069 }
2070 j_start += degree_square;
2071
2072 // Type 3
2073 for (unsigned int j = j_start;
2074 j < j_start + (degree - 1);
2075 j++)
2076 {
2077 constraints_matrix[i][j] =
2078 constraints_matrix_old[j + (degree - 1)];
2079 constraints_matrix[i][j + (degree - 1)] =
2080 constraints_matrix_old[j];
2081 }
2082 j_start += 2 * (degree - 1);
2083 }
2084 }
2085 }
2086
2087 {
2088 // faces:
2089 const unsigned int deg = degree - 1;
2090
2091 // copy the constraints
2092 std::vector<std::vector<double>> constraints_matrix_old(
2093 4 * fe.n_dofs_per_quad(),
2094 std::vector<double>(fe.n_dofs_per_quad(), 0));
2095 for (unsigned int i = 0;
2096 i < n_children_per_face * fe.n_dofs_per_quad();
2097 ++i)
2098 for (unsigned int j = 0; j < fe.n_dofs_per_quad(); ++j)
2099 constraints_matrix_old[i][j] = constraints_matrix
2100 [i + (n_lines_on_children * fe.n_dofs_per_line())]
2101 [j + (n_lines_on_mother * fe.n_dofs_per_line())];
2102
2103 // permute rows (on child)
2104 for (unsigned int child = 0; child < n_children_per_face;
2105 ++child)
2106 {
2107 if (!child_flip_xy[child])
2108 continue;
2109
2110 unsigned int i_start_new =
2111 n_lines_on_children * fe.n_dofs_per_line() +
2112 (child * fe.n_dofs_per_quad());
2113 unsigned int i_start_old = child * fe.n_dofs_per_quad();
2114
2115 unsigned int j_start =
2116 n_lines_on_mother * fe.n_dofs_per_line();
2117
2118 for (unsigned int block = 0; block < n_blocks; block++)
2119 {
2120 // Type 1:
2121 for (unsigned int ix = 0; ix < deg; ++ix)
2122 {
2123 for (unsigned int iy = 0; iy < deg; ++iy)
2124 {
2125 for (unsigned int j = 0;
2126 j < fe.n_dofs_per_quad();
2127 ++j)
2128 constraints_matrix[i_start_new + iy +
2129 (ix * deg)][j + j_start] =
2130 constraints_matrix_old[i_start_old + ix +
2131 (iy * deg)][j];
2132 }
2133 }
2134 i_start_new += deg * deg;
2135 i_start_old += deg * deg;
2136
2137 // Type 2:
2138 for (unsigned int ix = 0; ix < deg; ++ix)
2139 {
2140 for (unsigned int iy = 0; iy < deg; ++iy)
2141 {
2142 for (unsigned int j = 0;
2143 j < fe.n_dofs_per_quad();
2144 j++)
2145 constraints_matrix[i_start_new + iy +
2146 (ix * deg)][j + j_start] =
2147 -constraints_matrix_old[i_start_old + ix +
2148 (iy * deg)][j];
2149 }
2150 }
2151 i_start_new += deg * deg;
2152 i_start_old += deg * deg;
2153
2154 // Type 3:
2155 for (unsigned int ix = 0; ix < deg; ++ix)
2156 {
2157 for (unsigned int j = 0; j < fe.n_dofs_per_quad();
2158 ++j)
2159 constraints_matrix[i_start_new + ix][j +
2160 j_start] =
2161 constraints_matrix_old[i_start_old + ix + deg]
2162 [j];
2163 for (unsigned int j = 0; j < fe.n_dofs_per_quad();
2164 ++j)
2165 constraints_matrix[i_start_new + ix +
2166 deg][j + j_start] =
2167 constraints_matrix_old[i_start_old + ix][j];
2168 } // rof: ix
2169
2170 i_start_new += 2 * deg;
2171 i_start_old += 2 * deg;
2172 }
2173 }
2174
2175 // update the constraints_old
2176 for (unsigned int i = 0;
2177 i < n_children_per_face * fe.n_dofs_per_quad();
2178 i++)
2179 for (unsigned int j = 0; j < fe.n_dofs_per_quad(); j++)
2180 constraints_matrix_old[i][j] = constraints_matrix
2181 [i + (n_lines_on_children * fe.n_dofs_per_line())]
2182 [j + (n_lines_on_mother * fe.n_dofs_per_line())];
2183
2184 // Mother
2185 if (mother_flip_xy)
2186 {
2187 unsigned int i_start =
2188 n_lines_on_children * fe.n_dofs_per_line();
2189
2190 unsigned int j_start_new =
2191 n_lines_on_mother * fe.n_dofs_per_line();
2192 unsigned int j_start_old = 0;
2193
2194 for (unsigned int block = 0; block < n_blocks; ++block)
2195 {
2196 // Type 1:
2197 for (unsigned int jx = 0; jx < deg; ++jx)
2198 {
2199 for (unsigned int jy = 0; jy < deg; ++jy)
2200 {
2201 for (unsigned int i = 0;
2202 i <
2203 n_children_per_face * fe.n_dofs_per_quad();
2204 ++i)
2205 constraints_matrix[i + i_start][j_start_new +
2206 jy +
2207 (jx * deg)] =
2208 constraints_matrix_old[i][j_start_old + jx +
2209 (jy * deg)];
2210 }
2211 }
2212 j_start_new += deg * deg;
2213 j_start_old += deg * deg;
2214
2215 // Type 2:
2216 for (unsigned int jx = 0; jx < deg; ++jx)
2217 {
2218 for (unsigned int jy = 0; jy < deg; ++jy)
2219 {
2220 for (unsigned int i = 0;
2221 i <
2222 n_children_per_face * fe.n_dofs_per_quad();
2223 ++i)
2224 constraints_matrix[i + i_start][j_start_new +
2225 jy +
2226 (jx * deg)] =
2227 -constraints_matrix_old[i][j_start_old +
2228 jx + (jy * deg)];
2229 }
2230 }
2231 j_start_new += deg * deg;
2232 j_start_old += deg * deg;
2233
2234 // Type 3:
2235 for (unsigned int jx = 0; jx < deg; ++jx)
2236 {
2237 for (unsigned int i = 0;
2238 i < n_children_per_face * fe.n_dofs_per_quad();
2239 ++i)
2240 {
2241 constraints_matrix[i + i_start][j_start_new +
2242 jx] =
2243 constraints_matrix_old[i][j_start_old + jx +
2244 deg];
2245 constraints_matrix[i + i_start][j_start_new +
2246 jx + deg] =
2247 constraints_matrix_old[i][j_start_old + jx];
2248 }
2249 }
2250 j_start_new += 2 * deg;
2251 j_start_old += 2 * deg;
2252 }
2253 }
2254 }
2255
2256 // For each row in the AffineConstraints object for
2257 // this line, add the constraint. We split this into the different
2258 // cases.
2259
2260 // internal edges:
2261 for (unsigned int line = 0; line < n_internal_lines_on_children;
2262 ++line)
2263 {
2264 unsigned int row_start = line * fe.n_dofs_per_line();
2265
2266 for (unsigned int row = 0; row < fe.n_dofs_per_line(); ++row)
2267 {
2268 constraints.add_line(dofs_on_children[row_start + row]);
2269 for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
2270 {
2271 constraints.add_entry(
2272 dofs_on_children[row_start + row],
2273 dofs_on_mother[i],
2274 constraints_matrix[row_start + row][i]);
2275 }
2276 constraints.set_inhomogeneity(
2277 dofs_on_children[row_start + row], 0.);
2278 }
2279 }
2280
2281 // Exterior edges
2282 for (unsigned int line = 0; line < n_external_lines_on_children;
2283 ++line)
2284 {
2285 unsigned int row_start =
2286 (4 * fe.n_dofs_per_line()) + (line * fe.n_dofs_per_line());
2287
2288 for (unsigned int row = 0; row < fe.n_dofs_per_line(); ++row)
2289 {
2290 constraints.add_line(dofs_on_children[row_start + row]);
2291 for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
2292 {
2293 constraints.add_entry(
2294 dofs_on_children[row_start + row],
2295 dofs_on_mother[i],
2296 constraints_matrix[row_start + row][i]);
2297 }
2298 constraints.set_inhomogeneity(
2299 dofs_on_children[row_start + row], 0.);
2300 }
2301 }
2302
2303 // Faces:
2304 for (unsigned int f = 0; f < n_children_per_face; ++f)
2305 {
2306 unsigned int row_start =
2307 (n_lines_on_children * fe.n_dofs_per_line()) +
2308 (f * fe.n_dofs_per_quad());
2309
2310 for (unsigned int row = 0; row < fe.n_dofs_per_quad(); ++row)
2311 {
2312 constraints.add_line(dofs_on_children[row_start + row]);
2313
2314 for (unsigned int i = 0; i < dofs_on_mother.size(); ++i)
2315 {
2316 constraints.add_entry(
2317 dofs_on_children[row_start + row],
2318 dofs_on_mother[i],
2319 constraints_matrix[row_start + row][i]);
2320 }
2321
2322 constraints.set_inhomogeneity(
2323 dofs_on_children[row_start + row], 0.);
2324 }
2325 }
2326 }
2327 }
2328 }
2329
2330
2331 template <int dim, int spacedim, typename number>
2332 void
2333 make_hp_hanging_node_constraints(
2334 const DoFHandler<dim, spacedim> &dof_handler,
2335 AffineConstraints<number> &constraints)
2336 {
2337 // note: this function is going to be hard to understand if you haven't
2338 // read the hp-paper. however, we try to follow the notation laid out
2339 // there, so go read the paper before you try to understand what is going
2340 // on here
2341
2342
2343 // a matrix to be used for constraints below. declared here and simply
2344 // resized down below to avoid permanent re-allocation of memory
2345 FullMatrix<double> constraint_matrix;
2346
2347 // similarly have arrays that will hold primary and dependent dof numbers,
2348 // as well as a scratch array needed for the complicated case below
2349 std::vector<types::global_dof_index> primary_dofs;
2350 std::vector<types::global_dof_index> dependent_dofs;
2351 std::vector<types::global_dof_index> scratch_dofs;
2352
2353 // caches for the face and subface interpolation matrices between
2354 // different (or the same) finite elements. we compute them only once,
2355 // namely the first time they are needed, and then just reuse them
2356 Table<2, std::unique_ptr<FullMatrix<double>>> face_interpolation_matrices(
2357 n_finite_elements(dof_handler), n_finite_elements(dof_handler));
2358 Table<3, std::unique_ptr<FullMatrix<double>>>
2359 subface_interpolation_matrices(
2360 n_finite_elements(dof_handler),
2361 n_finite_elements(dof_handler),
2362 GeometryInfo<dim>::max_children_per_face);
2363
2364 // similarly have a cache for the matrices that are split into their
2365 // primary and dependent parts, and for which the primary part is
2366 // inverted. these two matrices are derived from the face interpolation
2367 // matrix
2368 // as described in the @ref hp_paper "hp-paper"
2369 Table<2,
2370 std::unique_ptr<std::pair<FullMatrix<double>, FullMatrix<double>>>>
2371 split_face_interpolation_matrices(n_finite_elements(dof_handler),
2372 n_finite_elements(dof_handler));
2373
2374 // finally, for each pair of finite elements, have a mask that states
2375 // which of the degrees of freedom on the coarse side of a refined face
2376 // will act as primary dofs.
2377 Table<2, std::unique_ptr<std::vector<bool>>> primary_dof_masks(
2378 n_finite_elements(dof_handler), n_finite_elements(dof_handler));
2379
2380 // loop over all faces
2381 //
2382 // note that even though we may visit a face twice if the neighboring
2383 // cells are equally refined, we can only visit each face with hanging
2384 // nodes once
2385 for (const auto &cell : dof_handler.active_cell_iterators())
2386 {
2387 // artificial cells can at best neighbor ghost cells, but we're not
2388 // interested in these interfaces
2389 if (cell->is_artificial())
2390 continue;
2391
2392 for (const unsigned int face : cell->face_indices())
2393 if (cell->face(face)->has_children())
2394 {
2395 // first of all, make sure that we treat a case which is
2396 // possible, i.e. either no dofs on the face at all or no
2397 // anisotropic refinement
2398 if (cell->get_fe().n_dofs_per_face(face) == 0)
2399 continue;
2400
2401 Assert(cell->face(face)->refinement_case() ==
2402 RefinementCase<dim - 1>::isotropic_refinement,
2403 ExcNotImplemented());
2404
2405 // so now we've found a face of an active cell that has
2406 // children. that means that there are hanging nodes here.
2407
2408 // in any case, faces can have at most two sets of active FE
2409 // indices, but here the face can have only one (namely the same
2410 // as that from the cell we're sitting on), and each of the
2411 // children can have only one as well. check this
2412 Assert(cell->face(face)->n_active_fe_indices() == 1,
2413 ExcInternalError());
2414 Assert(cell->face(face)->fe_index_is_active(
2415 cell->active_fe_index()) == true,
2416 ExcInternalError());
2417 for (unsigned int c = 0; c < cell->face(face)->n_children();
2418 ++c)
2419 if (!cell->neighbor_child_on_subface(face, c)
2420 ->is_artificial())
2421 Assert(cell->face(face)->child(c)->n_active_fe_indices() ==
2422 1,
2423 ExcInternalError());
2424
2425 // first find out whether we can constrain each of the subfaces
2426 // to the mother face. in the lingo of the hp-paper, this would
2427 // be the simple case. note that we can short-circuit this
2428 // decision if the dof_handler doesn't support hp at all
2429 //
2430 // ignore all interfaces with artificial cells
2431 FiniteElementDomination::Domination mother_face_dominates =
2432 FiniteElementDomination::either_element_can_dominate;
2433
2434 // auxiliary variable which holds FE indices of the mother face
2435 // and its subfaces. This knowledge will be needed in hp-case
2436 // with neither_element_dominates.
2437 std::set<unsigned int> fe_ind_face_subface;
2438
2439 if (dof_handler.has_hp_capabilities())
2440 {
2441 fe_ind_face_subface.insert(cell->active_fe_index());
2442 for (unsigned int c = 0;
2443 c < cell->face(face)->n_active_descendants();
2444 ++c)
2445 {
2446 const auto subcell =
2447 cell->neighbor_child_on_subface(face, c);
2448 if (!subcell->is_artificial())
2449 {
2450 mother_face_dominates =
2451 mother_face_dominates &
2452 (cell->get_fe().compare_for_domination(
2453 subcell->get_fe(), /*codim=*/1));
2454 fe_ind_face_subface.insert(
2455 subcell->active_fe_index());
2456 }
2457 }
2458 }
2459
2460 switch (mother_face_dominates)
2461 {
2462 case FiniteElementDomination::this_element_dominates:
2463 case FiniteElementDomination::either_element_can_dominate:
2464 {
2465 // Case 1 (the simple case and the only case that can
2466 // happen for non-hp-DoFHandlers): The coarse element
2467 // dominates the elements on the subfaces (or they are
2468 // all the same)
2469 //
2470 // so we are going to constrain the DoFs on the face
2471 // children against the DoFs on the face itself
2472 primary_dofs.resize(
2473 cell->get_fe().n_dofs_per_face(face));
2474
2475 cell->face(face)->get_dof_indices(
2476 primary_dofs, cell->active_fe_index());
2477
2478 // Now create constraints for the subfaces and
2479 // assemble it. ignore all interfaces with artificial
2480 // cells because we can only get to such interfaces if
2481 // the current cell is a ghost cell
2482 for (unsigned int c = 0;
2483 c < cell->face(face)->n_children();
2484 ++c)
2485 {
2486 if (cell->neighbor_child_on_subface(face, c)
2487 ->is_artificial())
2488 continue;
2489
2490 const typename DoFHandler<dim, spacedim>::
2491 active_face_iterator subface =
2492 cell->face(face)->child(c);
2493
2494 Assert(subface->n_active_fe_indices() == 1,
2495 ExcInternalError());
2496
2497 const types::fe_index subface_fe_index =
2498 subface->nth_active_fe_index(0);
2499
2500 // we sometime run into the situation where for
2501 // example on one big cell we have a FE_Q(1) and on
2502 // the subfaces we have a mixture of FE_Q(1) and
2503 // FE_Nothing. In that case, the face domination is
2504 // either_element_can_dominate for the whole
2505 // collection of subfaces, but on the particular
2506 // subface between FE_Q(1) and FE_Nothing, there are
2507 // no constraints that we need to take care of. in
2508 // that case, just continue
2509 if (cell->get_fe().compare_for_domination(
2510 subface->get_fe(subface_fe_index),
2511 /*codim=*/1) ==
2512 FiniteElementDomination::no_requirements)
2513 continue;
2514
2515 // Same procedure as for the mother cell. Extract
2516 // the face DoFs from the cell DoFs.
2517 dependent_dofs.resize(
2518 subface->get_fe(subface_fe_index)
2519 .n_dofs_per_face(face, c));
2520 subface->get_dof_indices(dependent_dofs,
2521 subface_fe_index);
2522
2523 for (const types::global_dof_index dependent_dof :
2524 dependent_dofs)
2525 {
2526 (void)dependent_dof;
2527 Assert(dependent_dof !=
2528 numbers::invalid_dof_index,
2529 ExcInternalError());
2530 }
2531
2532 // Now create the element constraint for this
2533 // subface.
2534 //
2535 // As a side remark, one may wonder the following:
2536 // neighbor_child is clearly computed correctly,
2537 // i.e. taking into account face_orientation (just
2538 // look at the implementation of that function).
2539 // however, we don't care about this here, when we
2540 // ask for subface_interpolation on subface c. the
2541 // question rather is: do we have to translate 'c'
2542 // here as well?
2543 //
2544 // the answer is in fact 'no'. if one does that,
2545 // results are wrong: constraints are added twice
2546 // for the same pair of nodes but with differing
2547 // weights. in addition, one can look at the
2548 // deal.II/project_*_03 tests that look at exactly
2549 // this case: there, we have a mesh with at least
2550 // one face_orientation==false and hanging nodes,
2551 // and the results of those tests show that the
2552 // result of projection verifies the approximation
2553 // properties of a finite element onto that mesh
2554 ensure_existence_of_subface_matrix(
2555 cell->get_fe(),
2556 subface->get_fe(subface_fe_index),
2557 c,
2558 subface_interpolation_matrices
2559 [cell->active_fe_index()][subface_fe_index][c]);
2560
2561 // Add constraints to global AffineConstraints
2562 // object.
2563 filter_constraints(primary_dofs,
2564 dependent_dofs,
2565 *(subface_interpolation_matrices
2566 [cell->active_fe_index()]
2567 [subface_fe_index][c]),
2568 constraints);
2569 } // loop over subfaces
2570
2571 break;
2572 } // Case 1
2573
2574 case FiniteElementDomination::other_element_dominates:
2575 case FiniteElementDomination::neither_element_dominates:
2576 {
2577 // Case 2 (the "complex" case): at least one (the
2578 // neither_... case) of the finer elements or all of
2579 // them (the other_... case) is dominating. See the hp-
2580 // paper for a way how to deal with this situation
2581 //
2582 // since this is something that can only happen for hp-
2583 // dof handlers, add a check here...
2584 Assert(dof_handler.has_hp_capabilities() == true,
2585 ExcInternalError());
2586
2587 const ::hp::FECollection<dim, spacedim>
2588 &fe_collection = dof_handler.get_fe_collection();
2589 // we first have to find the finite element that is able
2590 // to generate a space that all the other ones can be
2591 // constrained to. At this point we potentially have
2592 // different scenarios:
2593 //
2594 // 1) sub-faces dominate mother face and there is a
2595 // dominating FE among sub faces. We could loop over sub
2596 // faces to find the needed FE index. However, this will
2597 // not work in the case when ...
2598 //
2599 // 2) there is no dominating FE among sub faces (e.g.
2600 // Q1xQ2 vs Q2xQ1), but subfaces still dominate mother
2601 // face (e.g. Q2xQ2). To cover this case we would have
2602 // to find the least dominating element amongst all
2603 // finite elements on sub faces.
2604 //
2605 // 3) Finally, it could happen that we got here because
2606 // neither_element_dominates (e.g. Q1xQ1xQ2 and Q1xQ2xQ1
2607 // for subfaces and Q2xQ1xQ1 for mother face). This
2608 // requires finding the least dominating element amongst
2609 // all finite elements on sub faces and the mother face.
2610 //
2611 // Note that the last solution covers the first two
2612 // scenarios, thus we stick with it assuming that we
2613 // won't lose much time/efficiency.
2614
2615 // If we got here then the set must be nonempty
2616 Assert(fe_ind_face_subface.size() > 0,
2617 ExcInternalError());
2618 const types::fe_index dominating_fe_index =
2619 fe_collection.find_dominating_fe_extended(
2620 fe_ind_face_subface,
2621 /*codim=*/1);
2622
2623 AssertThrow(
2624 dominating_fe_index != numbers::invalid_fe_index,
2625 ExcMessage(
2626 "Could not find a least face dominating FE."));
2627
2628 const FiniteElement<dim, spacedim> &dominating_fe =
2629 dof_handler.get_fe(dominating_fe_index);
2630
2631 // first get the interpolation matrix from the mother to
2632 // the virtual dofs
2633 Assert(dominating_fe.n_dofs_per_face(face) <=
2634 cell->get_fe().n_dofs_per_face(face),
2635 ExcInternalError());
2636
2637 ensure_existence_of_face_matrix(
2638 dominating_fe,
2639 cell->get_fe(),
2640 face_interpolation_matrices[dominating_fe_index]
2641 [cell->active_fe_index()]);
2642
2643 // split this matrix into primary and dependent
2644 // components. invert the primary component
2645 ensure_existence_of_primary_dof_mask(
2646 cell->get_fe(),
2647 dominating_fe,
2648 (*face_interpolation_matrices
2649 [dominating_fe_index][cell->active_fe_index()]),
2650 primary_dof_masks[dominating_fe_index]
2651 [cell->active_fe_index()]);
2652
2653 ensure_existence_of_split_face_matrix(
2654 *face_interpolation_matrices[dominating_fe_index]
2655 [cell->active_fe_index()],
2656 (*primary_dof_masks[dominating_fe_index]
2657 [cell->active_fe_index()]),
2658 split_face_interpolation_matrices
2659 [dominating_fe_index][cell->active_fe_index()]);
2660
2661 const FullMatrix<double>
2662 &restrict_mother_to_virtual_primary_inv =
2663 (split_face_interpolation_matrices
2664 [dominating_fe_index][cell->active_fe_index()]
2665 ->first);
2666
2667 const FullMatrix<double>
2668 &restrict_mother_to_virtual_dependent =
2669 (split_face_interpolation_matrices
2670 [dominating_fe_index][cell->active_fe_index()]
2671 ->second);
2672
2673 // now compute the constraint matrix as the product
2674 // between the inverse matrix and the dependent part
2675 constraint_matrix.reinit(
2676 cell->get_fe().n_dofs_per_face(face) -
2677 dominating_fe.n_dofs_per_face(face),
2678 dominating_fe.n_dofs_per_face(face));
2679 restrict_mother_to_virtual_dependent.mmult(
2680 constraint_matrix,
2681 restrict_mother_to_virtual_primary_inv);
2682
2683 // then figure out the global numbers of primary and
2684 // dependent dofs and apply constraints
2685 scratch_dofs.resize(
2686 cell->get_fe().n_dofs_per_face(face));
2687 cell->face(face)->get_dof_indices(
2688 scratch_dofs, cell->active_fe_index());
2689
2690 // split dofs into primary and dependent components
2691 primary_dofs.clear();
2692 dependent_dofs.clear();
2693 for (unsigned int i = 0;
2694 i < cell->get_fe().n_dofs_per_face(face);
2695 ++i)
2696 if ((*primary_dof_masks[dominating_fe_index]
2697 [cell
2698 ->active_fe_index()])[i] ==
2699 true)
2700 primary_dofs.push_back(scratch_dofs[i]);
2701 else
2702 dependent_dofs.push_back(scratch_dofs[i]);
2703
2704 AssertDimension(primary_dofs.size(),
2705 dominating_fe.n_dofs_per_face(face));
2706 AssertDimension(dependent_dofs.size(),
2707 cell->get_fe().n_dofs_per_face(face) -
2708 dominating_fe.n_dofs_per_face(face));
2709
2710 filter_constraints(primary_dofs,
2711 dependent_dofs,
2712 constraint_matrix,
2713 constraints);
2714
2715
2716
2717 // next we have to deal with the subfaces. do as
2718 // discussed in the hp-paper
2719 for (unsigned int sf = 0;
2720 sf < cell->face(face)->n_children();
2721 ++sf)
2722 {
2723 // ignore interfaces with artificial cells as well
2724 // as interfaces between ghost cells in 2d
2725 if (cell->neighbor_child_on_subface(face, sf)
2726 ->is_artificial() ||
2727 (dim == 2 && cell->is_ghost() &&
2728 cell->neighbor_child_on_subface(face, sf)
2729 ->is_ghost()))
2730 continue;
2731
2732 Assert(cell->face(face)
2733 ->child(sf)
2734 ->n_active_fe_indices() == 1,
2735 ExcInternalError());
2736
2737 const types::fe_index subface_fe_index =
2738 cell->face(face)->child(sf)->nth_active_fe_index(
2739 0);
2740 const FiniteElement<dim, spacedim> &subface_fe =
2741 dof_handler.get_fe(subface_fe_index);
2742
2743 // first get the interpolation matrix from the
2744 // subface to the virtual dofs
2745 Assert(dominating_fe.n_dofs_per_face(face) <=
2746 subface_fe.n_dofs_per_face(face),
2747 ExcInternalError());
2748 ensure_existence_of_subface_matrix(
2749 dominating_fe,
2750 subface_fe,
2751 sf,
2752 subface_interpolation_matrices
2753 [dominating_fe_index][subface_fe_index][sf]);
2754
2755 const FullMatrix<double>
2756 &restrict_subface_to_virtual = *(
2757 subface_interpolation_matrices
2758 [dominating_fe_index][subface_fe_index][sf]);
2759
2760 constraint_matrix.reinit(
2761 subface_fe.n_dofs_per_face(face),
2762 dominating_fe.n_dofs_per_face(face));
2763
2764 restrict_subface_to_virtual.mmult(
2765 constraint_matrix,
2766 restrict_mother_to_virtual_primary_inv);
2767
2768 dependent_dofs.resize(
2769 subface_fe.n_dofs_per_face(face));
2770 cell->face(face)->child(sf)->get_dof_indices(
2771 dependent_dofs, subface_fe_index);
2772
2773 filter_constraints(primary_dofs,
2774 dependent_dofs,
2775 constraint_matrix,
2776 constraints);
2777 } // loop over subfaces
2778
2779 break;
2780 } // Case 2
2781
2782 case FiniteElementDomination::no_requirements:
2783 // there are no continuity requirements between the two
2784 // elements. record no constraints
2785 break;
2786
2787 default:
2788 // we shouldn't get here
2789 DEAL_II_ASSERT_UNREACHABLE();
2790 }
2791 }
2792 else
2793 {
2794 // this face has no children, but it could still be that it is
2795 // shared by two cells that use a different FE index
2796 Assert(cell->face(face)->fe_index_is_active(
2797 cell->active_fe_index()) == true,
2798 ExcInternalError());
2799
2800 // see if there is a neighbor that is an artificial cell. in
2801 // that case, we're not interested in this interface. we test
2802 // this case first since artificial cells may not have an
2803 // active FE index set, etc
2804 if (!cell->at_boundary(face) &&
2805 cell->neighbor(face)->is_artificial())
2806 continue;
2807
2808 // Only if there is a neighbor with a different active FE index
2809 // and the same h-level, some action has to be taken.
2810 if ((dof_handler.has_hp_capabilities()) &&
2811 !cell->face(face)->at_boundary() &&
2812 (cell->neighbor(face)->active_fe_index() !=
2813 cell->active_fe_index()) &&
2814 (!cell->face(face)->has_children() &&
2815 !cell->neighbor_is_coarser(face)))
2816 {
2817 const typename DoFHandler<dim,
2818 spacedim>::level_cell_iterator
2819 neighbor = cell->neighbor(face);
2820
2821 // see which side of the face we have to constrain
2822 switch (
2823 cell->get_fe().compare_for_domination(neighbor->get_fe(),
2824 /*codim=*/1))
2825 {
2826 case FiniteElementDomination::this_element_dominates:
2827 {
2828 // Get DoFs on dominating and dominated side of the
2829 // face
2830 primary_dofs.resize(
2831 cell->get_fe().n_dofs_per_face(face));
2832 cell->face(face)->get_dof_indices(
2833 primary_dofs, cell->active_fe_index());
2834
2835 // break if the n_primary_dofs == 0, because we are
2836 // attempting to constrain to an element that has no
2837 // face dofs
2838 if (primary_dofs.empty())
2839 break;
2840
2841 dependent_dofs.resize(
2842 neighbor->get_fe().n_dofs_per_face(face));
2843 cell->face(face)->get_dof_indices(
2844 dependent_dofs, neighbor->active_fe_index());
2845
2846 // make sure the element constraints for this face
2847 // are available
2848 ensure_existence_of_face_matrix(
2849 cell->get_fe(),
2850 neighbor->get_fe(),
2851 face_interpolation_matrices
2852 [cell->active_fe_index()]
2853 [neighbor->active_fe_index()]);
2854
2855 // Add constraints to global constraint matrix.
2856 filter_constraints(
2857 primary_dofs,
2858 dependent_dofs,
2859 *(face_interpolation_matrices
2860 [cell->active_fe_index()]
2861 [neighbor->active_fe_index()]),
2862 constraints);
2863
2864 break;
2865 }
2866
2867 case FiniteElementDomination::other_element_dominates:
2868 {
2869 // we don't do anything here since we will come back
2870 // to this face from the other cell, at which time
2871 // we will fall into the first case clause above
2872 break;
2873 }
2874
2875 case FiniteElementDomination::
2876 either_element_can_dominate:
2877 {
2878 // it appears as if neither element has any
2879 // constraints on its neighbor. this may be because
2880 // neither element has any DoFs on faces at all. or
2881 // that the two elements are actually the same,
2882 // although they happen to run under different
2883 // fe_indices (this is what happens in
2884 // hp/hp_hanging_nodes_01 for example).
2885 //
2886 // another possibility is what happens in crash_13.
2887 // there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs.
2888 // FESystem(FE_Q(1),FE_DGQ(1)). neither of them
2889 // dominates the other.
2890 //
2891 // a final possibility is that we have something
2892 // like FESystem(FE_Q(1),FE_Q(1)) vs
2893 // FESystem(FE_Q(1),FE_Nothing()), see
2894 // hp/fe_nothing_18/19.
2895 //
2896 // in any case, the point is that it doesn't matter.
2897 // there is nothing to do here.
2898 break;
2899 }
2900
2901 case FiniteElementDomination::neither_element_dominates:
2902 {
2903 // make sure we don't get here twice from each cell
2904 if (cell < neighbor)
2905 break;
2906
2907 // our best bet is to find the common space among
2908 // other FEs in FECollection and then constrain both
2909 // FEs to that one. More precisely, we follow the
2910 // strategy outlined on page 17 of the hp-paper:
2911 // First we find the dominant FE space S. Then we
2912 // divide our dofs in primary and dependent such
2913 // that I^{face,primary}_{S^{face}->S} is
2914 // invertible. And finally constrain dependent dofs
2915 // to primary dofs based on the interpolation
2916 // matrix.
2917
2918 const types::fe_index this_fe_index =
2919 cell->active_fe_index();
2920 const types::fe_index neighbor_fe_index =
2921 neighbor->active_fe_index();
2922 std::set<types::fe_index> fes;
2923 fes.insert(this_fe_index);
2924 fes.insert(neighbor_fe_index);
2925 const ::hp::FECollection<dim, spacedim>
2926 &fe_collection = dof_handler.get_fe_collection();
2927
2928 // TODO: Change set to types::fe_index
2929 const types::fe_index dominating_fe_index =
2930 fe_collection.find_dominating_fe_extended(
2931 {fes.begin(), fes.end()}, /*codim=*/1);
2932
2933 AssertThrow(
2934 dominating_fe_index != numbers::invalid_fe_index,
2935 ExcMessage(
2936 "Could not find the dominating FE for " +
2937 cell->get_fe().get_name() + " and " +
2938 neighbor->get_fe().get_name() +
2939 " inside FECollection."));
2940
2941 const FiniteElement<dim, spacedim> &dominating_fe =
2942 fe_collection[dominating_fe_index];
2943
2944 // TODO: until we hit the second face, the code is a
2945 // copy-paste from h-refinement case...
2946
2947 // first get the interpolation matrix from main FE
2948 // to the virtual dofs
2949 Assert(dominating_fe.n_dofs_per_face(face) <=
2950 cell->get_fe().n_dofs_per_face(face),
2951 ExcInternalError());
2952
2953 ensure_existence_of_face_matrix(
2954 dominating_fe,
2955 cell->get_fe(),
2956 face_interpolation_matrices
2957 [dominating_fe_index][cell->active_fe_index()]);
2958
2959 // split this matrix into primary and dependent
2960 // components. invert the primary component
2961 ensure_existence_of_primary_dof_mask(
2962 cell->get_fe(),
2963 dominating_fe,
2964 (*face_interpolation_matrices
2965 [dominating_fe_index]
2966 [cell->active_fe_index()]),
2967 primary_dof_masks[dominating_fe_index]
2968 [cell->active_fe_index()]);
2969
2970 ensure_existence_of_split_face_matrix(
2971 *face_interpolation_matrices
2972 [dominating_fe_index][cell->active_fe_index()],
2973 (*primary_dof_masks[dominating_fe_index]
2974 [cell->active_fe_index()]),
2975 split_face_interpolation_matrices
2976 [dominating_fe_index][cell->active_fe_index()]);
2977
2978 const FullMatrix<
2979 double> &restrict_mother_to_virtual_primary_inv =
2980 (split_face_interpolation_matrices
2981 [dominating_fe_index][cell->active_fe_index()]
2982 ->first);
2983
2984 const FullMatrix<
2985 double> &restrict_mother_to_virtual_dependent =
2986 (split_face_interpolation_matrices
2987 [dominating_fe_index][cell->active_fe_index()]
2988 ->second);
2989
2990 // now compute the constraint matrix as the product
2991 // between the inverse matrix and the dependent part
2992 constraint_matrix.reinit(
2993 cell->get_fe().n_dofs_per_face(face) -
2994 dominating_fe.n_dofs_per_face(face),
2995 dominating_fe.n_dofs_per_face(face));
2996 restrict_mother_to_virtual_dependent.mmult(
2997 constraint_matrix,
2998 restrict_mother_to_virtual_primary_inv);
2999
3000 // then figure out the global numbers of primary and
3001 // dependent dofs and apply constraints
3002 scratch_dofs.resize(
3003 cell->get_fe().n_dofs_per_face(face));
3004 cell->face(face)->get_dof_indices(
3005 scratch_dofs, cell->active_fe_index());
3006
3007 // split dofs into primary and dependent components
3008 primary_dofs.clear();
3009 dependent_dofs.clear();
3010 for (unsigned int i = 0;
3011 i < cell->get_fe().n_dofs_per_face(face);
3012 ++i)
3013 if ((*primary_dof_masks[dominating_fe_index]
3014 [cell->active_fe_index()])
3015 [i] == true)
3016 primary_dofs.push_back(scratch_dofs[i]);
3017 else
3018 dependent_dofs.push_back(scratch_dofs[i]);
3019
3020 AssertDimension(primary_dofs.size(),
3021 dominating_fe.n_dofs_per_face(
3022 face));
3023 AssertDimension(
3024 dependent_dofs.size(),
3025 cell->get_fe().n_dofs_per_face(face) -
3026 dominating_fe.n_dofs_per_face(face));
3027
3028 filter_constraints(primary_dofs,
3029 dependent_dofs,
3030 constraint_matrix,
3031 constraints);
3032
3033 // now do the same for another FE this is pretty
3034 // much the same we do above to resolve h-refinement
3035 // constraints
3036 Assert(dominating_fe.n_dofs_per_face(face) <=
3037 neighbor->get_fe().n_dofs_per_face(face),
3038 ExcInternalError());
3039
3040 ensure_existence_of_face_matrix(
3041 dominating_fe,
3042 neighbor->get_fe(),
3043 face_interpolation_matrices
3044 [dominating_fe_index]
3045 [neighbor->active_fe_index()]);
3046
3047 const FullMatrix<double>
3048 &restrict_secondface_to_virtual =
3049 *(face_interpolation_matrices
3050 [dominating_fe_index]
3051 [neighbor->active_fe_index()]);
3052
3053 constraint_matrix.reinit(
3054 neighbor->get_fe().n_dofs_per_face(face),
3055 dominating_fe.n_dofs_per_face(face));
3056
3057 restrict_secondface_to_virtual.mmult(
3058 constraint_matrix,
3059 restrict_mother_to_virtual_primary_inv);
3060
3061 dependent_dofs.resize(
3062 neighbor->get_fe().n_dofs_per_face(face));
3063 cell->face(face)->get_dof_indices(
3064 dependent_dofs, neighbor->active_fe_index());
3065
3066 filter_constraints(primary_dofs,
3067 dependent_dofs,
3068 constraint_matrix,
3069 constraints);
3070
3071 break;
3072 }
3073
3074 case FiniteElementDomination::no_requirements:
3075 {
3076 // nothing to do here
3077 break;
3078 }
3079
3080 default:
3081 // we shouldn't get here
3082 DEAL_II_ASSERT_UNREACHABLE();
3083 }
3084 }
3085 }
3086 }
3087 }
3088 } // namespace internal
3089
3090
3091
3092 template <int dim, int spacedim, typename number>
3093 void
3094 make_hanging_node_constraints(const DoFHandler<dim, spacedim> &dof_handler,
3095 AffineConstraints<number> &constraints)
3096 {
3097 Assert(dof_handler.has_active_dofs(),
3098 ExcMessage(
3099 "The given DoFHandler does not have any DoFs. Did you forget to "
3100 "call dof_handler.distribute_dofs()?"));
3101
3102 // Decide whether to use make_hanging_node_constraints_nedelec,
3103 // the new or old make_hanging_node_constraints
3104 // function. If all the FiniteElement or all elements in a FECollection
3105 // support the new face constraint matrix, the new code will be used.
3106 // Otherwise, the old implementation is used for the moment.
3107 if (dof_handler.get_fe().get_name().find("FE_NedelecSZ") !=
3108 std::string::npos)
3109 internal::make_hanging_node_constraints_nedelec(
3110 dof_handler, constraints, std::integral_constant<int, dim>());
3111 else if (dof_handler.get_fe_collection().hp_constraints_are_implemented())
3112 internal::make_hp_hanging_node_constraints(dof_handler, constraints);
3113 else
3114 internal::make_oldstyle_hanging_node_constraints(
3115 dof_handler, constraints, std::integral_constant<int, dim>());
3116 }
3117
3118
3119
3120 namespace internal
3121 {
3122 template <typename FaceIterator, typename number>
3123 void
3124 set_periodicity_constraints(
3125 const FaceIterator &face_1,
3126 const std_cxx20::type_identity_t<FaceIterator> &face_2,
3127 const FullMatrix<double> &transformation,
3128 AffineConstraints<number> &affine_constraints,
3129 const ComponentMask &component_mask,
3130 const types::geometric_orientation combined_orientation,
3131 const number periodicity_factor,
3132 const unsigned int level)
3133 {
3134 static const int dim = FaceIterator::AccessorType::dimension;
3135 static const int spacedim = FaceIterator::AccessorType::space_dimension;
3136
3137 // We need to inquire about some face things, but in this context we do
3138 // not have face numbers. Because of the assumptions above we can just
3139 // assume the face number is zero.
3140 const unsigned int face_no = 0;
3141
3142 const bool use_mg = (level != numbers::invalid_unsigned_int);
3143
3144 // If we don't use multigrid, we should be in the case where face_1 is
3145 // active, i.e. has no children. In the case of multigrid, constraints
3146 // between cells on the same level are set up.
3147 Assert(use_mg || (!face_1->has_children()), ExcInternalError());
3148 Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
3149
3150 AssertDimension(
3151 face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
3152 AssertDimension(
3153 face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
3154
3155 const types::fe_index face_1_index = face_1->nth_active_fe_index(0);
3156 const types::fe_index face_2_index = face_2->nth_active_fe_index(0);
3157 const FiniteElement<dim, spacedim> &fe = face_1->get_fe(face_1_index);
3158 Assert(face_1->get_fe(face_1_index) == face_2->get_fe(face_2_index),
3159 ExcMessage(
3160 "Matching periodic cells need to use the same finite element"));
3161 Assert(component_mask.represents_n_components(fe.n_components()),
3162 ExcMessage(
3163 "The number of components in the mask has to be either "
3164 "zero or equal to the number of components in the finite "
3165 "element."));
3166 const unsigned int dofs_per_face = fe.n_dofs_per_face(face_no);
3167
3168 // If we don't use multigrid and face_2 does have children,
3169 // then we need to iterate over these children and set periodic
3170 // constraints in the inverse direction. In the case of multigrid,
3171 // we don't need to do this, since constraints between cells on
3172 // the same level are set up.
3173
3174 if ((!use_mg) && face_2->has_children())
3175 {
3176 Assert(face_2->n_children() == face_2->reference_cell().n_children(),
3177 ExcNotImplemented());
3178
3179 // Skip further recursion if face_1 carries invalid dof indices,
3180 // i.e., it is on an artificial cell.
3181 std::vector<types::global_dof_index> dofs_1(dofs_per_face);
3182 face_1->get_dof_indices(dofs_1, face_1->nth_active_fe_index(0));
3183 for (unsigned int i = 0; i < dofs_per_face; ++i)
3184 if (dofs_1[i] == numbers::invalid_dof_index)
3185 {
3186 return;
3187 }
3188
3189 FullMatrix<double> child_transformation(dofs_per_face, dofs_per_face);
3190 FullMatrix<double> subface_interpolation(dofs_per_face,
3191 dofs_per_face);
3192
3193 for (unsigned int c = 0; c < face_2->n_children(); ++c)
3194 {
3195 // get the interpolation matrix recursively from the one that
3196 // interpolated from face_1 to face_2 by multiplying from the left
3197 // with the one that interpolates from face_2 to its child
3198 const auto &fe = face_1->get_fe(face_1->nth_active_fe_index(0));
3199 fe.get_subface_interpolation_matrix(fe,
3200 c,
3201 subface_interpolation,
3202 face_no);
3203 subface_interpolation.mmult(child_transformation, transformation);
3204
3205 set_periodicity_constraints(face_1,
3206 face_2->child(c),
3207 child_transformation,
3208 affine_constraints,
3209 component_mask,
3210 combined_orientation,
3211 periodicity_factor);
3212 }
3213 return;
3214 }
3215
3216 //
3217 // If we reached this point then both faces are active. Now all
3218 // that is left is to match the corresponding DoFs of both faces.
3219 //
3220
3221 std::vector<types::global_dof_index> dofs_1(dofs_per_face);
3222 std::vector<types::global_dof_index> dofs_2(dofs_per_face);
3223
3224 // Note that, in 3d, these functions take into account the line
3225 // orientations on each face: i.e., this function does not need to
3226 // consider line orientations.
3227 if (use_mg)
3228 face_1->get_mg_dof_indices(level, dofs_1, face_1_index);
3229 else
3230 face_1->get_dof_indices(dofs_1, face_1_index);
3231
3232 if (use_mg)
3233 face_2->get_mg_dof_indices(level, dofs_2, face_2_index);
3234 else
3235 face_2->get_dof_indices(dofs_2, face_2_index);
3236
3237 // If either of the two faces has an invalid dof index, stop. This is
3238 // so that there is no attempt to match artificial cells of parallel
3239 // distributed triangulations.
3240 //
3241 // While it seems like we ought to be able to avoid even calling
3242 // set_periodicity_constraints for artificial faces, this situation
3243 // can arise when a face that is being made periodic is only
3244 // partially touched by the local subdomain.
3245 // make_periodicity_constraints will be called recursively even for
3246 // the section of the face that is not touched by the local
3247 // subdomain.
3248 //
3249 // Until there is a better way to determine if the cells that
3250 // neighbor a face are artificial, we simply test to see if the face
3251 // does not have a valid dof initialization.
3252
3253 for (unsigned int i = 0; i < dofs_per_face; ++i)
3254 if (dofs_1[i] == numbers::invalid_dof_index ||
3255 dofs_2[i] == numbers::invalid_dof_index)
3256 {
3257 return;
3258 }
3259
3260 // In the case of shared Triangulation with artificial cells all
3261 // cells have valid DoF indices, i.e., the check above does not work.
3262 if (const auto tria = dynamic_cast<
3263 const parallel::shared::Triangulation<dim, spacedim> *>(
3264 &face_1->get_triangulation()))
3265 if (tria->with_artificial_cells() &&
3266 (affine_constraints.get_local_lines().size() != 0))
3267 for (unsigned int i = 0; i < dofs_per_face; ++i)
3268 if ((affine_constraints.get_local_lines().is_element(dofs_1[i]) ==
3269 false) ||
3270 (affine_constraints.get_local_lines().is_element(dofs_2[i]) ==
3271 false))
3272 {
3273 return;
3274 }
3275
3276 // To match DoFs we need to use combined_orientation to permute the face
3277 // DoFs. FiniteElement does not offer a function to do exactly that: i.e.,
3278 // adjust_line_dof_index_for_line_orientation() only considers the
3279 // orientation of a line and adjust_quad_dof_index_for_face_orientation()
3280 // is only for DoFs defined on quads. Hence we compute our own lookup
3281 // table with the necessary information:
3282 std::vector<unsigned int> cell_to_face_index(
3283 fe.dofs_per_cell, numbers::invalid_unsigned_int);
3284
3285 for (unsigned int face_dof = 0; face_dof < dofs_per_face; ++face_dof)
3286 cell_to_face_index[fe.face_to_cell_index(
3287 face_dof, face_no, numbers::default_geometric_orientation)] =
3288 face_dof;
3289
3290 // Build constraints in a vector of pairs that can be
3291 // arbitrarily large, but that holds up to 25 elements without
3292 // external memory allocation. This is good enough for hanging
3293 // node constraints of Q4 elements in 3d, so covers most
3294 // common cases.
3295 boost::container::small_vector<
3296 std::pair<typename AffineConstraints<number>::size_type, number>,
3297 25>
3298 constraint_entries;
3299
3300 //
3301 // Loop over all dofs on face 2 and constrain them against all
3302 // matching dofs on face 1:
3303 //
3304 for (unsigned int i = 0; i < dofs_per_face; ++i)
3305 {
3306 // Obey the component mask
3307 if ((component_mask.n_selected_components(fe.n_components()) !=
3308 fe.n_components()) &&
3309 !component_mask[fe.face_system_to_component_index(i, face_no)
3310 .first])
3311 continue;
3312
3313 // We have to be careful to treat so called "identity
3314 // constraints" special. These are constraints of the form
3315 // x1 == constraint_factor * x_2. In this case, if the constraint
3316 // x2 == 1./constraint_factor * x1 already exists we are in trouble.
3317 //
3318 // Consequently, we have to check that we have indeed such an
3319 // "identity constraint". We do this by looping over all entries
3320 // of the row of the transformation matrix and check whether we
3321 // find exactly one nonzero entry. If this is the case, set
3322 // "is_identity_constrained" to true and record the corresponding
3323 // index and constraint_factor.
3324
3325 bool is_identity_constrained = false;
3326 unsigned int target = numbers::invalid_unsigned_int;
3327 number constraint_factor = periodicity_factor;
3328
3329 constexpr double eps = 1.e-13;
3330 for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
3331 {
3332 const auto entry = transformation(i, jj);
3333 if (std::abs(entry) > eps)
3334 {
3335 if (is_identity_constrained)
3336 {
3337 // We did encounter more than one nonzero entry, so
3338 // the dof is not identity constrained. Set the
3339 // boolean to false and break out of the for loop.
3340 is_identity_constrained = false;
3341 break;
3342 }
3343 is_identity_constrained = true;
3344 target = jj;
3345 constraint_factor = entry * periodicity_factor;
3346 }
3347 }
3348
3349 // Next, we work on all constraints that are not identity
3350 // constraints, i.e., constraints that involve an interpolation
3351 // step that constrains the current dof (on face 2) to more than
3352 // one dof on face 1.
3353
3354 if (!is_identity_constrained)
3355 {
3356 // The current dof is already constrained. There is nothing
3357 // left to do.
3358 if (affine_constraints.is_constrained(dofs_2[i]))
3359 continue;
3360
3361 constraint_entries.clear();
3362 constraint_entries.reserve(dofs_per_face);
3363
3364 for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
3365 {
3366 // Get the correct dof index on face_1 respecting the
3367 // given orientation:
3368 const unsigned int j =
3369 cell_to_face_index[fe.face_to_cell_index(
3370 jj, face_no, combined_orientation)];
3371 Assert(j != numbers::invalid_unsigned_int,
3372 ExcInternalError());
3373
3374 if (std::abs(transformation(i, jj)) > eps)
3375 constraint_entries.emplace_back(dofs_1[j],
3376 transformation(i, jj));
3377 }
3378
3379 // Enter the constraint::
3380 affine_constraints.add_constraint(dofs_2[i],
3381 constraint_entries,
3382 0.);
3383
3384
3385 // Continue with next dof.
3386 continue;
3387 }
3388
3389 // We are left with an "identity constraint".
3390
3391 // Get the correct dof index on face_1 respecting the given
3392 // orientation:
3393 const unsigned int j = cell_to_face_index[fe.face_to_cell_index(
3394 target, face_no, combined_orientation)];
3395 Assert(j != numbers::invalid_unsigned_int, ExcInternalError());
3396
3397 auto dof_left = dofs_1[j];
3398 auto dof_right = dofs_2[i];
3399
3400 // If dof_left is already constrained, or dof_left < dof_right we
3401 // flip the order to ensure that dofs are constrained in a stable
3402 // manner on different MPI processes.
3403 if (affine_constraints.is_constrained(dof_left) ||
3404 (dof_left < dof_right &&
3405 !affine_constraints.is_constrained(dof_right)))
3406 {
3407 std::swap(dof_left, dof_right);
3408 constraint_factor = 1. / constraint_factor;
3409 }
3410
3411 // Next, we try to enter the constraint
3412 // dof_left = constraint_factor * dof_right;
3413
3414 // If both degrees of freedom are constrained, there is nothing we
3415 // can do. Simply continue with the next dof.
3416 if (affine_constraints.is_constrained(dof_left) &&
3417 affine_constraints.is_constrained(dof_right))
3418 continue;
3419
3420 // We have to be careful that adding the current identity
3421 // constraint does not create a constraint cycle. Thus, check for
3422 // a dependency cycle:
3423
3424 bool constraints_are_cyclic = true;
3425 number cycle_constraint_factor = constraint_factor;
3426
3427 for (auto test_dof = dof_right; test_dof != dof_left;)
3428 {
3429 if (!affine_constraints.is_constrained(test_dof))
3430 {
3431 constraints_are_cyclic = false;
3432 break;
3433 }
3434
3435 const auto &constraint_entries =
3436 *affine_constraints.get_constraint_entries(test_dof);
3437 if (constraint_entries.size() == 1)
3438 {
3439 test_dof = constraint_entries[0].first;
3440 cycle_constraint_factor *= constraint_entries[0].second;
3441 }
3442 else
3443 {
3444 constraints_are_cyclic = false;
3445 break;
3446 }
3447 }
3448
3449 // In case of a dependency cycle we, either
3450 // - do nothing if cycle_constraint_factor == 1. In this case all
3451 // degrees
3452 // of freedom are already periodically constrained,
3453 // - otherwise, force all dofs to zero (by setting dof_left to
3454 // zero). The reasoning behind this is the fact that
3455 // cycle_constraint_factor != 1 occurs in situations such as
3456 // x1 == x2 and x2 == -1. * x1. This system is only solved by
3457 // x_1 = x_2 = 0.
3458
3459 if (constraints_are_cyclic)
3460 {
3461 if (std::abs(cycle_constraint_factor - number(1.)) > eps)
3462 affine_constraints.constrain_dof_to_zero(dof_left);
3463 }
3464 else
3465 {
3466 affine_constraints.add_constraint(
3467 dof_left, {{dof_right, constraint_factor}}, 0.);
3468 // The number 1e10 in the assert below is arbitrary. If the
3469 // absolute value of constraint_factor is too large, then probably
3470 // the absolute value of periodicity_factor is too large or too
3471 // small. This would be equivalent to an evanescent wave that has
3472 // a very small wavelength. A quick calculation shows that if
3473 // |periodicity_factor| > 1e10 -> |np.exp(ikd)|> 1e10, therefore k
3474 // is imaginary (evanescent wave) and the evanescent wavelength is
3475 // 0.27 times smaller than the dimension of the structure,
3476 // lambda=((2*pi)/log(1e10))*d. Imaginary wavenumbers can be
3477 // interesting in some cases
3478 // (https://doi.org/10.1103/PhysRevA.94.033813).In order to
3479 // implement the case of in which the wavevector can be imaginary
3480 // it would be necessary to rewrite this function and the dof
3481 // ordering method should be modified.
3482 // Let's take the following constraint a*x1 + b*x2 = 0. You could
3483 // just always pick x1 = b/a*x2, but in practice this is not so
3484 // stable if a could be a small number -- intended to be zero, but
3485 // just very small due to roundoff. Of course, constraining x2 in
3486 // terms of x1 has the same problem. So one chooses x1 = b/a*x2 if
3487 // |b|<|a|, and x2 = a/b*x1 if |a|<|b|.
3488 Assert(std::abs(constraint_factor) < 1e10,
3489 ExcMessage("The periodicity constraint is too large. "
3490 "The parameter periodicity_factor might "
3491 "be too large or too small."));
3492 }
3493 } /* for dofs_per_face */
3494 }
3495 } // namespace internal
3496
3497
3498 namespace
3499 {
3500 // Internally used in make_periodicity_constraints.
3501 //
3502 // Build up a (possibly rotated) interpolation matrix that is used in
3503 // set_periodicity_constraints with the help of user supplied matrix and
3504 // first_vector_components.
3505 template <int dim, int spacedim>
3506 FullMatrix<double>
3507 compute_transformation(
3508 const FiniteElement<dim, spacedim> &fe,
3509 const FullMatrix<double> &matrix,
3510 const std::vector<unsigned int> &first_vector_components)
3511 {
3512 // TODO: the implementation makes the assumption that all faces have the
3513 // same number of dofs
3514 AssertDimension(fe.n_unique_faces(), 1);
3515 const unsigned int face_no = 0;
3516
3517 Assert(matrix.m() == matrix.n(), ExcInternalError());
3518
3519 const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no);
3520
3521 if (matrix.m() == n_dofs_per_face)
3522 {
3523 // In case of m == n == n_dofs_per_face the supplied matrix is already
3524 // an interpolation matrix, so we use it directly:
3525 return matrix;
3526 }
3527
3528 if (first_vector_components.empty() && matrix.m() == 0)
3529 {
3530 // Just the identity matrix in case no rotation is specified:
3531 return IdentityMatrix(n_dofs_per_face);
3532 }
3533
3534 // The matrix describes a rotation and we have to build a transformation
3535 // matrix, we assume that for a 0* rotation we would have to build the
3536 // identity matrix
3537
3538 Assert(matrix.m() == spacedim, ExcInternalError());
3539
3540 const Quadrature<dim - 1> quadrature(
3541 fe.get_unit_face_support_points(face_no));
3542
3543 // have an array that stores the location of each vector-dof tuple we want
3544 // to rotate.
3545 using DoFTuple = std::array<unsigned int, spacedim>;
3546
3547 // start with a pristine interpolation matrix...
3548 FullMatrix<double> transformation = IdentityMatrix(n_dofs_per_face);
3549
3550 for (unsigned int i = 0; i < n_dofs_per_face; ++i)
3551 {
3552 std::vector<unsigned int>::const_iterator comp_it =
3553 std::find(first_vector_components.begin(),
3554 first_vector_components.end(),
3555 fe.face_system_to_component_index(i, face_no).first);
3556 if (comp_it != first_vector_components.end())
3557 {
3558 const unsigned int first_vector_component = *comp_it;
3559
3560 // find corresponding other components of vector
3561 DoFTuple vector_dofs;
3562 vector_dofs[0] = i;
3563 unsigned int n_found = 1;
3564
3565 Assert(
3566 *comp_it + spacedim <= fe.n_components(),
3567 ExcMessage(
3568 "Error: the finite element does not have enough components "
3569 "to define rotated periodic boundaries."));
3570
3571 for (unsigned int k = 0; k < n_dofs_per_face; ++k)
3572 if ((k != i) && (quadrature.point(k) == quadrature.point(i)) &&
3573 (fe.face_system_to_component_index(k, face_no).first >=
3574 first_vector_component) &&
3575 (fe.face_system_to_component_index(k, face_no).first <
3576 first_vector_component + spacedim))
3577 {
3578 vector_dofs[fe.face_system_to_component_index(k, face_no)
3579 .first -
3580 first_vector_component] = k;
3581 ++n_found;
3582 if (n_found == dim)
3583 break;
3584 }
3585
3586 // ... and rotate all dofs belonging to vector valued components
3587 // that are selected by first_vector_components:
3588 for (unsigned int i = 0; i < spacedim; ++i)
3589 {
3590 transformation[vector_dofs[i]][vector_dofs[i]] = 0.;
3591 for (unsigned int j = 0; j < spacedim; ++j)
3592 transformation[vector_dofs[i]][vector_dofs[j]] =
3593 matrix[i][j];
3594 }
3595 }
3596 }
3597 return transformation;
3598 }
3599 } /*namespace*/
3600
3601
3602 // Low level interface:
3603
3604
3605 template <typename FaceIterator, typename number>
3606 void
3607 make_periodicity_constraints(
3608 const FaceIterator &face_1,
3609 const std_cxx20::type_identity_t<FaceIterator> &face_2,
3610 AffineConstraints<number> &affine_constraints,
3611 const ComponentMask &component_mask,
3612 const types::geometric_orientation combined_orientation,
3613 const FullMatrix<double> &matrix,
3614 const std::vector<unsigned int> &first_vector_components,
3615 const number periodicity_factor)
3616 {
3617 static const int dim = FaceIterator::AccessorType::dimension;
3618 static const int spacedim = FaceIterator::AccessorType::space_dimension;
3619
3620 if constexpr (running_in_debug_mode())
3621 {
3622 const auto [orientation, rotation, flip] =
3623 ::internal::split_face_orientation(combined_orientation);
3624
3625 Assert((dim != 1) ||
3626 (orientation == true && flip == false && rotation == false),
3627 ExcMessage(
3628 "The supplied orientation (orientation, rotation, flip) "
3629 "is invalid for 1d"));
3630
3631 Assert((dim != 2) || (flip == false && rotation == false),
3632 ExcMessage(
3633 "The supplied orientation (orientation, rotation, flip) "
3634 "is invalid for 2d"));
3635
3636 Assert(face_1 != face_2,
3637 ExcMessage("face_1 and face_2 are equal! Cannot constrain DoFs "
3638 "on the very same face"));
3639
3640 Assert(face_1->at_boundary() && face_2->at_boundary(),
3641 ExcMessage("Faces for periodicity constraints must be on the "
3642 "boundary"));
3643
3644 Assert(matrix.m() == matrix.n(),
3645 ExcMessage(
3646 "The supplied (rotation or interpolation) matrix must "
3647 "be a square matrix"));
3648
3649 Assert(first_vector_components.empty() || matrix.m() == spacedim,
3650 ExcMessage("first_vector_components is nonempty, so matrix must "
3651 "be a rotation matrix exactly of size spacedim"));
3652
3653 if (!face_1->has_children())
3654 {
3655 // TODO: the implementation makes the assumption that all faces have
3656 // the same number of dofs
3657 AssertDimension(
3658 face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(),
3659 1);
3660 const unsigned int face_no = 0;
3661
3662 Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
3663 const unsigned int n_dofs_per_face =
3664 face_1->get_fe(face_1->nth_active_fe_index(0))
3665 .n_dofs_per_face(face_no);
3666
3667 Assert(matrix.m() == 0 ||
3668 (first_vector_components.empty() &&
3669 matrix.m() == n_dofs_per_face) ||
3670 (!first_vector_components.empty() &&
3671 matrix.m() == spacedim),
3672 ExcMessage(
3673 "The matrix must have either size 0 or spacedim "
3674 "(if first_vector_components is nonempty) "
3675 "or the size must be equal to the # of DoFs on the face "
3676 "(if first_vector_components is empty)."));
3677 }
3678
3679 if (!face_2->has_children())
3680 {
3681 // TODO: the implementation makes the assumption that all faces have
3682 // the same number of dofs
3683 AssertDimension(
3684 face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(),
3685 1);
3686 const unsigned int face_no = 0;
3687
3688 Assert(face_2->n_active_fe_indices() == 1, ExcInternalError());
3689 const unsigned int n_dofs_per_face =
3690 face_2->get_fe(face_2->nth_active_fe_index(0))
3691 .n_dofs_per_face(face_no);
3692
3693 Assert(matrix.m() == 0 ||
3694 (first_vector_components.empty() &&
3695 matrix.m() == n_dofs_per_face) ||
3696 (!first_vector_components.empty() &&
3697 matrix.m() == spacedim),
3698 ExcMessage(
3699 "The matrix must have either size 0 or spacedim "
3700 "(if first_vector_components is nonempty) "
3701 "or the size must be equal to the # of DoFs on the face "
3702 "(if first_vector_components is empty)."));
3703 }
3704 }
3705
3706 if (face_1->has_children() && face_2->has_children())
3707 {
3708 // In the case that both faces have children, we loop over all children
3709 // and apply make_periodicity_constraints() recursively:
3710
3711 Assert(face_1->n_children() ==
3712 GeometryInfo<dim>::max_children_per_face &&
3713 face_2->n_children() ==
3714 GeometryInfo<dim>::max_children_per_face,
3715 ExcNotImplemented());
3716
3717 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face;
3718 ++i)
3719 {
3720 // We need to access the subface indices without knowing the face
3721 // number. Hence, we pick the lowest-value face: i.e., face 2 in 2D
3722 // has subfaces {0, 1} and face 4 in 3D has subfaces {0, 1, 2, 3}.
3723 const unsigned int face_no = dim == 2 ? 2 : 4;
3724
3725 // Lookup the index for the second face. Like the assertions above,
3726 // this is only presently valid for hypercube meshes.
3727 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
3728 const unsigned int j =
3729 reference_cell.child_cell_on_face(face_no,
3730 i,
3731 combined_orientation);
3732
3733 make_periodicity_constraints(face_1->child(i),
3734 face_2->child(j),
3735 affine_constraints,
3736 component_mask,
3737 combined_orientation,
3738 matrix,
3739 first_vector_components,
3740 periodicity_factor);
3741 }
3742 }
3743 else
3744 {
3745 // Otherwise at least one of the two faces is active and we need to do
3746 // some work and enter the constraints!
3747
3748 // The finite element that matters is the one on the active face:
3749 const FiniteElement<dim, spacedim> &fe =
3750 face_1->has_children() ?
3751 face_2->get_fe(face_2->nth_active_fe_index(0)) :
3752 face_1->get_fe(face_1->nth_active_fe_index(0));
3753
3754 // TODO: the implementation makes the assumption that all faces have the
3755 // same number of dofs
3756 AssertDimension(fe.n_unique_faces(), 1);
3757 const unsigned int face_no = 0;
3758
3759 const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no);
3760
3761 // Sometimes we just have nothing to do (for all finite elements, or
3762 // systems which accidentally don't have any dofs on the boundary).
3763 if (n_dofs_per_face == 0)
3764 return;
3765
3766 const FullMatrix<double> transformation =
3767 compute_transformation(fe, matrix, first_vector_components);
3768
3769 if (!face_2->has_children())
3770 {
3771 // Performance hack: We do not need to compute an inverse if the
3772 // matrix is the identity matrix.
3773 if (first_vector_components.empty() && matrix.m() == 0)
3774 {
3775 internal::set_periodicity_constraints(face_2,
3776 face_1,
3777 transformation,
3778 affine_constraints,
3779 component_mask,
3780 combined_orientation,
3781 periodicity_factor);
3782 }
3783 else
3784 {
3785 FullMatrix<double> inverse(transformation.m());
3786 inverse.invert(transformation);
3787
3788 internal::set_periodicity_constraints(face_2,
3789 face_1,
3790 inverse,
3791 affine_constraints,
3792 component_mask,
3793 combined_orientation,
3794 periodicity_factor);
3795 }
3796 }
3797 else
3798 {
3799 Assert(!face_1->has_children(), ExcInternalError());
3800
3801 // since the combined_orientation describes how we should rotate
3802 // face_1 to match face_2, we invert it to match the convention
3803 // expected by set_periodicity_constraints()
3804 const auto face_reference_cell = face_1->reference_cell();
3805 internal::set_periodicity_constraints(
3806 face_1,
3807 face_2,
3808 transformation,
3809 affine_constraints,
3810 component_mask,
3811 face_reference_cell.get_inverse_combined_orientation(
3812 combined_orientation),
3813 periodicity_factor);
3814 }
3815 }
3816 }
3817
3818
3819
3820 template <typename FaceIterator, typename number>
3821 void
3822 make_periodicity_constraints_on_level(
3823 const FaceIterator &face_1,
3824 const std_cxx20::type_identity_t<FaceIterator> &face_2,
3825 const unsigned int level,
3826 AffineConstraints<number> &affine_constraints,
3827 const ComponentMask &component_mask,
3828 const types::geometric_orientation combined_orientation,
3829 const FullMatrix<double> &matrix,
3830 const std::vector<unsigned int> &first_vector_components,
3831 const number periodicity_factor)
3832 {
3833 static const int dim = FaceIterator::AccessorType::dimension;
3834 static const int spacedim = FaceIterator::AccessorType::space_dimension;
3835
3836 if constexpr (running_in_debug_mode())
3837 {
3838 const auto [orientation, rotation, flip] =
3839 ::internal::split_face_orientation(combined_orientation);
3840 Assert((dim != 1) ||
3841 (orientation == true && flip == false && rotation == false),
3842 ExcMessage("The supplied face orientation is invalid for 1d."));
3843 Assert((dim != 2) || (flip == false && rotation == false),
3844 ExcMessage("The supplied face orientation is invalid for 2d."));
3845 }
3846
3847 Assert(face_1 != face_2,
3848 ExcMessage("face_1 and face_2 are equal! Cannot constrain DoFs "
3849 "on the very same face"));
3850 Assert(face_1->at_boundary() && face_2->at_boundary(),
3851 ExcMessage("Faces for periodicity constraints must be on the "
3852 "boundary"));
3853 Assert(&face_1->get_dof_handler() == &face_2->get_dof_handler(),
3854 ExcMessage("The two faces must belong to the same DoFHandler."));
3855 AssertIndexRange(
3856 level, face_1->get_dof_handler().get_triangulation().n_global_levels());
3857 Assert(face_1->get_dof_handler().has_hp_capabilities() == false,
3858 ExcNotImplemented());
3859 Assert(matrix.m() == matrix.n(),
3860 ExcMessage("The supplied rotation or interpolation matrix must "
3861 "be square."));
3862 Assert(first_vector_components.empty() || matrix.m() == spacedim,
3863 ExcMessage("If first_vector_components is nonempty, matrix must "
3864 "be a rotation matrix of size spacedim."));
3865
3866 const FiniteElement<dim, spacedim> &fe = face_1->get_dof_handler().get_fe();
3867 const FullMatrix<double> transformation =
3868 compute_transformation(fe, matrix, first_vector_components);
3869
3870 if (first_vector_components.empty() && matrix.m() == 0)
3871 internal::set_periodicity_constraints(face_2,
3872 face_1,
3873 transformation,
3874 affine_constraints,
3875 component_mask,
3876 combined_orientation,
3877 periodicity_factor,
3878 level);
3879 else
3880 {
3881 FullMatrix<double> inverse(transformation.m());
3882 inverse.invert(transformation);
3883
3884 internal::set_periodicity_constraints(face_2,
3885 face_1,
3886 inverse,
3887 affine_constraints,
3888 component_mask,
3889 combined_orientation,
3890 periodicity_factor,
3891 level);
3892 }
3893 }
3894
3895
3896
3897 template <int dim, int spacedim, typename number>
3898 void
3899 make_periodicity_constraints(
3900 const std::vector<GridTools::PeriodicFacePair<
3901 typename DoFHandler<dim, spacedim>::cell_iterator>> &periodic_faces,
3902 AffineConstraints<number> &constraints,
3903 const ComponentMask &component_mask,
3904 const std::vector<unsigned int> &first_vector_components,
3905 const number periodicity_factor)
3906 {
3907 // Loop over all periodic faces...
3908 for (auto &pair : periodic_faces)
3909 {
3910 using FaceIterator = typename DoFHandler<dim, spacedim>::face_iterator;
3911 const FaceIterator face_1 = pair.cell[0]->face(pair.face_idx[0]);
3912 const FaceIterator face_2 = pair.cell[1]->face(pair.face_idx[1]);
3913
3914 Assert(face_1->at_boundary() && face_2->at_boundary(),
3915 ExcInternalError());
3916
3917 Assert(face_1 != face_2, ExcInternalError());
3918
3919 // ... and apply the low level make_periodicity_constraints function to
3920 // every matching pair:
3921 make_periodicity_constraints(face_1,
3922 face_2,
3923 constraints,
3924 component_mask,
3925 pair.orientation,
3926 pair.matrix,
3927 first_vector_components,
3928 periodicity_factor);
3929 }
3930 }
3931
3932
3933 // High level interface variants:
3934
3935
3936 template <int dim, int spacedim, typename number>
3937 void
3938 make_periodicity_constraints(const DoFHandler<dim, spacedim> &dof_handler,
3939 const types::boundary_id b_id1,
3940 const types::boundary_id b_id2,
3941 const unsigned int direction,
3942 ::AffineConstraints<number> &constraints,
3943 const ComponentMask &component_mask,
3944 const number periodicity_factor)
3945 {
3946 AssertIndexRange(direction, spacedim);
3947
3948 Assert(b_id1 != b_id2,
3949 ExcMessage("The boundary indicators b_id1 and b_id2 must be "
3950 "different to denote different boundaries."));
3951
3952 std::vector<GridTools::PeriodicFacePair<
3953 typename DoFHandler<dim, spacedim>::cell_iterator>>
3954 matched_faces;
3955
3956 // Collect matching periodic cells on the coarsest level:
3957 GridTools::collect_periodic_faces(
3958 dof_handler, b_id1, b_id2, direction, matched_faces);
3959
3960 make_periodicity_constraints<dim, spacedim>(matched_faces,
3961 constraints,
3962 component_mask,
3963 std::vector<unsigned int>(),
3964 periodicity_factor);
3965 }
3966
3967
3968
3969 template <int dim, int spacedim, typename number>
3970 void
3971 make_periodicity_constraints(const DoFHandler<dim, spacedim> &dof_handler,
3972 const types::boundary_id b_id,
3973 const unsigned int direction,
3974 AffineConstraints<number> &constraints,
3975 const ComponentMask &component_mask,
3976 const number periodicity_factor)
3977 {
3978 AssertIndexRange(direction, spacedim);
3979
3980 Assert(dim == spacedim, ExcNotImplemented());
3981
3982 std::vector<GridTools::PeriodicFacePair<
3983 typename DoFHandler<dim, spacedim>::cell_iterator>>
3984 matched_faces;
3985
3986 // Collect matching periodic cells on the coarsest level:
3987 GridTools::collect_periodic_faces(dof_handler,
3988 b_id,
3989 direction,
3990 matched_faces);
3991
3992 make_periodicity_constraints<dim, spacedim>(matched_faces,
3993 constraints,
3994 component_mask,
3995 std::vector<unsigned int>(),
3996 periodicity_factor);
3997 }
3998
3999
4000
4001 namespace internal
4002 {
4003 namespace Assembler
4004 {
4005 // We don't actually need a scratch object, so use an empty class for it.
4006 struct Scratch
4007 {};
4008
4009
4010 template <int dim, int spacedim>
4011 struct CopyData
4012 {
4013 unsigned int dofs_per_cell;
4014 std::vector<types::global_dof_index> parameter_dof_indices;
4015#ifdef DEAL_II_WITH_MPI
4016 std::vector<::LinearAlgebra::distributed::Vector<double>>
4017 global_parameter_representation;
4018#else
4019 std::vector<::Vector<double>> global_parameter_representation;
4020#endif
4021 };
4022 } // namespace Assembler
4023
4024 namespace
4025 {
4031 template <int dim, int spacedim>
4032 void
4033 compute_intergrid_weights_3(
4034 const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
4035 Assembler::CopyData<dim, spacedim> &copy_data,
4036 const unsigned int coarse_component,
4037 const FiniteElement<dim, spacedim> &coarse_fe,
4038 const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map,
4039 const std::vector<::Vector<double>> &parameter_dofs)
4040 {
4041 // for each cell on the parameter grid: find out which degrees of
4042 // freedom on the fine grid correspond in which way to the degrees of
4043 // freedom on the parameter grid
4044 //
4045 // since for continuous FEs some dofs exist on more than one cell, we
4046 // have to track which ones were already visited. the problem is that if
4047 // we visit a dof first on one cell and compute its weight with respect
4048 // to some global dofs to be non-zero, and later visit the dof again on
4049 // another cell and (since we are on another cell) recompute the weights
4050 // with respect to the same dofs as above to be zero now, we have to
4051 // preserve them. we therefore overwrite all weights if they are nonzero
4052 // and do not enforce zero weights since that might be only due to the
4053 // fact that we are on another cell.
4054 //
4055 // example:
4056 // coarse grid
4057 // | | |
4058 // *-----*-----*
4059 // | cell|cell |
4060 // | 1 | 2 |
4061 // | | |
4062 // 0-----1-----*
4063 //
4064 // fine grid
4065 // | | | | |
4066 // *--*--*--*--*
4067 // | | | | |
4068 // *--*--*--*--*
4069 // | | | | |
4070 // *--x--y--*--*
4071 //
4072 // when on cell 1, we compute the weights of dof 'x' to be 1/2 from
4073 // parameter dofs 0 and 1, respectively. however, when later we are on
4074 // cell 2, we again compute the prolongation of shape function 1
4075 // restricted to cell 2 to the global grid and find that the weight of
4076 // global dof 'x' now is zero. however, we should not overwrite the old
4077 // value.
4078 //
4079 // we therefore always only set nonzero values. why adding up is not
4080 // useful: dof 'y' would get weight 1 from parameter dof 1 on both cells
4081 // 1 and 2, but the correct weight is nevertheless only 1.
4082
4083 // vector to hold the representation of a single degree of freedom on
4084 // the coarse grid (for the selected fe) on the fine grid
4085
4086 copy_data.dofs_per_cell = coarse_fe.n_dofs_per_cell();
4087 copy_data.parameter_dof_indices.resize(copy_data.dofs_per_cell);
4088
4089 // get the global indices of the parameter dofs on this parameter grid
4090 // cell
4091 cell->get_dof_indices(copy_data.parameter_dof_indices);
4092
4093 // loop over all dofs on this cell and check whether they are
4094 // interesting for us
4095 for (unsigned int local_dof = 0; local_dof < copy_data.dofs_per_cell;
4096 ++local_dof)
4097 if (coarse_fe.system_to_component_index(local_dof).first ==
4098 coarse_component)
4099 {
4100 // the how-many-th parameter is this on this cell?
4101 const unsigned int local_parameter_dof =
4102 coarse_fe.system_to_component_index(local_dof).second;
4103
4104 copy_data.global_parameter_representation[local_parameter_dof] =
4105 0.;
4106
4107 // distribute the representation of @p{local_parameter_dof} on the
4108 // parameter grid cell
4109 // @p{cell} to the global data space
4110 coarse_to_fine_grid_map[cell]->set_dof_values_by_interpolation(
4111 parameter_dofs[local_parameter_dof],
4112 copy_data.global_parameter_representation[local_parameter_dof]);
4113 }
4114 }
4115
4116
4117
4123 template <int dim, int spacedim>
4124 void
4125 copy_intergrid_weights_3(
4126 const Assembler::CopyData<dim, spacedim> &copy_data,
4127 const unsigned int coarse_component,
4128 const FiniteElement<dim, spacedim> &coarse_fe,
4129 const std::vector<types::global_dof_index> &weight_mapping,
4130 const bool is_called_in_parallel,
4131 std::vector<std::map<types::global_dof_index, float>> &weights)
4132 {
4133 unsigned int pos = 0;
4134 for (unsigned int local_dof = 0; local_dof < copy_data.dofs_per_cell;
4135 ++local_dof)
4136 if (coarse_fe.system_to_component_index(local_dof).first ==
4137 coarse_component)
4138 {
4139 // now that we've got the global representation of each parameter
4140 // dof, we've only got to clobber the non-zero entries in that
4141 // vector and store the result
4142 //
4143 // what we have learned: if entry @p{i} of the global vector holds
4144 // the value @p{v[i]}, then this is the weight with which the
4145 // present dof contributes to @p{i}. there may be several such
4146 // @p{i}s and their weights' sum should be one. Then, @p{v[i]}
4147 // should be equal to @p{\sum_j w_{ij} p[j]} with @p{p[j]} be the
4148 // values of the degrees of freedom on the coarse grid. we can
4149 // thus compute constraints which link the degrees of freedom
4150 // @p{v[i]} on the fine grid to those on the coarse grid,
4151 // @p{p[j]}. Now to use these as real constraints, rather than as
4152 // additional equations, we have to identify representants among
4153 // the @p{i} for each @p{j}. this will be done by simply taking
4154 // the first @p{i} for which @p{w_{ij}==1}.
4155 //
4156 // guard modification of the weights array by a Mutex. since it
4157 // should happen rather rarely that there are several threads
4158 // operating on different intergrid weights, have only one mutex
4159 // for all of them
4160 for (types::global_dof_index i = 0;
4161 i < copy_data.global_parameter_representation[pos].size();
4162 ++i)
4163 // set this weight if it belongs to a parameter dof.
4164 if (weight_mapping[i] != numbers::invalid_dof_index)
4165 {
4166 // only overwrite old value if not by zero
4167 if (copy_data.global_parameter_representation[pos](i) != 0)
4168 {
4169 const types::global_dof_index
4170 wi = copy_data.parameter_dof_indices[local_dof],
4171 wj = weight_mapping[i];
4172 weights[wi][wj] =
4173 copy_data.global_parameter_representation[pos](i);
4174 }
4175 }
4176 else if (!is_called_in_parallel)
4177 {
4178 // Note that when this function operates with distributed
4179 // fine grid, this assertion is switched off since the
4180 // condition does not necessarily hold
4181 Assert(copy_data.global_parameter_representation[pos](i) ==
4182 0,
4183 ExcInternalError());
4184 }
4185
4186 ++pos;
4187 }
4188 }
4189
4190
4191
4197 template <int dim, int spacedim>
4198 void
4199 compute_intergrid_weights_2(
4200 const DoFHandler<dim, spacedim> &coarse_grid,
4201 const unsigned int coarse_component,
4202 const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map,
4203 const std::vector<::Vector<double>> &parameter_dofs,
4204 const std::vector<types::global_dof_index> &weight_mapping,
4205 std::vector<std::map<types::global_dof_index, float>> &weights)
4206 {
4207 Assembler::CopyData<dim, spacedim> copy_data;
4208
4209 unsigned int n_interesting_dofs = 0;
4210 for (unsigned int local_dof = 0;
4211 local_dof < coarse_grid.get_fe().n_dofs_per_cell();
4212 ++local_dof)
4213 if (coarse_grid.get_fe().system_to_component_index(local_dof).first ==
4214 coarse_component)
4215 ++n_interesting_dofs;
4216
4217 copy_data.global_parameter_representation.resize(n_interesting_dofs);
4218
4219 bool is_called_in_parallel = false;
4220 for (std::size_t i = 0;
4221 i < copy_data.global_parameter_representation.size();
4222 ++i)
4223 {
4224#ifdef DEAL_II_WITH_MPI
4225 MPI_Comm communicator = MPI_COMM_SELF;
4226 try
4227 {
4228 const typename ::parallel::TriangulationBase<dim,
4229 spacedim>
4230 &tria = dynamic_cast<const typename ::parallel::
4231 TriangulationBase<dim, spacedim> &>(
4232 coarse_to_fine_grid_map.get_destination_grid()
4233 .get_triangulation());
4234 communicator = tria.get_mpi_communicator();
4235 is_called_in_parallel = true;
4236 }
4237 catch (std::bad_cast &)
4238 {
4239 // Nothing bad happened: the user used serial Triangulation
4240 }
4241
4242
4243 const IndexSet locally_relevant_dofs =
4244 DoFTools::extract_locally_relevant_dofs(
4245 coarse_to_fine_grid_map.get_destination_grid());
4246
4247 copy_data.global_parameter_representation[i].reinit(
4248 coarse_to_fine_grid_map.get_destination_grid()
4249 .locally_owned_dofs(),
4250 locally_relevant_dofs,
4251 communicator);
4252#else
4253 const types::global_dof_index n_fine_dofs = weight_mapping.size();
4254 copy_data.global_parameter_representation[i].reinit(n_fine_dofs);
4255#endif
4256 }
4257
4258 auto worker =
4259 [coarse_component,
4260 &coarse_grid,
4261 &coarse_to_fine_grid_map,
4262 &parameter_dofs](
4263 const typename DoFHandler<dim, spacedim>::active_cell_iterator
4264 &cell,
4265 const Assembler::Scratch &,
4266 Assembler::CopyData<dim, spacedim> &copy_data) {
4267 compute_intergrid_weights_3<dim, spacedim>(cell,
4268 copy_data,
4269 coarse_component,
4270 coarse_grid.get_fe(),
4271 coarse_to_fine_grid_map,
4272 parameter_dofs);
4273 };
4274
4275 auto copier =
4276 [coarse_component,
4277 &coarse_grid,
4278 &weight_mapping,
4279 is_called_in_parallel,
4280 &weights](const Assembler::CopyData<dim, spacedim> &copy_data) {
4281 copy_intergrid_weights_3<dim, spacedim>(copy_data,
4282 coarse_component,
4283 coarse_grid.get_fe(),
4284 weight_mapping,
4285 is_called_in_parallel,
4286 weights);
4287 };
4288
4289 WorkStream::run(coarse_grid.begin_active(),
4290 coarse_grid.end(),
4291 worker,
4292 copier,
4293 Assembler::Scratch(),
4294 copy_data);
4295
4296#ifdef DEAL_II_WITH_MPI
4297 for (std::size_t i = 0;
4298 i < copy_data.global_parameter_representation.size();
4299 ++i)
4300 copy_data.global_parameter_representation[i].update_ghost_values();
4301#endif
4302 }
4303
4304
4305
4311 template <int dim, int spacedim>
4312 unsigned int
4313 compute_intergrid_weights_1(
4314 const DoFHandler<dim, spacedim> &coarse_grid,
4315 const unsigned int coarse_component,
4316 const DoFHandler<dim, spacedim> &fine_grid,
4317 const unsigned int fine_component,
4318 const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map,
4319 std::vector<std::map<types::global_dof_index, float>> &weights,
4320 std::vector<types::global_dof_index> &weight_mapping)
4321 {
4322 // aliases to the finite elements used by the dof handlers:
4323 const FiniteElement<dim, spacedim> &coarse_fe = coarse_grid.get_fe(),
4324 &fine_fe = fine_grid.get_fe();
4325
4326 // global numbers of dofs
4327 const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(),
4328 n_fine_dofs = fine_grid.n_dofs();
4329
4330 // local numbers of dofs
4331 const unsigned int fine_dofs_per_cell = fine_fe.n_dofs_per_cell();
4332
4333 // alias the number of dofs per cell belonging to the coarse_component
4334 // which is to be the restriction of the fine grid:
4335 const unsigned int coarse_dofs_per_cell_component =
4336 coarse_fe
4337 .base_element(
4338 coarse_fe.component_to_base_index(coarse_component).first)
4339 .n_dofs_per_cell();
4340
4341
4342 // Try to find out whether the grids stem from the same coarse grid.
4343 // This is a rather crude test, but better than nothing
4344 Assert(coarse_grid.get_triangulation().n_cells(0) ==
4345 fine_grid.get_triangulation().n_cells(0),
4346 ExcGridsDontMatch());
4347
4348 // check whether the map correlates the right objects
4349 Assert(&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid,
4350 ExcGridsDontMatch());
4351 Assert(&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid,
4352 ExcGridsDontMatch());
4353
4354
4355 // check whether component numbers are valid
4356 AssertIndexRange(coarse_component, coarse_fe.n_components());
4357 AssertIndexRange(fine_component, fine_fe.n_components());
4358
4359 // check whether respective finite elements are equal
4360 Assert(coarse_fe.base_element(
4361 coarse_fe.component_to_base_index(coarse_component).first) ==
4362 fine_fe.base_element(
4363 fine_fe.component_to_base_index(fine_component).first),
4364 ExcFiniteElementsDontMatch());
4365
4366 if constexpr (running_in_debug_mode())
4367 {
4368 // if in debug mode, check whether the coarse grid is indeed coarser
4369 // everywhere than the fine grid
4370 for (const auto &cell : coarse_grid.active_cell_iterators())
4371 Assert(cell->level() <= coarse_to_fine_grid_map[cell]->level(),
4372 ExcGridNotCoarser());
4373 }
4374
4375 /*
4376 * From here on: the term `parameter' refers to the selected component
4377 * on the coarse grid and its analogon on the fine grid. The naming of
4378 * variables containing this term is due to the fact that
4379 * `selected_component' is longer, but also due to the fact that the
4380 * code of this function was initially written for a program where the
4381 * component which we wanted to match between grids was actually the
4382 * `parameter' variable.
4383 *
4384 * Likewise, the terms `parameter grid' and `state grid' refer to the
4385 * coarse and fine grids, respectively.
4386 *
4387 * Changing the names of variables would in principle be a good idea,
4388 * but would not make things simpler and would be another source of
4389 * errors. If anyone feels like doing so: patches would be welcome!
4390 */
4391
4392
4393
4394 // set up vectors of cell-local data; each vector represents one degree
4395 // of freedom of the coarse-grid variable in the fine-grid element
4396 std::vector<::Vector<double>> parameter_dofs(
4397 coarse_dofs_per_cell_component,
4398 ::Vector<double>(fine_dofs_per_cell));
4399 // for each coarse dof: find its position within the fine element and
4400 // set this value to one in the respective vector (all other values are
4401 // zero by construction)
4402 for (unsigned int local_coarse_dof = 0;
4403 local_coarse_dof < coarse_dofs_per_cell_component;
4404 ++local_coarse_dof)
4405 for (unsigned int fine_dof = 0; fine_dof < fine_fe.n_dofs_per_cell();
4406 ++fine_dof)
4407 if (fine_fe.system_to_component_index(fine_dof) ==
4408 std::make_pair(fine_component, local_coarse_dof))
4409 {
4410 parameter_dofs[local_coarse_dof](fine_dof) = 1.;
4411 break;
4412 }
4413
4414
4415 // find out how many DoFs there are on the grids belonging to the
4416 // components we want to match
4417 unsigned int n_parameters_on_fine_grid = 0;
4418 {
4419 // have a flag for each dof on the fine grid and set it to true if
4420 // this is an interesting dof. finally count how many true's there
4421 std::vector<bool> dof_is_interesting(fine_grid.n_dofs(), false);
4422 std::vector<types::global_dof_index> local_dof_indices(
4423 fine_fe.n_dofs_per_cell());
4424
4425 for (const auto &cell : fine_grid.active_cell_iterators() |
4426 IteratorFilters::LocallyOwnedCell())
4427 {
4428 cell->get_dof_indices(local_dof_indices);
4429 for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i)
4430 if (fine_fe.system_to_component_index(i).first ==
4431 fine_component)
4432 dof_is_interesting[local_dof_indices[i]] = true;
4433 }
4434
4435 n_parameters_on_fine_grid = std::count(dof_is_interesting.begin(),
4436 dof_is_interesting.end(),
4437 true);
4438 }
4439
4440
4441 // set up the weights mapping
4442 weights.clear();
4443 weights.resize(n_coarse_dofs);
4444
4445 weight_mapping.clear();
4446 weight_mapping.resize(n_fine_dofs, numbers::invalid_dof_index);
4447
4448 {
4449 std::vector<types::global_dof_index> local_dof_indices(
4450 fine_fe.n_dofs_per_cell());
4451 unsigned int next_free_index = 0;
4452 for (const auto &cell : fine_grid.active_cell_iterators() |
4453 IteratorFilters::LocallyOwnedCell())
4454 {
4455 cell->get_dof_indices(local_dof_indices);
4456 for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i)
4457 // if this DoF is a parameter dof and has not yet been
4458 // numbered, then do so
4459 if ((fine_fe.system_to_component_index(i).first ==
4460 fine_component) &&
4461 (weight_mapping[local_dof_indices[i]] ==
4462 numbers::invalid_dof_index))
4463 {
4464 weight_mapping[local_dof_indices[i]] = next_free_index;
4465 ++next_free_index;
4466 }
4467 }
4468
4469 Assert(next_free_index == n_parameters_on_fine_grid,
4470 ExcInternalError());
4471 }
4472
4473
4474 // for each cell on the parameter grid: find out which degrees of
4475 // freedom on the fine grid correspond in which way to the degrees of
4476 // freedom on the parameter grid
4477 //
4478 // do this in a separate function to allow for multithreading there. see
4479 // this function also if you want to read more information on the
4480 // algorithm used.
4481 compute_intergrid_weights_2(coarse_grid,
4482 coarse_component,
4483 coarse_to_fine_grid_map,
4484 parameter_dofs,
4485 weight_mapping,
4486 weights);
4487
4488
4489 // ok, now we have all weights for each dof on the fine grid. if in
4490 // debug mode lets see if everything went smooth, i.e. each dof has sum
4491 // of weights one
4492 //
4493 // in other words this means that if the sum of all shape functions on
4494 // the parameter grid is one (which is always the case), then the
4495 // representation on the state grid should be as well (division of
4496 // unity)
4497 //
4498 // if the parameter grid has more than one component, then the
4499 // respective dofs of the other components have sum of weights zero, of
4500 // course. we do not explicitly ask which component a dof belongs to,
4501 // but this at least tests some errors
4502 if constexpr (running_in_debug_mode())
4503 {
4504 for (unsigned int col = 0; col < n_parameters_on_fine_grid; ++col)
4505 {
4506 double sum = 0;
4507 for (types::global_dof_index row = 0; row < n_coarse_dofs;
4508 ++row)
4509 if (weights[row].find(col) != weights[row].end())
4510 sum += weights[row][col];
4511 Assert((std::fabs(sum - 1) < 1.e-12) ||
4512 ((coarse_fe.n_components() > 1) && (sum == 0)),
4513 ExcInternalError());
4514 }
4515 }
4516
4517
4518 return n_parameters_on_fine_grid;
4519 }
4520
4521
4522 } // namespace
4523 } // namespace internal
4524
4525
4526
4527 template <int dim, int spacedim>
4528 void
4529 compute_intergrid_constraints(
4530 const DoFHandler<dim, spacedim> &coarse_grid,
4531 const unsigned int coarse_component,
4532 const DoFHandler<dim, spacedim> &fine_grid,
4533 const unsigned int fine_component,
4534 const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map,
4535 AffineConstraints<double> &constraints)
4536 {
4537 Assert(coarse_grid.get_fe_collection().size() == 1 &&
4538 fine_grid.get_fe_collection().size() == 1,
4539 ExcMessage("This function is not yet implemented for DoFHandlers "
4540 "using hp-capabilities."));
4541 // store the weights with which a dof on the parameter grid contributes to a
4542 // dof on the fine grid. see the long doc below for more info
4543 //
4544 // allocate as many rows as there are parameter dofs on the coarse grid and
4545 // as many columns as there are parameter dofs on the fine grid.
4546 //
4547 // weight_mapping is used to map the global (fine grid) parameter dof
4548 // indices to the columns
4549 //
4550 // in the original implementation, the weights array was actually of
4551 // FullMatrix<double> type. this wasted huge amounts of memory, but was
4552 // fast. nonetheless, since the memory consumption was quadratic in the
4553 // number of degrees of freedom, this was not very practical, so we now use
4554 // a vector of rows of the matrix, and in each row a vector of pairs
4555 // (colnum,value). this seems like the best tradeoff between memory and
4556 // speed, as it is now linear in memory and still fast enough.
4557 //
4558 // to save some memory and since the weights are usually (negative) powers
4559 // of 2, we choose the value type of the matrix to be @p{float} rather than
4560 // @p{double}.
4561 std::vector<std::map<types::global_dof_index, float>> weights;
4562
4563 // this is this mapping. there is one entry for each dof on the fine grid;
4564 // if it is a parameter dof, then its value is the column in weights for
4565 // that parameter dof, if it is any other dof, then its value is -1,
4566 // indicating an error
4567 std::vector<types::global_dof_index> weight_mapping;
4568
4569 const unsigned int n_parameters_on_fine_grid =
4570 internal::compute_intergrid_weights_1(coarse_grid,
4571 coarse_component,
4572 fine_grid,
4573 fine_component,
4574 coarse_to_fine_grid_map,
4575 weights,
4576 weight_mapping);
4577 (void)n_parameters_on_fine_grid;
4578
4579 // global numbers of dofs
4580 const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(),
4581 n_fine_dofs = fine_grid.n_dofs();
4582
4583
4584 // get an array in which we store which dof on the coarse grid is a
4585 // parameter and which is not
4586 IndexSet coarse_dof_is_parameter;
4587 {
4588 std::vector<bool> mask(coarse_grid.get_fe(0).n_components(), false);
4589 mask[coarse_component] = true;
4590
4591 coarse_dof_is_parameter =
4592 extract_dofs<dim, spacedim>(coarse_grid, ComponentMask(mask));
4593 }
4594
4595 // now we know that the weights in each row constitute a constraint. enter
4596 // this into the constraints object
4597 //
4598 // first task: for each parameter dof on the parameter grid, find a
4599 // representant on the fine, global grid. this is possible since we use
4600 // conforming finite element. we take this representant to be the first
4601 // element in this row with weight identical to one. the representant will
4602 // become an unconstrained degree of freedom, while all others will be
4603 // constrained to this dof (and possibly others)
4604 std::vector<types::global_dof_index> representants(
4605 n_coarse_dofs, numbers::invalid_dof_index);
4606 for (types::global_dof_index parameter_dof = 0;
4607 parameter_dof < n_coarse_dofs;
4608 ++parameter_dof)
4609 if (coarse_dof_is_parameter.is_element(parameter_dof))
4610 {
4611 // if this is the line of a parameter dof on the coarse grid, then it
4612 // should have at least one dependent node on the fine grid
4613 Assert(weights[parameter_dof].size() > 0, ExcInternalError());
4614
4615 // find the column where the representant is mentioned
4616 std::map<types::global_dof_index, float>::const_iterator i =
4617 weights[parameter_dof].begin();
4618 for (; i != weights[parameter_dof].end(); ++i)
4619 if (i->second == 1)
4620 break;
4621 Assert(i != weights[parameter_dof].end(), ExcInternalError());
4622 const types::global_dof_index column = i->first;
4623
4624 // now we know in which column of weights the representant is, but we
4625 // don't know its global index. get it using the inverse operation of
4626 // the weight_mapping
4627 types::global_dof_index global_dof = 0;
4628 for (; global_dof < weight_mapping.size(); ++global_dof)
4629 if (weight_mapping[global_dof] ==
4630 static_cast<types::global_dof_index>(column))
4631 break;
4632 Assert(global_dof < weight_mapping.size(), ExcInternalError());
4633
4634 // now enter the representants global index into our list
4635 representants[parameter_dof] = global_dof;
4636 }
4637 else
4638 {
4639 // consistency check: if this is no parameter dof on the coarse grid,
4640 // then the respective row must be empty!
4641 Assert(weights[parameter_dof].empty(), ExcInternalError());
4642 }
4643
4644
4645
4646 // note for people that want to optimize this function: the largest part of
4647 // the computing time is spent in the following, rather innocent block of
4648 // code. basically, it must be the AffineConstraints::add_entry call which
4649 // takes the bulk of the time, but it is not known to the author how to make
4650 // it faster...
4651 std::vector<std::pair<types::global_dof_index, double>> constraint_line;
4652 for (types::global_dof_index global_dof = 0; global_dof < n_fine_dofs;
4653 ++global_dof)
4654 if (weight_mapping[global_dof] != numbers::invalid_dof_index)
4655 // this global dof is a parameter dof, so it may carry a constraint note
4656 // that for each global dof, the sum of weights shall be one, so we can
4657 // find out whether this dof is constrained in the following way: if the
4658 // only weight in this row is a one, and the representant for the
4659 // parameter dof of the line in which this one is is the present dof,
4660 // then we consider this dof to be unconstrained. otherwise, all other
4661 // dofs are constrained
4662 {
4663 const types::global_dof_index col = weight_mapping[global_dof];
4664 Assert(col < n_parameters_on_fine_grid, ExcInternalError());
4665
4666 types::global_dof_index first_used_row = 0;
4667
4668 {
4669 Assert(weights.size() > 0, ExcInternalError());
4670 std::map<types::global_dof_index, float>::const_iterator col_entry =
4671 weights[0].end();
4672 for (; first_used_row < n_coarse_dofs; ++first_used_row)
4673 {
4674 col_entry = weights[first_used_row].find(col);
4675 if (col_entry != weights[first_used_row].end())
4676 break;
4677 }
4678
4679 Assert(col_entry != weights[first_used_row].end(),
4680 ExcInternalError());
4681
4682 if ((col_entry->second == 1) &&
4683 (representants[first_used_row] == global_dof))
4684 // dof unconstrained or constrained to itself (in case this cell
4685 // is mapped to itself, rather than to children of itself)
4686 continue;
4687 }
4688
4689
4690 // otherwise enter all constraints
4691 constraint_line.clear();
4692 for (types::global_dof_index row = first_used_row;
4693 row < n_coarse_dofs;
4694 ++row)
4695 {
4696 const std::map<types::global_dof_index, float>::const_iterator j =
4697 weights[row].find(col);
4698 if ((j != weights[row].end()) && (j->second != 0))
4699 constraint_line.emplace_back(representants[row], j->second);
4700 }
4701
4702 constraints.add_constraint(global_dof, constraint_line, 0.);
4703 }
4704 }
4705
4706
4707
4708 template <int dim, int spacedim>
4709 void
4710 compute_intergrid_transfer_representation(
4711 const DoFHandler<dim, spacedim> &coarse_grid,
4712 const unsigned int coarse_component,
4713 const DoFHandler<dim, spacedim> &fine_grid,
4714 const unsigned int fine_component,
4715 const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map,
4716 std::vector<std::map<types::global_dof_index, float>>
4717 &transfer_representation)
4718 {
4719 Assert(coarse_grid.get_fe_collection().size() == 1 &&
4720 fine_grid.get_fe_collection().size() == 1,
4721 ExcMessage("This function is not yet implemented for DoFHandlers "
4722 "using hp-capabilities."));
4723 // store the weights with which a dof on the parameter grid contributes to a
4724 // dof on the fine grid. see the long doc below for more info
4725 //
4726 // allocate as many rows as there are parameter dofs on the coarse grid and
4727 // as many columns as there are parameter dofs on the fine grid.
4728 //
4729 // weight_mapping is used to map the global (fine grid) parameter dof
4730 // indices to the columns
4731 //
4732 // in the original implementation, the weights array was actually of
4733 // FullMatrix<double> type. this wasted huge amounts of memory, but was
4734 // fast. nonetheless, since the memory consumption was quadratic in the
4735 // number of degrees of freedom, this was not very practical, so we now use
4736 // a vector of rows of the matrix, and in each row a vector of pairs
4737 // (colnum,value). this seems like the best tradeoff between memory and
4738 // speed, as it is now linear in memory and still fast enough.
4739 //
4740 // to save some memory and since the weights are usually (negative) powers
4741 // of 2, we choose the value type of the matrix to be @p{float} rather than
4742 // @p{double}.
4743 std::vector<std::map<types::global_dof_index, float>> weights;
4744
4745 // this is this mapping. there is one entry for each dof on the fine grid;
4746 // if it is a parameter dof, then its value is the column in weights for
4747 // that parameter dof, if it is any other dof, then its value is -1,
4748 // indicating an error
4749 std::vector<types::global_dof_index> weight_mapping;
4750
4751 internal::compute_intergrid_weights_1(coarse_grid,
4752 coarse_component,
4753 fine_grid,
4754 fine_component,
4755 coarse_to_fine_grid_map,
4756 weights,
4757 weight_mapping);
4758
4759 // now compute the requested representation
4760 const types::global_dof_index n_global_parm_dofs =
4761 std::count_if(weight_mapping.begin(),
4762 weight_mapping.end(),
4763 [](const types::global_dof_index dof) {
4764 return dof != numbers::invalid_dof_index;
4765 });
4766
4767 // first construct the inverse mapping of weight_mapping
4768 std::vector<types::global_dof_index> inverse_weight_mapping(
4769 n_global_parm_dofs, numbers::invalid_dof_index);
4770 for (types::global_dof_index i = 0; i < weight_mapping.size(); ++i)
4771 {
4772 const types::global_dof_index parameter_dof = weight_mapping[i];
4773 // if this global dof is a parameter
4774 if (parameter_dof != numbers::invalid_dof_index)
4775 {
4776 Assert(parameter_dof < n_global_parm_dofs, ExcInternalError());
4777 Assert((inverse_weight_mapping[parameter_dof] ==
4778 numbers::invalid_dof_index),
4779 ExcInternalError());
4780
4781 inverse_weight_mapping[parameter_dof] = i;
4782 }
4783 }
4784
4785 // next copy over weights array and replace respective numbers
4786 const types::global_dof_index n_rows = weight_mapping.size();
4787
4788 transfer_representation.clear();
4789 transfer_representation.resize(n_rows);
4790
4791 const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs();
4792 for (types::global_dof_index i = 0; i < n_coarse_dofs; ++i)
4793 {
4794 std::map<types::global_dof_index, float>::const_iterator j =
4795 weights[i].begin();
4796 for (; j != weights[i].end(); ++j)
4797 {
4798 const types::global_dof_index p = inverse_weight_mapping[j->first];
4799 Assert(p < n_rows, ExcInternalError());
4800
4801 transfer_representation[p][i] = j->second;
4802 }
4803 }
4804 }
4805
4806
4807
4808 template <int dim, int spacedim, typename number>
4809 void
4810 make_zero_boundary_constraints(
4811 const DoFHandler<dim, spacedim> &dof,
4812 const types::boundary_id boundary_id,
4813 AffineConstraints<number> &zero_boundary_constraints,
4814 const ComponentMask &component_mask)
4815 {
4816 Assert(component_mask.represents_n_components(dof.get_fe(0).n_components()),
4817 ExcMessage("The number of components in the mask has to be either "
4818 "zero or equal to the number of components in the finite "
4819 "element."));
4820
4821 const unsigned int n_components = dof.get_fe_collection().n_components();
4822
4823 Assert(component_mask.n_selected_components(n_components) > 0,
4824 ComponentMask::ExcNoComponentSelected());
4825
4826 // a field to store the indices on the face
4827 std::vector<types::global_dof_index> face_dofs;
4828 face_dofs.reserve(dof.get_fe_collection().max_dofs_per_face());
4829 // a field to store the indices on the cell
4830 std::vector<types::global_dof_index> cell_dofs;
4831 cell_dofs.reserve(dof.get_fe_collection().max_dofs_per_cell());
4832
4833 // In looping over faces, we will encounter some DoFs multiple
4834 // times (namely, the ones on vertices and (in 3d) edges shared
4835 // between multiple boundary faces. Keep track of which DoFs we
4836 // have already encountered, so that we do not have to consider
4837 // them a second time.
4838 std::set<types::global_dof_index> dofs_already_treated;
4839
4840 for (const auto &cell : dof.active_cell_iterators())
4841 if (!cell->is_artificial() && cell->at_boundary())
4842 {
4843 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
4844
4845 // get global indices of dofs on the cell
4846 cell_dofs.resize(fe.n_dofs_per_cell());
4847 cell->get_dof_indices(cell_dofs);
4848
4849 for (const auto face_no : cell->face_indices())
4850 {
4851 const typename DoFHandler<dim, spacedim>::face_iterator face =
4852 cell->face(face_no);
4853
4854 // if face is on the boundary and satisfies the correct boundary
4855 // id property
4856 if (face->at_boundary() &&
4857 ((boundary_id == numbers::invalid_boundary_id) ||
4858 (face->boundary_id() == boundary_id)))
4859 {
4860 // get indices and physical location on this face
4861 face_dofs.resize(fe.n_dofs_per_face(face_no));
4862 face->get_dof_indices(face_dofs, cell->active_fe_index());
4863
4864 // enter those dofs into the list that match the component
4865 // signature.
4866 for (const types::global_dof_index face_dof : face_dofs)
4867 if (dofs_already_treated.find(face_dof) ==
4868 dofs_already_treated.end())
4869 {
4870 // Find out if a dof has a contribution in this
4871 // component, and if so, add it to the list
4872 const std::vector<types::global_dof_index>::iterator
4873 it_index_on_cell = std::find(cell_dofs.begin(),
4874 cell_dofs.end(),
4875 face_dof);
4876 Assert(it_index_on_cell != cell_dofs.end(),
4877 ExcInvalidIterator());
4878 const unsigned int index_on_cell =
4879 std::distance(cell_dofs.begin(), it_index_on_cell);
4880 const ComponentMask &nonzero_component_array =
4881 cell->get_fe().get_nonzero_components(index_on_cell);
4882
4883 bool nonzero = false;
4884 for (unsigned int c = 0; c < n_components; ++c)
4885 if (nonzero_component_array[c] && component_mask[c])
4886 {
4887 nonzero = true;
4888 break;
4889 }
4890
4891 if (nonzero)
4892 {
4893 // Check that either (i) the DoF is not
4894 // yet constrained, or (ii) if it is, its
4895 // inhomogeneity is zero:
4896 if (zero_boundary_constraints.is_constrained(
4897 face_dof) == false)
4898 zero_boundary_constraints.constrain_dof_to_zero(
4899 face_dof);
4900 else
4901 Assert(zero_boundary_constraints
4902 .is_inhomogeneously_constrained(
4903 face_dof) == false,
4904 ExcInternalError());
4905 }
4906
4907 // We already dealt with this DoF. Make sure we
4908 // don't touch it again.
4909 dofs_already_treated.insert(face_dof);
4910 }
4911 }
4912 }
4913 }
4914 }
4915
4916
4917
4918 template <int dim, int spacedim, typename number>
4919 void
4920 make_zero_boundary_constraints(
4921 const DoFHandler<dim, spacedim> &dof,
4922 AffineConstraints<number> &zero_boundary_constraints,
4923 const ComponentMask &component_mask)
4924 {
4925 make_zero_boundary_constraints(dof,
4926 numbers::invalid_boundary_id,
4927 zero_boundary_constraints,
4928 component_mask);
4929 }
4930
4931
4932} // end of namespace DoFTools
4933
4934
4935
4936// explicit instantiations
4937
4938#include "dofs/dof_tools_constraints.inst"
4939
4940
4941
4942DEAL_II_NAMESPACE_CLOSE
end
*  iterator end()
Definition ALE_Finite_Strain_Plasticity.h:10398
affine_constraints.h
AffineConstraints::add_line
void add_line(const size_type line_n)
Definition affine_constraints.h:2423
AffineConstraints::add_constraint
void add_constraint(const size_type constrained_dof, const ArrayView< const std::pair< size_type, number > > &dependencies, const number inhomogeneity=0)
AffineConstraints::add_entry
void add_entry(const size_type constrained_dof_index, const size_type column, const number weight)
Definition affine_constraints.h:2455
AffineConstraints::get_local_lines
const IndexSet & get_local_lines() const
Definition affine_constraints.h:2652
AffineConstraints::set_inhomogeneity
void set_inhomogeneity(const size_type constrained_dof_index, const number value)
Definition affine_constraints.h:2497
AffineConstraints::is_constrained
bool is_constrained(const size_type line_n) const
Definition affine_constraints.h:2558
AffineConstraints::get_constraint_entries
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
Definition affine_constraints.h:2588
AffineConstraints::constrain_dof_to_zero
void constrain_dof_to_zero(const size_type constrained_dof)
ComponentMask::represents_n_components
bool represents_n_components(const unsigned int n) const
ComponentMask::n_selected_components
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
DoFHandler::get_fe_collection
const hp::FECollection< dim, spacedim > & get_fe_collection() const
DoFHandler::has_active_dofs
bool has_active_dofs() const
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
DoFHandler::has_hp_capabilities
bool has_hp_capabilities() const
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
FiniteElementData::n_dofs_per_vertex
unsigned int n_dofs_per_vertex() const
FiniteElementData::degree
const unsigned int degree
Definition fe_data.h:450
FiniteElementData::n_dofs_per_cell
unsigned int n_dofs_per_cell() const
FiniteElementData::n_dofs_per_line
unsigned int n_dofs_per_line() const
FiniteElementData::n_dofs_per_face
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
FiniteElementData::n_components
unsigned int n_components() const
FiniteElementData::n_unique_faces
unsigned int n_unique_faces() const
FiniteElementData::n_dofs_per_quad
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
FiniteElementData::dofs_per_cell
const unsigned int dofs_per_cell
Definition fe_data.h:434
FiniteElement::get_name
virtual std::string get_name() const =0
FiniteElement::base_element
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
FiniteElement::component_to_base_index
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
FiniteElement::get_unit_face_support_points
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
FiniteElement::get_subface_interpolation_matrix
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const
FiniteElement::system_to_component_index
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
FiniteElement::constraints
const FullMatrix< double > & constraints(const ::internal::SubfaceCase< dim > &subface_case=::internal::SubfaceCase< dim >::case_isotropic) const
FiniteElement::face_to_cell_index
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const types::geometric_orientation combined_orientation=numbers::default_geometric_orientation) const
FiniteElement::get_face_interpolation_matrix
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const
FiniteElement::face_system_to_component_index
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
FullMatrix::mmult
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
FullMatrix::n
size_type n() const
FullMatrix::invert
void invert(const FullMatrix< number2 > &M)
FullMatrix::m
size_type m() const
IndexSet::size
size_type size() const
Definition index_set.h:1759
IndexSet::is_element
bool is_element(const size_type index) const
Definition index_set.h:1877
MPI_Comm
Triangulation::get_anisotropic_refinement_flag
bool get_anisotropic_refinement_flag() const
Triangulation::n_cells
unsigned int n_cells() const
Triangulation::compute_line_to_adjacent_cells_map
void compute_line_to_adjacent_cells_map()
hp::Collection::size
unsigned int size() const
Definition collection.h:314
hp::FECollection::find_dominating_fe_extended
unsigned int find_dominating_fe_extended(const std::set< unsigned int > &fes, const unsigned int codim=0) const
Definition fe_collection.cc:293
hp::FECollection::hp_constraints_are_implemented
bool hp_constraints_are_implemented() const
Definition fe_collection.h:1054
hp::FECollection::max_dofs_per_face
unsigned int max_dofs_per_face() const
Definition fe_collection.h:1025
hp::FECollection::n_components
unsigned int n_components() const
Definition fe_collection.h:903
hp::FECollection::max_dofs_per_cell
unsigned int max_dofs_per_cell() const
Definition fe_collection.h:1040
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
running_in_debug_mode
constexpr bool running_in_debug_mode()
Definition config.h:76
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
dof_accessor.h
dof_handler.h
dof_tools.h
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exception_macros.h:781
fe_values.h
fe_collection.h
fe_tools.h
filtered_iterator.h
level
unsigned int level
Definition grid_out.cc:4642
grid_tools.h
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
StandardExceptions::ExcInvalidIterator
static ::ExceptionBase & ExcInvalidIterator()
DoFTools::ExcGridNotCoarser
static ::ExceptionBase & ExcGridNotCoarser()
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
DoFTools::ExcFiniteElementsDontMatch
static ::ExceptionBase & ExcFiniteElementsDontMatch()
ComponentMask::ExcNoComponentSelected
static ::ExceptionBase & ExcNoComponentSelected()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
DoFTools::ExcGridsDontMatch
static ::ExceptionBase & ExcGridsDontMatch()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition dof_handler.h:475
DoFHandler::line_iterator
typename ActiveSelector::line_iterator line_iterator
Definition dof_handler.h:366
DoFHandler::face_iterator
typename ActiveSelector::face_iterator face_iterator
Definition dof_handler.h:493
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition dof_handler.h:447
DoFTools::compute_intergrid_transfer_representation
void compute_intergrid_transfer_representation(const DoFHandler< dim, spacedim > &coarse_grid, const unsigned int coarse_component, const DoFHandler< dim, spacedim > &fine_grid, const unsigned int fine_component, const InterGridMap< DoFHandler< dim, spacedim > > &coarse_to_fine_grid_map, std::vector< std::map< types::global_dof_index, float > > &transfer_representation)
Definition dof_tools_constraints.cc:4710
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
Definition dof_tools_constraints.cc:3094
DoFTools::compute_intergrid_constraints
void compute_intergrid_constraints(const DoFHandler< dim, spacedim > &coarse_grid, const unsigned int coarse_component, const DoFHandler< dim, spacedim > &fine_grid, const unsigned int fine_component, const InterGridMap< DoFHandler< dim, spacedim > > &coarse_to_fine_grid_map, AffineConstraints< double > &constraints)
Definition dof_tools_constraints.cc:4529
DoFTools::make_zero_boundary_constraints
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask={})
Definition dof_tools_constraints.cc:4810
fe_values.h
intergrid_map.h
la_parallel_vector.h
size
std::size_t size
Definition mpi.cc:733
Differentiation::SD::fabs
Expression fabs(const Expression &x)
Definition symengine_math.cc:272
DoFTools::internal::make_hp_hanging_node_constraints
void make_hp_hanging_node_constraints(const DoFHandler< 1, spacedim > &, AffineConstraints< number > &)
Definition dof_tools_constraints.cc:595
DoFTools::internal::set_periodicity_constraints
void set_periodicity_constraints(const FaceIterator &face_1, const std_cxx20::type_identity_t< FaceIterator > &face_2, const FullMatrix< double > &transformation, AffineConstraints< number > &affine_constraints, const ComponentMask &component_mask, const types::geometric_orientation combined_orientation, const number periodicity_factor, const unsigned int level=numbers::invalid_unsigned_int)
Definition dof_tools_constraints.cc:3124
DoFTools::internal::make_hanging_node_constraints_nedelec
void make_hanging_node_constraints_nedelec(const ::DoFHandler< 1, spacedim > &, AffineConstraints< number > &, std::integral_constant< int, 1 >)
Definition dof_tools_constraints.cc:606
DoFTools::internal::make_oldstyle_hanging_node_constraints
void make_oldstyle_hanging_node_constraints(const DoFHandler< 1, spacedim > &, AffineConstraints< number > &, std::integral_constant< int, 1 >)
Definition dof_tools_constraints.cc:618
DoFTools::extract_locally_relevant_dofs
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_tools.cc:1271
DoFTools::make_periodicity_constraints_on_level
void make_periodicity_constraints_on_level(const FaceIterator &face_1, const std_cxx20::type_identity_t< FaceIterator > &face_2, const unsigned int level, AffineConstraints< number > &constraints, const ComponentMask &component_mask={}, const types::geometric_orientation combined_orientation=numbers::default_geometric_orientation, const FullMatrix< double > &matrix=FullMatrix< double >(), const std::vector< unsigned int > &first_vector_components=std::vector< unsigned int >(), const number periodicity_factor=1.)
Definition dof_tools_constraints.cc:3822
DoFTools::make_periodicity_constraints
void make_periodicity_constraints(const FaceIterator &face_1, const std_cxx20::type_identity_t< FaceIterator > &face_2, AffineConstraints< number > &constraints, const ComponentMask &component_mask={}, const types::geometric_orientation combined_orientation=numbers::default_geometric_orientation, const FullMatrix< double > &matrix=FullMatrix< double >(), const std::vector< unsigned int > &first_vector_components=std::vector< unsigned int >(), const number periodicity_factor=1.)
Definition dof_tools_constraints.cc:3607
GridTools::collect_periodic_faces
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id1, const types::boundary_id b_id2, const unsigned int direction, std::vector< PeriodicFacePair< typename MeshType::cell_iterator > > &matched_pairs, const Tensor< 1, MeshType::space_dimension > &offset=::Tensor< 1, MeshType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >(), const double abs_tol=1e-10)
Definition grid_tools_dof_handlers.cc:2295
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:56
LAPACKSupport::N
constexpr char N
Definition lapack_support.h:155
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
PlasticityLab::if
*  *  if(update_pressure &update_flags) *  compute_pressure(constitutive_request
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm mpi_communicator)
WorkStream::run
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition work_stream.h:1709
internal::split_face_orientation
std::tuple< bool, bool, bool > split_face_orientation(const types::geometric_orientation combined_orientation)
Definition tria_orientation.h:44
numbers::invalid_dof_index
constexpr types::global_dof_index invalid_dof_index
Definition types.h:259
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
numbers::invalid_boundary_id
constexpr types::boundary_id invalid_boundary_id
Definition types.h:299
numbers::default_geometric_orientation
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:342
numbers::invalid_fe_index
constexpr types::fe_index invalid_fe_index
Definition types.h:250
std_cxx20::type_identity_t
typename type_identity< T >::type type_identity_t
Definition type_traits.h:93
std
STL namespace.
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6563
types::fe_index
unsigned short int fe_index
Definition types.h:70
types::geometric_orientation
std::uint8_t geometric_orientation
Definition types.h:38
DoFTools::internal::Assembler::CopyData::global_parameter_representation
std::vector<::LinearAlgebra::distributed::Vector< double > > global_parameter_representation
Definition dof_tools_constraints.cc:4017
DoFTools::internal::Assembler::CopyData::parameter_dof_indices
std::vector< types::global_dof_index > parameter_dof_indices
Definition dof_tools_constraints.cc:4014
GeometryInfo::n_children
static unsigned int n_children(const RefinementCase< dim > &refinement_case)
template_constraints.h
tria_iterator.h
work_stream.h
Generated by doxygen 1.9.8