Reference documentation for deal.II version Git 2618e0f 2017-11-23 17:25:26 +0100
dof_tools_constraints.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2017 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 
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);
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);
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> &,
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> &,
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> &,
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> &,
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> &,
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> &,
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,
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,
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,
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,
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,
871  Assert (cell->face(face)->child(c)->line(e)
872  ->fe_index_is_active(cell->active_fe_index()) == true,
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,
881  Assert (cell->face(face)->line(e)
882  ->fe_index_is_active(cell->active_fe_index()) == true,
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
1055  face_interpolation_matrices (n_finite_elements (dof_handler),
1056  n_finite_elements (dof_handler));
1058  subface_interpolation_matrices (n_finite_elements (dof_handler),
1059  n_finite_elements (dof_handler),
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"
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.
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
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  {
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  ==
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 
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 
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  {
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 
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 
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 
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 
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 (!constraint_matrix.is_constrained(dofs_2[i]))
1877  if ((component_mask.n_selected_components(fe.n_components())
1878  == fe.n_components())
1879  ||
1880  component_mask[fe.face_system_to_component_index(i).first])
1881  {
1882  // as mentioned in the comment above this function, we need
1883  // to be careful about treating identity constraints differently.
1884  // consequently, find out whether this dof 'i' will be
1885  // identity constrained
1886  //
1887  // to check whether this is the case, first see whether there are
1888  // any weights other than 0 and 1, then in a first stage make sure
1889  // that if so there is only one weight equal to 1
1890  //
1891  // afterwards do the same for constraints of type dof1=-dof2
1892  bool is_identity_constrained = true;
1893  const double eps = 1.e-13;
1894  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
1895  if (((std::abs(transformation(i,jj)) < eps) ||
1896  (std::abs(transformation(i,jj)-1) < eps)) == false)
1897  {
1898  is_identity_constrained = false;
1899  break;
1900  }
1901  unsigned int identity_constraint_target = numbers::invalid_unsigned_int;
1902  if (is_identity_constrained == true)
1903  {
1904  bool one_identity_found = false;
1905  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
1906  if (std::abs(transformation(i,jj)-1.) < eps)
1907  {
1908  if (one_identity_found == false)
1909  {
1910  one_identity_found = true;
1911  identity_constraint_target = jj;
1912  }
1913  else
1914  {
1915  is_identity_constrained = false;
1916  identity_constraint_target = numbers::invalid_unsigned_int;
1917  break;
1918  }
1919  }
1920  }
1921 
1922  bool is_inverse_constrained = !is_identity_constrained;
1923  unsigned int inverse_constraint_target = numbers::invalid_unsigned_int;
1924  if (is_inverse_constrained)
1925  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
1926  if (((std::abs(transformation(i,jj)) < eps) ||
1927  (std::abs(transformation(i,jj)+1) < eps)) == false)
1928  {
1929  is_inverse_constrained = false;
1930  break;
1931  }
1932  if (is_inverse_constrained)
1933  {
1934  bool one_identity_found = false;
1935  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
1936  if (std::abs(transformation(i,jj)+1) < eps)
1937  {
1938  if (one_identity_found == false)
1939  {
1940  one_identity_found = true;
1941  inverse_constraint_target = jj;
1942  }
1943  else
1944  {
1945  is_inverse_constrained = false;
1946  inverse_constraint_target = numbers::invalid_unsigned_int;
1947  break;
1948  }
1949  }
1950  }
1951 
1952  const unsigned int target = is_identity_constrained
1953  ? identity_constraint_target
1954  : inverse_constraint_target;
1955 
1956  // find out whether this dof also exists on face 1
1957  // if this is true and the constraint is no identity
1958  // constraint to itself, set it to zero
1959  bool constrained_set = false;
1960  for (unsigned int j=0; j<dofs_per_face; ++j)
1961  {
1962  if (dofs_2[i] == dofs_1[j])
1963  if (!(is_identity_constrained && target==i))
1964  {
1965  constraint_matrix.add_line(dofs_2[i]);
1966  constrained_set = true;
1967  }
1968  }
1969 
1970  if (!constrained_set)
1971  {
1972  // now treat constraints, either as an equality constraint or
1973  // as a sequence of constraints
1974  if (is_identity_constrained == true || is_inverse_constrained == true)
1975  {
1976  // Query the correct face_index on face_1 respecting the given
1977  // orientation:
1978  const unsigned int j
1979  = cell_to_rotated_face_index[fe.face_to_cell_index(target,
1980  0, /* It doesn't really matter, just assume
1981  * we're on the first face...
1982  */
1983  face_orientation, face_flip, face_rotation)];
1984 
1985  // if the two aren't already identity constrained (whichever way
1986  // around) or already identical (in case of rotated periodicity constraints),
1987  // then enter the constraint. otherwise there is nothing for us still to do
1988  bool enter_constraint = false;
1989  if (!constraint_matrix.is_constrained(dofs_1[j]))
1990  {
1991  if (dofs_2[i] != dofs_1[j])
1992  enter_constraint = true;
1993  }
1994  else //dofs_1[j] is constrained, is it identity or inverse constrained?
1995  {
1996  const std::vector<std::pair<types::global_dof_index, double > > *constraint_entries
1997  = constraint_matrix.get_constraint_entries(dofs_1[j]);
1998  if (constraint_entries->size()==1 && (*constraint_entries)[0].first == dofs_2[i])
1999  {
2000  if ((is_identity_constrained && std::abs((*constraint_entries)[0].second-1) > eps) ||
2001  (is_inverse_constrained && std::abs((*constraint_entries)[0].second+1) > eps))
2002  {
2003  //this pair of constraints means that both dofs have to be constrained to 0.
2004  constraint_matrix.add_line(dofs_2[i]);
2005  }
2006  }
2007  else
2008  enter_constraint = true;
2009  }
2010 
2011  if (enter_constraint)
2012  {
2013  constraint_matrix.add_line(dofs_2[i]);
2014  constraint_matrix.add_entry(dofs_2[i], dofs_1[j], is_identity_constrained?1.0:-1.0);
2015  }
2016  }
2017  else
2018  {
2019  // this is just a regular constraint. enter it piece by piece
2020  constraint_matrix.add_line(dofs_2[i]);
2021  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
2022  {
2023  // Query the correct face_index on face_1 respecting the given
2024  // orientation:
2025  const unsigned int j =
2026  cell_to_rotated_face_index[fe.face_to_cell_index
2027  (jj, 0, face_orientation, face_flip, face_rotation)];
2028 
2029  // And finally constrain the two DoFs respecting component_mask:
2030  if (transformation(i,jj) != 0)
2031  constraint_matrix.add_entry(dofs_2[i], dofs_1[j],
2032  transformation(i,jj));
2033  }
2034  }
2035  }
2036  }
2037  }
2038  }
2039 
2040 
2041  // Internally used in make_periodicity_constraints.
2042  //
2043  // Build up a (possibly rotated) interpolation matrix that is used in
2044  // set_periodicity_constraints with the help of user supplied matrix
2045  // and first_vector_components.
2046  template <int dim, int spacedim>
2047  FullMatrix<double> compute_transformation(
2048  const FiniteElement<dim, spacedim> &fe,
2049  const FullMatrix<double> &matrix,
2050  const std::vector<unsigned int> &first_vector_components)
2051  {
2052  Assert(matrix.m() == matrix.n(), ExcInternalError());
2053 
2054  const unsigned int n_dofs_per_face = fe.dofs_per_face;
2055 
2056  if (matrix.m() == n_dofs_per_face)
2057  {
2058  // In case of m == n == n_dofs_per_face the supplied matrix is already
2059  // an interpolation matrix, so we use it directly:
2060  return matrix;
2061  }
2062 
2063  if (first_vector_components.empty() && matrix.m() == 0)
2064  {
2065  // Just the identity matrix in case no rotation is specified:
2066  return IdentityMatrix(n_dofs_per_face);
2067  }
2068 
2069  // The matrix describes a rotation and we have to build a
2070  // transformation matrix, we assume that for a 0* rotation
2071  // we would have to build the identity matrix
2072 
2073  Assert(matrix.m() == (int)spacedim, ExcInternalError())
2074 
2076 
2077  // have an array that stores the location of each vector-dof tuple
2078  // we want to rotate.
2079  typedef std::array<unsigned int, spacedim> DoFTuple;
2080 
2081  // start with a pristine interpolation matrix...
2082  FullMatrix<double> transformation = IdentityMatrix(n_dofs_per_face);
2083 
2084  for (unsigned int i=0; i < n_dofs_per_face; ++i)
2085  {
2086  std::vector<unsigned int>::const_iterator comp_it
2087  = std::find (first_vector_components.begin(),
2088  first_vector_components.end(),
2089  fe.face_system_to_component_index(i).first);
2090  if (comp_it != first_vector_components.end())
2091  {
2092  const unsigned int first_vector_component = *comp_it;
2093 
2094  // find corresponding other components of vector
2095  DoFTuple vector_dofs;
2096  vector_dofs[0] = i;
2097  unsigned int n_found = 1;
2098 
2099  Assert(*comp_it + spacedim <= fe.n_components(),
2100  ExcMessage("Error: the finite element does not have enough components "
2101  "to define rotated periodic boundaries."));
2102 
2103  for (unsigned int k = 0; k < n_dofs_per_face; ++k)
2104  if ((k != i) &&
2105  (quadrature.point(k) == quadrature.point(i)) &&
2106  (fe.face_system_to_component_index(k).first >=
2107  first_vector_component) &&
2108  (fe.face_system_to_component_index(k).first <
2109  first_vector_component + spacedim))
2110  {
2111  vector_dofs[fe.face_system_to_component_index(k).first -
2112  first_vector_component]
2113  = k;
2114  n_found++;
2115  if (n_found==dim)
2116  break;
2117  }
2118 
2119  // ... and rotate all dofs belonging to vector valued
2120  // components that are selected by first_vector_components:
2121  for (int i = 0; i < spacedim; ++i)
2122  {
2123  transformation[vector_dofs[i]][vector_dofs[i]] = 0.;
2124  for (int j = 0; j < spacedim; ++j)
2125  transformation[vector_dofs[i]][vector_dofs[j]] = matrix[i][j];
2126  }
2127  }
2128  }
2129  return transformation;
2130  }
2131 
2132  } /*namespace*/
2133 
2134 
2135  // Low level interface:
2136 
2137 
2138  template <typename FaceIterator>
2139  void
2140  make_periodicity_constraints (const FaceIterator &face_1,
2141  const typename identity<FaceIterator>::type &face_2,
2142  ::ConstraintMatrix &constraint_matrix,
2143  const ComponentMask &component_mask,
2144  const bool face_orientation,
2145  const bool face_flip,
2146  const bool face_rotation,
2147  const FullMatrix<double> &matrix,
2148  const std::vector<unsigned int> &first_vector_components)
2149  {
2150  static const int dim = FaceIterator::AccessorType::dimension;
2151  static const int spacedim = FaceIterator::AccessorType::space_dimension;
2152 
2153  Assert( (dim != 1) ||
2154  (face_orientation == true &&
2155  face_flip == false &&
2156  face_rotation == false),
2157  ExcMessage ("The supplied orientation "
2158  "(face_orientation, face_flip, face_rotation) "
2159  "is invalid for 1D"));
2160 
2161  Assert( (dim != 2) ||
2162  (face_orientation == true &&
2163  face_rotation == false),
2164  ExcMessage ("The supplied orientation "
2165  "(face_orientation, face_flip, face_rotation) "
2166  "is invalid for 2D"));
2167 
2168  Assert(face_1 != face_2,
2169  ExcMessage ("face_1 and face_2 are equal! Cannot constrain DoFs "
2170  "on the very same face"));
2171 
2172  Assert(face_1->at_boundary() && face_2->at_boundary(),
2173  ExcMessage ("Faces for periodicity constraints must be on the "
2174  "boundary"));
2175 
2176  Assert(matrix.m() == matrix.n(),
2177  ExcMessage("The supplied (rotation or interpolation) matrix must "
2178  "be a square matrix"));
2179 
2180  Assert(first_vector_components.empty() || matrix.m() == (int)spacedim,
2181  ExcMessage ("first_vector_components is nonempty, so matrix must "
2182  "be a rotation matrix exactly of size spacedim"));
2183 
2184 #ifdef DEBUG
2185  if (!face_1->has_children())
2186  {
2187  Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
2188  const unsigned int n_dofs_per_face =
2189  face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face;
2190 
2191  Assert(matrix.m() == 0
2192  || (first_vector_components.empty() && matrix.m() == n_dofs_per_face)
2193  || (!first_vector_components.empty() && matrix.m() == (int)spacedim),
2194  ExcMessage ("The matrix must have either size 0 or spacedim "
2195  "(if first_vector_components is nonempty) "
2196  "or the size must be equal to the # of DoFs on the face "
2197  "(if first_vector_components is empty)."));
2198  }
2199 
2200  if (!face_2->has_children())
2201  {
2202  Assert(face_2->n_active_fe_indices() == 1, ExcInternalError());
2203  const unsigned int n_dofs_per_face =
2204  face_2->get_fe(face_2->nth_active_fe_index(0)).dofs_per_face;
2205 
2206  Assert(matrix.m() == 0
2207  || (first_vector_components.empty() && matrix.m() == n_dofs_per_face)
2208  || (!first_vector_components.empty() && matrix.m() == (int)spacedim),
2209  ExcMessage ("The matrix must have either size 0 or spacedim "
2210  "(if first_vector_components is nonempty) "
2211  "or the size must be equal to the # of DoFs on the face "
2212  "(if first_vector_components is empty)."));
2213  }
2214 #endif
2215 
2216  // A lookup table on how to go through the child faces depending on the
2217  // orientation:
2218 
2219  static const int lookup_table_2d[2][2] =
2220  {
2221  // flip:
2222  {0, 1}, // false
2223  {1, 0}, // true
2224  };
2225 
2226  static const int lookup_table_3d[2][2][2][4] =
2227  {
2228  // orientation flip rotation
2229  { { {0, 2, 1, 3}, // false false false
2230  {2, 3, 0, 1}, // false false true
2231  },
2232  { {3, 1, 2, 0}, // false true false
2233  {1, 0, 3, 2}, // false true true
2234  },
2235  },
2236  { { {0, 1, 2, 3}, // true false false
2237  {1, 3, 0, 2}, // true false true
2238  },
2239  { {3, 2, 1, 0}, // true true false
2240  {2, 0, 3, 1}, // true true true
2241  },
2242  },
2243  };
2244 
2245  if (face_1->has_children() && face_2->has_children())
2246  {
2247  // In the case that both faces have children, we loop over all
2248  // children and apply make_periodicty_constrains recursively:
2249 
2250  Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
2251  face_2->n_children() ==
2253  ExcNotImplemented());
2254 
2255  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face;
2256  ++i)
2257  {
2258  // Lookup the index for the second face
2259  unsigned int j;
2260  switch (dim)
2261  {
2262  case 2:
2263  j = lookup_table_2d[face_flip][i];
2264  break;
2265  case 3:
2266  j = lookup_table_3d[face_orientation][face_flip][face_rotation][i];
2267  break;
2268  default:
2269  AssertThrow(false, ExcNotImplemented());
2270  }
2271 
2272  make_periodicity_constraints (face_1->child(i),
2273  face_2->child(j),
2274  constraint_matrix,
2275  component_mask,
2276  face_orientation,
2277  face_flip,
2278  face_rotation,
2279  matrix,
2280  first_vector_components);
2281  }
2282  }
2283  else
2284  {
2285  // Otherwise at least one of the two faces is active and
2286  // we need to do some work and enter the constraints!
2287 
2288  // The finite element that matters is the one on the active face:
2289  const FiniteElement<dim,spacedim> &fe =
2290  face_1->has_children()
2291  ? face_2->get_fe(face_2->nth_active_fe_index(0))
2292  : face_1->get_fe(face_1->nth_active_fe_index(0));
2293 
2294  const unsigned int n_dofs_per_face = fe.dofs_per_face;
2295 
2296  // Sometimes we just have nothing to do (for all finite elements,
2297  // or systems which accidentally don't have any dofs on the
2298  // boundary).
2299  if (n_dofs_per_face == 0)
2300  return;
2301 
2302  const FullMatrix<double> transformation =
2303  compute_transformation(fe, matrix, first_vector_components);
2304 
2305  if (! face_2->has_children())
2306  {
2307  // Performance hack: We do not need to compute an inverse if
2308  // the matrix is the identity matrix.
2309  if (first_vector_components.empty() && matrix.m() == 0)
2310  {
2311  set_periodicity_constraints(face_2,
2312  face_1,
2313  transformation,
2314  constraint_matrix,
2315  component_mask,
2316  face_orientation,
2317  face_flip,
2318  face_rotation);
2319  }
2320  else
2321  {
2322  FullMatrix<double> inverse(transformation.m());
2323  inverse.invert(transformation);
2324 
2325  set_periodicity_constraints(face_2,
2326  face_1,
2327  inverse,
2328  constraint_matrix,
2329  component_mask,
2330  face_orientation,
2331  face_flip,
2332  face_rotation);
2333  }
2334  }
2335  else
2336  {
2337  Assert(!face_1->has_children(), ExcInternalError());
2338 
2339  // Important note:
2340  // In 3D we have to take care of the fact that face_rotation
2341  // gives the relative rotation of face_1 to face_2, i.e. we
2342  // have to invert the rotation when constraining face_2 to
2343  // face_1. Therefore face_flip has to be toggled if
2344  // face_rotation is true:
2345  // In case of inverted orientation, nothing has to be done.
2346  set_periodicity_constraints(face_1,
2347  face_2,
2348  transformation,
2349  constraint_matrix,
2350  component_mask,
2351  face_orientation,
2352  face_orientation
2353  ? face_rotation ^ face_flip
2354  : face_flip,
2355  face_rotation);
2356  }
2357  }
2358  }
2359 
2360 
2361 
2362  template <typename DoFHandlerType>
2363  void
2366  &periodic_faces,
2367  ::ConstraintMatrix &constraint_matrix,
2368  const ComponentMask &component_mask,
2369  const std::vector<unsigned int> &first_vector_components)
2370  {
2371  typedef std::vector<GridTools::PeriodicFacePair<typename DoFHandlerType::cell_iterator> >
2372  FaceVector;
2373  typename FaceVector::const_iterator it, end_periodic;
2374  it = periodic_faces.begin();
2375  end_periodic = periodic_faces.end();
2376 
2377  // Loop over all periodic faces...
2378  for (; it!=end_periodic; ++it)
2379  {
2380  typedef typename DoFHandlerType::face_iterator FaceIterator;
2381  const FaceIterator face_1 = it->cell[0]->face(it->face_idx[0]);
2382  const FaceIterator face_2 = it->cell[1]->face(it->face_idx[1]);
2383 
2384  Assert(face_1->at_boundary() && face_2->at_boundary(),
2385  ExcInternalError());
2386 
2387  Assert (face_1 != face_2,
2388  ExcInternalError());
2389 
2390  // ... and apply the low level make_periodicity_constraints function to
2391  // every matching pair:
2393  face_2,
2394  constraint_matrix,
2395  component_mask,
2396  it->orientation[0],
2397  it->orientation[1],
2398  it->orientation[2],
2399  it->matrix,
2400  first_vector_components);
2401  }
2402  }
2403 
2404 
2405  // High level interface variants:
2406 
2407 
2408  template <typename DoFHandlerType>
2409  void
2410  make_periodicity_constraints (const DoFHandlerType &dof_handler,
2411  const types::boundary_id b_id1,
2412  const types::boundary_id b_id2,
2413  const int direction,
2414  ::ConstraintMatrix &constraint_matrix,
2415  const ComponentMask &component_mask)
2416  {
2417  static const int space_dim = DoFHandlerType::space_dimension;
2418  (void)space_dim;
2419  Assert (0<=direction && direction<space_dim,
2420  ExcIndexRange (direction, 0, space_dim));
2421 
2422  Assert (b_id1 != b_id2,
2423  ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
2424  "different to denote different boundaries."));
2425 
2426  std::vector<GridTools::PeriodicFacePair
2427  <typename DoFHandlerType::cell_iterator> > matched_faces;
2428 
2429  // Collect matching periodic cells on the coarsest level:
2430  GridTools::collect_periodic_faces(dof_handler, b_id1, b_id2, direction,
2431  matched_faces);
2432 
2433  make_periodicity_constraints<DoFHandlerType>
2434  (matched_faces, constraint_matrix, component_mask);
2435  }
2436 
2437 
2438 
2439  template <typename DoFHandlerType>
2440  void
2441  make_periodicity_constraints (const DoFHandlerType &dof_handler,
2442  const types::boundary_id b_id,
2443  const int direction,
2444  ::ConstraintMatrix &constraint_matrix,
2445  const ComponentMask &component_mask)
2446  {
2447  static const int dim = DoFHandlerType::dimension;
2448  static const int space_dim = DoFHandlerType::space_dimension;
2449  (void)dim;
2450  (void)space_dim;
2451 
2452  Assert (0<=direction && direction<space_dim,
2453  ExcIndexRange (direction, 0, space_dim));
2454 
2455  Assert(dim == space_dim,
2456  ExcNotImplemented());
2457 
2458  std::vector<GridTools::PeriodicFacePair
2459  <typename DoFHandlerType::cell_iterator> > matched_faces;
2460 
2461  // Collect matching periodic cells on the coarsest level:
2462  GridTools::collect_periodic_faces(dof_handler, b_id, direction,
2463  matched_faces);
2464 
2465  make_periodicity_constraints<DoFHandlerType>
2466  (matched_faces, constraint_matrix, component_mask);
2467  }
2468 
2469 
2470 
2471  namespace internal
2472  {
2473  namespace Assembler
2474  {
2475  struct Scratch {};
2476 
2477 
2478  template <int dim,int spacedim>
2479  struct CopyData
2480  {
2481  unsigned int dofs_per_cell;
2482  std::vector<types::global_dof_index> parameter_dof_indices;
2483 #ifdef DEAL_II_WITH_MPI
2484  std::vector<::LinearAlgebra::distributed::Vector<double> > global_parameter_representation;
2485 #else
2486  std::vector<::Vector<double> > global_parameter_representation;
2487 #endif
2488  };
2489  }
2490 
2491  namespace
2492  {
2498  template <int dim, int spacedim>
2499  void compute_intergrid_weights_3 (
2500  const typename ::DoFHandler<dim,spacedim>::active_cell_iterator &cell,
2501  const Assembler::Scratch &,
2502  Assembler::CopyData<dim,spacedim> &copy_data,
2503  const unsigned int coarse_component,
2504  const FiniteElement<dim,spacedim> &coarse_fe,
2505  const InterGridMap<::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
2506  const std::vector<::Vector<double> > &parameter_dofs)
2507  {
2508  // for each cell on the parameter grid: find out which degrees of
2509  // freedom on the fine grid correspond in which way to the degrees
2510  // of freedom on the parameter grid
2511  //
2512  // since for continuous FEs some dofs exist on more than one cell,
2513  // we have to track which ones were already visited. the problem is
2514  // that if we visit a dof first on one cell and compute its weight
2515  // with respect to some global dofs to be non-zero, and later visit
2516  // the dof again on another cell and (since we are on another cell)
2517  // recompute the weights with respect to the same dofs as above to
2518  // be zero now, we have to preserve them. we therefore overwrite
2519  // all weights if they are nonzero and do not enforce zero weights
2520  // since that might be only due to the fact that we are on another
2521  // cell.
2522  //
2523  // example:
2524  // coarse grid
2525  // | | |
2526  // *-----*-----*
2527  // | cell|cell |
2528  // | 1 | 2 |
2529  // | | |
2530  // 0-----1-----*
2531  //
2532  // fine grid
2533  // | | | | |
2534  // *--*--*--*--*
2535  // | | | | |
2536  // *--*--*--*--*
2537  // | | | | |
2538  // *--x--y--*--*
2539  //
2540  // when on cell 1, we compute the weights of dof 'x' to be 1/2 from
2541  // parameter dofs 0 and 1, respectively. however, when later we are
2542  // on cell 2, we again compute the prolongation of shape function 1
2543  // restricted to cell 2 to the globla grid and find that the weight
2544  // of global dof 'x' now is zero. however, we should not overwrite
2545  // the old value.
2546  //
2547  // we therefore always only set nonzero values. why adding up is
2548  // not useful: dof 'y' would get weight 1 from parameter dof 1 on
2549  // both cells 1 and 2, but the correct weight is nevertheless only
2550  // 1.
2551 
2552  // vector to hold the representation of a single degree of freedom
2553  // on the coarse grid (for the selected fe) on the fine grid
2554 
2555  copy_data.dofs_per_cell = coarse_fe.dofs_per_cell;
2556  copy_data.parameter_dof_indices.resize(copy_data.dofs_per_cell);
2557 
2558  // get the global indices of the parameter dofs on this
2559  // parameter grid cell
2560  cell->get_dof_indices (copy_data.parameter_dof_indices);
2561 
2562  // loop over all dofs on this cell and check whether they are
2563  // interesting for us
2564  for (unsigned int local_dof=0; local_dof<copy_data.dofs_per_cell; ++local_dof)
2565  if (coarse_fe.system_to_component_index(local_dof).first
2566  ==
2567  coarse_component)
2568  {
2569  // the how-many-th parameter is this on this cell?
2570  const unsigned int local_parameter_dof
2571  = coarse_fe.system_to_component_index(local_dof).second;
2572 
2573  copy_data.global_parameter_representation[local_parameter_dof] = 0.;
2574 
2575  // distribute the representation of
2576  // @p{local_parameter_dof} on the parameter grid cell
2577  // @p{cell} to the global data space
2578  coarse_to_fine_grid_map[cell]->
2579  set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof],
2580  copy_data.global_parameter_representation[local_parameter_dof]);
2581  }
2582  }
2583 
2584 
2585 
2591  template <int dim,int spacedim>
2592  void copy_intergrid_weights_3(const Assembler::CopyData<dim,spacedim> &copy_data,
2593  const unsigned int coarse_component,
2594  const FiniteElement<dim,spacedim> &coarse_fe,
2595  const std::vector<types::global_dof_index> &weight_mapping,
2596  const bool is_called_in_parallel,
2597  std::vector<std::map<types::global_dof_index, float> > &weights)
2598  {
2599  unsigned int pos = 0;
2600  for (unsigned int local_dof=0; local_dof<copy_data.dofs_per_cell; ++local_dof)
2601  if (coarse_fe.system_to_component_index(local_dof).first
2602  ==
2603  coarse_component)
2604  {
2605  // now that we've got the global representation of each
2606  // parameter dof, we've only got to clobber the non-zero
2607  // entries in that vector and store the result
2608  //
2609  // what we have learned: if entry @p{i} of the global
2610  // vector holds the value @p{v[i]}, then this is the
2611  // weight with which the present dof contributes to
2612  // @p{i}. there may be several such @p{i}s and their
2613  // weights' sum should be one. Then, @p{v[i]} should be
2614  // equal to @p{\sum_j w_{ij} p[j]} with @p{p[j]} be the
2615  // values of the degrees of freedom on the coarse grid.
2616  // we can thus compute constraints which link the degrees
2617  // of freedom @p{v[i]} on the fine grid to those on the
2618  // coarse grid, @p{p[j]}. Now to use these as real
2619  // constraints, rather than as additional equations, we
2620  // have to identify representants among the @p{i} for
2621  // each @p{j}. this will be done by simply taking the
2622  // first @p{i} for which @p{w_{ij}==1}.
2623  //
2624  // guard modification of the weights array by a Mutex.
2625  // since it should happen rather rarely that there are
2626  // several threads operating on different intergrid
2627  // weights, have only one mutex for all of them
2628  for (types::global_dof_index i=0; i<copy_data.global_parameter_representation[pos].size();
2629  ++i)
2630  // set this weight if it belongs to a parameter dof.
2631  if (weight_mapping[i] != numbers::invalid_dof_index)
2632  {
2633  // only overwrite old value if not by zero
2634  if (copy_data.global_parameter_representation[pos](i) != 0)
2635  {
2636  const types::global_dof_index wi = copy_data.parameter_dof_indices[local_dof],
2637  wj = weight_mapping[i];
2638  weights[wi][wj] = copy_data.global_parameter_representation[pos](i);
2639  }
2640  }
2641  else if (!is_called_in_parallel)
2642  {
2643  // Note that when this function operates with distributed
2644  // fine grid, this assertion is switched off since the
2645  // condition does not necessarily hold
2646  Assert (copy_data.global_parameter_representation[pos](i) == 0,
2647  ExcInternalError());
2648  }
2649 
2650  ++pos;
2651  }
2652 
2653  }
2654 
2655 
2656 
2662  template <int dim, int spacedim>
2663  void
2664  compute_intergrid_weights_2 (
2665  const ::DoFHandler<dim,spacedim> &coarse_grid,
2666  const unsigned int coarse_component,
2667  const InterGridMap<::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
2668  const std::vector<::Vector<double> > &parameter_dofs,
2669  const std::vector<types::global_dof_index> &weight_mapping,
2670  std::vector<std::map<types::global_dof_index,float> > &weights)
2671  {
2672  Assembler::Scratch scratch;
2673  Assembler::CopyData<dim,spacedim> copy_data;
2674 
2675  unsigned int n_interesting_dofs = 0;
2676  for (unsigned int local_dof=0; local_dof<coarse_grid.get_fe().dofs_per_cell; ++local_dof)
2677  if (coarse_grid.get_fe().system_to_component_index(local_dof).first
2678  ==
2679  coarse_component)
2680  ++n_interesting_dofs;
2681 
2682  copy_data.global_parameter_representation.resize(n_interesting_dofs);
2683 
2684  bool is_called_in_parallel = false;
2685  for (size_t i=0; i<copy_data.global_parameter_representation.size(); ++i)
2686  {
2687 #ifdef DEAL_II_WITH_MPI
2688  MPI_Comm communicator = MPI_COMM_SELF;
2689  try
2690  {
2691  const typename ::parallel::Triangulation<dim, spacedim> &tria =
2692  dynamic_cast<const typename ::parallel::Triangulation<dim, spacedim>&>
2693  (coarse_to_fine_grid_map.get_destination_grid().get_triangulation());
2694  communicator = tria.get_communicator ();
2695  is_called_in_parallel = true;
2696  }
2697  catch (std::bad_cast &exp)
2698  {
2699  // Nothing bad happened: the user used serial Triangulation
2700  }
2701 
2702 
2703  IndexSet locally_relevant_dofs;
2705  (coarse_to_fine_grid_map.get_destination_grid (), locally_relevant_dofs);
2706 
2707  copy_data.global_parameter_representation[i].reinit
2708  (coarse_to_fine_grid_map.get_destination_grid().locally_owned_dofs(),
2709  locally_relevant_dofs, communicator);
2710 #else
2711  const types::global_dof_index n_fine_dofs = weight_mapping.size();
2712  copy_data.global_parameter_representation[i].reinit (n_fine_dofs);
2713 #endif
2714  }
2715 
2716  WorkStream::run(coarse_grid.begin_active(),
2717  coarse_grid.end(),
2718  std::bind(&compute_intergrid_weights_3<dim,spacedim>,
2719  std::placeholders::_1,
2720  std::placeholders::_2,
2721  std::placeholders::_3,
2722  coarse_component,
2723  std::cref(coarse_grid.get_fe()),
2724  std::cref(coarse_to_fine_grid_map),
2725  std::cref(parameter_dofs)),
2726  std::bind(&copy_intergrid_weights_3<dim,spacedim>,
2727  std::placeholders::_1,
2728  coarse_component,
2729  std::cref(coarse_grid.get_fe()),
2730  std::cref(weight_mapping),
2731  is_called_in_parallel,
2732  std::ref(weights)),
2733  scratch,
2734  copy_data);
2735 
2736 #ifdef DEAL_II_WITH_MPI
2737  for (size_t i=0; i<copy_data.global_parameter_representation.size(); ++i)
2738  copy_data.global_parameter_representation[i].update_ghost_values ();
2739 #endif
2740  }
2741 
2742 
2743 
2749  template <int dim, int spacedim>
2750  unsigned int
2751  compute_intergrid_weights_1 (
2752  const ::DoFHandler<dim,spacedim> &coarse_grid,
2753  const unsigned int coarse_component,
2754  const ::DoFHandler<dim,spacedim> &fine_grid,
2755  const unsigned int fine_component,
2756  const InterGridMap<::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
2757  std::vector<std::map<types::global_dof_index, float> > &weights,
2758  std::vector<types::global_dof_index> &weight_mapping)
2759  {
2760  // aliases to the finite elements used by the dof handlers:
2761  const FiniteElement<dim,spacedim> &coarse_fe = coarse_grid.get_fe(),
2762  &fine_fe = fine_grid.get_fe();
2763 
2764  // global numbers of dofs
2765  const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(),
2766  n_fine_dofs = fine_grid.n_dofs();
2767 
2768  // local numbers of dofs
2769  const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell;
2770 
2771  // alias the number of dofs per cell belonging to the
2772  // coarse_component which is to be the restriction of the fine
2773  // grid:
2774  const unsigned int coarse_dofs_per_cell_component
2775  = coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell;
2776 
2777 
2778  // Try to find out whether the grids stem from the same coarse
2779  // grid. This is a rather crude test, but better than nothing
2780  Assert (coarse_grid.get_triangulation().n_cells(0) == fine_grid.get_triangulation().n_cells(0),
2781  ExcGridsDontMatch());
2782 
2783  // check whether the map correlates the right objects
2784  Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid,
2785  ExcGridsDontMatch ());
2786  Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid,
2787  ExcGridsDontMatch ());
2788 
2789 
2790  // check whether component numbers are valid
2791  AssertIndexRange (coarse_component,coarse_fe.n_components());
2792  AssertIndexRange (fine_component, fine_fe.n_components());
2793 
2794  // check whether respective finite elements are equal
2795  Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first)
2796  ==
2797  fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first),
2799 
2800 #ifdef DEBUG
2801  // if in debug mode, check whether the coarse grid is indeed
2802  // coarser everywhere than the fine grid
2803  for (typename ::DoFHandler<dim,spacedim>::active_cell_iterator
2804  cell=coarse_grid.begin_active();
2805  cell != coarse_grid.end(); ++cell)
2806  Assert (cell->level() <= coarse_to_fine_grid_map[cell]->level(),
2807  ExcGridNotCoarser());
2808 #endif
2809 
2810  /*
2811  * From here on: the term `parameter' refers to the selected
2812  * component on the coarse grid and its analogon on the fine grid.
2813  * The naming of variables containing this term is due to the fact
2814  * that `selected_component' is longer, but also due to the fact
2815  * that the code of this function was initially written for a
2816  * program where the component which we wanted to match between
2817  * grids was actually the `parameter' variable.
2818  *
2819  * Likewise, the terms `parameter grid' and `state grid' refer to
2820  * the coarse and fine grids, respectively.
2821  *
2822  * Changing the names of variables would in principle be a good
2823  * idea, but would not make things simpler and would be another
2824  * source of errors. If anyone feels like doing so: patches would
2825  * be welcome!
2826  */
2827 
2828 
2829 
2830  // set up vectors of cell-local data; each vector represents one
2831  // degree of freedom of the coarse-grid variable in the fine-grid
2832  // element
2833  std::vector<::Vector<double> >
2834  parameter_dofs (coarse_dofs_per_cell_component,
2835  ::Vector<double>(fine_dofs_per_cell));
2836  // for each coarse dof: find its position within the fine element
2837  // and set this value to one in the respective vector (all other
2838  // values are zero by construction)
2839  for (unsigned int local_coarse_dof=0;
2840  local_coarse_dof<coarse_dofs_per_cell_component;
2841  ++local_coarse_dof)
2842  for (unsigned int fine_dof=0; fine_dof<fine_fe.dofs_per_cell; ++fine_dof)
2843  if (fine_fe.system_to_component_index(fine_dof)
2844  ==
2845  std::make_pair (fine_component, local_coarse_dof))
2846  {
2847  parameter_dofs[local_coarse_dof](fine_dof) = 1.;
2848  break;
2849  };
2850 
2851 
2852  // find out how many DoFs there are on the grids belonging to the
2853  // components we want to match
2854  unsigned int n_parameters_on_fine_grid=0;
2855  if (true)
2856  {
2857  // have a flag for each dof on the fine grid and set it to true
2858  // if this is an interesting dof. finally count how many true's
2859  // there
2860  std::vector<bool> dof_is_interesting (fine_grid.n_dofs(), false);
2861  std::vector<types::global_dof_index> local_dof_indices (fine_fe.dofs_per_cell);
2862 
2863  for (typename ::DoFHandler<dim,spacedim>::active_cell_iterator
2864  cell=fine_grid.begin_active();
2865  cell!=fine_grid.end(); ++cell)
2866  if (cell->is_locally_owned ())
2867  {
2868  cell->get_dof_indices (local_dof_indices);
2869  for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
2870  if (fine_fe.system_to_component_index(i).first == fine_component)
2871  dof_is_interesting[local_dof_indices[i]] = true;
2872  };
2873 
2874  n_parameters_on_fine_grid = std::count (dof_is_interesting.begin(),
2875  dof_is_interesting.end(),
2876  true);
2877  };
2878 
2879 
2880  // set up the weights mapping
2881  weights.clear ();
2882  weights.resize (n_coarse_dofs);
2883 
2884  weight_mapping.clear ();
2885  weight_mapping.resize (n_fine_dofs, numbers::invalid_dof_index);
2886 
2887  if (true)
2888  {
2889  std::vector<types::global_dof_index> local_dof_indices(fine_fe.dofs_per_cell);
2890  unsigned int next_free_index=0;
2891  for (typename ::DoFHandler<dim,spacedim>::active_cell_iterator
2892  cell=fine_grid.begin_active();
2893  cell != fine_grid.end(); ++cell)
2894  if (cell->is_locally_owned ())
2895  {
2896  cell->get_dof_indices (local_dof_indices);
2897  for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
2898  // if this DoF is a parameter dof and has not yet been
2899  // numbered, then do so
2900  if ((fine_fe.system_to_component_index(i).first == fine_component) &&
2901  (weight_mapping[local_dof_indices[i]] == numbers::invalid_dof_index))
2902  {
2903  weight_mapping[local_dof_indices[i]] = next_free_index;
2904  ++next_free_index;
2905  };
2906  };
2907 
2908  Assert (next_free_index == n_parameters_on_fine_grid,
2909  ExcInternalError());
2910  };
2911 
2912 
2913  // for each cell on the parameter grid: find out which degrees of
2914  // freedom on the fine grid correspond in which way to the degrees
2915  // of freedom on the parameter grid
2916  //
2917  // do this in a separate function to allow for multithreading
2918  // there. see this function also if you want to read more
2919  // information on the algorithm used.
2920  compute_intergrid_weights_2 (coarse_grid, coarse_component,
2921  coarse_to_fine_grid_map, parameter_dofs,
2922  weight_mapping, weights);
2923 
2924 
2925  // ok, now we have all weights for each dof on the fine grid. if in
2926  // debug mode lets see if everything went smooth, i.e. each dof has
2927  // sum of weights one
2928  //
2929  // in other words this means that if the sum of all shape functions
2930  // on the parameter grid is one (which is always the case), then
2931  // the representation on the state grid should be as well (division
2932  // of unity)
2933  //
2934  // if the parameter grid has more than one component, then the
2935  // respective dofs of the other components have sum of weights
2936  // zero, of course. we do not explicitly ask which component a dof
2937  // belongs to, but this at least tests some errors
2938 #ifdef DEBUG
2939  for (unsigned int col=0; col<n_parameters_on_fine_grid; ++col)
2940  {
2941  double sum=0;
2942  for (types::global_dof_index row=0; row<n_coarse_dofs; ++row)
2943  if (weights[row].find(col) != weights[row].end())
2944  sum += weights[row][col];
2945  Assert ((std::fabs(sum-1) < 1.e-12) ||
2946  ((coarse_fe.n_components()>1) && (sum==0)), ExcInternalError());
2947  };
2948 #endif
2949 
2950 
2951  return n_parameters_on_fine_grid;
2952  }
2953 
2954 
2955  }
2956  }
2957 
2958 
2959 
2960  template <int dim, int spacedim>
2961  void
2963  const DoFHandler<dim,spacedim> &coarse_grid,
2964  const unsigned int coarse_component,
2965  const DoFHandler<dim,spacedim> &fine_grid,
2966  const unsigned int fine_component,
2967  const InterGridMap<DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
2968  ConstraintMatrix &constraints)
2969  {
2970  // store the weights with which a dof on the parameter grid contributes
2971  // to a dof on the fine grid. see the long doc below for more info
2972  //
2973  // allocate as many rows as there are parameter dofs on the coarse grid
2974  // and as many columns as there are parameter dofs on the fine grid.
2975  //
2976  // weight_mapping is used to map the global (fine grid) parameter dof
2977  // indices to the columns
2978  //
2979  // in the original implementation, the weights array was actually of
2980  // FullMatrix<double> type. this wasted huge amounts of memory, but was
2981  // fast. nonetheless, since the memory consumption was quadratic in the
2982  // number of degrees of freedom, this was not very practical, so we now
2983  // use a vector of rows of the matrix, and in each row a vector of
2984  // pairs (colnum,value). this seems like the best tradeoff between
2985  // memory and speed, as it is now linear in memory and still fast
2986  // enough.
2987  //
2988  // to save some memory and since the weights are usually (negative)
2989  // powers of 2, we choose the value type of the matrix to be @p{float}
2990  // rather than @p{double}.
2991  std::vector<std::map<types::global_dof_index, float> > weights;
2992 
2993  // this is this mapping. there is one entry for each dof on the fine
2994  // grid; if it is a parameter dof, then its value is the column in
2995  // weights for that parameter dof, if it is any other dof, then its
2996  // value is -1, indicating an error
2997  std::vector<types::global_dof_index> weight_mapping;
2998 
2999  const unsigned int n_parameters_on_fine_grid
3000  = internal::compute_intergrid_weights_1 (coarse_grid, coarse_component,
3001  fine_grid, fine_component,
3002  coarse_to_fine_grid_map,
3003  weights, weight_mapping);
3004  (void)n_parameters_on_fine_grid;
3005 
3006  // global numbers of dofs
3007  const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(),
3008  n_fine_dofs = fine_grid.n_dofs();
3009 
3010 
3011  // get an array in which we store which dof on the coarse grid is a
3012  // parameter and which is not
3013  std::vector<bool> coarse_dof_is_parameter (coarse_grid.n_dofs());
3014  if (true)
3015  {
3016  std::vector<bool> mask (coarse_grid.get_fe(0).n_components(),
3017  false);
3018  mask[coarse_component] = true;
3019  extract_dofs (coarse_grid, ComponentMask(mask), coarse_dof_is_parameter);
3020  }
3021 
3022  // now we know that the weights in each row constitute a constraint.
3023  // enter this into the constraints object
3024  //
3025  // first task: for each parameter dof on the parameter grid, find a
3026  // representant on the fine, global grid. this is possible since we use
3027  // conforming finite element. we take this representant to be the first
3028  // element in this row with weight identical to one. the representant
3029  // will become an unconstrained degree of freedom, while all others
3030  // will be constrained to this dof (and possibly others)
3031  std::vector<types::global_dof_index> representants(n_coarse_dofs, numbers::invalid_dof_index);
3032  for (types::global_dof_index parameter_dof=0; parameter_dof<n_coarse_dofs;
3033  ++parameter_dof)
3034  if (coarse_dof_is_parameter[parameter_dof] == true)
3035  {
3036  // if this is the line of a parameter dof on the coarse grid,
3037  // then it should have at least one dependent node on the fine
3038  // grid
3039  Assert (weights[parameter_dof].size() > 0, ExcInternalError());
3040 
3041  // find the column where the representant is mentioned
3042  std::map<types::global_dof_index,float>::const_iterator i = weights[parameter_dof].begin();
3043  for (; i!=weights[parameter_dof].end(); ++i)
3044  if (i->second == 1)
3045  break;
3046  Assert (i!=weights[parameter_dof].end(), ExcInternalError());
3047  const types::global_dof_index column = i->first;
3048 
3049  // now we know in which column of weights the representant is,
3050  // but we don't know its global index. get it using the inverse
3051  // operation of the weight_mapping
3052  types::global_dof_index global_dof=0;
3053  for (; global_dof<weight_mapping.size(); ++global_dof)
3054  if (weight_mapping[global_dof] == static_cast<types::global_dof_index>(column))
3055  break;
3056  Assert (global_dof < weight_mapping.size(), ExcInternalError());
3057 
3058  // now enter the representants global index into our list
3059  representants[parameter_dof] = global_dof;
3060  }
3061  else
3062  {
3063  // consistency check: if this is no parameter dof on the coarse
3064  // grid, then the respective row must be empty!
3065  Assert (weights[parameter_dof].size() == 0, ExcInternalError());
3066  };
3067 
3068 
3069 
3070  // note for people that want to optimize this function: the largest
3071  // part of the computing time is spent in the following, rather
3072  // innocent block of code. basically, it must be the
3073  // ConstraintMatrix::add_entry call which takes the bulk of the time,
3074  // but it is not known to the author how to make it faster...
3075  std::vector<std::pair<types::global_dof_index,double> > constraint_line;
3076  for (types::global_dof_index global_dof=0; global_dof<n_fine_dofs; ++global_dof)
3077  if (weight_mapping[global_dof] != numbers::invalid_dof_index)
3078  // this global dof is a parameter dof, so it may carry a constraint
3079  // note that for each global dof, the sum of weights shall be one,
3080  // so we can find out whether this dof is constrained in the
3081  // following way: if the only weight in this row is a one, and the
3082  // representant for the parameter dof of the line in which this one
3083  // is is the present dof, then we consider this dof to be
3084  // unconstrained. otherwise, all other dofs are constrained
3085  {
3086  const types::global_dof_index col = weight_mapping[global_dof];
3087  Assert (col < n_parameters_on_fine_grid, ExcInternalError());
3088 
3089  types::global_dof_index first_used_row=0;
3090 
3091  {
3092  Assert (weights.size() > 0, ExcInternalError());
3093  std::map<types::global_dof_index,float>::const_iterator
3094  col_entry = weights[0].end();
3095  for (; first_used_row<n_coarse_dofs; ++first_used_row)
3096  {
3097  col_entry = weights[first_used_row].find(col);
3098  if (col_entry != weights[first_used_row].end())
3099  break;
3100  }
3101 
3102  Assert (col_entry != weights[first_used_row].end(), ExcInternalError());
3103 
3104  if ((col_entry->second == 1) &&
3105  (representants[first_used_row] == global_dof))
3106  // dof unconstrained or constrained to itself (in case this
3107  // cell is mapped to itself, rather than to children of
3108  // itself)
3109  continue;
3110  }
3111 
3112 
3113  // otherwise enter all constraints
3114  constraints.add_line (global_dof);
3115 
3116  constraint_line.clear ();
3117  for (types::global_dof_index row=first_used_row; row<n_coarse_dofs; ++row)
3118  {
3119  const std::map<types::global_dof_index,float>::const_iterator
3120  j = weights[row].find(col);
3121  if ((j != weights[row].end()) && (j->second != 0))
3122  constraint_line.emplace_back (representants[row], j->second);
3123  };
3124 
3125  constraints.add_entries (global_dof, constraint_line);
3126  };
3127  }
3128 
3129 
3130 
3131  template <int dim, int spacedim>
3132  void
3134  const DoFHandler<dim,spacedim> &coarse_grid,
3135  const unsigned int coarse_component,
3136  const DoFHandler<dim,spacedim> &fine_grid,
3137  const unsigned int fine_component,
3138  const InterGridMap<DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
3139  std::vector<std::map<types::global_dof_index, float> > &transfer_representation)
3140  {
3141  // store the weights with which a dof on the parameter grid contributes
3142  // to a dof on the fine grid. see the long doc below for more info
3143  //
3144  // allocate as many rows as there are parameter dofs on the coarse grid
3145  // and as many columns as there are parameter dofs on the fine grid.
3146  //
3147  // weight_mapping is used to map the global (fine grid) parameter dof
3148  // indices to the columns
3149  //
3150  // in the original implementation, the weights array was actually of
3151  // FullMatrix<double> type. this wasted huge amounts of memory, but was
3152  // fast. nonetheless, since the memory consumption was quadratic in the
3153  // number of degrees of freedom, this was not very practical, so we now
3154  // use a vector of rows of the matrix, and in each row a vector of
3155  // pairs (colnum,value). this seems like the best tradeoff between
3156  // memory and speed, as it is now linear in memory and still fast
3157  // enough.
3158  //
3159  // to save some memory and since the weights are usually (negative)
3160  // powers of 2, we choose the value type of the matrix to be @p{float}
3161  // rather than @p{double}.
3162  std::vector<std::map<types::global_dof_index, float> > weights;
3163 
3164  // this is this mapping. there is one entry for each dof on the fine
3165  // grid; if it is a parameter dof, then its value is the column in
3166  // weights for that parameter dof, if it is any other dof, then its
3167  // value is -1, indicating an error
3168  std::vector<types::global_dof_index> weight_mapping;
3169 
3170  internal::compute_intergrid_weights_1 (coarse_grid, coarse_component,
3171  fine_grid, fine_component,
3172  coarse_to_fine_grid_map,
3173  weights, weight_mapping);
3174 
3175  // now compute the requested representation
3176  const types::global_dof_index n_global_parm_dofs
3177  = std::count_if (weight_mapping.begin(), weight_mapping.end(),
3178  std::bind (std::not_equal_to<types::global_dof_index>(), std::placeholders::_1, numbers::invalid_dof_index));
3179 
3180  // first construct the inverse mapping of weight_mapping
3181  std::vector<types::global_dof_index> inverse_weight_mapping (n_global_parm_dofs,
3183  for (types::global_dof_index i=0; i<weight_mapping.size(); ++i)
3184  {
3185  const types::global_dof_index parameter_dof = weight_mapping[i];
3186  // if this global dof is a parameter
3187  if (parameter_dof != numbers::invalid_dof_index)
3188  {
3189  Assert (parameter_dof < n_global_parm_dofs, ExcInternalError());
3190  Assert ((inverse_weight_mapping[parameter_dof] == numbers::invalid_dof_index),
3191  ExcInternalError());
3192 
3193  inverse_weight_mapping[parameter_dof] = i;
3194  };
3195  };
3196 
3197  // next copy over weights array and replace respective numbers
3198  const types::global_dof_index n_rows = weight_mapping.size();
3199 
3200  transfer_representation.clear ();
3201  transfer_representation.resize (n_rows);
3202 
3203  const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs();
3204  for (types::global_dof_index i=0; i<n_coarse_dofs; ++i)
3205  {
3206  std::map<types::global_dof_index, float>::const_iterator j = weights[i].begin();
3207  for (; j!=weights[i].end(); ++j)
3208  {
3209  const types::global_dof_index p = inverse_weight_mapping[j->first];
3210  Assert (p<n_rows, ExcInternalError());
3211 
3212  transfer_representation[p][i] = j->second;
3213  };
3214  };
3215  }
3216 
3217 
3218 
3219  template <int dim, int spacedim, template <int, int> class DoFHandlerType>
3220  void
3221  make_zero_boundary_constraints (const DoFHandlerType<dim, spacedim> &dof,
3222  const types::boundary_id boundary_id,
3223  ConstraintMatrix &zero_boundary_constraints,
3224  const ComponentMask &component_mask)
3225  {
3226  Assert (component_mask.represents_n_components(dof.get_fe(0).n_components()),
3227  ExcMessage ("The number of components in the mask has to be either "
3228  "zero or equal to the number of components in the finite "
3229  "element."));
3230 
3231  const unsigned int n_components = DoFTools::n_components (dof);
3232 
3233  Assert (component_mask.n_selected_components(n_components) > 0,
3235 
3236  // a field to store the indices on the face
3237  std::vector<types::global_dof_index> face_dofs;
3238  face_dofs.reserve (max_dofs_per_face(dof));
3239  // a field to store the indices on the cell
3240  std::vector<types::global_dof_index> cell_dofs;
3241  cell_dofs.reserve (max_dofs_per_cell(dof));
3242 
3243  typename DoFHandlerType<dim,spacedim>::active_cell_iterator
3244  cell = dof.begin_active(),
3245  endc = dof.end();
3246  for (; cell!=endc; ++cell)
3247  if (!cell->is_artificial()
3248  &&
3249  cell->at_boundary ())
3250  {
3251  const FiniteElement<dim,spacedim> &fe = cell->get_fe();
3252 
3253  // get global indices of dofs on the cell
3254  cell_dofs.resize (fe.dofs_per_cell);
3255  cell->get_dof_indices (cell_dofs);
3256 
3257  for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
3258  ++face_no)
3259  {
3260  const typename DoFHandlerType<dim,spacedim>::face_iterator face = cell->face(face_no);
3261 
3262  // if face is on the boundary and satisfies the correct
3263  // boundary id property
3264  if (face->at_boundary ()
3265  &&
3266  ((boundary_id == numbers::invalid_boundary_id)
3267  ||
3268  (face->boundary_id() == boundary_id)))
3269  {
3270  // get indices and physical location on this face
3271  face_dofs.resize (fe.dofs_per_face);
3272  face->get_dof_indices (face_dofs, cell->active_fe_index());
3273 
3274  // enter those dofs into the list that match the component
3275  // signature.
3276  for (unsigned int i=0; i<face_dofs.size(); ++i)
3277  {
3278  // Find out if a dof has a contribution in this
3279  // component, and if so, add it to the list
3280  const std::vector<types::global_dof_index>::iterator it_index_on_cell
3281  = std::find (cell_dofs.begin(), cell_dofs.end(), face_dofs[i]);
3282  Assert (it_index_on_cell != cell_dofs.end(), ExcInvalidIterator());
3283  const unsigned int index_on_cell = std::distance(cell_dofs.begin(),
3284  it_index_on_cell);
3285  const ComponentMask &nonzero_component_array
3286  = cell->get_fe().get_nonzero_components (index_on_cell);
3287  bool nonzero = false;
3288  for (unsigned int c=0; c<n_components; ++c)
3289  if (nonzero_component_array[c] && component_mask[c])
3290  {
3291  nonzero = true;
3292  break;
3293  }
3294 
3295  if (nonzero)
3296  zero_boundary_constraints.add_line (face_dofs[i]);
3297  }
3298  }
3299  }
3300  }
3301  }
3302 
3303 
3304 
3305  template <int dim, int spacedim, template <int, int> class DoFHandlerType>
3306  void
3307  make_zero_boundary_constraints (const DoFHandlerType<dim, spacedim> &dof,
3308  ConstraintMatrix &zero_boundary_constraints,
3309  const ComponentMask &component_mask)
3310  {
3312  zero_boundary_constraints, component_mask);
3313  }
3314 
3315 
3316 } // end of namespace DoFTools
3317 
3318 
3319 
3320 // explicit instantiations
3321 
3322 #include "dof_tools_constraints.inst"
3323 
3324 
3325 
3326 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
Definition: types.h:173
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)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1183
static::ExceptionBase & ExcGridsDontMatch()
Contents is actually a matrix.
bool is_constrained(const size_type index) const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const std::vector< std::pair< size_type, double > > * get_constraint_entries(const size_type line) const
unsigned int boundary_id
Definition: types.h:110
#define AssertIndexRange(index, range)
Definition: exceptions.h:1207
const unsigned int dofs_per_quad
Definition: fe_base.h:250
void add_entries(const size_type line, const std::vector< std::pair< size_type, double > > &col_val_pairs)
static::ExceptionBase & ExcFiniteElementsDontMatch()
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:398
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const unsigned int dofs_per_line
Definition: fe_base.h:244
bool represents_n_components(const unsigned int n) const
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)
static::ExceptionBase & ExcInvalidIterator()
void add_line(const size_type line)
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, ConstraintMatrix &constraints)
static::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
static::ExceptionBase & ExcGridNotCoarser()
void add_entry(const size_type line, const size_type column, const double value)
void invert(const FullMatrix< number2 > &M)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:88
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:337
types::global_dof_index n_dofs() const
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const types::boundary_id invalid_boundary_id
Definition: types.h:204
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
Definition: fe.h:3124
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.cc:4542
static::ExceptionBase & ExcNoComponentSelected()
const FullMatrix< double > & constraints(const ::internal::SubfaceCase< dim > &subface_case=::internal::SubfaceCase< dim >::case_isotropic) const
Definition: fe.cc:793
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1080
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
unsigned int n_components() const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index) const
Definition: fe.h:3059
size_type m() const
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
Definition: fe.cc:862
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const
Definition: fe.cc:878
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 >())
const unsigned int dofs_per_cell
Definition: fe_base.h:295
void extract_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< bool > &selected_dofs)
Definition: dof_tools.cc:372
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3003
const unsigned int dofs_per_face
Definition: fe_base.h:288
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
static::ExceptionBase & ExcNotImplemented()
const unsigned int dofs_per_vertex
Definition: fe_base.h:238
Definition: table.h:33
const std::vector< Point< dim-1 > > & get_unit_face_support_points() const
Definition: fe.cc:1009
void set_inhomogeneity(const size_type line, const double value)
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())
const types::global_dof_index invalid_dof_index
Definition: types.h:187
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
Definition: fe.cc:1225
T max(const T &t, const MPI_Comm &mpi_communicator)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static::ExceptionBase & ExcInternalError()