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