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