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