Reference documentation for deal.II version Git 05e4468 2017-09-21 10:18:23 +0200
dof_tools_constraints.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2017 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/thread_management.h>
17 #include <deal.II/base/table.h>
18 #include <deal.II/base/template_constraints.h>
19 #include <deal.II/base/utilities.h>
20 #include <deal.II/base/work_stream.h>
21 #include <deal.II/lac/vector.h>
22 #include <deal.II/lac/constraint_matrix.h>
23 #include <deal.II/grid/tria.h>
24 #include <deal.II/grid/tria_iterator.h>
25 #include <deal.II/grid/intergrid_map.h>
26 #include <deal.II/grid/grid_tools.h>
27 #include <deal.II/dofs/dof_handler.h>
28 #include <deal.II/dofs/dof_accessor.h>
29 #include <deal.II/fe/fe.h>
30 #include <deal.II/fe/fe_values.h>
31 #include <deal.II/fe/fe_tools.h>
32 #include <deal.II/hp/fe_collection.h>
33 #include <deal.II/hp/fe_values.h>
34 #include <deal.II/dofs/dof_tools.h>
35 
36 #ifdef DEAL_II_WITH_MPI
37 #include <deal.II/lac/la_parallel_vector.h>
38 #endif
39 
40 #include <algorithm>
41 #include <array>
42 #include <numeric>
43 
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 
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::shared_ptr<std::vector<bool> > &master_dof_mask)
312  {
313  if (master_dof_mask == std::shared_ptr<std::vector<bool> >())
314  {
315  master_dof_mask = std::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::shared_ptr<FullMatrix<double> > &matrix)
337  {
338  if (matrix == std::shared_ptr<FullMatrix<double> >())
339  {
340  matrix = std::shared_ptr<FullMatrix<double> >
342  fe1.dofs_per_face));
344  *matrix);
345  }
346  }
347 
348 
349 
353  template <int dim, int spacedim>
354  void
355  ensure_existence_of_subface_matrix (const FiniteElement<dim,spacedim> &fe1,
356  const FiniteElement<dim,spacedim> &fe2,
357  const unsigned int subface,
358  std::shared_ptr<FullMatrix<double> > &matrix)
359  {
360  if (matrix == std::shared_ptr<FullMatrix<double> >())
361  {
362  matrix = std::shared_ptr<FullMatrix<double> >
364  fe1.dofs_per_face));
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::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::shared_ptr<std::pair<FullMatrix<double>,FullMatrix<double> > >())
390  {
391  split_matrix
392  = std::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_collection().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> &,
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> &,
568  std::integral_constant<int, 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  std::integral_constant<int, 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> &,
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> &,
610  std::integral_constant<int, 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> &,
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> &,
626  std::integral_constant<int, 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 // std::integral_constant<int, 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 // std::integral_constant<int, 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  std::integral_constant<int, 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,
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,
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  std::integral_constant<int, 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,
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,
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,
875  Assert (cell->face(face)->child(c)->line(e)
876  ->fe_index_is_active(cell->active_fe_index()) == true,
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,
885  Assert (cell->face(face)->line(e)
886  ->fe_index_is_active(cell->active_fe_index()) == true,
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 
1027 
1028  template <typename DoFHandlerType>
1029  void
1030  make_hp_hanging_node_constraints (const DoFHandlerType &dof_handler,
1031  ConstraintMatrix &constraints)
1032  {
1033  // note: this function is going to be hard to understand if you
1034  // haven't read the hp paper. however, we try to follow the notation
1035  // laid out there, so go read the paper before you try to understand
1036  // what is going on here
1037 
1038  const unsigned int dim = DoFHandlerType::dimension;
1039 
1040  const unsigned int spacedim = DoFHandlerType::space_dimension;
1041 
1042 
1043  // a matrix to be used for constraints below. declared here and
1044  // simply resized down below to avoid permanent re-allocation of
1045  // memory
1046  FullMatrix<double> constraint_matrix;
1047 
1048  // similarly have arrays that will hold master and slave dof numbers,
1049  // as well as a scratch array needed for the complicated case below
1050  std::vector<types::global_dof_index> master_dofs;
1051  std::vector<types::global_dof_index> slave_dofs;
1052  std::vector<types::global_dof_index> scratch_dofs;
1053 
1054  // caches for the face and subface interpolation matrices between
1055  // different (or the same) finite elements. we compute them only
1056  // once, namely the first time they are needed, and then just reuse
1057  // them
1059  face_interpolation_matrices (n_finite_elements (dof_handler),
1060  n_finite_elements (dof_handler));
1062  subface_interpolation_matrices (n_finite_elements (dof_handler),
1063  n_finite_elements (dof_handler),
1065 
1066  // similarly have a cache for the matrices that are split into their
1067  // master and slave parts, and for which the master part is inverted.
1068  // these two matrices are derived from the face interpolation matrix
1069  // as described in the @ref hp_paper "hp paper"
1071  split_face_interpolation_matrices (n_finite_elements (dof_handler),
1072  n_finite_elements (dof_handler));
1073 
1074  // finally, for each pair of finite elements, have a mask that states
1075  // which of the degrees of freedom on the coarse side of a refined
1076  // face will act as master dofs.
1078  master_dof_masks (n_finite_elements (dof_handler),
1079  n_finite_elements (dof_handler));
1080 
1081  // loop over all faces
1082  //
1083  // note that even though we may visit a face twice if the neighboring
1084  // cells are equally refined, we can only visit each face with
1085  // hanging nodes once
1086  typename DoFHandlerType::active_cell_iterator cell = dof_handler.begin_active(),
1087  endc = dof_handler.end();
1088  for (; cell!=endc; ++cell)
1089  {
1090  // artificial cells can at best neighbor ghost cells, but we're not
1091  // interested in these interfaces
1092  if (cell->is_artificial ())
1093  continue;
1094 
1095  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1096  if (cell->face(face)->has_children())
1097  {
1098  // first of all, make sure that we treat a case which is
1099  // possible, i.e. either no dofs on the face at all or no
1100  // anisotropic refinement
1101  if (cell->get_fe().dofs_per_face == 0)
1102  continue;
1103 
1104  Assert(cell->face(face)->refinement_case()==RefinementCase<dim-1>::isotropic_refinement,
1105  ExcNotImplemented());
1106 
1107  // so now we've found a face of an active cell that has
1108  // children. that means that there are hanging nodes here.
1109 
1110  // in any case, faces can have at most two sets of active
1111  // fe indices, but here the face can have only one (namely
1112  // the same as that from the cell we're sitting on), and
1113  // each of the children can have only one as well. check
1114  // this
1115  Assert (cell->face(face)->n_active_fe_indices() == 1,
1116  ExcInternalError());
1117  Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
1118  == true,
1119  ExcInternalError());
1120  for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
1121  Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
1122  ExcInternalError());
1123 
1124  // first find out whether we can constrain each of the
1125  // subfaces to the mother face. in the lingo of the hp
1126  // paper, this would be the simple case. note that we can
1127  // short-circuit this decision if the dof_handler doesn't
1128  // support hp at all
1129  //
1130  // ignore all interfaces with artificial cells
1133 
1134  // auxiliary variable which holds FE indices of the mother face
1135  // and its subfaces. This knowledge will be needed in hp-case
1136  // with neither_element_dominates.
1137  std::set<unsigned int> fe_ind_face_subface;
1138  fe_ind_face_subface.insert(cell->active_fe_index());
1139 
1140  if (DoFHandlerSupportsDifferentFEs<DoFHandlerType>::value == true)
1141  for (unsigned int c=0; c<cell->face(face)->number_of_children(); ++c)
1142  if (!cell->neighbor_child_on_subface (face, c)->is_artificial())
1143  {
1144  mother_face_dominates = mother_face_dominates &
1145  (cell->get_fe().compare_for_face_domination
1146  (cell->neighbor_child_on_subface (face, c)->get_fe()));
1147  fe_ind_face_subface.insert(cell->neighbor_child_on_subface (face, c)->active_fe_index());
1148  }
1149 
1150  switch (mother_face_dominates)
1151  {
1154  {
1155  // Case 1 (the simple case and the only case that can
1156  // happen for non-hp DoFHandlers): The coarse element
1157  // dominates the elements on the subfaces (or they are
1158  // all the same)
1159  //
1160  // so we are going to constrain the DoFs on the face
1161  // children against the DoFs on the face itself
1162  master_dofs.resize (cell->get_fe().dofs_per_face);
1163 
1164  cell->face(face)->get_dof_indices (master_dofs,
1165  cell->active_fe_index ());
1166 
1167  // Now create constraint matrix for the subfaces and
1168  // assemble it. ignore all interfaces with artificial
1169  // cells because we can only get to such interfaces if
1170  // the current cell is a ghost cell
1171  for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
1172  {
1173  if (cell->neighbor_child_on_subface (face, c)->is_artificial())
1174  continue;
1175 
1176  const typename DoFHandlerType::active_face_iterator
1177  subface = cell->face(face)->child(c);
1178 
1179  Assert (subface->n_active_fe_indices() == 1,
1180  ExcInternalError());
1181 
1182  const unsigned int
1183  subface_fe_index = subface->nth_active_fe_index(0);
1184 
1185  // we sometime run into the situation where for
1186  // example on one big cell we have a FE_Q(1) and on
1187  // the subfaces we have a mixture of FE_Q(1) and
1188  // FE_Nothing. In that case, the face domination is
1189  // either_element_can_dominate for the whole
1190  // collection of subfaces, but on the particular
1191  // subface between FE_Q(1) and FE_Nothing, there
1192  // are no constraints that we need to take care of.
1193  // in that case, just continue
1194  if (cell->get_fe().compare_for_face_domination
1195  (subface->get_fe(subface_fe_index))
1196  ==
1198  continue;
1199 
1200  // Same procedure as for the mother cell. Extract
1201  // the face DoFs from the cell DoFs.
1202  slave_dofs.resize (subface->get_fe(subface_fe_index)
1203  .dofs_per_face);
1204  subface->get_dof_indices (slave_dofs, subface_fe_index);
1205 
1206  for (unsigned int i=0; i<slave_dofs.size(); ++i)
1207  Assert (slave_dofs[i] != numbers::invalid_dof_index,
1208  ExcInternalError());
1209 
1210  // Now create the element constraint for this
1211  // subface.
1212  //
1213  // As a side remark, one may wonder the following:
1214  // neighbor_child is clearly computed correctly,
1215  // i.e. taking into account face_orientation (just
1216  // look at the implementation of that function).
1217  // however, we don't care about this here, when we
1218  // ask for subface_interpolation on subface c. the
1219  // question rather is: do we have to translate 'c'
1220  // here as well?
1221  //
1222  // the answer is in fact 'no'. if one does that,
1223  // results are wrong: constraints are added twice
1224  // for the same pair of nodes but with differing
1225  // weights. in addition, one can look at the
1226  // deal.II/project_*_03 tests that look at exactly
1227  // this case: there, we have a mesh with at least
1228  // one face_orientation==false and hanging nodes,
1229  // and the results of those tests show that the
1230  // result of projection verifies the approximation
1231  // properties of a finite element onto that mesh
1232  ensure_existence_of_subface_matrix
1233  (cell->get_fe(),
1234  subface->get_fe(subface_fe_index),
1235  c,
1236  subface_interpolation_matrices
1237  [cell->active_fe_index()][subface_fe_index][c]);
1238 
1239  // Add constraints to global constraint matrix.
1240  filter_constraints (master_dofs,
1241  slave_dofs,
1242  *(subface_interpolation_matrices
1243  [cell->active_fe_index()][subface_fe_index][c]),
1244  constraints);
1245  }
1246 
1247  break;
1248  }
1249 
1252  {
1253  // Case 2 (the "complex" case): at least one (the
1254  // neither_... case) of the finer elements or all of
1255  // them (the other_... case) is dominating. See the hp
1256  // paper for a way how to deal with this situation
1257  //
1258  // since this is something that can only happen for hp
1259  // dof handlers, add a check here...
1260  Assert (DoFHandlerSupportsDifferentFEs<DoFHandlerType>::value == true,
1261  ExcInternalError());
1262 
1263  const ::hp::FECollection<dim,spacedim> &fe_collection =
1264  dof_handler.get_fe_collection ();
1265  // we first have to find the finite element that is
1266  // able to generate a space that all the other ones can
1267  // be constrained to.
1268  // At this point we potentially have different scenarios:
1269  // 1) sub-faces dominate mother face and there is a
1270  // dominating FE among sub faces. We could loop over sub
1271  // faces to find the needed FE index. However, this will not
1272  // work in the case when
1273  // 2) there is no dominating FE among sub faces (e.g. Q1xQ2 vs Q2xQ1),
1274  // but subfaces still dominate mother face (e.g. Q2xQ2).
1275  // To cover this case we would have to use find_least_face_dominating_fe()
1276  // of FECollection with fe_indices of sub faces.
1277  // 3) Finally, it could happen that we got here because
1278  // neither_element_dominates (e.g. Q1xQ1xQ2 and Q1xQ2xQ1 for
1279  // subfaces and Q2xQ1xQ1 for mother face).
1280  // This requires usage of find_least_face_dominating_fe()
1281  // with fe_indices of sub-faces and the mother face.
1282  // Note that the last solution covers the first two scenarios,
1283  // thus we stick with it assuming that we won't lose much time/efficiency.
1284  const unsigned int dominating_fe_index = fe_collection.find_least_face_dominating_fe(fe_ind_face_subface);
1285  AssertThrow(dominating_fe_index != numbers::invalid_unsigned_int,
1286  ExcMessage("Could not find a least face dominating FE."));
1287 
1288  const FiniteElement<dim,spacedim> &dominating_fe
1289  = dof_handler.get_fe(dominating_fe_index);
1290 
1291  // first get the interpolation matrix from the mother
1292  // to the virtual dofs
1293  Assert (dominating_fe.dofs_per_face <=
1294  cell->get_fe().dofs_per_face,
1295  ExcInternalError());
1296 
1297  ensure_existence_of_face_matrix
1298  (dominating_fe,
1299  cell->get_fe(),
1300  face_interpolation_matrices
1301  [dominating_fe_index][cell->active_fe_index()]);
1302 
1303  // split this matrix into master and slave components.
1304  // invert the master component
1305  ensure_existence_of_master_dof_mask
1306  (cell->get_fe(),
1307  dominating_fe,
1308  (*face_interpolation_matrices
1309  [dominating_fe_index]
1310  [cell->active_fe_index()]),
1311  master_dof_masks
1312  [dominating_fe_index]
1313  [cell->active_fe_index()]);
1314 
1315  ensure_existence_of_split_face_matrix
1316  (*face_interpolation_matrices
1317  [dominating_fe_index][cell->active_fe_index()],
1318  (*master_dof_masks
1319  [dominating_fe_index][cell->active_fe_index()]),
1320  split_face_interpolation_matrices
1321  [dominating_fe_index][cell->active_fe_index()]);
1322 
1323  const FullMatrix<double> &restrict_mother_to_virtual_master_inv
1324  = (split_face_interpolation_matrices
1325  [dominating_fe_index][cell->active_fe_index()]->first);
1326 
1327  const FullMatrix<double> &restrict_mother_to_virtual_slave
1328  = (split_face_interpolation_matrices
1329  [dominating_fe_index][cell->active_fe_index()]->second);
1330 
1331  // now compute the constraint matrix as the product
1332  // between the inverse matrix and the slave part
1333  constraint_matrix.reinit (cell->get_fe().dofs_per_face -
1334  dominating_fe.dofs_per_face,
1335  dominating_fe.dofs_per_face);
1336  restrict_mother_to_virtual_slave
1337  .mmult (constraint_matrix,
1338  restrict_mother_to_virtual_master_inv);
1339 
1340  // then figure out the global numbers of master and
1341  // slave dofs and apply constraints
1342  scratch_dofs.resize (cell->get_fe().dofs_per_face);
1343  cell->face(face)->get_dof_indices (scratch_dofs,
1344  cell->active_fe_index ());
1345 
1346  // split dofs into master and slave components
1347  master_dofs.clear ();
1348  slave_dofs.clear ();
1349  for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
1350  if ((*master_dof_masks
1351  [dominating_fe_index][cell->active_fe_index()])[i] == true)
1352  master_dofs.push_back (scratch_dofs[i]);
1353  else
1354  slave_dofs.push_back (scratch_dofs[i]);
1355 
1356  AssertDimension (master_dofs.size(), dominating_fe.dofs_per_face);
1357  AssertDimension (slave_dofs.size(),
1358  cell->get_fe().dofs_per_face - dominating_fe.dofs_per_face);
1359 
1360  filter_constraints (master_dofs,
1361  slave_dofs,
1362  constraint_matrix,
1363  constraints);
1364 
1365 
1366 
1367  // next we have to deal with the subfaces. do as
1368  // discussed in the hp paper
1369  for (unsigned int sf=0;
1370  sf<cell->face(face)->n_children(); ++sf)
1371  {
1372  // ignore interfaces with artificial cells as well
1373  // as interfaces between ghost cells in 2d
1374  if (cell->neighbor_child_on_subface (face, sf)->is_artificial()
1375  ||
1376  (dim==2 && cell->is_ghost()
1377  &&
1378  cell->neighbor_child_on_subface (face, sf)->is_ghost()))
1379  continue;
1380 
1381  Assert (cell->face(face)->child(sf)
1382  ->n_active_fe_indices() == 1,
1383  ExcInternalError());
1384 
1385  const unsigned int subface_fe_index
1386  = cell->face(face)->child(sf)->nth_active_fe_index(0);
1387  const FiniteElement<dim,spacedim> &subface_fe
1388  = dof_handler.get_fe(subface_fe_index);
1389 
1390  // first get the interpolation matrix from the
1391  // subface to the virtual dofs
1392  Assert (dominating_fe.dofs_per_face <=
1393  subface_fe.dofs_per_face,
1394  ExcInternalError());
1395  ensure_existence_of_subface_matrix
1396  (dominating_fe,
1397  subface_fe,
1398  sf,
1399  subface_interpolation_matrices
1400  [dominating_fe_index][subface_fe_index][sf]);
1401 
1402  const FullMatrix<double> &restrict_subface_to_virtual
1403  = *(subface_interpolation_matrices
1404  [dominating_fe_index][subface_fe_index][sf]);
1405 
1406  constraint_matrix.reinit (subface_fe.dofs_per_face,
1407  dominating_fe.dofs_per_face);
1408 
1409  restrict_subface_to_virtual
1410  .mmult (constraint_matrix,
1411  restrict_mother_to_virtual_master_inv);
1412 
1413  slave_dofs.resize (subface_fe.dofs_per_face);
1414  cell->face(face)->child(sf)->get_dof_indices (slave_dofs,
1415  subface_fe_index);
1416 
1417  filter_constraints (master_dofs,
1418  slave_dofs,
1419  constraint_matrix,
1420  constraints);
1421  }
1422 
1423  break;
1424  }
1425 
1427  // there are no continuity requirements between the two
1428  // elements. record no constraints
1429  break;
1430 
1431  default:
1432  // we shouldn't get here
1433  Assert (false, ExcInternalError());
1434  }
1435  }
1436  else
1437  {
1438  // this face has no children, but it could still be that it
1439  // is shared by two cells that use a different fe index
1440  Assert (cell->face(face)
1441  ->fe_index_is_active(cell->active_fe_index()) == true,
1442  ExcInternalError());
1443 
1444  // see if there is a neighbor that is an artificial cell.
1445  // in that case, we're not interested in this interface. we
1446  // test this case first since artificial cells may not have
1447  // an active_fe_index set, etc
1448  if (!cell->at_boundary(face)
1449  &&
1450  cell->neighbor(face)->is_artificial())
1451  continue;
1452 
1453  // Only if there is a neighbor with a different
1454  // active_fe_index and the same h-level, some action has to
1455  // be taken.
1456  if ((DoFHandlerSupportsDifferentFEs<DoFHandlerType>::value == true)
1457  &&
1458  !cell->face(face)->at_boundary ()
1459  &&
1460  (cell->neighbor(face)->active_fe_index () !=
1461  cell->active_fe_index ())
1462  &&
1463  (!cell->face(face)->has_children() &&
1464  !cell->neighbor_is_coarser(face) ))
1465  {
1466  const typename DoFHandlerType::level_cell_iterator neighbor = cell->neighbor (face);
1467 
1468  // see which side of the face we have to constrain
1469  switch (cell->get_fe().compare_for_face_domination (neighbor->get_fe ()))
1470  {
1472  {
1473  // Get DoFs on dominating and dominated side of the
1474  // face
1475  master_dofs.resize (cell->get_fe().dofs_per_face);
1476  cell->face(face)->get_dof_indices (master_dofs,
1477  cell->active_fe_index ());
1478 
1479  slave_dofs.resize (neighbor->get_fe().dofs_per_face);
1480  cell->face(face)->get_dof_indices (slave_dofs,
1481  neighbor->active_fe_index ());
1482 
1483  // break if the n_master_dofs == 0, because we are
1484  // attempting to constrain to an element that has
1485  // no face dofs
1486  if (master_dofs.size() == 0) break;
1487 
1488  // make sure the element constraints for this face
1489  // are available
1490  ensure_existence_of_face_matrix
1491  (cell->get_fe(),
1492  neighbor->get_fe(),
1493  face_interpolation_matrices
1494  [cell->active_fe_index()][neighbor->active_fe_index()]);
1495 
1496  // Add constraints to global constraint matrix.
1497  filter_constraints (master_dofs,
1498  slave_dofs,
1499  *(face_interpolation_matrices
1500  [cell->active_fe_index()]
1501  [neighbor->active_fe_index()]),
1502  constraints);
1503 
1504  break;
1505  }
1506 
1508  {
1509  // we don't do anything here since we will come
1510  // back to this face from the other cell, at which
1511  // time we will fall into the first case clause
1512  // above
1513  break;
1514  }
1515 
1517  {
1518  // it appears as if neither element has any
1519  // constraints on its neighbor. this may be because
1520  // neither element has any DoFs on faces at all. or
1521  // that the two elements are actually the same,
1522  // although they happen to run under different
1523  // fe_indices (this is what happens in
1524  // hp/hp_hanging_nodes_01 for example).
1525  //
1526  // another possibility is what happens in crash_13.
1527  // there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs.
1528  // FESystem(FE_Q(1),FE_DGQ(1)). neither of them
1529  // dominates the other.
1530  //
1531  // a final possibility is that we have something like
1532  // FESystem(FE_Q(1),FE_Q(1)) vs
1533  // FESystem(FE_Q(1),FE_Nothing()), see
1534  // hp/fe_nothing_18/19.
1535  //
1536  // in any case, the point is that it doesn't
1537  // matter. there is nothing to do here.
1538  break;
1539  }
1540 
1542  {
1543  // make sure we don't get here twice from each cell
1544  if (cell < neighbor)
1545  break;
1546 
1547  // our best bet is to find the common space among other
1548  // FEs in FECollection and then constrain both FEs
1549  // to that one.
1550  // More precisely, we follow the strategy outlined on
1551  // page 17 of the hp paper:
1552  // First we find the dominant FE space S.
1553  // Then we divide our dofs in master and slave such that
1554  // I^{face,master}_{S^{face}->S} is invertible.
1555  // And finally constrain slave dofs to master dofs based
1556  // on the interpolation matrix.
1557 
1558  const unsigned int this_fe_index = cell->active_fe_index();
1559  const unsigned int neighbor_fe_index = neighbor->active_fe_index();
1560  std::set<unsigned int> fes;
1561  fes.insert(this_fe_index);
1562  fes.insert(neighbor_fe_index);
1563  const ::hp::FECollection<dim,spacedim> &fe_collection =
1564  dof_handler.get_fe_collection ();
1565  const unsigned int dominating_fe_index = fe_collection.find_least_face_dominating_fe(fes);
1566 
1567  AssertThrow(dominating_fe_index != numbers::invalid_unsigned_int,
1568  ExcMessage("Could not find the dominating FE for "
1569  +cell->get_fe().get_name()
1570  +" and "
1571  +neighbor->get_fe().get_name()
1572  +" inside FECollection."));
1573 
1574  const FiniteElement<dim,spacedim> &dominating_fe = fe_collection[dominating_fe_index];
1575 
1576  // TODO: until we hit the second face, the code is
1577  // a copy-paste from h-refinement case...
1578 
1579  // first get the interpolation matrix from main FE
1580  // to the virtual dofs
1581  Assert (dominating_fe.dofs_per_face <=
1582  cell->get_fe().dofs_per_face,
1583  ExcInternalError());
1584 
1585  ensure_existence_of_face_matrix
1586  (dominating_fe,
1587  cell->get_fe(),
1588  face_interpolation_matrices
1589  [dominating_fe_index][cell->active_fe_index()]);
1590 
1591  // split this matrix into master and slave components.
1592  // invert the master component
1593  ensure_existence_of_master_dof_mask
1594  (cell->get_fe(),
1595  dominating_fe,
1596  (*face_interpolation_matrices
1597  [dominating_fe_index]
1598  [cell->active_fe_index()]),
1599  master_dof_masks
1600  [dominating_fe_index]
1601  [cell->active_fe_index()]);
1602 
1603  ensure_existence_of_split_face_matrix
1604  (*face_interpolation_matrices
1605  [dominating_fe_index][cell->active_fe_index()],
1606  (*master_dof_masks
1607  [dominating_fe_index][cell->active_fe_index()]),
1608  split_face_interpolation_matrices
1609  [dominating_fe_index][cell->active_fe_index()]);
1610 
1611  const FullMatrix<double> &restrict_mother_to_virtual_master_inv
1612  = (split_face_interpolation_matrices
1613  [dominating_fe_index][cell->active_fe_index()]->first);
1614 
1615  const FullMatrix<double> &restrict_mother_to_virtual_slave
1616  = (split_face_interpolation_matrices
1617  [dominating_fe_index][cell->active_fe_index()]->second);
1618 
1619  // now compute the constraint matrix as the product
1620  // between the inverse matrix and the slave part
1621  constraint_matrix.reinit (cell->get_fe().dofs_per_face -
1622  dominating_fe.dofs_per_face,
1623  dominating_fe.dofs_per_face);
1624  restrict_mother_to_virtual_slave
1625  .mmult (constraint_matrix,
1626  restrict_mother_to_virtual_master_inv);
1627 
1628  // then figure out the global numbers of master and
1629  // slave dofs and apply constraints
1630  scratch_dofs.resize (cell->get_fe().dofs_per_face);
1631  cell->face(face)->get_dof_indices (scratch_dofs,
1632  cell->active_fe_index ());
1633 
1634  // split dofs into master and slave components
1635  master_dofs.clear ();
1636  slave_dofs.clear ();
1637  for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
1638  if ((*master_dof_masks
1639  [dominating_fe_index][cell->active_fe_index()])[i] == true)
1640  master_dofs.push_back (scratch_dofs[i]);
1641  else
1642  slave_dofs.push_back (scratch_dofs[i]);
1643 
1644  AssertDimension (master_dofs.size(), dominating_fe.dofs_per_face);
1645  AssertDimension (slave_dofs.size(),
1646  cell->get_fe().dofs_per_face - dominating_fe.dofs_per_face);
1647 
1648  filter_constraints (master_dofs,
1649  slave_dofs,
1650  constraint_matrix,
1651  constraints);
1652 
1653  // now do the same for another FE
1654  // this is pretty much the same we do above to
1655  // resolve h-refinement constraints
1656  Assert (dominating_fe.dofs_per_face <=
1657  neighbor->get_fe().dofs_per_face,
1658  ExcInternalError());
1659 
1660  ensure_existence_of_face_matrix
1661  (dominating_fe,
1662  neighbor->get_fe(),
1663  face_interpolation_matrices
1664  [dominating_fe_index][neighbor->active_fe_index()]);
1665 
1666  const FullMatrix<double> &restrict_secondface_to_virtual
1667  = *(face_interpolation_matrices
1668  [dominating_fe_index][neighbor->active_fe_index()]);
1669 
1670  constraint_matrix.reinit (neighbor->get_fe().dofs_per_face,
1671  dominating_fe.dofs_per_face);
1672 
1673  restrict_secondface_to_virtual
1674  .mmult (constraint_matrix,
1675  restrict_mother_to_virtual_master_inv);
1676 
1677  slave_dofs.resize (neighbor->get_fe().dofs_per_face);
1678  cell->face(face)->get_dof_indices (slave_dofs,
1679  neighbor->active_fe_index ());
1680 
1681  filter_constraints (master_dofs,
1682  slave_dofs,
1683  constraint_matrix,
1684  constraints);
1685 
1686  break;
1687  }
1688 
1690  {
1691  // nothing to do here
1692  break;
1693  }
1694 
1695  default:
1696  // we shouldn't get here
1697  Assert (false, ExcInternalError());
1698  }
1699  }
1700  }
1701  }
1702  }
1703  }
1704 
1705 
1706 
1707 
1708  template <typename DoFHandlerType>
1709  void
1710  make_hanging_node_constraints (const DoFHandlerType &dof_handler,
1711  ConstraintMatrix &constraints)
1712  {
1713  // Decide whether to use the new or old make_hanging_node_constraints
1714  // function. If all the FiniteElement or all elements in a FECollection
1715  // support the new face constraint matrix, the new code will be used.
1716  // Otherwise, the old implementation is used for the moment.
1717  if (dof_handler.get_fe_collection().hp_constraints_are_implemented ())
1718  internal::
1719  make_hp_hanging_node_constraints (dof_handler,
1720  constraints);
1721  else
1722  internal::
1723  make_oldstyle_hanging_node_constraints (dof_handler,
1724  constraints,
1725  std::integral_constant<int, DoFHandlerType::dimension>());
1726  }
1727 
1728 
1729 
1730  namespace
1731  {
1755  template <typename FaceIterator>
1756  void
1757  set_periodicity_constraints (const FaceIterator &face_1,
1758  const typename identity<FaceIterator>::type &face_2,
1759  const FullMatrix<double> &transformation,
1760  ::ConstraintMatrix &constraint_matrix,
1761  const ComponentMask &component_mask,
1762  const bool face_orientation,
1763  const bool face_flip,
1764  const bool face_rotation)
1765  {
1766  static const int dim = FaceIterator::AccessorType::dimension;
1767  static const int spacedim = FaceIterator::AccessorType::space_dimension;
1768 
1769  // we should be in the case where face_1 is active, i.e. has no children:
1770  Assert (!face_1->has_children(),
1771  ExcInternalError());
1772 
1773  Assert (face_1->n_active_fe_indices() == 1,
1774  ExcInternalError());
1775 
1776  // if face_2 does have children, then we need to iterate over them
1777  if (face_2->has_children())
1778  {
1779  Assert (face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
1780  ExcNotImplemented());
1781  const unsigned int dofs_per_face
1782  = face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face;
1783  FullMatrix<double> child_transformation (dofs_per_face, dofs_per_face);
1784  FullMatrix<double> subface_interpolation (dofs_per_face, dofs_per_face);
1785  for (unsigned int c=0; c<face_2->n_children(); ++c)
1786  {
1787  // get the interpolation matrix recursively from the one that
1788  // interpolated from face_1 to face_2 by multiplying from the
1789  // left with the one that interpolates from face_2 to
1790  // its child
1791  face_1->get_fe(face_1->nth_active_fe_index(0))
1792  .get_subface_interpolation_matrix (face_1->get_fe(face_1->nth_active_fe_index(0)),
1793  c,
1794  subface_interpolation);
1795  subface_interpolation.mmult (child_transformation, transformation);
1796  set_periodicity_constraints(face_1, face_2->child(c),
1797  child_transformation,
1798  constraint_matrix, component_mask,
1799  face_orientation, face_flip, face_rotation);
1800  }
1801  }
1802  else
1803  // both faces are active. we need to match the corresponding DoFs of both faces
1804  {
1805  const unsigned int face_1_index = face_1->nth_active_fe_index(0);
1806  const unsigned int face_2_index = face_2->nth_active_fe_index(0);
1807  Assert(face_1->get_fe(face_1_index) == face_2->get_fe(face_2_index),
1808  ExcMessage ("Matching periodic cells need to use the same finite element"));
1809 
1810  const FiniteElement<dim, spacedim> &fe = face_1->get_fe(face_1_index);
1811 
1812  Assert(component_mask.represents_n_components(fe.n_components()),
1813  ExcMessage ("The number of components in the mask has to be either "
1814  "zero or equal to the number of components in the finite " "element."));
1815 
1816  const unsigned int dofs_per_face = fe.dofs_per_face;
1817 
1818  std::vector<types::global_dof_index> dofs_1(dofs_per_face);
1819  std::vector<types::global_dof_index> dofs_2(dofs_per_face);
1820 
1821  face_1->get_dof_indices(dofs_1, face_1_index);
1822  face_2->get_dof_indices(dofs_2, face_2_index);
1823 
1824  for (unsigned int i=0; i < dofs_per_face; i++)
1825  {
1826  if (dofs_1[i] == numbers::invalid_dof_index ||
1827  dofs_2[i] == numbers::invalid_dof_index)
1828  {
1829  /* If either of these faces have no indices, stop. This is so
1830  * that there is no attempt to match artificial cells of
1831  * parallel distributed triangulations.
1832  *
1833  * While it seems like we ought to be able to avoid even calling
1834  * set_periodicity_constraints for artificial faces, this
1835  * situation can arise when a face that is being made periodic
1836  * is only partially touched by the local subdomain.
1837  * make_periodicity_constraints will be called recursively even
1838  * for the section of the face that is not touched by the local
1839  * subdomain.
1840  *
1841  * Until there is a better way to determine if the cells that
1842  * neighbor a face are artificial, we simply test to see if the
1843  * face does not have a valid dof initialization.
1844  */
1845  return;
1846  }
1847  }
1848 
1849  // Well, this is a hack:
1850  //
1851  // There is no
1852  // face_to_face_index(face_index,
1853  // face_orientation,
1854  // face_flip,
1855  // face_rotation)
1856  // function in FiniteElementData, so we have to use
1857  // face_to_cell_index(face_index, face
1858  // face_orientation,
1859  // face_flip,
1860  // face_rotation)
1861  // But this will give us an index on a cell - something we cannot work
1862  // with directly. But luckily we can match them back :-]
1863 
1864  std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
1865 
1866  // Build up a cell to face index for face_2:
1867  for (unsigned int i = 0; i < dofs_per_face; ++i)
1868  {
1869  const unsigned int cell_index = fe.face_to_cell_index(i, 0, /* It doesn't really matter, just assume
1870  * we're on the first face...
1871  */
1872  true, false, false // default orientation
1873  );
1874  cell_to_rotated_face_index[cell_index] = i;
1875  }
1876 
1877  // loop over all dofs on face 2 and constrain them against the ones on face 1
1878  for (unsigned int i=0; i<dofs_per_face; ++i)
1879  if (!constraint_matrix.is_constrained(dofs_2[i]))
1880  if ((component_mask.n_selected_components(fe.n_components())
1881  == fe.n_components())
1882  ||
1883  component_mask[fe.face_system_to_component_index(i).first])
1884  {
1885  // as mentioned in the comment above this function, we need
1886  // to be careful about treating identity constraints differently.
1887  // consequently, find out whether this dof 'i' will be
1888  // identity constrained
1889  //
1890  // to check whether this is the case, first see whether there are
1891  // any weights other than 0 and 1, then in a first stage make sure
1892  // that if so there is only one weight equal to 1
1893  //
1894  // afterwards do the same for constraints of type dof1=-dof2
1895  bool is_identity_constrained = true;
1896  const double eps = 1.e-13;
1897  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
1898  if (((std::abs(transformation(i,jj)) < eps) ||
1899  (std::abs(transformation(i,jj)-1) < eps)) == false)
1900  {
1901  is_identity_constrained = false;
1902  break;
1903  }
1904  unsigned int identity_constraint_target = numbers::invalid_unsigned_int;
1905  if (is_identity_constrained == true)
1906  {
1907  bool one_identity_found = false;
1908  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
1909  if (std::abs(transformation(i,jj)-1.) < eps)
1910  {
1911  if (one_identity_found == false)
1912  {
1913  one_identity_found = true;
1914  identity_constraint_target = jj;
1915  }
1916  else
1917  {
1918  is_identity_constrained = false;
1919  identity_constraint_target = numbers::invalid_unsigned_int;
1920  break;
1921  }
1922  }
1923  }
1924 
1925  bool is_inverse_constrained = !is_identity_constrained;
1926  unsigned int inverse_constraint_target = numbers::invalid_unsigned_int;
1927  if (is_inverse_constrained)
1928  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
1929  if (((std::abs(transformation(i,jj)) < eps) ||
1930  (std::abs(transformation(i,jj)+1) < eps)) == false)
1931  {
1932  is_inverse_constrained = false;
1933  break;
1934  }
1935  if (is_inverse_constrained)
1936  {
1937  bool one_identity_found = false;
1938  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
1939  if (std::abs(transformation(i,jj)+1) < eps)
1940  {
1941  if (one_identity_found == false)
1942  {
1943  one_identity_found = true;
1944  inverse_constraint_target = jj;
1945  }
1946  else
1947  {
1948  is_inverse_constrained = false;
1949  inverse_constraint_target = numbers::invalid_unsigned_int;
1950  break;
1951  }
1952  }
1953  }
1954 
1955  const unsigned int target = is_identity_constrained
1956  ? identity_constraint_target
1957  : inverse_constraint_target;
1958 
1959  // find out whether this dof also exists on face 1
1960  // if this is true and the constraint is no identity
1961  // constraint to itself, set it to zero
1962  bool constrained_set = false;
1963  for (unsigned int j=0; j<dofs_per_face; ++j)
1964  {
1965  if (dofs_2[i] == dofs_1[j])
1966  if (!(is_identity_constrained && target==i))
1967  {
1968  constraint_matrix.add_line(dofs_2[i]);
1969  constrained_set = true;
1970  }
1971  }
1972 
1973  if (!constrained_set)
1974  {
1975  // now treat constraints, either as an equality constraint or
1976  // as a sequence of constraints
1977  if (is_identity_constrained == true || is_inverse_constrained == true)
1978  {
1979  // Query the correct face_index on face_1 respecting the given
1980  // orientation:
1981  const unsigned int j
1982  = cell_to_rotated_face_index[fe.face_to_cell_index(target,
1983  0, /* It doesn't really matter, just assume
1984  * we're on the first face...
1985  */
1986  face_orientation, face_flip, face_rotation)];
1987 
1988  // if the two aren't already identity constrained (whichever way
1989  // around) or already identical (in case of rotated periodicity constraints),
1990  // then enter the constraint. otherwise there is nothing for us still to do
1991  bool enter_constraint = false;
1992  if (!constraint_matrix.is_constrained(dofs_1[j]))
1993  {
1994  if (dofs_2[i] != dofs_1[j])
1995  enter_constraint = true;
1996  }
1997  else //dofs_1[j] is constrained, is it identity or inverse constrained?
1998  {
1999  const std::vector<std::pair<types::global_dof_index, double > > *constraint_entries
2000  = constraint_matrix.get_constraint_entries(dofs_1[j]);
2001  if (constraint_entries->size()==1 && (*constraint_entries)[0].first == dofs_2[i])
2002  {
2003  if ((is_identity_constrained && std::abs((*constraint_entries)[0].second-1) > eps) ||
2004  (is_inverse_constrained && std::abs((*constraint_entries)[0].second+1) > eps))
2005  {
2006  //this pair of constraints means that both dofs have to be constrained to 0.
2007  constraint_matrix.add_line(dofs_2[i]);
2008  }
2009  }
2010  else
2011  enter_constraint = true;
2012  }
2013 
2014  if (enter_constraint)
2015  {
2016  constraint_matrix.add_line(dofs_2[i]);
2017  constraint_matrix.add_entry(dofs_2[i], dofs_1[j], is_identity_constrained?1.0:-1.0);
2018  }
2019  }
2020  else
2021  {
2022  // this is just a regular constraint. enter it piece by piece
2023  constraint_matrix.add_line(dofs_2[i]);
2024  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
2025  {
2026  // Query the correct face_index on face_1 respecting the given
2027  // orientation:
2028  const unsigned int j =
2029  cell_to_rotated_face_index[fe.face_to_cell_index
2030  (jj, 0, face_orientation, face_flip, face_rotation)];
2031 
2032  // And finally constrain the two DoFs respecting component_mask:
2033  if (transformation(i,jj) != 0)
2034  constraint_matrix.add_entry(dofs_2[i], dofs_1[j],
2035  transformation(i,jj));
2036  }
2037  }
2038  }
2039  }
2040  }
2041  }
2042 
2043 
2044  // Internally used in make_periodicity_constraints.
2045  //
2046  // Build up a (possibly rotated) interpolation matrix that is used in
2047  // set_periodicity_constraints with the help of user supplied matrix
2048  // and first_vector_components.
2049  template <int dim, int spacedim>
2050  FullMatrix<double> compute_transformation(
2051  const FiniteElement<dim, spacedim> &fe,
2052  const FullMatrix<double> &matrix,
2053  const std::vector<unsigned int> &first_vector_components)
2054  {
2055  Assert(matrix.m() == matrix.n(), ExcInternalError());
2056 
2057  const unsigned int n_dofs_per_face = fe.dofs_per_face;
2058 
2059  if (matrix.m() == n_dofs_per_face)
2060  {
2061  // In case of m == n == n_dofs_per_face the supplied matrix is already
2062  // an interpolation matrix, so we use it directly:
2063  return matrix;
2064  }
2065 
2066  if (first_vector_components.empty() && matrix.m() == 0)
2067  {
2068  // Just the identity matrix in case no rotation is specified:
2069  return IdentityMatrix(n_dofs_per_face);
2070  }
2071 
2072  // The matrix describes a rotation and we have to build a
2073  // transformation matrix, we assume that for a 0* rotation
2074  // we would have to build the identity matrix
2075 
2076  Assert(matrix.m() == (int)spacedim, ExcInternalError())
2077 
2079 
2080  // have an array that stores the location of each vector-dof tuple
2081  // we want to rotate.
2082  typedef std::array<unsigned int, spacedim> DoFTuple;
2083 
2084  // start with a pristine interpolation matrix...
2085  FullMatrix<double> transformation = IdentityMatrix(n_dofs_per_face);
2086 
2087  for (unsigned int i=0; i < n_dofs_per_face; ++i)
2088  {
2089  std::vector<unsigned int>::const_iterator comp_it
2090  = std::find (first_vector_components.begin(),
2091  first_vector_components.end(),
2092  fe.face_system_to_component_index(i).first);
2093  if (comp_it != first_vector_components.end())
2094  {
2095  const unsigned int first_vector_component = *comp_it;
2096 
2097  // find corresponding other components of vector
2098  DoFTuple vector_dofs;
2099  vector_dofs[0] = i;
2100  unsigned int n_found = 1;
2101 
2102  Assert(*comp_it + spacedim <= fe.n_components(),
2103  ExcMessage("Error: the finite element does not have enough components "
2104  "to define rotated periodic boundaries."));
2105 
2106  for (unsigned int k = 0; k < n_dofs_per_face; ++k)
2107  if ((k != i) &&
2108  (quadrature.point(k) == quadrature.point(i)) &&
2109  (fe.face_system_to_component_index(k).first >=
2110  first_vector_component) &&
2111  (fe.face_system_to_component_index(k).first <
2112  first_vector_component + spacedim))
2113  {
2114  vector_dofs[fe.face_system_to_component_index(k).first -
2115  first_vector_component]
2116  = k;
2117  n_found++;
2118  if (n_found==dim)
2119  break;
2120  }
2121 
2122  // ... and rotate all dofs belonging to vector valued
2123  // components that are selected by first_vector_components:
2124  for (int i = 0; i < spacedim; ++i)
2125  {
2126  transformation[vector_dofs[i]][vector_dofs[i]] = 0.;
2127  for (int j = 0; j < spacedim; ++j)
2128  transformation[vector_dofs[i]][vector_dofs[j]] = matrix[i][j];
2129  }
2130  }
2131  }
2132  return transformation;
2133  }
2134 
2135  } /*namespace*/
2136 
2137 
2138  // Low level interface:
2139 
2140 
2141  template <typename FaceIterator>
2142  void
2143  make_periodicity_constraints (const FaceIterator &face_1,
2144  const typename identity<FaceIterator>::type &face_2,
2145  ::ConstraintMatrix &constraint_matrix,
2146  const ComponentMask &component_mask,
2147  const bool face_orientation,
2148  const bool face_flip,
2149  const bool face_rotation,
2150  const FullMatrix<double> &matrix,
2151  const std::vector<unsigned int> &first_vector_components)
2152  {
2153  static const int dim = FaceIterator::AccessorType::dimension;
2154  static const int spacedim = FaceIterator::AccessorType::space_dimension;
2155 
2156  Assert( (dim != 1) ||
2157  (face_orientation == true &&
2158  face_flip == false &&
2159  face_rotation == false),
2160  ExcMessage ("The supplied orientation "
2161  "(face_orientation, face_flip, face_rotation) "
2162  "is invalid for 1D"));
2163 
2164  Assert( (dim != 2) ||
2165  (face_orientation == true &&
2166  face_rotation == false),
2167  ExcMessage ("The supplied orientation "
2168  "(face_orientation, face_flip, face_rotation) "
2169  "is invalid for 2D"));
2170 
2171  Assert(face_1 != face_2,
2172  ExcMessage ("face_1 and face_2 are equal! Cannot constrain DoFs "
2173  "on the very same face"));
2174 
2175  Assert(face_1->at_boundary() && face_2->at_boundary(),
2176  ExcMessage ("Faces for periodicity constraints must be on the "
2177  "boundary"));
2178 
2179  Assert(matrix.m() == matrix.n(),
2180  ExcMessage("The supplied (rotation or interpolation) matrix must "
2181  "be a square matrix"));
2182 
2183  Assert(first_vector_components.empty() || matrix.m() == (int)spacedim,
2184  ExcMessage ("first_vector_components is nonempty, so matrix must "
2185  "be a rotation matrix exactly of size spacedim"));
2186 
2187 #ifdef DEBUG
2188  if (!face_1->has_children())
2189  {
2190  Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
2191  const unsigned int n_dofs_per_face =
2192  face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face;
2193 
2194  Assert(matrix.m() == 0
2195  || (first_vector_components.empty() && matrix.m() == n_dofs_per_face)
2196  || (!first_vector_components.empty() && matrix.m() == (int)spacedim),
2197  ExcMessage ("The matrix must have either size 0 or spacedim "
2198  "(if first_vector_components is nonempty) "
2199  "or the size must be equal to the # of DoFs on the face "
2200  "(if first_vector_components is empty)."));
2201  }
2202 
2203  if (!face_2->has_children())
2204  {
2205  Assert(face_2->n_active_fe_indices() == 1, ExcInternalError());
2206  const unsigned int n_dofs_per_face =
2207  face_2->get_fe(face_2->nth_active_fe_index(0)).dofs_per_face;
2208 
2209  Assert(matrix.m() == 0
2210  || (first_vector_components.empty() && matrix.m() == n_dofs_per_face)
2211  || (!first_vector_components.empty() && matrix.m() == (int)spacedim),
2212  ExcMessage ("The matrix must have either size 0 or spacedim "
2213  "(if first_vector_components is nonempty) "
2214  "or the size must be equal to the # of DoFs on the face "
2215  "(if first_vector_components is empty)."));
2216  }
2217 #endif
2218 
2219  // A lookup table on how to go through the child faces depending on the
2220  // orientation:
2221 
2222  static const int lookup_table_2d[2][2] =
2223  {
2224  // flip:
2225  {0, 1}, // false
2226  {1, 0}, // true
2227  };
2228 
2229  static const int lookup_table_3d[2][2][2][4] =
2230  {
2231  // orientation flip rotation
2232  { { {0, 2, 1, 3}, // false false false
2233  {2, 3, 0, 1}, // false false true
2234  },
2235  { {3, 1, 2, 0}, // false true false
2236  {1, 0, 3, 2}, // false true true
2237  },
2238  },
2239  { { {0, 1, 2, 3}, // true false false
2240  {1, 3, 0, 2}, // true false true
2241  },
2242  { {3, 2, 1, 0}, // true true false
2243  {2, 0, 3, 1}, // true true true
2244  },
2245  },
2246  };
2247 
2248  if (face_1->has_children() && face_2->has_children())
2249  {
2250  // In the case that both faces have children, we loop over all
2251  // children and apply make_periodicty_constrains recursively:
2252 
2253  Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
2254  face_2->n_children() ==
2256  ExcNotImplemented());
2257 
2258  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face;
2259  ++i)
2260  {
2261  // Lookup the index for the second face
2262  unsigned int j;
2263  switch (dim)
2264  {
2265  case 2:
2266  j = lookup_table_2d[face_flip][i];
2267  break;
2268  case 3:
2269  j = lookup_table_3d[face_orientation][face_flip][face_rotation][i];
2270  break;
2271  default:
2272  AssertThrow(false, ExcNotImplemented());
2273  }
2274 
2275  make_periodicity_constraints (face_1->child(i),
2276  face_2->child(j),
2277  constraint_matrix,
2278  component_mask,
2279  face_orientation,
2280  face_flip,
2281  face_rotation,
2282  matrix,
2283  first_vector_components);
2284  }
2285  }
2286  else
2287  {
2288  // Otherwise at least one of the two faces is active and
2289  // we need to do some work and enter the constraints!
2290 
2291  // The finite element that matters is the one on the active face:
2292  const FiniteElement<dim,spacedim> &fe =
2293  face_1->has_children()
2294  ? face_2->get_fe(face_2->nth_active_fe_index(0))
2295  : face_1->get_fe(face_1->nth_active_fe_index(0));
2296 
2297  const unsigned int n_dofs_per_face = fe.dofs_per_face;
2298 
2299  // Sometimes we just have nothing to do (for all finite elements,
2300  // or systems which accidentally don't have any dofs on the
2301  // boundary).
2302  if (n_dofs_per_face == 0)
2303  return;
2304 
2305  const FullMatrix<double> transformation =
2306  compute_transformation(fe, matrix, first_vector_components);
2307 
2308  if (! face_2->has_children())
2309  {
2310  // Performance hack: We do not need to compute an inverse if
2311  // the matrix is the identity matrix.
2312  if (first_vector_components.empty() && matrix.m() == 0)
2313  {
2314  set_periodicity_constraints(face_2,
2315  face_1,
2316  transformation,
2317  constraint_matrix,
2318  component_mask,
2319  face_orientation,
2320  face_flip,
2321  face_rotation);
2322  }
2323  else
2324  {
2325  FullMatrix<double> inverse(transformation.m());
2326  inverse.invert(transformation);
2327 
2328  set_periodicity_constraints(face_2,
2329  face_1,
2330  inverse,
2331  constraint_matrix,
2332  component_mask,
2333  face_orientation,
2334  face_flip,
2335  face_rotation);
2336  }
2337  }
2338  else
2339  {
2340  Assert(!face_1->has_children(), ExcInternalError());
2341 
2342  // Important note:
2343  // In 3D we have to take care of the fact that face_rotation
2344  // gives the relative rotation of face_1 to face_2, i.e. we
2345  // have to invert the rotation when constraining face_2 to
2346  // face_1. Therefore face_flip has to be toggled if
2347  // face_rotation is true:
2348  // In case of inverted orientation, nothing has to be done.
2349  set_periodicity_constraints(face_1,
2350  face_2,
2351  transformation,
2352  constraint_matrix,
2353  component_mask,
2354  face_orientation,
2355  face_orientation
2356  ? face_rotation ^ face_flip
2357  : face_flip,
2358  face_rotation);
2359  }
2360  }
2361  }
2362 
2363 
2364 
2365  template <typename DoFHandlerType>
2366  void
2369  &periodic_faces,
2370  ::ConstraintMatrix &constraint_matrix,
2371  const ComponentMask &component_mask,
2372  const std::vector<unsigned int> &first_vector_components)
2373  {
2374  typedef std::vector<GridTools::PeriodicFacePair<typename DoFHandlerType::cell_iterator> >
2375  FaceVector;
2376  typename FaceVector::const_iterator it, end_periodic;
2377  it = periodic_faces.begin();
2378  end_periodic = periodic_faces.end();
2379 
2380  // Loop over all periodic faces...
2381  for (; it!=end_periodic; ++it)
2382  {
2383  typedef typename DoFHandlerType::face_iterator FaceIterator;
2384  const FaceIterator face_1 = it->cell[0]->face(it->face_idx[0]);
2385  const FaceIterator face_2 = it->cell[1]->face(it->face_idx[1]);
2386 
2387  Assert(face_1->at_boundary() && face_2->at_boundary(),
2388  ExcInternalError());
2389 
2390  Assert (face_1 != face_2,
2391  ExcInternalError());
2392 
2393  // ... and apply the low level make_periodicity_constraints function to
2394  // every matching pair:
2396  face_2,
2397  constraint_matrix,
2398  component_mask,
2399  it->orientation[0],
2400  it->orientation[1],
2401  it->orientation[2],
2402  it->matrix,
2403  first_vector_components);
2404  }
2405  }
2406 
2407 
2408  // High level interface variants:
2409 
2410 
2411  template <typename DoFHandlerType>
2412  void
2413  make_periodicity_constraints (const DoFHandlerType &dof_handler,
2414  const types::boundary_id b_id1,
2415  const types::boundary_id b_id2,
2416  const int direction,
2417  ::ConstraintMatrix &constraint_matrix,
2418  const ComponentMask &component_mask)
2419  {
2420  static const int space_dim = DoFHandlerType::space_dimension;
2421  (void)space_dim;
2422  Assert (0<=direction && direction<space_dim,
2423  ExcIndexRange (direction, 0, space_dim));
2424 
2425  Assert (b_id1 != b_id2,
2426  ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
2427  "different to denote different boundaries."));
2428 
2429  std::vector<GridTools::PeriodicFacePair
2430  <typename DoFHandlerType::cell_iterator> > matched_faces;
2431 
2432  // Collect matching periodic cells on the coarsest level:
2433  GridTools::collect_periodic_faces(dof_handler, b_id1, b_id2, direction,
2434  matched_faces);
2435 
2436  make_periodicity_constraints<DoFHandlerType>
2437  (matched_faces, constraint_matrix, component_mask);
2438  }
2439 
2440 
2441 
2442  template <typename DoFHandlerType>
2443  void
2444  make_periodicity_constraints (const DoFHandlerType &dof_handler,
2445  const types::boundary_id b_id,
2446  const int direction,
2447  ::ConstraintMatrix &constraint_matrix,
2448  const ComponentMask &component_mask)
2449  {
2450  static const int dim = DoFHandlerType::dimension;
2451  static const int space_dim = DoFHandlerType::space_dimension;
2452  (void)dim;
2453  (void)space_dim;
2454 
2455  Assert (0<=direction && direction<space_dim,
2456  ExcIndexRange (direction, 0, space_dim));
2457 
2458  Assert(dim == space_dim,
2459  ExcNotImplemented());
2460 
2461  std::vector<GridTools::PeriodicFacePair
2462  <typename DoFHandlerType::cell_iterator> > matched_faces;
2463 
2464  // Collect matching periodic cells on the coarsest level:
2465  GridTools::collect_periodic_faces(dof_handler, b_id, direction,
2466  matched_faces);
2467 
2468  make_periodicity_constraints<DoFHandlerType>
2469  (matched_faces, constraint_matrix, component_mask);
2470  }
2471 
2472 
2473 
2474  namespace internal
2475  {
2476  namespace Assembler
2477  {
2478  struct Scratch {};
2479 
2480 
2481  template <int dim,int spacedim>
2482  struct CopyData
2483  {
2484  unsigned int dofs_per_cell;
2485  std::vector<types::global_dof_index> parameter_dof_indices;
2486 #ifdef DEAL_II_WITH_MPI
2487  std::vector<::LinearAlgebra::distributed::Vector<double> > global_parameter_representation;
2488 #else
2489  std::vector<::Vector<double> > global_parameter_representation;
2490 #endif
2491  };
2492  }
2493 
2494  namespace
2495  {
2501  template <int dim, int spacedim>
2502  void compute_intergrid_weights_3 (
2503  const typename ::DoFHandler<dim,spacedim>::active_cell_iterator &cell,
2504  const Assembler::Scratch &,
2505  Assembler::CopyData<dim,spacedim> &copy_data,
2506  const unsigned int coarse_component,
2507  const FiniteElement<dim,spacedim> &coarse_fe,
2508  const InterGridMap<::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
2509  const std::vector<::Vector<double> > &parameter_dofs)
2510  {
2511  // for each cell on the parameter grid: find out which degrees of
2512  // freedom on the fine grid correspond in which way to the degrees
2513  // of freedom on the parameter grid
2514  //
2515  // since for continuous FEs some dofs exist on more than one cell,
2516  // we have to track which ones were already visited. the problem is
2517  // that if we visit a dof first on one cell and compute its weight
2518  // with respect to some global dofs to be non-zero, and later visit
2519  // the dof again on another cell and (since we are on another cell)
2520  // recompute the weights with respect to the same dofs as above to
2521  // be zero now, we have to preserve them. we therefore overwrite
2522  // all weights if they are nonzero and do not enforce zero weights
2523  // since that might be only due to the fact that we are on another
2524  // cell.
2525  //
2526  // example:
2527  // coarse grid
2528  // | | |
2529  // *-----*-----*
2530  // | cell|cell |
2531  // | 1 | 2 |
2532  // | | |
2533  // 0-----1-----*
2534  //
2535  // fine grid
2536  // | | | | |
2537  // *--*--*--*--*
2538  // | | | | |
2539  // *--*--*--*--*
2540  // | | | | |
2541  // *--x--y--*--*
2542  //
2543  // when on cell 1, we compute the weights of dof 'x' to be 1/2 from
2544  // parameter dofs 0 and 1, respectively. however, when later we are
2545  // on cell 2, we again compute the prolongation of shape function 1
2546  // restricted to cell 2 to the globla grid and find that the weight
2547  // of global dof 'x' now is zero. however, we should not overwrite
2548  // the old value.
2549  //
2550  // we therefore always only set nonzero values. why adding up is
2551  // not useful: dof 'y' would get weight 1 from parameter dof 1 on
2552  // both cells 1 and 2, but the correct weight is nevertheless only
2553  // 1.
2554 
2555  // vector to hold the representation of a single degree of freedom
2556  // on the coarse grid (for the selected fe) on the fine grid
2557 
2558  copy_data.dofs_per_cell = coarse_fe.dofs_per_cell;
2559  copy_data.parameter_dof_indices.resize(copy_data.dofs_per_cell);
2560 
2561  // get the global indices of the parameter dofs on this
2562  // parameter grid cell
2563  cell->get_dof_indices (copy_data.parameter_dof_indices);
2564 
2565  // loop over all dofs on this cell and check whether they are
2566  // interesting for us
2567  for (unsigned int local_dof=0; local_dof<copy_data.dofs_per_cell; ++local_dof)
2568  if (coarse_fe.system_to_component_index(local_dof).first
2569  ==
2570  coarse_component)
2571  {
2572  // the how-many-th parameter is this on this cell?
2573  const unsigned int local_parameter_dof
2574  = coarse_fe.system_to_component_index(local_dof).second;
2575 
2576  copy_data.global_parameter_representation[local_parameter_dof] = 0.;
2577 
2578  // distribute the representation of
2579  // @p{local_parameter_dof} on the parameter grid cell
2580  // @p{cell} to the global data space
2581  coarse_to_fine_grid_map[cell]->
2582  set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof],
2583  copy_data.global_parameter_representation[local_parameter_dof]);
2584  }
2585  }
2586 
2587 
2588 
2594  template <int dim,int spacedim>
2595  void copy_intergrid_weights_3(const Assembler::CopyData<dim,spacedim> &copy_data,
2596  const unsigned int coarse_component,
2597  const FiniteElement<dim,spacedim> &coarse_fe,
2598  const std::vector<types::global_dof_index> &weight_mapping,
2599  const bool is_called_in_parallel,
2600  std::vector<std::map<types::global_dof_index, float> > &weights)
2601  {
2602  unsigned int pos = 0;
2603  for (unsigned int local_dof=0; local_dof<copy_data.dofs_per_cell; ++local_dof)
2604  if (coarse_fe.system_to_component_index(local_dof).first
2605  ==
2606  coarse_component)
2607  {
2608  // now that we've got the global representation of each
2609  // parameter dof, we've only got to clobber the non-zero
2610  // entries in that vector and store the result
2611  //
2612  // what we have learned: if entry @p{i} of the global
2613  // vector holds the value @p{v[i]}, then this is the
2614  // weight with which the present dof contributes to
2615  // @p{i}. there may be several such @p{i}s and their
2616  // weights' sum should be one. Then, @p{v[i]} should be
2617  // equal to @p{\sum_j w_{ij} p[j]} with @p{p[j]} be the
2618  // values of the degrees of freedom on the coarse grid.
2619  // we can thus compute constraints which link the degrees
2620  // of freedom @p{v[i]} on the fine grid to those on the
2621  // coarse grid, @p{p[j]}. Now to use these as real
2622  // constraints, rather than as additional equations, we
2623  // have to identify representants among the @p{i} for
2624  // each @p{j}. this will be done by simply taking the
2625  // first @p{i} for which @p{w_{ij}==1}.
2626  //
2627  // guard modification of the weights array by a Mutex.
2628  // since it should happen rather rarely that there are
2629  // several threads operating on different intergrid
2630  // weights, have only one mutex for all of them
2631  for (types::global_dof_index i=0; i<copy_data.global_parameter_representation[pos].size();
2632  ++i)
2633  // set this weight if it belongs to a parameter dof.
2634  if (weight_mapping[i] != numbers::invalid_dof_index)
2635  {
2636  // only overwrite old value if not by zero
2637  if (copy_data.global_parameter_representation[pos](i) != 0)
2638  {
2639  const types::global_dof_index wi = copy_data.parameter_dof_indices[local_dof],
2640  wj = weight_mapping[i];
2641  weights[wi][wj] = copy_data.global_parameter_representation[pos](i);
2642  }
2643  }
2644  else if (!is_called_in_parallel)
2645  {
2646  // Note that when this function operates with distributed
2647  // fine grid, this assertion is switched off since the
2648  // condition does not necessarily hold
2649  Assert (copy_data.global_parameter_representation[pos](i) == 0,
2650  ExcInternalError());
2651  }
2652 
2653  ++pos;
2654  }
2655 
2656  }
2657 
2658 
2659 
2665  template <int dim, int spacedim>
2666  void
2667  compute_intergrid_weights_2 (
2668  const ::DoFHandler<dim,spacedim> &coarse_grid,
2669  const unsigned int coarse_component,
2670  const InterGridMap<::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
2671  const std::vector<::Vector<double> > &parameter_dofs,
2672  const std::vector<types::global_dof_index> &weight_mapping,
2673  std::vector<std::map<types::global_dof_index,float> > &weights)
2674  {
2675  Assembler::Scratch scratch;
2676  Assembler::CopyData<dim,spacedim> copy_data;
2677 
2678  unsigned int n_interesting_dofs = 0;
2679  for (unsigned int local_dof=0; local_dof<coarse_grid.get_fe().dofs_per_cell; ++local_dof)
2680  if (coarse_grid.get_fe().system_to_component_index(local_dof).first
2681  ==
2682  coarse_component)
2683  ++n_interesting_dofs;
2684 
2685  copy_data.global_parameter_representation.resize(n_interesting_dofs);
2686 
2687  bool is_called_in_parallel = false;
2688  for (size_t i=0; i<copy_data.global_parameter_representation.size(); ++i)
2689  {
2690 #ifdef DEAL_II_WITH_MPI
2691  MPI_Comm communicator = MPI_COMM_SELF;
2692  try
2693  {
2694  const typename ::parallel::Triangulation<dim, spacedim> &tria =
2695  dynamic_cast<const typename ::parallel::Triangulation<dim, spacedim>&>
2696  (coarse_to_fine_grid_map.get_destination_grid().get_triangulation());
2697  communicator = tria.get_communicator ();
2698  is_called_in_parallel = true;
2699  }
2700  catch (std::bad_cast &exp)
2701  {
2702  // Nothing bad happened: the user used serial Triangulation
2703  }
2704 
2705  IndexSet locally_owned_dofs, locally_relevant_dofs;
2707  (coarse_to_fine_grid_map.get_destination_grid (), locally_owned_dofs);
2709  (coarse_to_fine_grid_map.get_destination_grid (), locally_relevant_dofs);
2710 
2711  copy_data.global_parameter_representation[i].reinit
2712  (locally_owned_dofs, locally_relevant_dofs, communicator);
2713 #else
2714  const types::global_dof_index n_fine_dofs = weight_mapping.size();
2715  copy_data.global_parameter_representation[i].reinit (n_fine_dofs);
2716 #endif
2717  }
2718 
2719  WorkStream::run(coarse_grid.begin_active(),
2720  coarse_grid.end(),
2721  std::bind(&compute_intergrid_weights_3<dim,spacedim>,
2722  std::placeholders::_1,
2723  std::placeholders::_2,
2724  std::placeholders::_3,
2725  coarse_component,
2726  std::cref(coarse_grid.get_fe()),
2727  std::cref(coarse_to_fine_grid_map),
2728  std::cref(parameter_dofs)),
2729  std::bind(&copy_intergrid_weights_3<dim,spacedim>,
2730  std::placeholders::_1,
2731  coarse_component,
2732  std::cref(coarse_grid.get_fe()),
2733  std::cref(weight_mapping),
2734  is_called_in_parallel,
2735  std::ref(weights)),
2736  scratch,
2737  copy_data);
2738 
2739 #ifdef DEAL_II_WITH_MPI
2740  for (size_t i=0; i<copy_data.global_parameter_representation.size(); ++i)
2741  copy_data.global_parameter_representation[i].update_ghost_values ();
2742 #endif
2743  }
2744 
2745 
2746 
2752  template <int dim, int spacedim>
2753  unsigned int
2754  compute_intergrid_weights_1 (
2755  const ::DoFHandler<dim,spacedim> &coarse_grid,
2756  const unsigned int coarse_component,
2757  const ::DoFHandler<dim,spacedim> &fine_grid,
2758  const unsigned int fine_component,
2759  const InterGridMap<::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
2760  std::vector<std::map<types::global_dof_index, float> > &weights,
2761  std::vector<types::global_dof_index> &weight_mapping)
2762  {
2763  // aliases to the finite elements used by the dof handlers:
2764  const FiniteElement<dim,spacedim> &coarse_fe = coarse_grid.get_fe(),
2765  &fine_fe = fine_grid.get_fe();
2766 
2767  // global numbers of dofs
2768  const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(),
2769  n_fine_dofs = fine_grid.n_dofs();
2770 
2771  // local numbers of dofs
2772  const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell;
2773 
2774  // alias the number of dofs per cell belonging to the
2775  // coarse_component which is to be the restriction of the fine
2776  // grid:
2777  const unsigned int coarse_dofs_per_cell_component
2778  = coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell;
2779 
2780 
2781  // Try to find out whether the grids stem from the same coarse
2782  // grid. This is a rather crude test, but better than nothing
2783  Assert (coarse_grid.get_triangulation().n_cells(0) == fine_grid.get_triangulation().n_cells(0),
2784  ExcGridsDontMatch());
2785 
2786  // check whether the map correlates the right objects
2787  Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid,
2788  ExcGridsDontMatch ());
2789  Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid,
2790  ExcGridsDontMatch ());
2791 
2792 
2793  // check whether component numbers are valid
2794  AssertIndexRange (coarse_component,coarse_fe.n_components());
2795  AssertIndexRange (fine_component, fine_fe.n_components());
2796 
2797  // check whether respective finite elements are equal
2798  Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first)
2799  ==
2800  fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first),
2802 
2803 #ifdef DEBUG
2804  // if in debug mode, check whether the coarse grid is indeed
2805  // coarser everywhere than the fine grid
2806  for (typename ::DoFHandler<dim,spacedim>::active_cell_iterator
2807  cell=coarse_grid.begin_active();
2808  cell != coarse_grid.end(); ++cell)
2809  Assert (cell->level() <= coarse_to_fine_grid_map[cell]->level(),
2810  ExcGridNotCoarser());
2811 #endif
2812 
2813  /*
2814  * From here on: the term `parameter' refers to the selected
2815  * component on the coarse grid and its analogon on the fine grid.
2816  * The naming of variables containing this term is due to the fact
2817  * that `selected_component' is longer, but also due to the fact
2818  * that the code of this function was initially written for a
2819  * program where the component which we wanted to match between
2820  * grids was actually the `parameter' variable.
2821  *
2822  * Likewise, the terms `parameter grid' and `state grid' refer to
2823  * the coarse and fine grids, respectively.
2824  *
2825  * Changing the names of variables would in principle be a good
2826  * idea, but would not make things simpler and would be another
2827  * source of errors. If anyone feels like doing so: patches would
2828  * be welcome!
2829  */
2830 
2831 
2832 
2833  // set up vectors of cell-local data; each vector represents one
2834  // degree of freedom of the coarse-grid variable in the fine-grid
2835  // element
2836  std::vector<::Vector<double> >
2837  parameter_dofs (coarse_dofs_per_cell_component,
2838  ::Vector<double>(fine_dofs_per_cell));
2839  // for each coarse dof: find its position within the fine element
2840  // and set this value to one in the respective vector (all other
2841  // values are zero by construction)
2842  for (unsigned int local_coarse_dof=0;
2843  local_coarse_dof<coarse_dofs_per_cell_component;
2844  ++local_coarse_dof)
2845  for (unsigned int fine_dof=0; fine_dof<fine_fe.dofs_per_cell; ++fine_dof)
2846  if (fine_fe.system_to_component_index(fine_dof)
2847  ==
2848  std::make_pair (fine_component, local_coarse_dof))
2849  {
2850  parameter_dofs[local_coarse_dof](fine_dof) = 1.;
2851  break;
2852  };
2853 
2854 
2855  // find out how many DoFs there are on the grids belonging to the
2856  // components we want to match
2857  unsigned int n_parameters_on_fine_grid=0;
2858  if (true)
2859  {
2860  // have a flag for each dof on the fine grid and set it to true
2861  // if this is an interesting dof. finally count how many true's
2862  // there
2863  std::vector<bool> dof_is_interesting (fine_grid.n_dofs(), false);
2864  std::vector<types::global_dof_index> local_dof_indices (fine_fe.dofs_per_cell);
2865 
2866  for (typename ::DoFHandler<dim,spacedim>::active_cell_iterator
2867  cell=fine_grid.begin_active();
2868  cell!=fine_grid.end(); ++cell)
2869  if (cell->is_locally_owned ())
2870  {
2871  cell->get_dof_indices (local_dof_indices);
2872  for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
2873  if (fine_fe.system_to_component_index(i).first == fine_component)
2874  dof_is_interesting[local_dof_indices[i]] = true;
2875  };
2876 
2877  n_parameters_on_fine_grid = std::count (dof_is_interesting.begin(),
2878  dof_is_interesting.end(),
2879  true);
2880  };
2881 
2882 
2883  // set up the weights mapping
2884  weights.clear ();
2885  weights.resize (n_coarse_dofs);
2886 
2887  weight_mapping.clear ();
2888  weight_mapping.resize (n_fine_dofs, numbers::invalid_dof_index);
2889 
2890  if (true)
2891  {
2892  std::vector<types::global_dof_index> local_dof_indices(fine_fe.dofs_per_cell);
2893  unsigned int next_free_index=0;
2894  for (typename ::DoFHandler<dim,spacedim>::active_cell_iterator
2895  cell=fine_grid.begin_active();
2896  cell != fine_grid.end(); ++cell)
2897  if (cell->is_locally_owned ())
2898  {
2899  cell->get_dof_indices (local_dof_indices);
2900  for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
2901  // if this DoF is a parameter dof and has not yet been
2902  // numbered, then do so
2903  if ((fine_fe.system_to_component_index(i).first == fine_component) &&
2904  (weight_mapping[local_dof_indices[i]] == numbers::invalid_dof_index))
2905  {
2906  weight_mapping[local_dof_indices[i]] = next_free_index;
2907  ++next_free_index;
2908  };
2909  };
2910 
2911  Assert (next_free_index == n_parameters_on_fine_grid,
2912  ExcInternalError());
2913  };
2914 
2915 
2916  // for each cell on the parameter grid: find out which degrees of
2917  // freedom on the fine grid correspond in which way to the degrees
2918  // of freedom on the parameter grid
2919  //
2920  // do this in a separate function to allow for multithreading
2921  // there. see this function also if you want to read more
2922  // information on the algorithm used.
2923  compute_intergrid_weights_2 (coarse_grid, coarse_component,
2924  coarse_to_fine_grid_map, parameter_dofs,
2925  weight_mapping, weights);
2926 
2927 
2928  // ok, now we have all weights for each dof on the fine grid. if in
2929  // debug mode lets see if everything went smooth, i.e. each dof has
2930  // sum of weights one
2931  //
2932  // in other words this means that if the sum of all shape functions
2933  // on the parameter grid is one (which is always the case), then
2934  // the representation on the state grid should be as well (division
2935  // of unity)
2936  //
2937  // if the parameter grid has more than one component, then the
2938  // respective dofs of the other components have sum of weights
2939  // zero, of course. we do not explicitly ask which component a dof
2940  // belongs to, but this at least tests some errors
2941 #ifdef DEBUG
2942  for (unsigned int col=0; col<n_parameters_on_fine_grid; ++col)
2943  {
2944  double sum=0;
2945  for (types::global_dof_index row=0; row<n_coarse_dofs; ++row)
2946  if (weights[row].find(col) != weights[row].end())
2947  sum += weights[row][col];
2948  Assert ((std::fabs(sum-1) < 1.e-12) ||
2949  ((coarse_fe.n_components()>1) && (sum==0)), ExcInternalError());
2950  };
2951 #endif
2952 
2953 
2954  return n_parameters_on_fine_grid;
2955  }
2956 
2957 
2958  }
2959  }
2960 
2961 
2962 
2963  template <int dim, int spacedim>
2964  void
2966  const DoFHandler<dim,spacedim> &coarse_grid,
2967  const unsigned int coarse_component,
2968  const DoFHandler<dim,spacedim> &fine_grid,
2969  const unsigned int fine_component,
2970  const InterGridMap<DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
2971  ConstraintMatrix &constraints)
2972  {
2973  // store the weights with which a dof on the parameter grid contributes
2974  // to a dof on the fine grid. see the long doc below for more info
2975  //
2976  // allocate as many rows as there are parameter dofs on the coarse grid
2977  // and as many columns as there are parameter dofs on the fine grid.
2978  //
2979  // weight_mapping is used to map the global (fine grid) parameter dof
2980  // indices to the columns
2981  //
2982  // in the original implementation, the weights array was actually of
2983  // FullMatrix<double> type. this wasted huge amounts of memory, but was
2984  // fast. nonetheless, since the memory consumption was quadratic in the
2985  // number of degrees of freedom, this was not very practical, so we now
2986  // use a vector of rows of the matrix, and in each row a vector of
2987  // pairs (colnum,value). this seems like the best tradeoff between
2988  // memory and speed, as it is now linear in memory and still fast
2989  // enough.
2990  //
2991  // to save some memory and since the weights are usually (negative)
2992  // powers of 2, we choose the value type of the matrix to be @p{float}
2993  // rather than @p{double}.
2994  std::vector<std::map<types::global_dof_index, float> > weights;
2995 
2996  // this is this mapping. there is one entry for each dof on the fine
2997  // grid; if it is a parameter dof, then its value is the column in
2998  // weights for that parameter dof, if it is any other dof, then its
2999  // value is -1, indicating an error
3000  std::vector<types::global_dof_index> weight_mapping;
3001 
3002  const unsigned int n_parameters_on_fine_grid
3003  = internal::compute_intergrid_weights_1 (coarse_grid, coarse_component,
3004  fine_grid, fine_component,
3005  coarse_to_fine_grid_map,
3006  weights, weight_mapping);
3007  (void)n_parameters_on_fine_grid;
3008 
3009  // global numbers of dofs
3010  const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(),
3011  n_fine_dofs = fine_grid.n_dofs();
3012 
3013 
3014  // get an array in which we store which dof on the coarse grid is a
3015  // parameter and which is not
3016  std::vector<bool> coarse_dof_is_parameter (coarse_grid.n_dofs());
3017  if (true)
3018  {
3019  std::vector<bool> mask (coarse_grid.get_fe(0).n_components(),
3020  false);
3021  mask[coarse_component] = true;
3022  extract_dofs (coarse_grid, ComponentMask(mask), coarse_dof_is_parameter);
3023  }
3024 
3025  // now we know that the weights in each row constitute a constraint.
3026  // enter this into the constraints object
3027  //
3028  // first task: for each parameter dof on the parameter grid, find a
3029  // representant on the fine, global grid. this is possible since we use
3030  // conforming finite element. we take this representant to be the first
3031  // element in this row with weight identical to one. the representant
3032  // will become an unconstrained degree of freedom, while all others
3033  // will be constrained to this dof (and possibly others)
3034  std::vector<types::global_dof_index> representants(n_coarse_dofs, numbers::invalid_dof_index);
3035  for (types::global_dof_index parameter_dof=0; parameter_dof<n_coarse_dofs;
3036  ++parameter_dof)
3037  if (coarse_dof_is_parameter[parameter_dof] == true)
3038  {
3039  // if this is the line of a parameter dof on the coarse grid,
3040  // then it should have at least one dependent node on the fine
3041  // grid
3042  Assert (weights[parameter_dof].size() > 0, ExcInternalError());
3043 
3044  // find the column where the representant is mentioned
3045  std::map<types::global_dof_index,float>::const_iterator i = weights[parameter_dof].begin();
3046  for (; i!=weights[parameter_dof].end(); ++i)
3047  if (i->second == 1)
3048  break;
3049  Assert (i!=weights[parameter_dof].end(), ExcInternalError());
3050  const types::global_dof_index column = i->first;
3051 
3052  // now we know in which column of weights the representant is,
3053  // but we don't know its global index. get it using the inverse
3054  // operation of the weight_mapping
3055  types::global_dof_index global_dof=0;
3056  for (; global_dof<weight_mapping.size(); ++global_dof)
3057  if (weight_mapping[global_dof] == static_cast<types::global_dof_index>(column))
3058  break;
3059  Assert (global_dof < weight_mapping.size(), ExcInternalError());
3060 
3061  // now enter the representants global index into our list
3062  representants[parameter_dof] = global_dof;
3063  }
3064  else
3065  {
3066  // consistency check: if this is no parameter dof on the coarse
3067  // grid, then the respective row must be empty!
3068  Assert (weights[parameter_dof].size() == 0, ExcInternalError());
3069  };
3070 
3071 
3072 
3073  // note for people that want to optimize this function: the largest
3074  // part of the computing time is spent in the following, rather
3075  // innocent block of code. basically, it must be the
3076  // ConstraintMatrix::add_entry call which takes the bulk of the time,
3077  // but it is not known to the author how to make it faster...
3078  std::vector<std::pair<types::global_dof_index,double> > constraint_line;
3079  for (types::global_dof_index global_dof=0; global_dof<n_fine_dofs; ++global_dof)
3080  if (weight_mapping[global_dof] != numbers::invalid_dof_index)
3081  // this global dof is a parameter dof, so it may carry a constraint
3082  // note that for each global dof, the sum of weights shall be one,
3083  // so we can find out whether this dof is constrained in the
3084  // following way: if the only weight in this row is a one, and the
3085  // representant for the parameter dof of the line in which this one
3086  // is is the present dof, then we consider this dof to be
3087  // unconstrained. otherwise, all other dofs are constrained
3088  {
3089  const types::global_dof_index col = weight_mapping[global_dof];
3090  Assert (col < n_parameters_on_fine_grid, ExcInternalError());
3091 
3092  types::global_dof_index first_used_row=0;
3093 
3094  {
3095  Assert (weights.size() > 0, ExcInternalError());
3096  std::map<types::global_dof_index,float>::const_iterator
3097  col_entry = weights[0].end();
3098  for (; first_used_row<n_coarse_dofs; ++first_used_row)
3099  {
3100  col_entry = weights[first_used_row].find(col);
3101  if (col_entry != weights[first_used_row].end())
3102  break;
3103  }
3104 
3105  Assert (col_entry != weights[first_used_row].end(), ExcInternalError());
3106 
3107  if ((col_entry->second == 1) &&
3108  (representants[first_used_row] == global_dof))
3109  // dof unconstrained or constrained to itself (in case this
3110  // cell is mapped to itself, rather than to children of
3111  // itself)
3112  continue;
3113  }
3114 
3115 
3116  // otherwise enter all constraints
3117  constraints.add_line (global_dof);
3118 
3119  constraint_line.clear ();
3120  for (types::global_dof_index row=first_used_row; row<n_coarse_dofs; ++row)
3121  {
3122  const std::map<types::global_dof_index,float>::const_iterator
3123  j = weights[row].find(col);
3124  if ((j != weights[row].end()) && (j->second != 0))
3125  constraint_line.emplace_back (representants[row], j->second);
3126  };
3127 
3128  constraints.add_entries (global_dof, constraint_line);
3129  };
3130  }
3131 
3132 
3133 
3134  template <int dim, int spacedim>
3135  void
3137  const DoFHandler<dim,spacedim> &coarse_grid,
3138  const unsigned int coarse_component,
3139  const DoFHandler<dim,spacedim> &fine_grid,
3140  const unsigned int fine_component,
3141  const InterGridMap<DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
3142  std::vector<std::map<types::global_dof_index, float> > &transfer_representation)
3143  {
3144  // store the weights with which a dof on the parameter grid contributes
3145  // to a dof on the fine grid. see the long doc below for more info
3146  //
3147  // allocate as many rows as there are parameter dofs on the coarse grid
3148  // and as many columns as there are parameter dofs on the fine grid.
3149  //
3150  // weight_mapping is used to map the global (fine grid) parameter dof
3151  // indices to the columns
3152  //
3153  // in the original implementation, the weights array was actually of
3154  // FullMatrix<double> type. this wasted huge amounts of memory, but was
3155  // fast. nonetheless, since the memory consumption was quadratic in the
3156  // number of degrees of freedom, this was not very practical, so we now
3157  // use a vector of rows of the matrix, and in each row a vector of
3158  // pairs (colnum,value). this seems like the best tradeoff between
3159  // memory and speed, as it is now linear in memory and still fast
3160  // enough.
3161  //
3162  // to save some memory and since the weights are usually (negative)
3163  // powers of 2, we choose the value type of the matrix to be @p{float}
3164  // rather than @p{double}.
3165  std::vector<std::map<types::global_dof_index, float> > weights;
3166 
3167  // this is this mapping. there is one entry for each dof on the fine
3168  // grid; if it is a parameter dof, then its value is the column in
3169  // weights for that parameter dof, if it is any other dof, then its
3170  // value is -1, indicating an error
3171  std::vector<types::global_dof_index> weight_mapping;
3172 
3173  internal::compute_intergrid_weights_1 (coarse_grid, coarse_component,
3174  fine_grid, fine_component,
3175  coarse_to_fine_grid_map,
3176  weights, weight_mapping);
3177 
3178  // now compute the requested representation
3179  const types::global_dof_index n_global_parm_dofs
3180  = std::count_if (weight_mapping.begin(), weight_mapping.end(),
3181  std::bind (std::not_equal_to<types::global_dof_index>(), std::placeholders::_1, numbers::invalid_dof_index));
3182 
3183  // first construct the inverse mapping of weight_mapping
3184  std::vector<types::global_dof_index> inverse_weight_mapping (n_global_parm_dofs,
3186  for (types::global_dof_index i=0; i<weight_mapping.size(); ++i)
3187  {
3188  const types::global_dof_index parameter_dof = weight_mapping[i];
3189  // if this global dof is a parameter
3190  if (parameter_dof != numbers::invalid_dof_index)
3191  {
3192  Assert (parameter_dof < n_global_parm_dofs, ExcInternalError());
3193  Assert ((inverse_weight_mapping[parameter_dof] == numbers::invalid_dof_index),
3194  ExcInternalError());
3195 
3196  inverse_weight_mapping[parameter_dof] = i;
3197  };
3198  };
3199 
3200  // next copy over weights array and replace respective numbers
3201  const types::global_dof_index n_rows = weight_mapping.size();
3202 
3203  transfer_representation.clear ();
3204  transfer_representation.resize (n_rows);
3205 
3206  const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs();
3207  for (types::global_dof_index i=0; i<n_coarse_dofs; ++i)
3208  {
3209  std::map<types::global_dof_index, float>::const_iterator j = weights[i].begin();
3210  for (; j!=weights[i].end(); ++j)
3211  {
3212  const types::global_dof_index p = inverse_weight_mapping[j->first];
3213  Assert (p<n_rows, ExcInternalError());
3214 
3215  transfer_representation[p][i] = j->second;
3216  };
3217  };
3218  }
3219 
3220 
3221 
3222  template <int dim, int spacedim, template <int, int> class DoFHandlerType>
3223  void
3224  make_zero_boundary_constraints (const DoFHandlerType<dim, spacedim> &dof,
3225  const types::boundary_id boundary_id,
3226  ConstraintMatrix &zero_boundary_constraints,
3227  const ComponentMask &component_mask)
3228  {
3229  Assert (component_mask.represents_n_components(dof.get_fe(0).n_components()),
3230  ExcMessage ("The number of components in the mask has to be either "
3231  "zero or equal to the number of components in the finite "
3232  "element."));
3233 
3234  const unsigned int n_components = DoFTools::n_components (dof);
3235 
3236  Assert (component_mask.n_selected_components(n_components) > 0,
3238 
3239  // a field to store the indices on the face
3240  std::vector<types::global_dof_index> face_dofs;
3241  face_dofs.reserve (max_dofs_per_face(dof));
3242  // a field to store the indices on the cell
3243  std::vector<types::global_dof_index> cell_dofs;
3244  cell_dofs.reserve (max_dofs_per_cell(dof));
3245 
3246  typename DoFHandlerType<dim,spacedim>::active_cell_iterator
3247  cell = dof.begin_active(),
3248  endc = dof.end();
3249  for (; cell!=endc; ++cell)
3250  if (!cell->is_artificial()
3251  &&
3252  cell->at_boundary ())
3253  {
3254  const FiniteElement<dim,spacedim> &fe = cell->get_fe();
3255 
3256  // get global indices of dofs on the cell
3257  cell_dofs.resize (fe.dofs_per_cell);
3258  cell->get_dof_indices (cell_dofs);
3259 
3260  for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
3261  ++face_no)
3262  {
3263  const typename DoFHandlerType<dim,spacedim>::face_iterator face = cell->face(face_no);
3264 
3265  // if face is on the boundary and satisfies the correct
3266  // boundary id property
3267  if (face->at_boundary ()
3268  &&
3269  ((boundary_id == numbers::invalid_boundary_id)
3270  ||
3271  (face->boundary_id() == boundary_id)))
3272  {
3273  // get indices and physical location on this face
3274  face_dofs.resize (fe.dofs_per_face);
3275  face->get_dof_indices (face_dofs, cell->active_fe_index());
3276 
3277  // enter those dofs into the list that match the component
3278  // signature.
3279  for (unsigned int i=0; i<face_dofs.size(); ++i)
3280  {
3281  // Find out if a dof has a contribution in this
3282  // component, and if so, add it to the list
3283  const std::vector<types::global_dof_index>::iterator it_index_on_cell
3284  = std::find (cell_dofs.begin(), cell_dofs.end(), face_dofs[i]);
3285  Assert (it_index_on_cell != cell_dofs.end(), ExcInvalidIterator());
3286  const unsigned int index_on_cell = std::distance(cell_dofs.begin(),
3287  it_index_on_cell);
3288  const ComponentMask &nonzero_component_array
3289  = cell->get_fe().get_nonzero_components (index_on_cell);
3290  bool nonzero = false;
3291  for (unsigned int c=0; c<n_components; ++c)
3292  if (nonzero_component_array[c] && component_mask[c])
3293  {
3294  nonzero = true;
3295  break;
3296  }
3297 
3298  if (nonzero)
3299  zero_boundary_constraints.add_line (face_dofs[i]);
3300  }
3301  }
3302  }
3303  }
3304  }
3305 
3306 
3307 
3308  template <int dim, int spacedim, template <int, int> class DoFHandlerType>
3309  void
3310  make_zero_boundary_constraints (const DoFHandlerType<dim, spacedim> &dof,
3311  ConstraintMatrix &zero_boundary_constraints,
3312  const ComponentMask &component_mask)
3313  {
3315  zero_boundary_constraints, component_mask);
3316  }
3317 
3318 
3319 } // end of namespace DoFTools
3320 
3321 
3322 
3323 // explicit instantiations
3324 
3325 #include "dof_tools_constraints.inst"
3326 
3327 
3328 
3329 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
Definition: types.h:173
void compute_intergrid_transfer_representation(const DoFHandler< dim, spacedim > &coarse_grid, const unsigned int coarse_component, const DoFHandler< dim, spacedim > &fine_grid, const unsigned int fine_component, const InterGridMap< DoFHandler< dim, spacedim > > &coarse_to_fine_grid_map, std::vector< std::map< types::global_dof_index, float > > &transfer_representation)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1175
static::ExceptionBase & ExcGridsDontMatch()
bool is_constrained(const size_type index) const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const std::vector< std::pair< size_type, double > > * get_constraint_entries(const size_type line) const
unsigned int boundary_id
Definition: types.h:110
#define AssertIndexRange(index, range)
Definition: exceptions.h:1199
const unsigned int dofs_per_quad
Definition: fe_base.h:250
void add_entries(const size_type line, const std::vector< std::pair< size_type, double > > &col_val_pairs)
static::ExceptionBase & ExcFiniteElementsDontMatch()
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:398
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const unsigned int dofs_per_line
Definition: fe_base.h:244
bool represents_n_components(const unsigned int n) const
void compute_intergrid_constraints(const DoFHandler< dim, spacedim > &coarse_grid, const unsigned int coarse_component, const DoFHandler< dim, spacedim > &fine_grid, const unsigned int fine_component, const InterGridMap< DoFHandler< dim, spacedim > > &coarse_to_fine_grid_map, ConstraintMatrix &constraints)
static::ExceptionBase & ExcInvalidIterator()
void add_line(const size_type line)
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, ConstraintMatrix &constraints)
static::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
static::ExceptionBase & ExcGridNotCoarser()
void add_entry(const size_type line, const size_type column, const double value)
void invert(const FullMatrix< number2 > &M)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:88
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:337
types::global_dof_index n_dofs() const
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void extract_locally_owned_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1032
const types::boundary_id invalid_boundary_id
Definition: types.h:204
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
Definition: fe.h:3030
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:4326
static::ExceptionBase & ExcNoComponentSelected()
const FullMatrix< double > & constraints(const ::internal::SubfaceCase< dim > &subface_case=::internal::SubfaceCase< dim >::case_isotropic) const
Definition: fe.cc:793
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1080
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
unsigned int n_components() const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index) const
Definition: fe.h:2966
size_type m() const
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const
Definition: fe.cc:863
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const
Definition: fe.cc:880
void make_periodicity_constraints(const FaceIterator &face_1, const typename identity< FaceIterator >::type &face_2,::ConstraintMatrix &constraint_matrix, const ComponentMask &component_mask=ComponentMask(), const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false, const FullMatrix< double > &matrix=FullMatrix< double >(), const std::vector< unsigned int > &first_vector_components=std::vector< unsigned int >())
const unsigned int dofs_per_cell
Definition: fe_base.h:295
void extract_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< bool > &selected_dofs)
Definition: dof_tools.cc:372
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:2910
const unsigned int dofs_per_face
Definition: fe_base.h:288
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1106
static::ExceptionBase & ExcNotImplemented()
const unsigned int dofs_per_vertex
Definition: fe_base.h:238
Definition: table.h:33
const std::vector< Point< dim-1 > > & get_unit_face_support_points() const
Definition: fe.cc:1011
void set_inhomogeneity(const size_type line, const double value)
void make_zero_boundary_constraints(const DoFHandlerType< dim, spacedim > &dof, const types::boundary_id boundary_id, ConstraintMatrix &zero_boundary_constraints, const ComponentMask &component_mask=ComponentMask())
const types::global_dof_index invalid_dof_index
Definition: types.h:187
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
Definition: fe.cc:1176
T max(const T &t, const MPI_Comm &mpi_communicator)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static::ExceptionBase & ExcInternalError()