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