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