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