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