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