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