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