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