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