Reference documentation for deal.II version 9.3.3
\(\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\}}\)
dof_tools_constraints.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#include <deal.II/base/table.h>
17#include <deal.II/base/template_constraints.h>
18#include <deal.II/base/utilities.h>
19#include <deal.II/base/work_stream.h>
20
21#include <deal.II/dofs/dof_accessor.h>
22#include <deal.II/dofs/dof_handler.h>
23#include <deal.II/dofs/dof_tools.h>
24
25#include <deal.II/fe/fe.h>
26#include <deal.II/fe/fe_tools.h>
27#include <deal.II/fe/fe_values.h>
28
29#include <deal.II/grid/grid_tools.h>
30#include <deal.II/grid/intergrid_map.h>
31#include <deal.II/grid/tria.h>
32#include <deal.II/grid/tria_iterator.h>
33
34#include <deal.II/hp/fe_collection.h>
35#include <deal.II/hp/fe_values.h>
36
37#include <deal.II/lac/affine_constraints.h>
38#include <deal.II/lac/vector.h>
39
40#ifdef DEAL_II_WITH_MPI
41# include <deal.II/lac/la_parallel_vector.h>
42#endif
43
44#include <algorithm>
45#include <array>
46#include <memory>
47#include <numeric>
48
49DEAL_II_NAMESPACE_OPEN
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(
173 const FiniteElement<dim, spacedim> &fe1,
174 const FiniteElement<dim, spacedim> &fe2,
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
180 AssertDimension(fe1.n_unique_faces(), 1);
181 AssertDimension(fe2.n_unique_faces(), 1);
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),
186 ExcInternalError());
187 AssertDimension(primary_dof_mask.size(), fe1.n_dofs_per_face(face_no));
188
189 Assert(fe2.n_dofs_per_vertex() <= fe1.n_dofs_per_vertex(),
190 ExcInternalError());
191 Assert(fe2.n_dofs_per_line() <= fe1.n_dofs_per_line(),
192 ExcInternalError());
193 Assert((dim < 3) ||
194 (fe2.n_dofs_per_quad(face_no) <= fe1.n_dofs_per_quad(face_no)),
195 ExcInternalError());
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
229 Assert(i < fe1.n_dofs_per_vertex(), ExcInternalError());
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 {
258 Assert(i < fe1.n_dofs_per_line(), ExcInternalError());
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 {
281 Assert(i < fe1.n_dofs_per_quad(q), ExcInternalError());
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(
313 const FiniteElement<dim, spacedim> &fe1,
314 const FiniteElement<dim, spacedim> &fe2,
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
320 AssertDimension(fe1.n_unique_faces(), 1);
321 AssertDimension(fe2.n_unique_faces(), 1);
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(
345 const FiniteElement<dim, spacedim> & fe1,
346 const FiniteElement<dim, spacedim> & fe2,
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
351 AssertDimension(fe1.n_unique_faces(), 1);
352 AssertDimension(fe2.n_unique_faces(), 1);
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(
371 const FiniteElement<dim, spacedim> & fe1,
372 const FiniteElement<dim, spacedim> & fe2,
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
378 AssertDimension(fe1.n_unique_faces(), 1);
379 AssertDimension(fe2.n_unique_faces(), 1);
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));
386 fe1.get_subface_interpolation_matrix(fe2,
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()),
412 ExcInternalError());
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,
506 ExcInternalError());
507 for (unsigned int col = 0; col != n_primary_dofs; ++col)
508 Assert(primary_dofs[col] != numbers::invalid_dof_index,
509 ExcInternalError());
510
511
512 for (unsigned int row = 0; row != n_dependent_dofs; ++row)
513 if (constraints.is_constrained(dependent_dofs[row]) == false)
514 {
515 // Check if we have an identity constraint, i.e.,
516 // something of the form
517 // U(dependent_dof[row])==U(primary_dof[row]),
518 // where
519 // dependent_dof[row] == primary_dof[row].
520 // This can happen in the hp context where we have previously
521 // unified DoF indices, for example, the middle node on the
522 // face of a Q4 element will have gotten the same index
523 // as the middle node of the Q2 element on the neighbor
524 // cell. But because the other Q4 nodes will still have to be
525 // constrained, so the middle node shows up again here.
526 //
527 // If we find such a constraint, then it is trivially
528 // satisfied, and we can move on to the next dependent
529 // DoF (row). The only thing we should make sure is that the
530 // row of the matrix really just contains this one entry.
531 {
532 bool is_trivial_constraint = false;
533
534 for (unsigned int i = 0; i < n_primary_dofs; ++i)
535 if (face_constraints(row, i) == 1.0)
536 if (dependent_dofs[row] == primary_dofs[i])
537 {
538 is_trivial_constraint = true;
539
540 for (unsigned int ii = 0; ii < n_primary_dofs; ++ii)
541 if (ii != i)
542 Assert(face_constraints(row, ii) == 0.0,
543 ExcInternalError());
544
545 break;
546 }
547
548 if (is_trivial_constraint == true)
549 continue;
550 }
551 // add up the absolute values of all constraints in this line
552 // to get a measure of their absolute size
553 number1 abs_sum = 0;
554 for (unsigned int i = 0; i < n_primary_dofs; ++i)
555 abs_sum += std::abs(face_constraints(row, i));
556
557 // then enter those constraints that are larger than
558 // 1e-14*abs_sum. everything else probably originated from
559 // inexact inversion of matrices and similar effects. having
560 // those constraints in here will only lead to problems
561 // because it makes sparsity patterns fuller than necessary
562 // without producing any significant effect
563 constraints.add_line(dependent_dofs[row]);
564 for (unsigned int i = 0; i < n_primary_dofs; ++i)
565 if ((face_constraints(row, i) != 0) &&
566 (std::fabs(face_constraints(row, i)) >= 1e-14 * abs_sum))
567 constraints.add_entry(dependent_dofs[row],
568 primary_dofs[i],
569 face_constraints(row, i));
570 constraints.set_inhomogeneity(dependent_dofs[row], 0.);
571 }
572 }
573
574 } // namespace
575
576
577 template <typename number>
578 void
579 make_hp_hanging_node_constraints(const ::DoFHandler<1> &,
580 AffineConstraints<number> &)
581 {
582 // nothing to do for regular dof handlers in 1d
583 }
584
585
586 template <typename number>
587 void
588 make_oldstyle_hanging_node_constraints(const ::DoFHandler<1> &,
589 AffineConstraints<number> &,
590 std::integral_constant<int, 1>)
591 {
592 // nothing to do for regular dof handlers in 1d
593 }
594
595
596 template <typename number>
597 void
598 make_hp_hanging_node_constraints(const ::DoFHandler<1, 2> &,
599 AffineConstraints<number> &)
600 {
601 // nothing to do for regular dof handlers in 1d
602 }
603
604
605 template <typename number>
606 void
607 make_oldstyle_hanging_node_constraints(const ::DoFHandler<1, 2> &,
608 AffineConstraints<number> &,
609 std::integral_constant<int, 1>)
610 {
611 // nothing to do for regular dof handlers in 1d
612 }
613
614
615 template <typename number, int spacedim>
616 void
617 make_hp_hanging_node_constraints(
618 const ::DoFHandler<1, spacedim> & /*dof_handler*/,
619 AffineConstraints<number> & /*constraints*/)
620 {
621 // nothing to do for dof handlers in 1d
622 }
623
624
625 template <typename number, int spacedim>
626 void
627 make_oldstyle_hanging_node_constraints(
628 const ::DoFHandler<1, spacedim> & /*dof_handler*/,
629 AffineConstraints<number> & /*constraints*/,
630 std::integral_constant<int, 1>)
631 {
632 // nothing to do for dof handlers in 1d
633 }
634
635 template <int dim_, int spacedim, typename number>
636 void
637 make_oldstyle_hanging_node_constraints(
638 const DoFHandler<dim_, spacedim> &dof_handler,
639 AffineConstraints<number> & constraints,
640 std::integral_constant<int, 2>)
641 {
642 const unsigned int dim = 2;
643
644 std::vector<types::global_dof_index> dofs_on_mother;
645 std::vector<types::global_dof_index> dofs_on_children;
646
647 // loop over all lines; only on lines there can be constraints. We do so
648 // by looping over all active cells and checking whether any of the faces
649 // are refined which can only be from the neighboring cell because this
650 // one is active. In that case, the face is subject to constraints
651 //
652 // note that even though we may visit a face twice if the neighboring
653 // cells are equally refined, we can only visit each face with hanging
654 // nodes once
655 typename DoFHandler<dim_, spacedim>::active_cell_iterator
656 cell = dof_handler.begin_active(),
657 endc = dof_handler.end();
658 for (; cell != endc; ++cell)
659 {
660 // artificial cells can at best neighbor ghost cells, but we're not
661 // interested in these interfaces
662 if (cell->is_artificial())
663 continue;
664
665 for (const unsigned int face : cell->face_indices())
666 if (cell->face(face)->has_children())
667 {
668 // in any case, faces can have at most two active FE indices,
669 // but here the face can have only one (namely the same as that
670 // from the cell we're sitting on), and each of the children can
671 // have only one as well. check this
672 Assert(cell->face(face)->n_active_fe_indices() == 1,
673 ExcInternalError());
674 Assert(cell->face(face)->fe_index_is_active(
675 cell->active_fe_index()) == true,
676 ExcInternalError());
677 for (unsigned int c = 0; c < cell->face(face)->n_children();
678 ++c)
679 if (!cell->neighbor_child_on_subface(face, c)
680 ->is_artificial())
681 Assert(cell->face(face)->child(c)->n_active_fe_indices() ==
682 1,
683 ExcInternalError());
684
685 // right now, all that is implemented is the case that both
686 // sides use the same FE
687 for (unsigned int c = 0; c < cell->face(face)->n_children();
688 ++c)
689 if (!cell->neighbor_child_on_subface(face, c)
690 ->is_artificial())
691 Assert(cell->face(face)->child(c)->fe_index_is_active(
692 cell->active_fe_index()) == true,
693 ExcNotImplemented());
694
695 // ok, start up the work
696 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
697 const unsigned int fe_index = cell->active_fe_index();
698
699 const unsigned int n_dofs_on_mother =
700 2 * fe.n_dofs_per_vertex() +
701 fe.n_dofs_per_line(),
702 n_dofs_on_children =
703 fe.n_dofs_per_vertex() +
704 2 * fe.n_dofs_per_line();
705
706 dofs_on_mother.resize(n_dofs_on_mother);
707 // we might not use all of those in case of artificial cells, so
708 // do not resize(), but reserve() and use push_back later.
709 dofs_on_children.clear();
710 dofs_on_children.reserve(n_dofs_on_children);
711
712 Assert(n_dofs_on_mother == fe.constraints().n(),
713 ExcDimensionMismatch(n_dofs_on_mother,
714 fe.constraints().n()));
715 Assert(n_dofs_on_children == fe.constraints().m(),
716 ExcDimensionMismatch(n_dofs_on_children,
717 fe.constraints().m()));
718
719 const typename DoFHandler<dim_, spacedim>::line_iterator
720 this_face = cell->face(face);
721
722 // fill the dofs indices. Use same enumeration scheme as in
723 // @p{FiniteElement::constraints()}
724 unsigned int next_index = 0;
725 for (unsigned int vertex = 0; vertex < 2; ++vertex)
726 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
727 ++dof)
728 dofs_on_mother[next_index++] =
729 this_face->vertex_dof_index(vertex, dof, fe_index);
730 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof)
731 dofs_on_mother[next_index++] =
732 this_face->dof_index(dof, fe_index);
733 AssertDimension(next_index, dofs_on_mother.size());
734
735 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); ++dof)
736 dofs_on_children.push_back(
737 this_face->child(0)->vertex_dof_index(1, dof, fe_index));
738 for (unsigned int child = 0; child < 2; ++child)
739 {
740 // skip artificial cells
741 if (cell->neighbor_child_on_subface(face, child)
742 ->is_artificial())
743 continue;
744 for (unsigned int dof = 0; dof != fe.n_dofs_per_line();
745 ++dof)
746 dofs_on_children.push_back(
747 this_face->child(child)->dof_index(dof, fe_index));
748 }
749 // note: can get fewer DoFs when we have artificial cells
750 Assert(dofs_on_children.size() <= n_dofs_on_children,
751 ExcInternalError());
752
753 // for each row in the AffineConstraints object for this line:
754 for (unsigned int row = 0; row != dofs_on_children.size();
755 ++row)
756 {
757 constraints.add_line(dofs_on_children[row]);
758 for (unsigned int i = 0; i != dofs_on_mother.size(); ++i)
759 constraints.add_entry(dofs_on_children[row],
760 dofs_on_mother[i],
761 fe.constraints()(row, i));
762
763 constraints.set_inhomogeneity(dofs_on_children[row], 0.);
764 }
765 }
766 else
767 {
768 // this face has no children, but it could still be that it is
769 // shared by two cells that use a different FE index. check a
770 // couple of things, but ignore the case that the neighbor is an
771 // artificial cell
772 if (!cell->at_boundary(face) &&
773 !cell->neighbor(face)->is_artificial())
774 {
775 Assert(cell->face(face)->n_active_fe_indices() == 1,
776 ExcNotImplemented());
777 Assert(cell->face(face)->fe_index_is_active(
778 cell->active_fe_index()) == true,
779 ExcInternalError());
780 }
781 }
782 }
783 }
784
785
786
787 template <int dim_, int spacedim, typename number>
788 void
789 make_oldstyle_hanging_node_constraints(
790 const DoFHandler<dim_, spacedim> &dof_handler,
791 AffineConstraints<number> & constraints,
792 std::integral_constant<int, 3>)
793 {
794 const unsigned int dim = 3;
795
796 std::vector<types::global_dof_index> dofs_on_mother;
797 std::vector<types::global_dof_index> dofs_on_children;
798
799 // loop over all quads; only on quads there can be constraints. We do so
800 // by looping over all active cells and checking whether any of the faces
801 // are refined which can only be from the neighboring cell because this
802 // one is active. In that case, the face is subject to constraints
803 //
804 // note that even though we may visit a face twice if the neighboring
805 // cells are equally refined, we can only visit each face with hanging
806 // nodes once
807 typename DoFHandler<dim_, spacedim>::active_cell_iterator
808 cell = dof_handler.begin_active(),
809 endc = dof_handler.end();
810 for (; cell != endc; ++cell)
811 {
812 // artificial cells can at best neighbor ghost cells, but we're not
813 // interested in these interfaces
814 if (cell->is_artificial())
815 continue;
816
817 for (const unsigned int face : cell->face_indices())
818 if (cell->face(face)->has_children())
819 {
820 // first of all, make sure that we treat a case which is
821 // possible, i.e. either no dofs on the face at all or no
822 // anisotropic refinement
823 if (cell->get_fe().n_dofs_per_face(face) == 0)
824 continue;
825
826 Assert(cell->face(face)->refinement_case() ==
827 RefinementCase<dim - 1>::isotropic_refinement,
828 ExcNotImplemented());
829
830 // in any case, faces can have at most two active FE indices,
831 // but here the face can have only one (namely the same as that
832 // from the cell we're sitting on), and each of the children can
833 // have only one as well. check this
834 AssertDimension(cell->face(face)->n_active_fe_indices(), 1);
835 Assert(cell->face(face)->fe_index_is_active(
836 cell->active_fe_index()) == true,
837 ExcInternalError());
838 for (unsigned int c = 0; c < cell->face(face)->n_children();
839 ++c)
840 if (!cell->neighbor_child_on_subface(face, c)
841 ->is_artificial())
842 AssertDimension(
843 cell->face(face)->child(c)->n_active_fe_indices(), 1);
844
845 // right now, all that is implemented is the case that both
846 // sides use the same fe, and not only that but also that all
847 // lines bounding this face and the children have the same FE
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())
852 {
853 Assert(cell->face(face)->child(c)->fe_index_is_active(
854 cell->active_fe_index()) == true,
855 ExcNotImplemented());
856 for (unsigned int e = 0; e < 4; ++e)
857 {
858 Assert(cell->face(face)
859 ->child(c)
860 ->line(e)
861 ->n_active_fe_indices() == 1,
862 ExcNotImplemented());
863 Assert(cell->face(face)
864 ->child(c)
865 ->line(e)
866 ->fe_index_is_active(
867 cell->active_fe_index()) == true,
868 ExcNotImplemented());
869 }
870 }
871 for (unsigned int e = 0; e < 4; ++e)
872 {
873 Assert(cell->face(face)->line(e)->n_active_fe_indices() ==
874 1,
875 ExcNotImplemented());
876 Assert(cell->face(face)->line(e)->fe_index_is_active(
877 cell->active_fe_index()) == true,
878 ExcNotImplemented());
879 }
880
881 // ok, start up the work
882 const FiniteElement<dim> &fe = cell->get_fe();
883 const unsigned int fe_index = cell->active_fe_index();
884
885 const unsigned int n_dofs_on_mother = fe.n_dofs_per_face(face);
886 const unsigned int n_dofs_on_children =
887 (5 * fe.n_dofs_per_vertex() + 12 * fe.n_dofs_per_line() +
888 4 * fe.n_dofs_per_quad(face));
889
890 // TODO[TL]: think about this and the following in case of
891 // anisotropic refinement
892
893 dofs_on_mother.resize(n_dofs_on_mother);
894 // we might not use all of those in case of artificial cells, so
895 // do not resize(), but reserve() and use push_back later.
896 dofs_on_children.clear();
897 dofs_on_children.reserve(n_dofs_on_children);
898
899 Assert(n_dofs_on_mother == fe.constraints().n(),
900 ExcDimensionMismatch(n_dofs_on_mother,
901 fe.constraints().n()));
902 Assert(n_dofs_on_children == fe.constraints().m(),
903 ExcDimensionMismatch(n_dofs_on_children,
904 fe.constraints().m()));
905
906 const typename DoFHandler<dim_, spacedim>::face_iterator
907 this_face = cell->face(face);
908
909 // fill the dofs indices. Use same enumeration scheme as in
910 // @p{FiniteElement::constraints()}
911 unsigned int next_index = 0;
912 for (unsigned int vertex = 0; vertex < 4; ++vertex)
913 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
914 ++dof)
915 dofs_on_mother[next_index++] =
916 this_face->vertex_dof_index(vertex, dof, fe_index);
917 for (unsigned int line = 0; line < 4; ++line)
918 for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof)
919 dofs_on_mother[next_index++] =
920 this_face->line(line)->dof_index(dof, fe_index);
921 for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face);
922 ++dof)
923 dofs_on_mother[next_index++] =
924 this_face->dof_index(dof, fe_index);
925 AssertDimension(next_index, dofs_on_mother.size());
926
927 // TODO: assert some consistency assumptions
928
929 // TODO[TL]: think about this in case of anisotropic refinement
930
931 Assert(dof_handler.get_triangulation()
932 .get_anisotropic_refinement_flag() ||
933 ((this_face->child(0)->vertex_index(3) ==
934 this_face->child(1)->vertex_index(2)) &&
935 (this_face->child(0)->vertex_index(3) ==
936 this_face->child(2)->vertex_index(1)) &&
937 (this_face->child(0)->vertex_index(3) ==
938 this_face->child(3)->vertex_index(0))),
939 ExcInternalError());
940
941 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex(); ++dof)
942 dofs_on_children.push_back(
943 this_face->child(0)->vertex_dof_index(3, dof));
944
945 // dof numbers on the centers of the lines bounding this face
946 for (unsigned int line = 0; line < 4; ++line)
947 for (unsigned int dof = 0; dof != fe.n_dofs_per_vertex();
948 ++dof)
949 dofs_on_children.push_back(
950 this_face->line(line)->child(0)->vertex_dof_index(
951 1, dof, fe_index));
952
953 // next the dofs on the lines interior to the face; the order of
954 // these lines is laid down in the FiniteElement class
955 // documentation
956 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
957 dofs_on_children.push_back(
958 this_face->child(0)->line(1)->dof_index(dof, fe_index));
959 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
960 dofs_on_children.push_back(
961 this_face->child(2)->line(1)->dof_index(dof, fe_index));
962 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
963 dofs_on_children.push_back(
964 this_face->child(0)->line(3)->dof_index(dof, fe_index));
965 for (unsigned int dof = 0; dof < fe.n_dofs_per_line(); ++dof)
966 dofs_on_children.push_back(
967 this_face->child(1)->line(3)->dof_index(dof, fe_index));
968
969 // dofs on the bordering lines
970 for (unsigned int line = 0; line < 4; ++line)
971 for (unsigned int child = 0; child < 2; ++child)
972 {
973 for (unsigned int dof = 0; dof != fe.n_dofs_per_line();
974 ++dof)
975 dofs_on_children.push_back(
976 this_face->line(line)->child(child)->dof_index(
977 dof, fe_index));
978 }
979
980 // finally, for the dofs interior to the four child faces
981 for (unsigned int child = 0; child < 4; ++child)
982 {
983 // skip artificial cells
984 if (cell->neighbor_child_on_subface(face, child)
985 ->is_artificial())
986 continue;
987 for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face);
988 ++dof)
989 dofs_on_children.push_back(
990 this_face->child(child)->dof_index(dof, fe_index));
991 }
992
993 // note: can get fewer DoFs when we have artificial cells:
994 Assert(dofs_on_children.size() <= n_dofs_on_children,
995 ExcInternalError());
996
997 // for each row in the AffineConstraints object for this line:
998 for (unsigned int row = 0; row != dofs_on_children.size();
999 ++row)
1000 {
1001 constraints.add_line(dofs_on_children[row]);
1002 for (unsigned int i = 0; i != dofs_on_mother.size(); ++i)
1003 constraints.add_entry(dofs_on_children[row],
1004 dofs_on_mother[i],
1005 fe.constraints()(row, i));
1006
1007 constraints.set_inhomogeneity(dofs_on_children[row], 0.);
1008 }
1009 }
1010 else
1011 {
1012 // this face has no children, but it could still be that it is
1013 // shared by two cells that use a different FE index. check a
1014 // couple of things, but ignore the case that the neighbor is an
1015 // artificial cell
1016 if (!cell->at_boundary(face) &&
1017 !cell->neighbor(face)->is_artificial())
1018 {
1019 Assert(cell->face(face)->n_active_fe_indices() == 1,
1020 ExcNotImplemented());
1021 Assert(cell->face(face)->fe_index_is_active(
1022 cell->active_fe_index()) == true,
1023 ExcInternalError());
1024 }
1025 }
1026 }
1027 }
1028
1029
1030
1031 template <int dim, int spacedim, typename number>
1032 void
1033 make_hp_hanging_node_constraints(
1034 const DoFHandler<dim, spacedim> &dof_handler,
1035 AffineConstraints<number> & constraints)
1036 {
1037 // note: this function is going to be hard to understand if you haven't
1038 // read the hp-paper. however, we try to follow the notation laid out
1039 // there, so go read the paper before you try to understand what is going
1040 // on here
1041
1042
1043 // a matrix to be used for constraints below. declared here and simply
1044 // resized down below to avoid permanent re-allocation of memory
1045 FullMatrix<double> constraint_matrix;
1046
1047 // similarly have arrays that will hold primary and dependent dof numbers,
1048 // as well as a scratch array needed for the complicated case below
1049 std::vector<types::global_dof_index> primary_dofs;
1050 std::vector<types::global_dof_index> dependent_dofs;
1051 std::vector<types::global_dof_index> scratch_dofs;
1052
1053 // caches for the face and subface interpolation matrices between
1054 // different (or the same) finite elements. we compute them only once,
1055 // namely the first time they are needed, and then just reuse them
1056 Table<2, std::unique_ptr<FullMatrix<double>>> face_interpolation_matrices(
1057 n_finite_elements(dof_handler), n_finite_elements(dof_handler));
1058 Table<3, std::unique_ptr<FullMatrix<double>>>
1059 subface_interpolation_matrices(
1060 n_finite_elements(dof_handler),
1061 n_finite_elements(dof_handler),
1062 GeometryInfo<dim>::max_children_per_face);
1063
1064 // similarly have a cache for the matrices that are split into their
1065 // primary and dependent parts, and for which the primary part is
1066 // inverted. these two matrices are derived from the face interpolation
1067 // matrix
1068 // as described in the @ref hp_paper "hp-paper"
1069 Table<2,
1070 std::unique_ptr<std::pair<FullMatrix<double>, FullMatrix<double>>>>
1071 split_face_interpolation_matrices(n_finite_elements(dof_handler),
1072 n_finite_elements(dof_handler));
1073
1074 // finally, for each pair of finite elements, have a mask that states
1075 // which of the degrees of freedom on the coarse side of a refined face
1076 // will act as primary dofs.
1077 Table<2, std::unique_ptr<std::vector<bool>>> primary_dof_masks(
1078 n_finite_elements(dof_handler), n_finite_elements(dof_handler));
1079
1080 // loop over all faces
1081 //
1082 // note that even though we may visit a face twice if the neighboring
1083 // cells are equally refined, we can only visit each face with hanging
1084 // nodes once
1085 for (const auto &cell : dof_handler.active_cell_iterators())
1086 {
1087 // artificial cells can at best neighbor ghost cells, but we're not
1088 // interested in these interfaces
1089 if (cell->is_artificial())
1090 continue;
1091
1092 for (const unsigned int face : cell->face_indices())
1093 if (cell->face(face)->has_children())
1094 {
1095 // first of all, make sure that we treat a case which is
1096 // possible, i.e. either no dofs on the face at all or no
1097 // anisotropic refinement
1098 if (cell->get_fe().n_dofs_per_face(face) == 0)
1099 continue;
1100
1101 Assert(cell->face(face)->refinement_case() ==
1102 RefinementCase<dim - 1>::isotropic_refinement,
1103 ExcNotImplemented());
1104
1105 // so now we've found a face of an active cell that has
1106 // children. that means that there are hanging nodes here.
1107
1108 // in any case, faces can have at most two sets of active FE
1109 // indices, but here the face can have only one (namely the same
1110 // as that from the cell we're sitting on), and each of the
1111 // children can have only one as well. check this
1112 Assert(cell->face(face)->n_active_fe_indices() == 1,
1113 ExcInternalError());
1114 Assert(cell->face(face)->fe_index_is_active(
1115 cell->active_fe_index()) == true,
1116 ExcInternalError());
1117 for (unsigned int c = 0; c < cell->face(face)->n_children();
1118 ++c)
1119 if (!cell->neighbor_child_on_subface(face, c)
1120 ->is_artificial())
1121 Assert(cell->face(face)->child(c)->n_active_fe_indices() ==
1122 1,
1123 ExcInternalError());
1124
1125 // first find out whether we can constrain each of the subfaces
1126 // to the mother face. in the lingo of the hp-paper, this would
1127 // be the simple case. note that we can short-circuit this
1128 // decision if the dof_handler doesn't support hp at all
1129 //
1130 // ignore all interfaces with artificial cells
1131 FiniteElementDomination::Domination mother_face_dominates =
1132 FiniteElementDomination::either_element_can_dominate;
1133
1134 // auxiliary variable which holds FE indices of the mother face
1135 // and its subfaces. This knowledge will be needed in hp-case
1136 // with neither_element_dominates.
1137 std::set<unsigned int> fe_ind_face_subface;
1138 fe_ind_face_subface.insert(cell->active_fe_index());
1139
1140 if (dof_handler.has_hp_capabilities())
1141 for (unsigned int c = 0;
1142 c < cell->face(face)->n_active_descendants();
1143 ++c)
1144 {
1145 const auto subcell =
1146 cell->neighbor_child_on_subface(face, c);
1147 if (!subcell->is_artificial())
1148 {
1149 mother_face_dominates =
1150 mother_face_dominates &
1151 (cell->get_fe().compare_for_domination(
1152 subcell->get_fe(), /*codim=*/1));
1153 fe_ind_face_subface.insert(
1154 subcell->active_fe_index());
1155 }
1156 }
1157
1158 switch (mother_face_dominates)
1159 {
1160 case FiniteElementDomination::this_element_dominates:
1161 case FiniteElementDomination::either_element_can_dominate:
1162 {
1163 // Case 1 (the simple case and the only case that can
1164 // happen for non-hp-DoFHandlers): The coarse element
1165 // dominates the elements on the subfaces (or they are
1166 // all the same)
1167 //
1168 // so we are going to constrain the DoFs on the face
1169 // children against the DoFs on the face itself
1170 primary_dofs.resize(
1171 cell->get_fe().n_dofs_per_face(face));
1172
1173 cell->face(face)->get_dof_indices(
1174 primary_dofs, cell->active_fe_index());
1175
1176 // Now create constraints for the subfaces and
1177 // assemble it. ignore all interfaces with artificial
1178 // cells because we can only get to such interfaces if
1179 // the current cell is a ghost cell
1180 for (unsigned int c = 0;
1181 c < cell->face(face)->n_children();
1182 ++c)
1183 {
1184 if (cell->neighbor_child_on_subface(face, c)
1185 ->is_artificial())
1186 continue;
1187
1188 const typename DoFHandler<dim, spacedim>::
1189 active_face_iterator subface =
1190 cell->face(face)->child(c);
1191
1192 Assert(subface->n_active_fe_indices() == 1,
1193 ExcInternalError());
1194
1195 const unsigned int subface_fe_index =
1196 subface->nth_active_fe_index(0);
1197
1198 // we sometime run into the situation where for
1199 // example on one big cell we have a FE_Q(1) and on
1200 // the subfaces we have a mixture of FE_Q(1) and
1201 // FE_Nothing. In that case, the face domination is
1202 // either_element_can_dominate for the whole
1203 // collection of subfaces, but on the particular
1204 // subface between FE_Q(1) and FE_Nothing, there are
1205 // no constraints that we need to take care of. in
1206 // that case, just continue
1207 if (cell->get_fe().compare_for_domination(
1208 subface->get_fe(subface_fe_index),
1209 /*codim=*/1) ==
1210 FiniteElementDomination::no_requirements)
1211 continue;
1212
1213 // Same procedure as for the mother cell. Extract
1214 // the face DoFs from the cell DoFs.
1215 dependent_dofs.resize(
1216 subface->get_fe(subface_fe_index)
1217 .n_dofs_per_face(face, c));
1218 subface->get_dof_indices(dependent_dofs,
1219 subface_fe_index);
1220
1221 for (const types::global_dof_index dependent_dof :
1222 dependent_dofs)
1223 {
1224 (void)dependent_dof;
1225 Assert(dependent_dof !=
1226 numbers::invalid_dof_index,
1227 ExcInternalError());
1228 }
1229
1230 // Now create the element constraint for this
1231 // subface.
1232 //
1233 // As a side remark, one may wonder the following:
1234 // neighbor_child is clearly computed correctly,
1235 // i.e. taking into account face_orientation (just
1236 // look at the implementation of that function).
1237 // however, we don't care about this here, when we
1238 // ask for subface_interpolation on subface c. the
1239 // question rather is: do we have to translate 'c'
1240 // here as well?
1241 //
1242 // the answer is in fact 'no'. if one does that,
1243 // results are wrong: constraints are added twice
1244 // for the same pair of nodes but with differing
1245 // weights. in addition, one can look at the
1246 // deal.II/project_*_03 tests that look at exactly
1247 // this case: there, we have a mesh with at least
1248 // one face_orientation==false and hanging nodes,
1249 // and the results of those tests show that the
1250 // result of projection verifies the approximation
1251 // properties of a finite element onto that mesh
1252 ensure_existence_of_subface_matrix(
1253 cell->get_fe(),
1254 subface->get_fe(subface_fe_index),
1255 c,
1256 subface_interpolation_matrices
1257 [cell->active_fe_index()][subface_fe_index][c]);
1258
1259 // Add constraints to global AffineConstraints
1260 // object.
1261 filter_constraints(primary_dofs,
1262 dependent_dofs,
1263 *(subface_interpolation_matrices
1264 [cell->active_fe_index()]
1265 [subface_fe_index][c]),
1266 constraints);
1267 } // loop over subfaces
1268
1269 break;
1270 } // Case 1
1271
1272 case FiniteElementDomination::other_element_dominates:
1273 case FiniteElementDomination::neither_element_dominates:
1274 {
1275 // Case 2 (the "complex" case): at least one (the
1276 // neither_... case) of the finer elements or all of
1277 // them (the other_... case) is dominating. See the hp-
1278 // paper for a way how to deal with this situation
1279 //
1280 // since this is something that can only happen for hp-
1281 // dof handlers, add a check here...
1282 Assert(dof_handler.has_hp_capabilities() == true,
1283 ExcInternalError());
1284
1285 const ::hp::FECollection<dim, spacedim>
1286 &fe_collection = dof_handler.get_fe_collection();
1287 // we first have to find the finite element that is able
1288 // to generate a space that all the other ones can be
1289 // constrained to. At this point we potentially have
1290 // different scenarios:
1291 //
1292 // 1) sub-faces dominate mother face and there is a
1293 // dominating FE among sub faces. We could loop over sub
1294 // faces to find the needed FE index. However, this will
1295 // not work in the case when ...
1296 //
1297 // 2) there is no dominating FE among sub faces (e.g.
1298 // Q1xQ2 vs Q2xQ1), but subfaces still dominate mother
1299 // face (e.g. Q2xQ2). To cover this case we would have
1300 // to find the least dominating element amongst all
1301 // finite elements on sub faces.
1302 //
1303 // 3) Finally, it could happen that we got here because
1304 // neither_element_dominates (e.g. Q1xQ1xQ2 and Q1xQ2xQ1
1305 // for subfaces and Q2xQ1xQ1 for mother face). This
1306 // requires finding the least dominating element amongst
1307 // all finite elements on sub faces and the mother face.
1308 //
1309 // Note that the last solution covers the first two
1310 // scenarios, thus we stick with it assuming that we
1311 // won't lose much time/efficiency.
1312 const unsigned int dominating_fe_index =
1313 fe_collection.find_dominating_fe_extended(
1314 fe_ind_face_subface, /*codim=*/1);
1315
1316 AssertThrow(
1317 dominating_fe_index != numbers::invalid_unsigned_int,
1318 ExcMessage(
1319 "Could not find a least face dominating FE."));
1320
1321 const FiniteElement<dim, spacedim> &dominating_fe =
1322 dof_handler.get_fe(dominating_fe_index);
1323
1324 // first get the interpolation matrix from the mother to
1325 // the virtual dofs
1326 Assert(dominating_fe.n_dofs_per_face(face) <=
1327 cell->get_fe().n_dofs_per_face(face),
1328 ExcInternalError());
1329
1330 ensure_existence_of_face_matrix(
1331 dominating_fe,
1332 cell->get_fe(),
1333 face_interpolation_matrices[dominating_fe_index]
1334 [cell->active_fe_index()]);
1335
1336 // split this matrix into primary and dependent
1337 // components. invert the primary component
1338 ensure_existence_of_primary_dof_mask(
1339 cell->get_fe(),
1340 dominating_fe,
1341 (*face_interpolation_matrices
1342 [dominating_fe_index][cell->active_fe_index()]),
1343 primary_dof_masks[dominating_fe_index]
1344 [cell->active_fe_index()]);
1345
1346 ensure_existence_of_split_face_matrix(
1347 *face_interpolation_matrices[dominating_fe_index]
1348 [cell->active_fe_index()],
1349 (*primary_dof_masks[dominating_fe_index]
1350 [cell->active_fe_index()]),
1351 split_face_interpolation_matrices
1352 [dominating_fe_index][cell->active_fe_index()]);
1353
1354 const FullMatrix<double>
1355 &restrict_mother_to_virtual_primary_inv =
1356 (split_face_interpolation_matrices
1357 [dominating_fe_index][cell->active_fe_index()]
1358 ->first);
1359
1360 const FullMatrix<double>
1361 &restrict_mother_to_virtual_dependent =
1362 (split_face_interpolation_matrices
1363 [dominating_fe_index][cell->active_fe_index()]
1364 ->second);
1365
1366 // now compute the constraint matrix as the product
1367 // between the inverse matrix and the dependent part
1368 constraint_matrix.reinit(
1369 cell->get_fe().n_dofs_per_face(face) -
1370 dominating_fe.n_dofs_per_face(face),
1371 dominating_fe.n_dofs_per_face(face));
1372 restrict_mother_to_virtual_dependent.mmult(
1373 constraint_matrix,
1374 restrict_mother_to_virtual_primary_inv);
1375
1376 // then figure out the global numbers of primary and
1377 // dependent dofs and apply constraints
1378 scratch_dofs.resize(
1379 cell->get_fe().n_dofs_per_face(face));
1380 cell->face(face)->get_dof_indices(
1381 scratch_dofs, cell->active_fe_index());
1382
1383 // split dofs into primary and dependent components
1384 primary_dofs.clear();
1385 dependent_dofs.clear();
1386 for (unsigned int i = 0;
1387 i < cell->get_fe().n_dofs_per_face(face);
1388 ++i)
1389 if ((*primary_dof_masks[dominating_fe_index]
1390 [cell
1391 ->active_fe_index()])[i] ==
1392 true)
1393 primary_dofs.push_back(scratch_dofs[i]);
1394 else
1395 dependent_dofs.push_back(scratch_dofs[i]);
1396
1397 AssertDimension(primary_dofs.size(),
1398 dominating_fe.n_dofs_per_face(face));
1399 AssertDimension(dependent_dofs.size(),
1400 cell->get_fe().n_dofs_per_face(face) -
1401 dominating_fe.n_dofs_per_face(face));
1402
1403 filter_constraints(primary_dofs,
1404 dependent_dofs,
1405 constraint_matrix,
1406 constraints);
1407
1408
1409
1410 // next we have to deal with the subfaces. do as
1411 // discussed in the hp-paper
1412 for (unsigned int sf = 0;
1413 sf < cell->face(face)->n_children();
1414 ++sf)
1415 {
1416 // ignore interfaces with artificial cells as well
1417 // as interfaces between ghost cells in 2d
1418 if (cell->neighbor_child_on_subface(face, sf)
1419 ->is_artificial() ||
1420 (dim == 2 && cell->is_ghost() &&
1421 cell->neighbor_child_on_subface(face, sf)
1422 ->is_ghost()))
1423 continue;
1424
1425 Assert(cell->face(face)
1426 ->child(sf)
1427 ->n_active_fe_indices() == 1,
1428 ExcInternalError());
1429
1430 const unsigned int subface_fe_index =
1431 cell->face(face)->child(sf)->nth_active_fe_index(
1432 0);
1433 const FiniteElement<dim, spacedim> &subface_fe =
1434 dof_handler.get_fe(subface_fe_index);
1435
1436 // first get the interpolation matrix from the
1437 // subface to the virtual dofs
1438 Assert(dominating_fe.n_dofs_per_face(face) <=
1439 subface_fe.n_dofs_per_face(face),
1440 ExcInternalError());
1441 ensure_existence_of_subface_matrix(
1442 dominating_fe,
1443 subface_fe,
1444 sf,
1445 subface_interpolation_matrices
1446 [dominating_fe_index][subface_fe_index][sf]);
1447
1448 const FullMatrix<double>
1449 &restrict_subface_to_virtual = *(
1450 subface_interpolation_matrices
1451 [dominating_fe_index][subface_fe_index][sf]);
1452
1453 constraint_matrix.reinit(
1454 subface_fe.n_dofs_per_face(face),
1455 dominating_fe.n_dofs_per_face(face));
1456
1457 restrict_subface_to_virtual.mmult(
1458 constraint_matrix,
1459 restrict_mother_to_virtual_primary_inv);
1460
1461 dependent_dofs.resize(
1462 subface_fe.n_dofs_per_face(face));
1463 cell->face(face)->child(sf)->get_dof_indices(
1464 dependent_dofs, subface_fe_index);
1465
1466 filter_constraints(primary_dofs,
1467 dependent_dofs,
1468 constraint_matrix,
1469 constraints);
1470 } // loop over subfaces
1471
1472 break;
1473 } // Case 2
1474
1475 case FiniteElementDomination::no_requirements:
1476 // there are no continuity requirements between the two
1477 // elements. record no constraints
1478 break;
1479
1480 default:
1481 // we shouldn't get here
1482 Assert(false, ExcInternalError());
1483 }
1484 }
1485 else
1486 {
1487 // this face has no children, but it could still be that it is
1488 // shared by two cells that use a different FE index
1489 Assert(cell->face(face)->fe_index_is_active(
1490 cell->active_fe_index()) == true,
1491 ExcInternalError());
1492
1493 // see if there is a neighbor that is an artificial cell. in
1494 // that case, we're not interested in this interface. we test
1495 // this case first since artificial cells may not have an
1496 // active FE index set, etc
1497 if (!cell->at_boundary(face) &&
1498 cell->neighbor(face)->is_artificial())
1499 continue;
1500
1501 // Only if there is a neighbor with a different active FE index
1502 // and the same h-level, some action has to be taken.
1503 if ((dof_handler.has_hp_capabilities()) &&
1504 !cell->face(face)->at_boundary() &&
1505 (cell->neighbor(face)->active_fe_index() !=
1506 cell->active_fe_index()) &&
1507 (!cell->face(face)->has_children() &&
1508 !cell->neighbor_is_coarser(face)))
1509 {
1510 const typename DoFHandler<dim,
1511 spacedim>::level_cell_iterator
1512 neighbor = cell->neighbor(face);
1513
1514 // see which side of the face we have to constrain
1515 switch (
1516 cell->get_fe().compare_for_domination(neighbor->get_fe(),
1517 /*codim=*/1))
1518 {
1519 case FiniteElementDomination::this_element_dominates:
1520 {
1521 // Get DoFs on dominating and dominated side of the
1522 // face
1523 primary_dofs.resize(
1524 cell->get_fe().n_dofs_per_face(face));
1525 cell->face(face)->get_dof_indices(
1526 primary_dofs, cell->active_fe_index());
1527
1528 // break if the n_primary_dofs == 0, because we are
1529 // attempting to constrain to an element that has no
1530 // face dofs
1531 if (primary_dofs.size() == 0)
1532 break;
1533
1534 dependent_dofs.resize(
1535 neighbor->get_fe().n_dofs_per_face(face));
1536 cell->face(face)->get_dof_indices(
1537 dependent_dofs, neighbor->active_fe_index());
1538
1539 // make sure the element constraints for this face
1540 // are available
1541 ensure_existence_of_face_matrix(
1542 cell->get_fe(),
1543 neighbor->get_fe(),
1544 face_interpolation_matrices
1545 [cell->active_fe_index()]
1546 [neighbor->active_fe_index()]);
1547
1548 // Add constraints to global constraint matrix.
1549 filter_constraints(
1550 primary_dofs,
1551 dependent_dofs,
1552 *(face_interpolation_matrices
1553 [cell->active_fe_index()]
1554 [neighbor->active_fe_index()]),
1555 constraints);
1556
1557 break;
1558 }
1559
1560 case FiniteElementDomination::other_element_dominates:
1561 {
1562 // we don't do anything here since we will come back
1563 // to this face from the other cell, at which time
1564 // we will fall into the first case clause above
1565 break;
1566 }
1567
1568 case FiniteElementDomination::
1569 either_element_can_dominate:
1570 {
1571 // it appears as if neither element has any
1572 // constraints on its neighbor. this may be because
1573 // neither element has any DoFs on faces at all. or
1574 // that the two elements are actually the same,
1575 // although they happen to run under different
1576 // fe_indices (this is what happens in
1577 // hp/hp_hanging_nodes_01 for example).
1578 //
1579 // another possibility is what happens in crash_13.
1580 // there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs.
1581 // FESystem(FE_Q(1),FE_DGQ(1)). neither of them
1582 // dominates the other.
1583 //
1584 // a final possibility is that we have something
1585 // like FESystem(FE_Q(1),FE_Q(1)) vs
1586 // FESystem(FE_Q(1),FE_Nothing()), see
1587 // hp/fe_nothing_18/19.
1588 //
1589 // in any case, the point is that it doesn't matter.
1590 // there is nothing to do here.
1591 break;
1592 }
1593
1594 case FiniteElementDomination::neither_element_dominates:
1595 {
1596 // make sure we don't get here twice from each cell
1597 if (cell < neighbor)
1598 break;
1599
1600 // our best bet is to find the common space among
1601 // other FEs in FECollection and then constrain both
1602 // FEs to that one. More precisely, we follow the
1603 // strategy outlined on page 17 of the hp-paper:
1604 // First we find the dominant FE space S. Then we
1605 // divide our dofs in primary and dependent such
1606 // that I^{face,primary}_{S^{face}->S} is
1607 // invertible. And finally constrain dependent dofs
1608 // to primary dofs based on the interpolation
1609 // matrix.
1610
1611 const unsigned int this_fe_index =
1612 cell->active_fe_index();
1613 const unsigned int neighbor_fe_index =
1614 neighbor->active_fe_index();
1615 std::set<unsigned int> fes;
1616 fes.insert(this_fe_index);
1617 fes.insert(neighbor_fe_index);
1618 const ::hp::FECollection<dim, spacedim>
1619 &fe_collection = dof_handler.get_fe_collection();
1620
1621 const unsigned int dominating_fe_index =
1622 fe_collection.find_dominating_fe_extended(
1623 fes, /*codim=*/1);
1624
1625 AssertThrow(
1626 dominating_fe_index !=
1627 numbers::invalid_unsigned_int,
1628 ExcMessage(
1629 "Could not find the dominating FE for " +
1630 cell->get_fe().get_name() + " and " +
1631 neighbor->get_fe().get_name() +
1632 " inside FECollection."));
1633
1634 const FiniteElement<dim, spacedim> &dominating_fe =
1635 fe_collection[dominating_fe_index];
1636
1637 // TODO: until we hit the second face, the code is a
1638 // copy-paste from h-refinement case...
1639
1640 // first get the interpolation matrix from main FE
1641 // to the virtual dofs
1642 Assert(dominating_fe.n_dofs_per_face(face) <=
1643 cell->get_fe().n_dofs_per_face(face),
1644 ExcInternalError());
1645
1646 ensure_existence_of_face_matrix(
1647 dominating_fe,
1648 cell->get_fe(),
1649 face_interpolation_matrices
1650 [dominating_fe_index][cell->active_fe_index()]);
1651
1652 // split this matrix into primary and dependent
1653 // components. invert the primary component
1654 ensure_existence_of_primary_dof_mask(
1655 cell->get_fe(),
1656 dominating_fe,
1657 (*face_interpolation_matrices
1658 [dominating_fe_index]
1659 [cell->active_fe_index()]),
1660 primary_dof_masks[dominating_fe_index]
1661 [cell->active_fe_index()]);
1662
1663 ensure_existence_of_split_face_matrix(
1664 *face_interpolation_matrices
1665 [dominating_fe_index][cell->active_fe_index()],
1666 (*primary_dof_masks[dominating_fe_index]
1667 [cell->active_fe_index()]),
1668 split_face_interpolation_matrices
1669 [dominating_fe_index][cell->active_fe_index()]);
1670
1671 const FullMatrix<
1672 double> &restrict_mother_to_virtual_primary_inv =
1673 (split_face_interpolation_matrices
1674 [dominating_fe_index][cell->active_fe_index()]
1675 ->first);
1676
1677 const FullMatrix<
1678 double> &restrict_mother_to_virtual_dependent =
1679 (split_face_interpolation_matrices
1680 [dominating_fe_index][cell->active_fe_index()]
1681 ->second);
1682
1683 // now compute the constraint matrix as the product
1684 // between the inverse matrix and the dependent part
1685 constraint_matrix.reinit(
1686 cell->get_fe().n_dofs_per_face(face) -
1687 dominating_fe.n_dofs_per_face(face),
1688 dominating_fe.n_dofs_per_face(face));
1689 restrict_mother_to_virtual_dependent.mmult(
1690 constraint_matrix,
1691 restrict_mother_to_virtual_primary_inv);
1692
1693 // then figure out the global numbers of primary and
1694 // dependent dofs and apply constraints
1695 scratch_dofs.resize(
1696 cell->get_fe().n_dofs_per_face(face));
1697 cell->face(face)->get_dof_indices(
1698 scratch_dofs, cell->active_fe_index());
1699
1700 // split dofs into primary and dependent components
1701 primary_dofs.clear();
1702 dependent_dofs.clear();
1703 for (unsigned int i = 0;
1704 i < cell->get_fe().n_dofs_per_face(face);
1705 ++i)
1706 if ((*primary_dof_masks[dominating_fe_index]
1707 [cell->active_fe_index()])
1708 [i] == true)
1709 primary_dofs.push_back(scratch_dofs[i]);
1710 else
1711 dependent_dofs.push_back(scratch_dofs[i]);
1712
1713 AssertDimension(primary_dofs.size(),
1714 dominating_fe.n_dofs_per_face(
1715 face));
1716 AssertDimension(
1717 dependent_dofs.size(),
1718 cell->get_fe().n_dofs_per_face(face) -
1719 dominating_fe.n_dofs_per_face(face));
1720
1721 filter_constraints(primary_dofs,
1722 dependent_dofs,
1723 constraint_matrix,
1724 constraints);
1725
1726 // now do the same for another FE this is pretty
1727 // much the same we do above to resolve h-refinement
1728 // constraints
1729 Assert(dominating_fe.n_dofs_per_face(face) <=
1730 neighbor->get_fe().n_dofs_per_face(face),
1731 ExcInternalError());
1732
1733 ensure_existence_of_face_matrix(
1734 dominating_fe,
1735 neighbor->get_fe(),
1736 face_interpolation_matrices
1737 [dominating_fe_index]
1738 [neighbor->active_fe_index()]);
1739
1740 const FullMatrix<double>
1741 &restrict_secondface_to_virtual =
1742 *(face_interpolation_matrices
1743 [dominating_fe_index]
1744 [neighbor->active_fe_index()]);
1745
1746 constraint_matrix.reinit(
1747 neighbor->get_fe().n_dofs_per_face(face),
1748 dominating_fe.n_dofs_per_face(face));
1749
1750 restrict_secondface_to_virtual.mmult(
1751 constraint_matrix,
1752 restrict_mother_to_virtual_primary_inv);
1753
1754 dependent_dofs.resize(
1755 neighbor->get_fe().n_dofs_per_face(face));
1756 cell->face(face)->get_dof_indices(
1757 dependent_dofs, neighbor->active_fe_index());
1758
1759 filter_constraints(primary_dofs,
1760 dependent_dofs,
1761 constraint_matrix,
1762 constraints);
1763
1764 break;
1765 }
1766
1767 case FiniteElementDomination::no_requirements:
1768 {
1769 // nothing to do here
1770 break;
1771 }
1772
1773 default:
1774 // we shouldn't get here
1775 Assert(false, ExcInternalError());
1776 }
1777 }
1778 }
1779 }
1780 }
1781 } // namespace internal
1782
1783
1784
1785 template <int dim, int spacedim, typename number>
1786 void
1787 make_hanging_node_constraints(const DoFHandler<dim, spacedim> &dof_handler,
1788 AffineConstraints<number> & constraints)
1789 {
1790 Assert(dof_handler.has_active_dofs(),
1791 ExcMessage(
1792 "The given DoFHandler does not have any DoFs. Did you forget to "
1793 "call dof_handler.distribute_dofs()?"));
1794
1795 // Decide whether to use the new or old make_hanging_node_constraints
1796 // function. If all the FiniteElement or all elements in a FECollection
1797 // support the new face constraint matrix, the new code will be used.
1798 // Otherwise, the old implementation is used for the moment.
1799 if (dof_handler.get_fe_collection().hp_constraints_are_implemented())
1800 internal::make_hp_hanging_node_constraints(dof_handler, constraints);
1801 else
1802 internal::make_oldstyle_hanging_node_constraints(
1803 dof_handler, constraints, std::integral_constant<int, dim>());
1804 }
1805
1806
1807
1808 namespace
1809 {
1835 template <typename FaceIterator, typename number>
1836 void
1837 set_periodicity_constraints(
1838 const FaceIterator & face_1,
1839 const typename identity<FaceIterator>::type &face_2,
1840 const FullMatrix<double> & transformation,
1841 AffineConstraints<number> & affine_constraints,
1842 const ComponentMask & component_mask,
1843 const bool face_orientation,
1844 const bool face_flip,
1845 const bool face_rotation,
1846 const number periodicity_factor)
1847 {
1848 static const int dim = FaceIterator::AccessorType::dimension;
1849 static const int spacedim = FaceIterator::AccessorType::space_dimension;
1850
1851 // we should be in the case where face_1 is active, i.e. has no children:
1852 Assert(!face_1->has_children(), ExcInternalError());
1853
1854 Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
1855
1856 // TODO: the implementation makes the assumption that all faces have the
1857 // same number of dofs
1858 AssertDimension(
1859 face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
1860 AssertDimension(
1861 face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
1862 const unsigned int face_no = 0;
1863
1864 // If face_2 does have children, then we need to iterate over these
1865 // children and set periodic constraints in the inverse direction:
1866
1867 if (face_2->has_children())
1868 {
1869 Assert(face_2->n_children() ==
1870 GeometryInfo<dim>::max_children_per_face,
1871 ExcNotImplemented());
1872
1873 const unsigned int dofs_per_face =
1874 face_1->get_fe(face_1->nth_active_fe_index(0))
1875 .n_dofs_per_face(face_no);
1876 FullMatrix<double> child_transformation(dofs_per_face, dofs_per_face);
1877 FullMatrix<double> subface_interpolation(dofs_per_face,
1878 dofs_per_face);
1879
1880 for (unsigned int c = 0; c < face_2->n_children(); ++c)
1881 {
1882 // get the interpolation matrix recursively from the one that
1883 // interpolated from face_1 to face_2 by multiplying from the left
1884 // with the one that interpolates from face_2 to its child
1885 const auto &fe = face_1->get_fe(face_1->nth_active_fe_index(0));
1886 fe.get_subface_interpolation_matrix(fe,
1887 c,
1888 subface_interpolation,
1889 face_no);
1890 subface_interpolation.mmult(child_transformation, transformation);
1891
1892 set_periodicity_constraints(face_1,
1893 face_2->child(c),
1894 child_transformation,
1895 affine_constraints,
1896 component_mask,
1897 face_orientation,
1898 face_flip,
1899 face_rotation,
1900 periodicity_factor);
1901 }
1902 return;
1903 }
1904
1905 //
1906 // If we reached this point then both faces are active. Now all
1907 // that is left is to match the corresponding DoFs of both faces.
1908 //
1909
1910 const unsigned int face_1_index = face_1->nth_active_fe_index(0);
1911 const unsigned int face_2_index = face_2->nth_active_fe_index(0);
1912 Assert(face_1->get_fe(face_1_index) == face_2->get_fe(face_2_index),
1913 ExcMessage(
1914 "Matching periodic cells need to use the same finite element"));
1915
1916 const FiniteElement<dim, spacedim> &fe = face_1->get_fe(face_1_index);
1917
1918 Assert(component_mask.represents_n_components(fe.n_components()),
1919 ExcMessage(
1920 "The number of components in the mask has to be either "
1921 "zero or equal to the number of components in the finite "
1922 "element."));
1923
1924 const unsigned int dofs_per_face = fe.n_dofs_per_face(face_no);
1925
1926 std::vector<types::global_dof_index> dofs_1(dofs_per_face);
1927 std::vector<types::global_dof_index> dofs_2(dofs_per_face);
1928
1929 face_1->get_dof_indices(dofs_1, face_1_index);
1930 face_2->get_dof_indices(dofs_2, face_2_index);
1931
1932 // If either of the two faces has an invalid dof index, stop. This is
1933 // so that there is no attempt to match artificial cells of parallel
1934 // distributed triangulations.
1935 //
1936 // While it seems like we ought to be able to avoid even calling
1937 // set_periodicity_constraints for artificial faces, this situation
1938 // can arise when a face that is being made periodic is only
1939 // partially touched by the local subdomain.
1940 // make_periodicity_constraints will be called recursively even for
1941 // the section of the face that is not touched by the local
1942 // subdomain.
1943 //
1944 // Until there is a better way to determine if the cells that
1945 // neighbor a face are artificial, we simply test to see if the face
1946 // does not have a valid dof initialization.
1947
1948 for (unsigned int i = 0; i < dofs_per_face; i++)
1949 if (dofs_1[i] == numbers::invalid_dof_index ||
1950 dofs_2[i] == numbers::invalid_dof_index)
1951 {
1952 return;
1953 }
1954
1955 // Well, this is a hack:
1956 //
1957 // There is no
1958 // face_to_face_index(face_index,
1959 // face_orientation,
1960 // face_flip,
1961 // face_rotation)
1962 // function in FiniteElementData, so we have to use
1963 // face_to_cell_index(face_index, face
1964 // face_orientation,
1965 // face_flip,
1966 // face_rotation)
1967 // But this will give us an index on a cell - something we cannot work
1968 // with directly. But luckily we can match them back :-]
1969
1970 std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
1971
1972 // Build up a cell to face index for face_2:
1973 for (unsigned int i = 0; i < dofs_per_face; ++i)
1974 {
1975 const unsigned int cell_index =
1976 fe.face_to_cell_index(i,
1977 0, /* It doesn't really matter, just
1978 * assume we're on the first face...
1979 */
1980 true,
1981 false,
1982 false // default orientation
1983 );
1984 cell_to_rotated_face_index[cell_index] = i;
1985 }
1986
1987 //
1988 // Loop over all dofs on face 2 and constrain them against all
1989 // matching dofs on face 1:
1990 //
1991
1992 for (unsigned int i = 0; i < dofs_per_face; ++i)
1993 {
1994 // Obey the component mask
1995 if ((component_mask.n_selected_components(fe.n_components()) !=
1996 fe.n_components()) &&
1997 !component_mask[fe.face_system_to_component_index(i, face_no)
1998 .first])
1999 continue;
2000
2001 // We have to be careful to treat so called "identity
2002 // constraints" special. These are constraints of the form
2003 // x1 == constraint_factor * x_2. In this case, if the constraint
2004 // x2 == 1./constraint_factor * x1 already exists we are in trouble.
2005 //
2006 // Consequently, we have to check that we have indeed such an
2007 // "identity constraint". We do this by looping over all entries
2008 // of the row of the transformation matrix and check whether we
2009 // find exactly one nonzero entry. If this is the case, set
2010 // "is_identity_constrained" to true and record the corresponding
2011 // index and constraint_factor.
2012
2013 bool is_identity_constrained = false;
2014 unsigned int target = numbers::invalid_unsigned_int;
2015 number constraint_factor = periodicity_factor;
2016
2017 constexpr double eps = 1.e-13;
2018 for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
2019 {
2020 const auto entry = transformation(i, jj);
2021 if (std::abs(entry) > eps)
2022 {
2023 if (is_identity_constrained)
2024 {
2025 // We did encounter more than one nonzero entry, so
2026 // the dof is not identity constrained. Set the
2027 // boolean to false and break out of the for loop.
2028 is_identity_constrained = false;
2029 break;
2030 }
2031 is_identity_constrained = true;
2032 target = jj;
2033 constraint_factor = entry * periodicity_factor;
2034 }
2035 }
2036
2037 // Next, we work on all constraints that are not identity
2038 // constraints, i.e., constraints that involve an interpolation
2039 // step that constrains the current dof (on face 2) to more than
2040 // one dof on face 1.
2041
2042 if (!is_identity_constrained)
2043 {
2044 // The current dof is already constrained. There is nothing
2045 // left to do.
2046 if (affine_constraints.is_constrained(dofs_2[i]))
2047 continue;
2048
2049 // Enter the constraint piece by piece:
2050 affine_constraints.add_line(dofs_2[i]);
2051
2052 for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
2053 {
2054 // Get the correct dof index on face_1 respecting the
2055 // given orientation:
2056 const unsigned int j =
2057 cell_to_rotated_face_index[fe.face_to_cell_index(
2058 jj, 0, face_orientation, face_flip, face_rotation)];
2059
2060 if (std::abs(transformation(i, jj)) > eps)
2061 affine_constraints.add_entry(dofs_2[i],
2062 dofs_1[j],
2063 transformation(i, jj));
2064 }
2065
2066 // Continue with next dof.
2067 continue;
2068 }
2069
2070 // We are left with an "identity constraint".
2071
2072 // Get the correct dof index on face_1 respecting the given
2073 // orientation:
2074 const unsigned int j =
2075 cell_to_rotated_face_index[fe.face_to_cell_index(
2076 target, 0, face_orientation, face_flip, face_rotation)];
2077
2078 auto dof_left = dofs_1[j];
2079 auto dof_right = dofs_2[i];
2080
2081 // If dof_left is already constrained, or dof_left < dof_right we
2082 // flip the order to ensure that dofs are constrained in a stable
2083 // manner on different MPI processes.
2084 if (affine_constraints.is_constrained(dof_left) ||
2085 (dof_left < dof_right &&
2086 !affine_constraints.is_constrained(dof_right)))
2087 {
2088 std::swap(dof_left, dof_right);
2089 constraint_factor = 1. / constraint_factor;
2090 }
2091
2092 // Next, we try to enter the constraint
2093 // dof_left = constraint_factor * dof_right;
2094
2095 // If both degrees of freedom are constrained, there is nothing we
2096 // can do. Simply continue with the next dof.
2097 if (affine_constraints.is_constrained(dof_left) &&
2098 affine_constraints.is_constrained(dof_right))
2099 continue;
2100
2101 // We have to be careful that adding the current identity
2102 // constraint does not create a constraint cycle. Thus, check for
2103 // a dependency cycle:
2104
2105 bool constraints_are_cyclic = true;
2106 number cycle_constraint_factor = constraint_factor;
2107
2108 for (auto test_dof = dof_right; test_dof != dof_left;)
2109 {
2110 if (!affine_constraints.is_constrained(test_dof))
2111 {
2112 constraints_are_cyclic = false;
2113 break;
2114 }
2115
2116 const auto &constraint_entries =
2117 *affine_constraints.get_constraint_entries(test_dof);
2118 if (constraint_entries.size() == 1)
2119 {
2120 test_dof = constraint_entries[0].first;
2121 cycle_constraint_factor *= constraint_entries[0].second;
2122 }
2123 else
2124 {
2125 constraints_are_cyclic = false;
2126 break;
2127 }
2128 }
2129
2130 // In case of a dependency cycle we, either
2131 // - do nothing if cycle_constraint_factor == 1. In this case all
2132 // degrees
2133 // of freedom are already periodically constrained,
2134 // - otherwise, force all dofs to zero (by setting dof_left to
2135 // zero). The reasoning behind this is the fact that
2136 // cycle_constraint_factor != 1 occurs in situations such as
2137 // x1 == x2 and x2 == -1. * x1. This system is only solved by
2138 // x_1 = x_2 = 0.
2139
2140 if (constraints_are_cyclic)
2141 {
2142 if (std::abs(cycle_constraint_factor - 1.) > eps)
2143 affine_constraints.add_line(dof_left);
2144 }
2145 else
2146 {
2147 affine_constraints.add_line(dof_left);
2148 affine_constraints.add_entry(dof_left,
2149 dof_right,
2150 constraint_factor);
2151 // The number 1e10 in the assert below is arbitrary. If the
2152 // absolute value of constraint_factor is too large, then probably
2153 // the absolute value of periodicity_factor is too large or too
2154 // small. This would be equivalent to an evanescent wave that has
2155 // a very small wavelength. A quick calculation shows that if
2156 // |periodicity_factor| > 1e10 -> |np.exp(ikd)|> 1e10, therefore k
2157 // is imaginary (evanescent wave) and the evanescent wavelength is
2158 // 0.27 times smaller than the dimension of the structure,
2159 // lambda=((2*pi)/log(1e10))*d. Imaginary wavenumbers can be
2160 // interesting in some cases
2161 // (https://doi.org/10.1103/PhysRevA.94.033813).In order to
2162 // implement the case of in which the wavevector can be imaginary
2163 // it would be necessary to rewrite this function and the dof
2164 // ordering method should be modified.
2165 // Let's take the following constraint a*x1 + b*x2 = 0. You could
2166 // just always pick x1 = b/a*x2, but in practice this is not so
2167 // stable if a could be a small number -- intended to be zero, but
2168 // just very small due to roundoff. Of course, constraining x2 in
2169 // terms of x1 has the same problem. So one chooses x1 = b/a*x2 if
2170 // |b|<|a|, and x2 = a/b*x1 if |a|<|b|.
2171 Assert(
2172 std::abs(constraint_factor) < 1e10,
2173 ExcMessage(
2174 "The periodicity constraint is too large. The parameter periodicity_factor might be too large or too small."));
2175 }
2176 } /* for dofs_per_face */
2177 }
2178
2179
2180
2181 // Internally used in make_periodicity_constraints.
2182 //
2183 // Build up a (possibly rotated) interpolation matrix that is used in
2184 // set_periodicity_constraints with the help of user supplied matrix and
2185 // first_vector_components.
2186 template <int dim, int spacedim>
2187 FullMatrix<double>
2188 compute_transformation(
2189 const FiniteElement<dim, spacedim> &fe,
2190 const FullMatrix<double> & matrix,
2191 const std::vector<unsigned int> & first_vector_components)
2192 {
2193 // TODO: the implementation makes the assumption that all faces have the
2194 // same number of dofs
2195 AssertDimension(fe.n_unique_faces(), 1);
2196 const unsigned int face_no = 0;
2197
2198 Assert(matrix.m() == matrix.n(), ExcInternalError());
2199
2200 const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no);
2201
2202 if (matrix.m() == n_dofs_per_face)
2203 {
2204 // In case of m == n == n_dofs_per_face the supplied matrix is already
2205 // an interpolation matrix, so we use it directly:
2206 return matrix;
2207 }
2208
2209 if (first_vector_components.empty() && matrix.m() == 0)
2210 {
2211 // Just the identity matrix in case no rotation is specified:
2212 return IdentityMatrix(n_dofs_per_face);
2213 }
2214
2215 // The matrix describes a rotation and we have to build a transformation
2216 // matrix, we assume that for a 0* rotation we would have to build the
2217 // identity matrix
2218
2219 Assert(matrix.m() == spacedim, ExcInternalError())
2220
2221 Quadrature<dim - 1>
2222 quadrature(fe.get_unit_face_support_points(face_no));
2223
2224 // have an array that stores the location of each vector-dof tuple we want
2225 // to rotate.
2226 using DoFTuple = std::array<unsigned int, spacedim>;
2227
2228 // start with a pristine interpolation matrix...
2229 FullMatrix<double> transformation = IdentityMatrix(n_dofs_per_face);
2230
2231 for (unsigned int i = 0; i < n_dofs_per_face; ++i)
2232 {
2233 std::vector<unsigned int>::const_iterator comp_it =
2234 std::find(first_vector_components.begin(),
2235 first_vector_components.end(),
2236 fe.face_system_to_component_index(i, face_no).first);
2237 if (comp_it != first_vector_components.end())
2238 {
2239 const unsigned int first_vector_component = *comp_it;
2240
2241 // find corresponding other components of vector
2242 DoFTuple vector_dofs;
2243 vector_dofs[0] = i;
2244 unsigned int n_found = 1;
2245
2246 Assert(
2247 *comp_it + spacedim <= fe.n_components(),
2248 ExcMessage(
2249 "Error: the finite element does not have enough components "
2250 "to define rotated periodic boundaries."));
2251
2252 for (unsigned int k = 0; k < n_dofs_per_face; ++k)
2253 if ((k != i) && (quadrature.point(k) == quadrature.point(i)) &&
2254 (fe.face_system_to_component_index(k, face_no).first >=
2255 first_vector_component) &&
2256 (fe.face_system_to_component_index(k, face_no).first <
2257 first_vector_component + spacedim))
2258 {
2259 vector_dofs[fe.face_system_to_component_index(k, face_no)
2260 .first -
2261 first_vector_component] = k;
2262 n_found++;
2263 if (n_found == dim)
2264 break;
2265 }
2266
2267 // ... and rotate all dofs belonging to vector valued components
2268 // that are selected by first_vector_components:
2269 for (int i = 0; i < spacedim; ++i)
2270 {
2271 transformation[vector_dofs[i]][vector_dofs[i]] = 0.;
2272 for (int j = 0; j < spacedim; ++j)
2273 transformation[vector_dofs[i]][vector_dofs[j]] =
2274 matrix[i][j];
2275 }
2276 }
2277 }
2278 return transformation;
2279 }
2280 } /*namespace*/
2281
2282
2283 // Low level interface:
2284
2285
2286 template <typename FaceIterator, typename number>
2287 void
2288 make_periodicity_constraints(
2289 const FaceIterator & face_1,
2290 const typename identity<FaceIterator>::type &face_2,
2291 AffineConstraints<number> & affine_constraints,
2292 const ComponentMask & component_mask,
2293 const bool face_orientation,
2294 const bool face_flip,
2295 const bool face_rotation,
2296 const FullMatrix<double> & matrix,
2297 const std::vector<unsigned int> & first_vector_components,
2298 const number periodicity_factor)
2299 {
2300 // TODO: the implementation makes the assumption that all faces have the
2301 // same number of dofs
2302 AssertDimension(
2303 face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
2304 AssertDimension(
2305 face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
2306 const unsigned int face_no = 0;
2307
2308 static const int dim = FaceIterator::AccessorType::dimension;
2309 static const int spacedim = FaceIterator::AccessorType::space_dimension;
2310
2311 Assert((dim != 1) || (face_orientation == true && face_flip == false &&
2312 face_rotation == false),
2313 ExcMessage("The supplied orientation "
2314 "(face_orientation, face_flip, face_rotation) "
2315 "is invalid for 1D"));
2316
2317 Assert((dim != 2) || (face_orientation == true && face_rotation == false),
2318 ExcMessage("The supplied orientation "
2319 "(face_orientation, face_flip, face_rotation) "
2320 "is invalid for 2D"));
2321
2322 Assert(face_1 != face_2,
2323 ExcMessage("face_1 and face_2 are equal! Cannot constrain DoFs "
2324 "on the very same face"));
2325
2326 Assert(face_1->at_boundary() && face_2->at_boundary(),
2327 ExcMessage("Faces for periodicity constraints must be on the "
2328 "boundary"));
2329
2330 Assert(matrix.m() == matrix.n(),
2331 ExcMessage("The supplied (rotation or interpolation) matrix must "
2332 "be a square matrix"));
2333
2334 Assert(first_vector_components.empty() || matrix.m() == spacedim,
2335 ExcMessage("first_vector_components is nonempty, so matrix must "
2336 "be a rotation matrix exactly of size spacedim"));
2337
2338#ifdef DEBUG
2339 if (!face_1->has_children())
2340 {
2341 Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
2342 const unsigned int n_dofs_per_face =
2343 face_1->get_fe(face_1->nth_active_fe_index(0))
2344 .n_dofs_per_face(face_no);
2345
2346 Assert(matrix.m() == 0 ||
2347 (first_vector_components.empty() &&
2348 matrix.m() == n_dofs_per_face) ||
2349 (!first_vector_components.empty() && matrix.m() == spacedim),
2350 ExcMessage(
2351 "The matrix must have either size 0 or spacedim "
2352 "(if first_vector_components is nonempty) "
2353 "or the size must be equal to the # of DoFs on the face "
2354 "(if first_vector_components is empty)."));
2355 }
2356
2357 if (!face_2->has_children())
2358 {
2359 Assert(face_2->n_active_fe_indices() == 1, ExcInternalError());
2360 const unsigned int n_dofs_per_face =
2361 face_2->get_fe(face_2->nth_active_fe_index(0))
2362 .n_dofs_per_face(face_no);
2363
2364 Assert(matrix.m() == 0 ||
2365 (first_vector_components.empty() &&
2366 matrix.m() == n_dofs_per_face) ||
2367 (!first_vector_components.empty() && matrix.m() == spacedim),
2368 ExcMessage(
2369 "The matrix must have either size 0 or spacedim "
2370 "(if first_vector_components is nonempty) "
2371 "or the size must be equal to the # of DoFs on the face "
2372 "(if first_vector_components is empty)."));
2373 }
2374#endif
2375
2376 // A lookup table on how to go through the child faces depending on the
2377 // orientation:
2378
2379 static const int lookup_table_2d[2][2] = {
2380 // flip:
2381 {0, 1}, // false
2382 {1, 0}, // true
2383 };
2384
2385 static const int lookup_table_3d[2][2][2][4] = {
2386 // orientation flip rotation
2387 {
2388 {
2389 {0, 2, 1, 3}, // false false false
2390 {2, 3, 0, 1}, // false false true
2391 },
2392 {
2393 {3, 1, 2, 0}, // false true false
2394 {1, 0, 3, 2}, // false true true
2395 },
2396 },
2397 {
2398 {
2399 {0, 1, 2, 3}, // true false false
2400 {1, 3, 0, 2}, // true false true
2401 },
2402 {
2403 {3, 2, 1, 0}, // true true false
2404 {2, 0, 3, 1}, // true true true
2405 },
2406 },
2407 };
2408
2409 if (face_1->has_children() && face_2->has_children())
2410 {
2411 // In the case that both faces have children, we loop over all children
2412 // and apply make_periodicty_constrains recursively:
2413
2414 Assert(face_1->n_children() ==
2415 GeometryInfo<dim>::max_children_per_face &&
2416 face_2->n_children() ==
2417 GeometryInfo<dim>::max_children_per_face,
2418 ExcNotImplemented());
2419
2420 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face;
2421 ++i)
2422 {
2423 // Lookup the index for the second face
2424 unsigned int j;
2425 switch (dim)
2426 {
2427 case 2:
2428 j = lookup_table_2d[face_flip][i];
2429 break;
2430 case 3:
2431 j = lookup_table_3d[face_orientation][face_flip]
2432 [face_rotation][i];
2433 break;
2434 default:
2435 AssertThrow(false, ExcNotImplemented());
2436 }
2437
2438 make_periodicity_constraints(face_1->child(i),
2439 face_2->child(j),
2440 affine_constraints,
2441 component_mask,
2442 face_orientation,
2443 face_flip,
2444 face_rotation,
2445 matrix,
2446 first_vector_components,
2447 periodicity_factor);
2448 }
2449 }
2450 else
2451 {
2452 // Otherwise at least one of the two faces is active and we need to do
2453 // some work and enter the constraints!
2454
2455 // The finite element that matters is the one on the active face:
2456 const FiniteElement<dim, spacedim> &fe =
2457 face_1->has_children() ?
2458 face_2->get_fe(face_2->nth_active_fe_index(0)) :
2459 face_1->get_fe(face_1->nth_active_fe_index(0));
2460
2461 const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no);
2462
2463 // Sometimes we just have nothing to do (for all finite elements, or
2464 // systems which accidentally don't have any dofs on the boundary).
2465 if (n_dofs_per_face == 0)
2466 return;
2467
2468 const FullMatrix<double> transformation =
2469 compute_transformation(fe, matrix, first_vector_components);
2470
2471 if (!face_2->has_children())
2472 {
2473 // Performance hack: We do not need to compute an inverse if the
2474 // matrix is the identity matrix.
2475 if (first_vector_components.empty() && matrix.m() == 0)
2476 {
2477 set_periodicity_constraints(face_2,
2478 face_1,
2479 transformation,
2480 affine_constraints,
2481 component_mask,
2482 face_orientation,
2483 face_flip,
2484 face_rotation,
2485 periodicity_factor);
2486 }
2487 else
2488 {
2489 FullMatrix<double> inverse(transformation.m());
2490 inverse.invert(transformation);
2491
2492 set_periodicity_constraints(face_2,
2493 face_1,
2494 inverse,
2495 affine_constraints,
2496 component_mask,
2497 face_orientation,
2498 face_flip,
2499 face_rotation,
2500 periodicity_factor);
2501 }
2502 }
2503 else
2504 {
2505 Assert(!face_1->has_children(), ExcInternalError());
2506
2507 // Important note:
2508 // In 3D we have to take care of the fact that face_rotation gives
2509 // the relative rotation of face_1 to face_2, i.e. we have to invert
2510 // the rotation when constraining face_2 to face_1. Therefore
2511 // face_flip has to be toggled if face_rotation is true: In case of
2512 // inverted orientation, nothing has to be done.
2513 set_periodicity_constraints(face_1,
2514 face_2,
2515 transformation,
2516 affine_constraints,
2517 component_mask,
2518 face_orientation,
2519 face_orientation ?
2520 face_rotation ^ face_flip :
2521 face_flip,
2522 face_rotation,
2523 periodicity_factor);
2524 }
2525 }
2526 }
2527
2528
2529
2530 template <int dim, int spacedim, typename number>
2531 void
2532 make_periodicity_constraints(
2533 const std::vector<GridTools::PeriodicFacePair<
2534 typename DoFHandler<dim, spacedim>::cell_iterator>> &periodic_faces,
2535 AffineConstraints<number> & constraints,
2536 const ComponentMask & component_mask,
2537 const std::vector<unsigned int> &first_vector_components,
2538 const number periodicity_factor)
2539 {
2540 // Loop over all periodic faces...
2541 for (auto &pair : periodic_faces)
2542 {
2543 using FaceIterator = typename DoFHandler<dim, spacedim>::face_iterator;
2544 const FaceIterator face_1 = pair.cell[0]->face(pair.face_idx[0]);
2545 const FaceIterator face_2 = pair.cell[1]->face(pair.face_idx[1]);
2546
2547 Assert(face_1->at_boundary() && face_2->at_boundary(),
2548 ExcInternalError());
2549
2550 Assert(face_1 != face_2, ExcInternalError());
2551
2552 // ... and apply the low level make_periodicity_constraints function to
2553 // every matching pair:
2554 make_periodicity_constraints(face_1,
2555 face_2,
2556 constraints,
2557 component_mask,
2558 pair.orientation[0],
2559 pair.orientation[1],
2560 pair.orientation[2],
2561 pair.matrix,
2562 first_vector_components,
2563 periodicity_factor);
2564 }
2565 }
2566
2567
2568 // High level interface variants:
2569
2570
2571 template <int dim, int spacedim, typename number>
2572 void
2573 make_periodicity_constraints(const DoFHandler<dim, spacedim> & dof_handler,
2574 const types::boundary_id b_id1,
2575 const types::boundary_id b_id2,
2576 const unsigned int direction,
2577 ::AffineConstraints<number> &constraints,
2578 const ComponentMask &component_mask,
2579 const number periodicity_factor)
2580 {
2581 AssertIndexRange(direction, spacedim);
2582
2583 Assert(b_id1 != b_id2,
2584 ExcMessage("The boundary indicators b_id1 and b_id2 must be "
2585 "different to denote different boundaries."));
2586
2587 std::vector<GridTools::PeriodicFacePair<
2588 typename DoFHandler<dim, spacedim>::cell_iterator>>
2589 matched_faces;
2590
2591 // Collect matching periodic cells on the coarsest level:
2592 GridTools::collect_periodic_faces(
2593 dof_handler, b_id1, b_id2, direction, matched_faces);
2594
2595 make_periodicity_constraints<dim, spacedim>(matched_faces,
2596 constraints,
2597 component_mask,
2598 std::vector<unsigned int>(),
2599 periodicity_factor);
2600 }
2601
2602
2603
2604 template <int dim, int spacedim, typename number>
2605 void
2606 make_periodicity_constraints(const DoFHandler<dim, spacedim> &dof_handler,
2607 const types::boundary_id b_id,
2608 const unsigned int direction,
2609 AffineConstraints<number> & constraints,
2610 const ComponentMask & component_mask,
2611 const number periodicity_factor)
2612 {
2613 AssertIndexRange(direction, spacedim);
2614
2615 Assert(dim == spacedim, ExcNotImplemented());
2616
2617 std::vector<GridTools::PeriodicFacePair<
2618 typename DoFHandler<dim, spacedim>::cell_iterator>>
2619 matched_faces;
2620
2621 // Collect matching periodic cells on the coarsest level:
2622 GridTools::collect_periodic_faces(dof_handler,
2623 b_id,
2624 direction,
2625 matched_faces);
2626
2627 make_periodicity_constraints<dim, spacedim>(matched_faces,
2628 constraints,
2629 component_mask,
2630 std::vector<unsigned int>(),
2631 periodicity_factor);
2632 }
2633
2634
2635
2636 namespace internal
2637 {
2638 namespace Assembler
2639 {
2640 struct Scratch
2641 {};
2642
2643
2644 template <int dim, int spacedim>
2645 struct CopyData
2646 {
2647 unsigned int dofs_per_cell;
2648 std::vector<types::global_dof_index> parameter_dof_indices;
2649#ifdef DEAL_II_WITH_MPI
2650 std::vector<::LinearAlgebra::distributed::Vector<double>>
2651 global_parameter_representation;
2652#else
2653 std::vector<::Vector<double>> global_parameter_representation;
2654#endif
2655 };
2656 } // namespace Assembler
2657
2658 namespace
2659 {
2665 template <int dim, int spacedim>
2666 void
2667 compute_intergrid_weights_3(
2668 const typename ::DoFHandler<dim, spacedim>::active_cell_iterator
2669 &cell,
2670 const Assembler::Scratch &,
2671 Assembler::CopyData<dim, spacedim> &copy_data,
2672 const unsigned int coarse_component,
2673 const FiniteElement<dim, spacedim> &coarse_fe,
2674 const InterGridMap<::DoFHandler<dim, spacedim>>
2675 & coarse_to_fine_grid_map,
2676 const std::vector<::Vector<double>> &parameter_dofs)
2677 {
2678 // for each cell on the parameter grid: find out which degrees of
2679 // freedom on the fine grid correspond in which way to the degrees of
2680 // freedom on the parameter grid
2681 //
2682 // since for continuous FEs some dofs exist on more than one cell, we
2683 // have to track which ones were already visited. the problem is that if
2684 // we visit a dof first on one cell and compute its weight with respect
2685 // to some global dofs to be non-zero, and later visit the dof again on
2686 // another cell and (since we are on another cell) recompute the weights
2687 // with respect to the same dofs as above to be zero now, we have to
2688 // preserve them. we therefore overwrite all weights if they are nonzero
2689 // and do not enforce zero weights since that might be only due to the
2690 // fact that we are on another cell.
2691 //
2692 // example:
2693 // coarse grid
2694 // | | |
2695 // *-----*-----*
2696 // | cell|cell |
2697 // | 1 | 2 |
2698 // | | |
2699 // 0-----1-----*
2700 //
2701 // fine grid
2702 // | | | | |
2703 // *--*--*--*--*
2704 // | | | | |
2705 // *--*--*--*--*
2706 // | | | | |
2707 // *--x--y--*--*
2708 //
2709 // when on cell 1, we compute the weights of dof 'x' to be 1/2 from
2710 // parameter dofs 0 and 1, respectively. however, when later we are on
2711 // cell 2, we again compute the prolongation of shape function 1
2712 // restricted to cell 2 to the globla grid and find that the weight of
2713 // global dof 'x' now is zero. however, we should not overwrite the old
2714 // value.
2715 //
2716 // we therefore always only set nonzero values. why adding up is not
2717 // useful: dof 'y' would get weight 1 from parameter dof 1 on both cells
2718 // 1 and 2, but the correct weight is nevertheless only 1.
2719
2720 // vector to hold the representation of a single degree of freedom on
2721 // the coarse grid (for the selected fe) on the fine grid
2722
2723 copy_data.dofs_per_cell = coarse_fe.n_dofs_per_cell();
2724 copy_data.parameter_dof_indices.resize(copy_data.dofs_per_cell);
2725
2726 // get the global indices of the parameter dofs on this parameter grid
2727 // cell
2728 cell->get_dof_indices(copy_data.parameter_dof_indices);
2729
2730 // loop over all dofs on this cell and check whether they are
2731 // interesting for us
2732 for (unsigned int local_dof = 0; local_dof < copy_data.dofs_per_cell;
2733 ++local_dof)
2734 if (coarse_fe.system_to_component_index(local_dof).first ==
2735 coarse_component)
2736 {
2737 // the how-many-th parameter is this on this cell?
2738 const unsigned int local_parameter_dof =
2739 coarse_fe.system_to_component_index(local_dof).second;
2740
2741 copy_data.global_parameter_representation[local_parameter_dof] =
2742 0.;
2743
2744 // distribute the representation of @p{local_parameter_dof} on the
2745 // parameter grid cell
2746 // @p{cell} to the global data space
2747 coarse_to_fine_grid_map[cell]->set_dof_values_by_interpolation(
2748 parameter_dofs[local_parameter_dof],
2749 copy_data.global_parameter_representation[local_parameter_dof]);
2750 }
2751 }
2752
2753
2754
2760 template <int dim, int spacedim>
2761 void
2762 copy_intergrid_weights_3(
2763 const Assembler::CopyData<dim, spacedim> & copy_data,
2764 const unsigned int coarse_component,
2765 const FiniteElement<dim, spacedim> & coarse_fe,
2766 const std::vector<types::global_dof_index> &weight_mapping,
2767 const bool is_called_in_parallel,
2768 std::vector<std::map<types::global_dof_index, float>> &weights)
2769 {
2770 unsigned int pos = 0;
2771 for (unsigned int local_dof = 0; local_dof < copy_data.dofs_per_cell;
2772 ++local_dof)
2773 if (coarse_fe.system_to_component_index(local_dof).first ==
2774 coarse_component)
2775 {
2776 // now that we've got the global representation of each parameter
2777 // dof, we've only got to clobber the non-zero entries in that
2778 // vector and store the result
2779 //
2780 // what we have learned: if entry @p{i} of the global vector holds
2781 // the value @p{v[i]}, then this is the weight with which the
2782 // present dof contributes to @p{i}. there may be several such
2783 // @p{i}s and their weights' sum should be one. Then, @p{v[i]}
2784 // should be equal to @p{\sum_j w_{ij} p[j]} with @p{p[j]} be the
2785 // values of the degrees of freedom on the coarse grid. we can
2786 // thus compute constraints which link the degrees of freedom
2787 // @p{v[i]} on the fine grid to those on the coarse grid,
2788 // @p{p[j]}. Now to use these as real constraints, rather than as
2789 // additional equations, we have to identify representants among
2790 // the @p{i} for each @p{j}. this will be done by simply taking
2791 // the first @p{i} for which @p{w_{ij}==1}.
2792 //
2793 // guard modification of the weights array by a Mutex. since it
2794 // should happen rather rarely that there are several threads
2795 // operating on different intergrid weights, have only one mutex
2796 // for all of them
2797 for (types::global_dof_index i = 0;
2798 i < copy_data.global_parameter_representation[pos].size();
2799 ++i)
2800 // set this weight if it belongs to a parameter dof.
2801 if (weight_mapping[i] != numbers::invalid_dof_index)
2802 {
2803 // only overwrite old value if not by zero
2804 if (copy_data.global_parameter_representation[pos](i) != 0)
2805 {
2806 const types::global_dof_index
2807 wi = copy_data.parameter_dof_indices[local_dof],
2808 wj = weight_mapping[i];
2809 weights[wi][wj] =
2810 copy_data.global_parameter_representation[pos](i);
2811 }
2812 }
2813 else if (!is_called_in_parallel)
2814 {
2815 // Note that when this function operates with distributed
2816 // fine grid, this assertion is switched off since the
2817 // condition does not necessarily hold
2818 Assert(copy_data.global_parameter_representation[pos](i) ==
2819 0,
2820 ExcInternalError());
2821 }
2822
2823 ++pos;
2824 }
2825 }
2826
2827
2828
2834 template <int dim, int spacedim>
2835 void
2836 compute_intergrid_weights_2(
2837 const ::DoFHandler<dim, spacedim> &coarse_grid,
2838 const unsigned int coarse_component,
2839 const InterGridMap<::DoFHandler<dim, spacedim>>
2840 & coarse_to_fine_grid_map,
2841 const std::vector<::Vector<double>> & parameter_dofs,
2842 const std::vector<types::global_dof_index> &weight_mapping,
2843 std::vector<std::map<types::global_dof_index, float>> &weights)
2844 {
2845 Assembler::Scratch scratch;
2846 Assembler::CopyData<dim, spacedim> copy_data;
2847
2848 unsigned int n_interesting_dofs = 0;
2849 for (unsigned int local_dof = 0;
2850 local_dof < coarse_grid.get_fe().n_dofs_per_cell();
2851 ++local_dof)
2852 if (coarse_grid.get_fe().system_to_component_index(local_dof).first ==
2853 coarse_component)
2854 ++n_interesting_dofs;
2855
2856 copy_data.global_parameter_representation.resize(n_interesting_dofs);
2857
2858 bool is_called_in_parallel = false;
2859 for (std::size_t i = 0;
2860 i < copy_data.global_parameter_representation.size();
2861 ++i)
2862 {
2863#ifdef DEAL_II_WITH_MPI
2864 MPI_Comm communicator = MPI_COMM_SELF;
2865 try
2866 {
2867 const typename ::parallel::TriangulationBase<dim,
2868 spacedim>
2869 &tria = dynamic_cast<const typename ::parallel::
2870 TriangulationBase<dim, spacedim> &>(
2871 coarse_to_fine_grid_map.get_destination_grid()
2872 .get_triangulation());
2873 communicator = tria.get_communicator();
2874 is_called_in_parallel = true;
2875 }
2876 catch (std::bad_cast &)
2877 {
2878 // Nothing bad happened: the user used serial Triangulation
2879 }
2880
2881
2882 IndexSet locally_relevant_dofs;
2883 DoFTools::extract_locally_relevant_dofs(
2884 coarse_to_fine_grid_map.get_destination_grid(),
2885 locally_relevant_dofs);
2886
2887 copy_data.global_parameter_representation[i].reinit(
2888 coarse_to_fine_grid_map.get_destination_grid()
2889 .locally_owned_dofs(),
2890 locally_relevant_dofs,
2891 communicator);
2892#else
2893 const types::global_dof_index n_fine_dofs = weight_mapping.size();
2894 copy_data.global_parameter_representation[i].reinit(n_fine_dofs);
2895#endif
2896 }
2897
2898 auto worker =
2899 [coarse_component,
2900 &coarse_grid,
2901 &coarse_to_fine_grid_map,
2902 &parameter_dofs](const typename ::DoFHandler<dim, spacedim>::
2903 active_cell_iterator & cell,
2904 const Assembler::Scratch & scratch_data,
2905 Assembler::CopyData<dim, spacedim> &copy_data) {
2906 compute_intergrid_weights_3<dim, spacedim>(cell,
2907 scratch_data,
2908 copy_data,
2909 coarse_component,
2910 coarse_grid.get_fe(),
2911 coarse_to_fine_grid_map,
2912 parameter_dofs);
2913 };
2914
2915 auto copier =
2916 [coarse_component,
2917 &coarse_grid,
2918 &weight_mapping,
2919 is_called_in_parallel,
2920 &weights](const Assembler::CopyData<dim, spacedim> &copy_data) {
2921 copy_intergrid_weights_3<dim, spacedim>(copy_data,
2922 coarse_component,
2923 coarse_grid.get_fe(),
2924 weight_mapping,
2925 is_called_in_parallel,
2926 weights);
2927 };
2928
2929 WorkStream::run(coarse_grid.begin_active(),
2930 coarse_grid.end(),
2931 worker,
2932 copier,
2933 scratch,
2934 copy_data);
2935
2936#ifdef DEAL_II_WITH_MPI
2937 for (std::size_t i = 0;
2938 i < copy_data.global_parameter_representation.size();
2939 ++i)
2940 copy_data.global_parameter_representation[i].update_ghost_values();
2941#endif
2942 }
2943
2944
2945
2951 template <int dim, int spacedim>
2952 unsigned int
2953 compute_intergrid_weights_1(
2954 const ::DoFHandler<dim, spacedim> &coarse_grid,
2955 const unsigned int coarse_component,
2956 const ::DoFHandler<dim, spacedim> &fine_grid,
2957 const unsigned int fine_component,
2958 const InterGridMap<::DoFHandler<dim, spacedim>>
2959 &coarse_to_fine_grid_map,
2960 std::vector<std::map<types::global_dof_index, float>> &weights,
2961 std::vector<types::global_dof_index> & weight_mapping)
2962 {
2963 // aliases to the finite elements used by the dof handlers:
2964 const FiniteElement<dim, spacedim> &coarse_fe = coarse_grid.get_fe(),
2965 &fine_fe = fine_grid.get_fe();
2966
2967 // global numbers of dofs
2968 const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(),
2969 n_fine_dofs = fine_grid.n_dofs();
2970
2971 // local numbers of dofs
2972 const unsigned int fine_dofs_per_cell = fine_fe.n_dofs_per_cell();
2973
2974 // alias the number of dofs per cell belonging to the coarse_component
2975 // which is to be the restriction of the fine grid:
2976 const unsigned int coarse_dofs_per_cell_component =
2977 coarse_fe
2978 .base_element(
2979 coarse_fe.component_to_base_index(coarse_component).first)
2980 .n_dofs_per_cell();
2981
2982
2983 // Try to find out whether the grids stem from the same coarse grid.
2984 // This is a rather crude test, but better than nothing
2985 Assert(coarse_grid.get_triangulation().n_cells(0) ==
2986 fine_grid.get_triangulation().n_cells(0),
2987 ExcGridsDontMatch());
2988
2989 // check whether the map correlates the right objects
2990 Assert(&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid,
2991 ExcGridsDontMatch());
2992 Assert(&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid,
2993 ExcGridsDontMatch());
2994
2995
2996 // check whether component numbers are valid
2997 AssertIndexRange(coarse_component, coarse_fe.n_components());
2998 AssertIndexRange(fine_component, fine_fe.n_components());
2999
3000 // check whether respective finite elements are equal
3001 Assert(coarse_fe.base_element(
3002 coarse_fe.component_to_base_index(coarse_component).first) ==
3003 fine_fe.base_element(
3004 fine_fe.component_to_base_index(fine_component).first),
3005 ExcFiniteElementsDontMatch());
3006
3007#ifdef DEBUG
3008 // if in debug mode, check whether the coarse grid is indeed coarser
3009 // everywhere than the fine grid
3010 for (const auto &cell : coarse_grid.active_cell_iterators())
3011 Assert(cell->level() <= coarse_to_fine_grid_map[cell]->level(),
3012 ExcGridNotCoarser());
3013#endif
3014
3015 /*
3016 * From here on: the term `parameter' refers to the selected component
3017 * on the coarse grid and its analogon on the fine grid. The naming of
3018 * variables containing this term is due to the fact that
3019 * `selected_component' is longer, but also due to the fact that the
3020 * code of this function was initially written for a program where the
3021 * component which we wanted to match between grids was actually the
3022 * `parameter' variable.
3023 *
3024 * Likewise, the terms `parameter grid' and `state grid' refer to the
3025 * coarse and fine grids, respectively.
3026 *
3027 * Changing the names of variables would in principle be a good idea,
3028 * but would not make things simpler and would be another source of
3029 * errors. If anyone feels like doing so: patches would be welcome!
3030 */
3031
3032
3033
3034 // set up vectors of cell-local data; each vector represents one degree
3035 // of freedom of the coarse-grid variable in the fine-grid element
3036 std::vector<::Vector<double>> parameter_dofs(
3037 coarse_dofs_per_cell_component,
3038 ::Vector<double>(fine_dofs_per_cell));
3039 // for each coarse dof: find its position within the fine element and
3040 // set this value to one in the respective vector (all other values are
3041 // zero by construction)
3042 for (unsigned int local_coarse_dof = 0;
3043 local_coarse_dof < coarse_dofs_per_cell_component;
3044 ++local_coarse_dof)
3045 for (unsigned int fine_dof = 0; fine_dof < fine_fe.n_dofs_per_cell();
3046 ++fine_dof)
3047 if (fine_fe.system_to_component_index(fine_dof) ==
3048 std::make_pair(fine_component, local_coarse_dof))
3049 {
3050 parameter_dofs[local_coarse_dof](fine_dof) = 1.;
3051 break;
3052 }
3053
3054
3055 // find out how many DoFs there are on the grids belonging to the
3056 // components we want to match
3057 unsigned int n_parameters_on_fine_grid = 0;
3058 {
3059 // have a flag for each dof on the fine grid and set it to true if
3060 // this is an interesting dof. finally count how many true's there
3061 std::vector<bool> dof_is_interesting(fine_grid.n_dofs(), false);
3062 std::vector<types::global_dof_index> local_dof_indices(
3063 fine_fe.n_dofs_per_cell());
3064
3065 for (const auto &cell : fine_grid.active_cell_iterators())
3066 if (cell->is_locally_owned())
3067 {
3068 cell->get_dof_indices(local_dof_indices);
3069 for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i)
3070 if (fine_fe.system_to_component_index(i).first ==
3071 fine_component)
3072 dof_is_interesting[local_dof_indices[i]] = true;
3073 }
3074
3075 n_parameters_on_fine_grid = std::count(dof_is_interesting.begin(),
3076 dof_is_interesting.end(),
3077 true);
3078 }
3079
3080
3081 // set up the weights mapping
3082 weights.clear();
3083 weights.resize(n_coarse_dofs);
3084
3085 weight_mapping.clear();
3086 weight_mapping.resize(n_fine_dofs, numbers::invalid_dof_index);
3087
3088 {
3089 std::vector<types::global_dof_index> local_dof_indices(
3090 fine_fe.n_dofs_per_cell());
3091 unsigned int next_free_index = 0;
3092 for (const auto &cell : fine_grid.active_cell_iterators())
3093 if (cell->is_locally_owned())
3094 {
3095 cell->get_dof_indices(local_dof_indices);
3096 for (unsigned int i = 0; i < fine_fe.n_dofs_per_cell(); ++i)
3097 // if this DoF is a parameter dof and has not yet been
3098 // numbered, then do so
3099 if ((fine_fe.system_to_component_index(i).first ==
3100 fine_component) &&
3101 (weight_mapping[local_dof_indices[i]] ==
3102 numbers::invalid_dof_index))
3103 {
3104 weight_mapping[local_dof_indices[i]] = next_free_index;
3105 ++next_free_index;
3106 }
3107 }
3108
3109 Assert(next_free_index == n_parameters_on_fine_grid,
3110 ExcInternalError());
3111 }
3112
3113
3114 // for each cell on the parameter grid: find out which degrees of
3115 // freedom on the fine grid correspond in which way to the degrees of
3116 // freedom on the parameter grid
3117 //
3118 // do this in a separate function to allow for multithreading there. see
3119 // this function also if you want to read more information on the
3120 // algorithm used.
3121 compute_intergrid_weights_2(coarse_grid,
3122 coarse_component,
3123 coarse_to_fine_grid_map,
3124 parameter_dofs,
3125 weight_mapping,
3126 weights);
3127
3128
3129 // ok, now we have all weights for each dof on the fine grid. if in
3130 // debug mode lets see if everything went smooth, i.e. each dof has sum
3131 // of weights one
3132 //
3133 // in other words this means that if the sum of all shape functions on
3134 // the parameter grid is one (which is always the case), then the
3135 // representation on the state grid should be as well (division of
3136 // unity)
3137 //
3138 // if the parameter grid has more than one component, then the
3139 // respective dofs of the other components have sum of weights zero, of
3140 // course. we do not explicitly ask which component a dof belongs to,
3141 // but this at least tests some errors
3142#ifdef DEBUG
3143 for (unsigned int col = 0; col < n_parameters_on_fine_grid; ++col)
3144 {
3145 double sum = 0;
3146 for (types::global_dof_index row = 0; row < n_coarse_dofs; ++row)
3147 if (weights[row].find(col) != weights[row].end())
3148 sum += weights[row][col];
3149 Assert((std::fabs(sum - 1) < 1.e-12) ||
3150 ((coarse_fe.n_components() > 1) && (sum == 0)),
3151 ExcInternalError());
3152 }
3153#endif
3154
3155
3156 return n_parameters_on_fine_grid;
3157 }
3158
3159
3160 } // namespace
3161 } // namespace internal
3162
3163
3164
3165 template <int dim, int spacedim>
3166 void
3167 compute_intergrid_constraints(
3168 const DoFHandler<dim, spacedim> & coarse_grid,
3169 const unsigned int coarse_component,
3170 const DoFHandler<dim, spacedim> & fine_grid,
3171 const unsigned int fine_component,
3172 const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map,
3173 AffineConstraints<double> & constraints)
3174 {
3175 Assert(coarse_grid.get_fe_collection().size() == 1 &&
3176 fine_grid.get_fe_collection().size() == 1,
3177 ExcMessage("This function is not yet implemented for DoFHandlers "
3178 "using hp-capabilities."));
3179 // store the weights with which a dof on the parameter grid contributes to a
3180 // dof on the fine grid. see the long doc below for more info
3181 //
3182 // allocate as many rows as there are parameter dofs on the coarse grid and
3183 // as many columns as there are parameter dofs on the fine grid.
3184 //
3185 // weight_mapping is used to map the global (fine grid) parameter dof
3186 // indices to the columns
3187 //
3188 // in the original implementation, the weights array was actually of
3189 // FullMatrix<double> type. this wasted huge amounts of memory, but was
3190 // fast. nonetheless, since the memory consumption was quadratic in the
3191 // number of degrees of freedom, this was not very practical, so we now use
3192 // a vector of rows of the matrix, and in each row a vector of pairs
3193 // (colnum,value). this seems like the best tradeoff between memory and
3194 // speed, as it is now linear in memory and still fast enough.
3195 //
3196 // to save some memory and since the weights are usually (negative) powers
3197 // of 2, we choose the value type of the matrix to be @p{float} rather than
3198 // @p{double}.
3199 std::vector<std::map<types::global_dof_index, float>> weights;
3200
3201 // this is this mapping. there is one entry for each dof on the fine grid;
3202 // if it is a parameter dof, then its value is the column in weights for
3203 // that parameter dof, if it is any other dof, then its value is -1,
3204 // indicating an error
3205 std::vector<types::global_dof_index> weight_mapping;
3206
3207 const unsigned int n_parameters_on_fine_grid =
3208 internal::compute_intergrid_weights_1(coarse_grid,
3209 coarse_component,
3210 fine_grid,
3211 fine_component,
3212 coarse_to_fine_grid_map,
3213 weights,
3214 weight_mapping);
3215 (void)n_parameters_on_fine_grid;
3216
3217 // global numbers of dofs
3218 const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(),
3219 n_fine_dofs = fine_grid.n_dofs();
3220
3221
3222 // get an array in which we store which dof on the coarse grid is a
3223 // parameter and which is not
3224 IndexSet coarse_dof_is_parameter;
3225 {
3226 std::vector<bool> mask(coarse_grid.get_fe(0).n_components(), false);
3227 mask[coarse_component] = true;
3228
3229 coarse_dof_is_parameter =
3230 extract_dofs<dim, spacedim>(coarse_grid, ComponentMask(mask));
3231 }
3232
3233 // now we know that the weights in each row constitute a constraint. enter
3234 // this into the constraints object
3235 //
3236 // first task: for each parameter dof on the parameter grid, find a
3237 // representant on the fine, global grid. this is possible since we use
3238 // conforming finite element. we take this representant to be the first
3239 // element in this row with weight identical to one. the representant will
3240 // become an unconstrained degree of freedom, while all others will be
3241 // constrained to this dof (and possibly others)
3242 std::vector<types::global_dof_index> representants(
3243 n_coarse_dofs, numbers::invalid_dof_index);
3244 for (types::global_dof_index parameter_dof = 0;
3245 parameter_dof < n_coarse_dofs;
3246 ++parameter_dof)
3247 if (coarse_dof_is_parameter.is_element(parameter_dof))
3248 {
3249 // if this is the line of a parameter dof on the coarse grid, then it
3250 // should have at least one dependent node on the fine grid
3251 Assert(weights[parameter_dof].size() > 0, ExcInternalError());
3252
3253 // find the column where the representant is mentioned
3254 std::map<types::global_dof_index, float>::const_iterator i =
3255 weights[parameter_dof].begin();
3256 for (; i != weights[parameter_dof].end(); ++i)
3257 if (i->second == 1)
3258 break;
3259 Assert(i != weights[parameter_dof].end(), ExcInternalError());
3260 const types::global_dof_index column = i->first;
3261
3262 // now we know in which column of weights the representant is, but we
3263 // don't know its global index. get it using the inverse operation of
3264 // the weight_mapping
3265 types::global_dof_index global_dof = 0;
3266 for (; global_dof < weight_mapping.size(); ++global_dof)
3267 if (weight_mapping[global_dof] ==
3268 static_cast<types::global_dof_index>(column))
3269 break;
3270 Assert(global_dof < weight_mapping.size(), ExcInternalError());
3271
3272 // now enter the representants global index into our list
3273 representants[parameter_dof] = global_dof;
3274 }
3275 else
3276 {
3277 // consistency check: if this is no parameter dof on the coarse grid,
3278 // then the respective row must be empty!
3279 Assert(weights[parameter_dof].size() == 0, ExcInternalError());
3280 }
3281
3282
3283
3284 // note for people that want to optimize this function: the largest part of
3285 // the computing time is spent in the following, rather innocent block of
3286 // code. basically, it must be the AffineConstraints::add_entry call which
3287 // takes the bulk of the time, but it is not known to the author how to make
3288 // it faster...
3289 std::vector<std::pair<types::global_dof_index, double>> constraint_line;
3290 for (types::global_dof_index global_dof = 0; global_dof < n_fine_dofs;
3291 ++global_dof)
3292 if (weight_mapping[global_dof] != numbers::invalid_dof_index)
3293 // this global dof is a parameter dof, so it may carry a constraint note
3294 // that for each global dof, the sum of weights shall be one, so we can
3295 // find out whether this dof is constrained in the following way: if the
3296 // only weight in this row is a one, and the representant for the
3297 // parameter dof of the line in which this one is is the present dof,
3298 // then we consider this dof to be unconstrained. otherwise, all other
3299 // dofs are constrained
3300 {
3301 const types::global_dof_index col = weight_mapping[global_dof];
3302 Assert(col < n_parameters_on_fine_grid, ExcInternalError());
3303
3304 types::global_dof_index first_used_row = 0;
3305
3306 {
3307 Assert(weights.size() > 0, ExcInternalError());
3308 std::map<types::global_dof_index, float>::const_iterator col_entry =
3309 weights[0].end();
3310 for (; first_used_row < n_coarse_dofs; ++first_used_row)
3311 {
3312 col_entry = weights[first_used_row].find(col);
3313 if (col_entry != weights[first_used_row].end())
3314 break;
3315 }
3316
3317 Assert(col_entry != weights[first_used_row].end(),
3318 ExcInternalError());
3319
3320 if ((col_entry->second == 1) &&
3321 (representants[first_used_row] == global_dof))
3322 // dof unconstrained or constrained to itself (in case this cell
3323 // is mapped to itself, rather than to children of itself)
3324 continue;
3325 }
3326
3327
3328 // otherwise enter all constraints
3329 constraints.add_line(global_dof);
3330
3331 constraint_line.clear();
3332 for (types::global_dof_index row = first_used_row;
3333 row < n_coarse_dofs;
3334 ++row)
3335 {
3336 const std::map<types::global_dof_index, float>::const_iterator j =
3337 weights[row].find(col);
3338 if ((j != weights[row].end()) && (j->second != 0))
3339 constraint_line.emplace_back(representants[row], j->second);
3340 }
3341
3342 constraints.add_entries(global_dof, constraint_line);
3343 }
3344 }
3345
3346
3347
3348 template <int dim, int spacedim>
3349 void
3350 compute_intergrid_transfer_representation(
3351 const DoFHandler<dim, spacedim> & coarse_grid,
3352 const unsigned int coarse_component,
3353 const DoFHandler<dim, spacedim> & fine_grid,
3354 const unsigned int fine_component,
3355 const InterGridMap<DoFHandler<dim, spacedim>> &coarse_to_fine_grid_map,
3356 std::vector<std::map<types::global_dof_index, float>>
3357 &transfer_representation)
3358 {
3359 Assert(coarse_grid.get_fe_collection().size() == 1 &&
3360 fine_grid.get_fe_collection().size() == 1,
3361 ExcMessage("This function is not yet implemented for DoFHandlers "
3362 "using hp-capabilities."));
3363 // store the weights with which a dof on the parameter grid contributes to a
3364 // dof on the fine grid. see the long doc below for more info
3365 //
3366 // allocate as many rows as there are parameter dofs on the coarse grid and
3367 // as many columns as there are parameter dofs on the fine grid.
3368 //
3369 // weight_mapping is used to map the global (fine grid) parameter dof
3370 // indices to the columns
3371 //
3372 // in the original implementation, the weights array was actually of
3373 // FullMatrix<double> type. this wasted huge amounts of memory, but was
3374 // fast. nonetheless, since the memory consumption was quadratic in the
3375 // number of degrees of freedom, this was not very practical, so we now use
3376 // a vector of rows of the matrix, and in each row a vector of pairs
3377 // (colnum,value). this seems like the best tradeoff between memory and
3378 // speed, as it is now linear in memory and still fast enough.
3379 //
3380 // to save some memory and since the weights are usually (negative) powers
3381 // of 2, we choose the value type of the matrix to be @p{float} rather than
3382 // @p{double}.
3383 std::vector<std::map<types::global_dof_index, float>> weights;
3384
3385 // this is this mapping. there is one entry for each dof on the fine grid;
3386 // if it is a parameter dof, then its value is the column in weights for
3387 // that parameter dof, if it is any other dof, then its value is -1,
3388 // indicating an error
3389 std::vector<types::global_dof_index> weight_mapping;
3390
3391 internal::compute_intergrid_weights_1(coarse_grid,
3392 coarse_component,
3393 fine_grid,
3394 fine_component,
3395 coarse_to_fine_grid_map,
3396 weights,
3397 weight_mapping);
3398
3399 // now compute the requested representation
3400 const types::global_dof_index n_global_parm_dofs =
3401 std::count_if(weight_mapping.begin(),
3402 weight_mapping.end(),
3403 [](const types::global_dof_index dof) {
3404 return dof != numbers::invalid_dof_index;
3405 });
3406
3407 // first construct the inverse mapping of weight_mapping
3408 std::vector<types::global_dof_index> inverse_weight_mapping(
3409 n_global_parm_dofs, numbers::invalid_dof_index);
3410 for (types::global_dof_index i = 0; i < weight_mapping.size(); ++i)
3411 {
3412 const types::global_dof_index parameter_dof = weight_mapping[i];
3413 // if this global dof is a parameter
3414 if (parameter_dof != numbers::invalid_dof_index)
3415 {
3416 Assert(parameter_dof < n_global_parm_dofs, ExcInternalError());
3417 Assert((inverse_weight_mapping[parameter_dof] ==
3418 numbers::invalid_dof_index),
3419 ExcInternalError());
3420
3421 inverse_weight_mapping[parameter_dof] = i;
3422 }
3423 }
3424
3425 // next copy over weights array and replace respective numbers
3426 const types::global_dof_index n_rows = weight_mapping.size();
3427
3428 transfer_representation.clear();
3429 transfer_representation.resize(n_rows);
3430
3431 const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs();
3432 for (types::global_dof_index i = 0; i < n_coarse_dofs; ++i)
3433 {
3434 std::map<types::global_dof_index, float>::const_iterator j =
3435 weights[i].begin();
3436 for (; j != weights[i].end(); ++j)
3437 {
3438 const types::global_dof_index p = inverse_weight_mapping[j->first];
3439 Assert(p < n_rows, ExcInternalError());
3440
3441 transfer_representation[p][i] = j->second;
3442 }
3443 }
3444 }
3445
3446
3447
3448 template <int dim, int spacedim, typename number>
3449 void
3450 make_zero_boundary_constraints(
3451 const DoFHandler<dim, spacedim> &dof,
3452 const types::boundary_id boundary_id,
3453 AffineConstraints<number> & zero_boundary_constraints,
3454 const ComponentMask & component_mask)
3455 {
3456 Assert(component_mask.represents_n_components(dof.get_fe(0).n_components()),
3457 ExcMessage("The number of components in the mask has to be either "
3458 "zero or equal to the number of components in the finite "
3459 "element."));
3460
3461 const unsigned int n_components = dof.get_fe_collection().n_components();
3462
3463 Assert(component_mask.n_selected_components(n_components) > 0,
3464 ComponentMask::ExcNoComponentSelected());
3465
3466 // a field to store the indices on the face
3467 std::vector<types::global_dof_index> face_dofs;
3468 face_dofs.reserve(dof.get_fe_collection().max_dofs_per_face());
3469 // a field to store the indices on the cell
3470 std::vector<types::global_dof_index> cell_dofs;
3471 cell_dofs.reserve(dof.get_fe_collection().max_dofs_per_cell());
3472
3473 typename DoFHandler<dim, spacedim>::active_cell_iterator
3474 cell = dof.begin_active(),
3475 endc = dof.end();
3476 for (; cell != endc; ++cell)
3477 if (!cell->is_artificial() && cell->at_boundary())
3478 {
3479 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
3480
3481 // get global indices of dofs on the cell
3482 cell_dofs.resize(fe.n_dofs_per_cell());
3483 cell->get_dof_indices(cell_dofs);
3484
3485 for (const auto face_no : cell->face_indices())
3486 {
3487 const typename DoFHandler<dim, spacedim>::face_iterator face =
3488 cell->face(face_no);
3489
3490 // if face is on the boundary and satisfies the correct boundary
3491 // id property
3492 if (face->at_boundary() &&
3493 ((boundary_id == numbers::invalid_boundary_id) ||
3494 (face->boundary_id() == boundary_id)))
3495 {
3496 // get indices and physical location on this face
3497 face_dofs.resize(fe.n_dofs_per_face(face_no));
3498 face->get_dof_indices(face_dofs, cell->active_fe_index());
3499
3500 // enter those dofs into the list that match the component
3501 // signature.
3502 for (const types::global_dof_index face_dof : face_dofs)
3503 {
3504 // Find out if a dof has a contribution in this component,
3505 // and if so, add it to the list
3506 const std::vector<types::global_dof_index>::iterator
3507 it_index_on_cell = std::find(cell_dofs.begin(),
3508 cell_dofs.end(),
3509 face_dof);
3510 Assert(it_index_on_cell != cell_dofs.end(),
3511 ExcInvalidIterator());
3512 const unsigned int index_on_cell =
3513 std::distance(cell_dofs.begin(), it_index_on_cell);
3514 const ComponentMask &nonzero_component_array =
3515 cell->get_fe().get_nonzero_components(index_on_cell);
3516 bool nonzero = false;
3517 for (unsigned int c = 0; c < n_components; ++c)
3518 if (nonzero_component_array[c] && component_mask[c])
3519 {
3520 nonzero = true;
3521 break;
3522 }
3523
3524 if (nonzero)
3525 zero_boundary_constraints.add_line(face_dof);
3526 }
3527 }
3528 }
3529 }
3530 }
3531
3532
3533
3534 template <int dim, int spacedim, typename number>
3535 void
3536 make_zero_boundary_constraints(
3537 const DoFHandler<dim, spacedim> &dof,
3538 AffineConstraints<number> & zero_boundary_constraints,
3539 const ComponentMask & component_mask)
3540 {
3541 make_zero_boundary_constraints(dof,
3542 numbers::invalid_boundary_id,
3543 zero_boundary_constraints,
3544 component_mask);
3545 }
3546
3547
3548} // end of namespace DoFTools
3549
3550
3551
3552// explicit instantiations
3553
3554#include "dof_tools_constraints.inst"
3555
3556
3557
3558DEAL_II_NAMESPACE_CLOSE
affine_constraints.h
AffineConstraints::add_entries
void add_entries(const size_type constrained_dof_index, const std::vector< std::pair< size_type, number > > &col_weight_pairs)
AffineConstraints::add_line
void add_line(const size_type line_n)
Definition: affine_constraints.h:2010
AffineConstraints::add_entry
void add_entry(const size_type constrained_dof_index, const size_type column, const number weight)
Definition: affine_constraints.h:2042
AffineConstraints::set_inhomogeneity
void set_inhomogeneity(const size_type constrained_dof_index, const number value)
Definition: affine_constraints.h:2084
AffineConstraints::is_constrained
bool is_constrained(const size_type line_n) const
Definition: affine_constraints.h:2122
AffineConstraints::get_constraint_entries
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
Definition: affine_constraints.h:2149
ComponentMask::represents_n_components
bool represents_n_components(const unsigned int n) const
ComponentMask::n_selected_components
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
DoFHandler::end
cell_iterator end() const
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
DoFHandler::get_fe_collection
const hp::FECollection< dim, spacedim > & get_fe_collection() const
DoFHandler::has_active_dofs
bool has_active_dofs() const
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
DoFHandler::has_hp_capabilities
bool has_hp_capabilities() const
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
FiniteElementData::n_dofs_per_vertex
unsigned int n_dofs_per_vertex() const
FiniteElementData::n_dofs_per_cell
unsigned int n_dofs_per_cell() const
FiniteElementData::n_dofs_per_line
unsigned int n_dofs_per_line() const
FiniteElementData::n_dofs_per_face
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
FiniteElementData::n_components
unsigned int n_components() const
FiniteElementData::n_unique_faces
unsigned int n_unique_faces() const
FiniteElementData::n_dofs_per_quad
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
FiniteElement::get_nonzero_components
const ComponentMask & get_nonzero_components(const unsigned int i) const
FiniteElement::base_element
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
FiniteElement::component_to_base_index
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
FiniteElement::face_to_cell_index
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
FiniteElement::get_subface_interpolation_matrix
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const
FiniteElement::system_to_component_index
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
FiniteElement::constraints
const FullMatrix< double > & constraints(const ::internal::SubfaceCase< dim > &subface_case=::internal::SubfaceCase< dim >::case_isotropic) const
FiniteElement::get_face_interpolation_matrix
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const
FiniteElement::face_system_to_component_index
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
FullMatrix::mmult
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
FullMatrix::n
size_type n() const
FullMatrix::invert
void invert(const FullMatrix< number2 > &M)
FullMatrix::m
size_type m() const
IndexSet::is_element
bool is_element(const size_type index) const
Definition: index_set.h:1765
MPI_Comm
Triangulation::get_anisotropic_refinement_flag
bool get_anisotropic_refinement_flag() const
Vector
Definition: vector.h:110
hp::Collection::size
unsigned int size() const
Definition: collection.h:109
hp::FECollection::find_dominating_fe_extended
unsigned int find_dominating_fe_extended(const std::set< unsigned int > &fes, const unsigned int codim=0) const
Definition: fe_collection.cc:201
hp::FECollection::hp_constraints_are_implemented
bool hp_constraints_are_implemented() const
Definition: fe_collection.h:939
hp::FECollection::max_dofs_per_face
unsigned int max_dofs_per_face() const
Definition: fe_collection.h:910
hp::FECollection::n_components
unsigned int n_components() const
Definition: fe_collection.h:788
hp::FECollection::max_dofs_per_cell
unsigned int max_dofs_per_cell() const
Definition: fe_collection.h:925
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
dof_accessor.h
dof_tools.h
dof_handler.h
fe_values.h
fe_collection.h
fe_tools.h
cell_index
unsigned int cell_index
Definition: grid_tools.cc:1092
grid_tools.h
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
StandardExceptions::ExcInvalidIterator
static ::ExceptionBase & ExcInvalidIterator()
DoFTools::ExcGridNotCoarser
static ::ExceptionBase & ExcGridNotCoarser()
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
DataOutBase::eps
@ eps
Definition: data_out_base.h:1580
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
DoFTools::ExcFiniteElementsDontMatch
static ::ExceptionBase & ExcFiniteElementsDontMatch()
ComponentMask::ExcNoComponentSelected
static ::ExceptionBase & ExcNoComponentSelected()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
DoFTools::ExcGridsDontMatch
static ::ExceptionBase & ExcGridsDontMatch()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
DoFHandler::line_iterator
typename ActiveSelector::line_iterator line_iterator
Definition: dof_handler.h:357
DoFHandler::face_iterator
typename ActiveSelector::face_iterator face_iterator
Definition: dof_handler.h:484
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
DoFTools::compute_intergrid_transfer_representation
void compute_intergrid_transfer_representation(const DoFHandler< dim, spacedim > &coarse_grid, const unsigned int coarse_component, const DoFHandler< dim, spacedim > &fine_grid, const unsigned int fine_component, const InterGridMap< DoFHandler< dim, spacedim > > &coarse_to_fine_grid_map, std::vector< std::map< types::global_dof_index, float > > &transfer_representation)
Definition: dof_tools_constraints.cc:3350
DoFTools::make_zero_boundary_constraints
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask=ComponentMask())
Definition: dof_tools_constraints.cc:3450
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
Definition: dof_tools_constraints.cc:1787
DoFTools::compute_intergrid_constraints
void compute_intergrid_constraints(const DoFHandler< dim, spacedim > &coarse_grid, const unsigned int coarse_component, const DoFHandler< dim, spacedim > &fine_grid, const unsigned int fine_component, const InterGridMap< DoFHandler< dim, spacedim > > &coarse_to_fine_grid_map, AffineConstraints< double > &constraints)
Definition: dof_tools_constraints.cc:3167
fe_values.h
intergrid_map.h
la_parallel_vector.h
Differentiation::SD::fabs
Expression fabs(const Expression &x)
Definition: symengine_math.cc:273
DoFTools::internal::make_hp_hanging_node_constraints
void make_hp_hanging_node_constraints(const ::DoFHandler< 1 > &, AffineConstraints< number > &)
Definition: dof_tools_constraints.cc:579
DoFTools::internal::make_oldstyle_hanging_node_constraints
void make_oldstyle_hanging_node_constraints(const ::DoFHandler< 1 > &, AffineConstraints< number > &, std::integral_constant< int, 1 >)
Definition: dof_tools_constraints.cc:588
DoFTools::make_periodicity_constraints
void make_periodicity_constraints(const FaceIterator &face_1, const typename identity< FaceIterator >::type &face_2, AffineConstraints< number > &constraints, const ComponentMask &component_mask=ComponentMask(), 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.)
Definition: dof_tools_constraints.cc:2288
DoFTools::extract_locally_relevant_dofs
void extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1210
GridTools::collect_periodic_faces
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id1, const types::boundary_id b_id2, const 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 >())
Definition: grid_tools_dof_handlers.cc:2189
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:58
LAPACKSupport::N
static const char N
Definition: lapack_support.h:157
MemorySpace::swap
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
Definition: memory_space_data.h:92
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:67
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
WorkStream::run
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1252
numbers::invalid_boundary_id
const types::boundary_id invalid_boundary_id
Definition: types.h:239
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:196
numbers::invalid_dof_index
const types::global_dof_index invalid_dof_index
Definition: types.h:211
std
STL namespace.
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5468
types::boundary_id
unsigned int boundary_id
Definition: types.h:129
DoFTools::internal::Assembler::CopyData::global_parameter_representation
std::vector<::LinearAlgebra::distributed::Vector< double > > global_parameter_representation
Definition: dof_tools_constraints.cc:2651
DoFTools::internal::Assembler::CopyData::parameter_dof_indices
std::vector< types::global_dof_index > parameter_dof_indices
Definition: dof_tools_constraints.cc:2648
template_constraints.h
tria_iterator.h
work_stream.h
Generated by doxygen 1.9.4