Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_handler.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
20 
21 #include <deal.II/distributed/cell_data_transfer.templates.h>
24 
27 
28 #include <deal.II/fe/fe.h>
29 
31 #include <deal.II/grid/tria.h>
35 
36 #include <deal.II/hp/dof_faces.h>
37 #include <deal.II/hp/dof_handler.h>
38 #include <deal.II/hp/dof_level.h>
39 
40 #include <boost/serialization/array.hpp>
41 
42 #include <algorithm>
43 #include <functional>
44 #include <set>
45 #include <unordered_set>
46 
48 
49 // The following is necessary for compilation under Visual Studio which is
50 // unable to correctly distinguish between ::DoFHandler and
51 // ::hp::DoFHandler. Plus it makes code in dof_handler.cc easier to read.
52 #if defined(_MSC_VER) && (_MSC_VER >= 1800)
53 template <int dim, int spacedim>
55 #else
56 // When using older Visual Studio or a different compiler just fall back.
57 # define HpDoFHandler DoFHandler
58 #endif
59 
60 namespace parallel
61 {
62  namespace distributed
63  {
64  template <int, int>
65  class Triangulation;
66  }
67 } // namespace parallel
68 
69 
70 
71 namespace internal
72 {
73  namespace hp
74  {
75  namespace DoFHandlerImplementation
76  {
77  // access class ::hp::DoFHandler instead of namespace
78  // internal::hp::DoFHandler, etc
79  using ::hp::DoFHandler;
80 
86  {
92  template <int dim, int spacedim>
93  static void
95  DoFHandler<dim, spacedim> &dof_handler)
96  {
97  (void)dof_handler;
98  for (const auto &cell : dof_handler.active_cell_iterators())
99  if (cell->is_locally_owned())
100  Assert(
101  !cell->future_fe_index_set(),
102  ExcMessage(
103  "There shouldn't be any cells flagged for p-adaptation when partitioning."));
104  }
105 
106 
107 
112  template <int dim, int spacedim>
113  static void
115  {
116  // Release all space except the fields for active_fe_indices and
117  // refinement flags which we have to back up before
118  {
119  std::vector<std::vector<DoFLevel::active_fe_index_type>>
120  active_fe_backup(dof_handler.levels.size()),
121  future_fe_backup(dof_handler.levels.size());
122  for (unsigned int level = 0; level < dof_handler.levels.size();
123  ++level)
124  {
125  active_fe_backup[level] =
126  std::move(dof_handler.levels[level]->active_fe_indices);
127  future_fe_backup[level] =
128  std::move(dof_handler.levels[level]->future_fe_indices);
129  }
130 
131  // delete all levels and set them up newly, since vectors
132  // are troublesome if you want to change their size
133  dof_handler.clear_space();
134 
135  for (unsigned int level = 0; level < dof_handler.tria->n_levels();
136  ++level)
137  {
138  dof_handler.levels.emplace_back(new internal::hp::DoFLevel);
139  // recover backups
140  dof_handler.levels[level]->active_fe_indices =
141  std::move(active_fe_backup[level]);
142  dof_handler.levels[level]->future_fe_indices =
143  std::move(future_fe_backup[level]);
144  }
145 
146  if (dim > 1)
147  dof_handler.faces =
148  std_cxx14::make_unique<internal::hp::DoFIndicesOnFaces<dim>>();
149  }
150  }
151 
152 
153 
158  template <int dim, int spacedim>
159  static void
161  {
162  // The final step in all of the reserve_space() functions is to set
163  // up vertex dof information. since vertices are sequentially
164  // numbered, what we do first is to set up an array in which
165  // we record whether a vertex is associated with any of the
166  // given fe's, by setting a bit. in a later step, we then
167  // actually allocate memory for the required dofs
168  //
169  // in the following, we only need to consider vertices that are
170  // adjacent to either a locally owned or a ghost cell; we never
171  // store anything on vertices that are only surrounded by
172  // artificial cells. so figure out that subset of vertices
173  // first
174  std::vector<bool> locally_used_vertices(
175  dof_handler.tria->n_vertices(), false);
176  for (const auto &cell : dof_handler.active_cell_iterators())
177  if (!cell->is_artificial())
178  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
179  locally_used_vertices[cell->vertex_index(v)] = true;
180 
181  std::vector<std::vector<bool>> vertex_fe_association(
182  dof_handler.fe_collection.size(),
183  std::vector<bool>(dof_handler.tria->n_vertices(), false));
184 
185  for (const auto &cell : dof_handler.active_cell_iterators())
186  if (!cell->is_artificial())
187  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
188  vertex_fe_association[cell->active_fe_index()]
189  [cell->vertex_index(v)] = true;
190 
191  // in debug mode, make sure that each vertex is associated
192  // with at least one fe (note that except for unused
193  // vertices, all vertices are actually active). this is of
194  // course only true for vertices that are part of either
195  // ghost or locally owned cells
196 #ifdef DEBUG
197  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
198  if (locally_used_vertices[v] == true)
199  if (dof_handler.tria->vertex_used(v) == true)
200  {
201  unsigned int fe = 0;
202  for (; fe < dof_handler.fe_collection.size(); ++fe)
203  if (vertex_fe_association[fe][v] == true)
204  break;
205  Assert(fe != dof_handler.fe_collection.size(),
206  ExcInternalError());
207  }
208 #endif
209 
210  // next count how much memory we actually need. for each
211  // vertex, we need one slot per fe to store the fe_index,
212  // plus dofs_per_vertex for this fe. in addition, we need
213  // one slot as the end marker for the fe_indices. at the
214  // same time already fill the vertex_dof_offsets field
215  dof_handler.vertex_dof_offsets.resize(dof_handler.tria->n_vertices(),
217 
218  unsigned int vertex_slots_needed = 0;
219  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
220  if (dof_handler.tria->vertex_used(v) == true)
221  if (locally_used_vertices[v] == true)
222  {
223  dof_handler.vertex_dof_offsets[v] = vertex_slots_needed;
224 
225  for (unsigned int fe = 0;
226  fe < dof_handler.fe_collection.size();
227  ++fe)
228  if (vertex_fe_association[fe][v] == true)
229  vertex_slots_needed +=
230  dof_handler.get_fe(fe).dofs_per_vertex + 1;
231 
232  // don't forget the end_marker:
233  ++vertex_slots_needed;
234  }
235 
236  // now allocate the space we have determined we need, and
237  // set up the linked lists for each of the vertices
238  dof_handler.vertex_dofs.resize(vertex_slots_needed,
240  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
241  if (dof_handler.tria->vertex_used(v) == true)
242  if (locally_used_vertices[v] == true)
243  {
244  unsigned int current_index =
245  dof_handler.vertex_dof_offsets[v];
246  for (unsigned int fe = 0;
247  fe < dof_handler.fe_collection.size();
248  ++fe)
249  if (vertex_fe_association[fe][v] == true)
250  {
251  // if this vertex uses this fe, then set the
252  // fe_index and move the pointer ahead
253  dof_handler.vertex_dofs[current_index] = fe;
254  current_index +=
255  dof_handler.get_fe(fe).dofs_per_vertex + 1;
256  }
257  // finally place the end marker
258  dof_handler.vertex_dofs[current_index] =
260  }
261  }
262 
263 
264 
269  template <int dim, int spacedim>
270  static void
272  {
273  // count how much space we need on each level for the cell
274  // dofs and set the dof_*_offsets data. initially set the
275  // latter to an invalid index, and only later set it to
276  // something reasonable for active dof_handler.cells
277  //
278  // note that for dof_handler.cells, the situation is simpler
279  // than for other (lower dimensional) objects since exactly
280  // one finite element is used for it
281  for (unsigned int level = 0; level < dof_handler.tria->n_levels();
282  ++level)
283  {
284  dof_handler.levels[level]->dof_offsets =
285  std::vector<DoFLevel::offset_type>(
286  dof_handler.tria->n_raw_cells(level),
287  static_cast<DoFLevel::offset_type>(-1));
288  dof_handler.levels[level]->cell_cache_offsets =
289  std::vector<DoFLevel::offset_type>(
290  dof_handler.tria->n_raw_cells(level),
291  static_cast<DoFLevel::offset_type>(-1));
292 
293  types::global_dof_index next_free_dof = 0;
294  types::global_dof_index cache_size = 0;
295  typename HpDoFHandler<dim, spacedim>::active_cell_iterator
296  cell = dof_handler.begin_active(level),
297  endc = dof_handler.end_active(level);
298  for (; cell != endc; ++cell)
299  if (cell->is_active() && !cell->is_artificial())
300  {
301  dof_handler.levels[level]->dof_offsets[cell->index()] =
302  next_free_dof;
303  next_free_dof +=
304  cell->get_fe().template n_dofs_per_object<dim>();
305 
306  dof_handler.levels[level]
307  ->cell_cache_offsets[cell->index()] = cache_size;
308  cache_size += cell->get_fe().dofs_per_cell;
309  }
310 
311  dof_handler.levels[level]->dof_indices =
312  std::vector<types::global_dof_index>(
313  next_free_dof, numbers::invalid_dof_index);
314  dof_handler.levels[level]->cell_dof_indices_cache =
315  std::vector<types::global_dof_index>(
316  cache_size, numbers::invalid_dof_index);
317  }
318 
319  // safety check: make sure that the number of DoFs we
320  // allocated is actually correct (above we have also set the
321  // dof_*_offsets field, so we couldn't use this simpler
322  // algorithm)
323 #ifdef DEBUG
324  for (unsigned int level = 0; level < dof_handler.tria->n_levels();
325  ++level)
326  {
327  types::global_dof_index counter = 0;
328  for (const auto &cell :
330  if (cell->is_active() && !cell->is_artificial())
331  counter += cell->get_fe().template n_dofs_per_object<dim>();
332 
333  Assert(dof_handler.levels[level]->dof_indices.size() == counter,
334  ExcInternalError());
335 
336  // also check that the number of unassigned slots in the
337  // dof_offsets equals the number of cells on that level minus the
338  // number of active, non-artificial cells (because these are
339  // exactly the cells on which we do something)
340  unsigned int n_active_non_artificial_cells = 0;
341  for (const auto &cell :
343  if (cell->is_active() && !cell->is_artificial())
344  ++n_active_non_artificial_cells;
345 
346  Assert(static_cast<unsigned int>(std::count(
347  dof_handler.levels[level]->dof_offsets.begin(),
348  dof_handler.levels[level]->dof_offsets.end(),
349  static_cast<DoFLevel::offset_type>(-1))) ==
350  dof_handler.tria->n_raw_cells(level) -
351  n_active_non_artificial_cells,
352  ExcInternalError());
353  }
354 #endif
355  }
356 
357 
358 
363  template <int dim, int spacedim>
364  static void
366  {
367  // make the code generic between lines and quads
368  std::vector<unsigned int> &face_dof_offsets =
369  (dim == 2 ?
370  reinterpret_cast<::internal::hp::DoFIndicesOnFaces<2> &>(
371  *dof_handler.faces)
372  .lines.dof_offsets :
373  reinterpret_cast<::internal::hp::DoFIndicesOnFaces<3> &>(
374  *dof_handler.faces)
375  .quads.dof_offsets);
376 
377  std::vector<types::global_dof_index> &face_dof_indices =
378  (dim == 2 ?
379  reinterpret_cast<::internal::hp::DoFIndicesOnFaces<2> &>(
380  *dof_handler.faces)
381  .lines.dofs :
382  reinterpret_cast<::internal::hp::DoFIndicesOnFaces<3> &>(
383  *dof_handler.faces)
384  .quads.dofs);
385 
386  // FACE DOFS
387  //
388  // Count face dofs, then allocate as much space
389  // as we need and prime the linked list for faces (see the
390  // description in hp::DoFLevel) with the indices we will
391  // need. Note that our task is more complicated than for the
392  // cell case above since two adjacent cells may have different
393  // active_fe_indices, in which case we need to allocate
394  // *two* sets of face dofs for the same face. But they don't
395  // *have* to be different, and so we need to prepare for this
396  // as well.
397  //
398  // The way we do things is that we loop over all active
399  // cells (these are the only ones that have DoFs
400  // anyway) and all their faces. We note in the
401  // user flags whether we have previously visited a face and
402  // if so skip it (consequently, we have to save and later
403  // restore the face flags)
404  {
405  std::vector<bool> saved_face_user_flags;
406  switch (dim)
407  {
408  case 2:
409  {
410  const_cast<::Triangulation<dim, spacedim> &>(
411  *dof_handler.tria)
412  .save_user_flags_line(saved_face_user_flags);
413  const_cast<::Triangulation<dim, spacedim> &>(
414  *dof_handler.tria)
415  .clear_user_flags_line();
416 
417  break;
418  }
419 
420  case 3:
421  {
422  const_cast<::Triangulation<dim, spacedim> &>(
423  *dof_handler.tria)
424  .save_user_flags_quad(saved_face_user_flags);
425  const_cast<::Triangulation<dim, spacedim> &>(
426  *dof_handler.tria)
427  .clear_user_flags_quad();
428 
429  break;
430  }
431 
432  default:
433  Assert(false, ExcNotImplemented());
434  }
435 
436  // An array to hold how many slots (see the hp::DoFLevel
437  // class) we will have to store on each level
438  unsigned int n_face_slots = 0;
439 
440  for (const auto &cell : dof_handler.active_cell_iterators())
441  if (!cell->is_artificial())
442  for (const unsigned int face :
444  if (cell->face(face)->user_flag_set() == false)
445  {
446  // Ok, face has not been visited. So we need to
447  // allocate space for it. Let's see how much we
448  // need: we need one set if a) there is no
449  // neighbor behind this face, or b) the neighbor
450  // is either coarser or finer than we are, or c)
451  // the neighbor is artificial, or d) the neighbor
452  // is neither coarser nor finer, nor is artificial,
453  // and just so happens to have the same active_fe_index :
454  if (cell->at_boundary(face) ||
455  cell->face(face)->has_children() ||
456  cell->neighbor_is_coarser(face) ||
457  (!cell->at_boundary(face) &&
458  cell->neighbor(face)->is_artificial()) ||
459  (!cell->at_boundary(face) &&
460  !cell->neighbor(face)->is_artificial() &&
461  (cell->active_fe_index() ==
462  cell->neighbor(face)->active_fe_index())))
463  // Ok, one set of dofs. that makes one active_fe_index,
464  // 1 times dofs_per_face dofs, and one stop index
465  n_face_slots +=
466  1 +
467  dof_handler.get_fe(cell->active_fe_index())
468  .template n_dofs_per_object<dim - 1>() +
469  1;
470 
471  // Otherwise we do indeed need two sets, i.e. two
472  // active_fe_indices, two sets of dofs, and one stop
473  // index:
474  else
475  n_face_slots +=
476  (1 + // the active_fe_index
477  dof_handler.get_fe(cell->active_fe_index())
478  .template n_dofs_per_object<
479  dim - 1>() // actual DoF indices
480  + 1 // the second active_fe_index
481  + dof_handler
482  .get_fe(cell->neighbor(face)->active_fe_index())
483  .template n_dofs_per_object<
484  dim - 1>() // actual DoF indices
485  + 1); // stop marker
486 
487  // mark this face as visited
488  cell->face(face)->set_user_flag();
489  }
490 
491  // Now that we know how many face dofs we will have to
492  // have, allocate the memory. Note that we
493  // allocate offsets for all faces, though only the active
494  // ones will have a non-invalid value later on
495  face_dof_offsets =
496  std::vector<unsigned int>(dof_handler.tria->n_raw_faces(),
498  face_dof_indices =
499  std::vector<types::global_dof_index>(n_face_slots,
501 
502  // With the memory now allocated, loop over the
503  // dof_handler cells again and prime the _offset values as
504  // well as the fe_index fields
505  switch (dim)
506  {
507  case 2:
508  {
509  const_cast<::Triangulation<dim, spacedim> &>(
510  *dof_handler.tria)
511  .clear_user_flags_line();
512 
513  break;
514  }
515 
516  case 3:
517  {
518  const_cast<::Triangulation<dim, spacedim> &>(
519  *dof_handler.tria)
520  .clear_user_flags_quad();
521 
522  break;
523  }
524 
525  default:
526  Assert(false, ExcNotImplemented());
527  }
528 
529  unsigned int next_free_face_slot = 0;
530 
531  for (const auto &cell : dof_handler.active_cell_iterators())
532  if (!cell->is_artificial())
533  for (const unsigned int face :
535  if (!cell->face(face)->user_flag_set())
536  {
537  // Same decision tree as before
538  if (cell->at_boundary(face) ||
539  cell->face(face)->has_children() ||
540  cell->neighbor_is_coarser(face) ||
541  (!cell->at_boundary(face) &&
542  cell->neighbor(face)->is_artificial()) ||
543  (!cell->at_boundary(face) &&
544  !cell->neighbor(face)->is_artificial() &&
545  (cell->active_fe_index() ==
546  cell->neighbor(face)->active_fe_index())))
547  {
548  // Only one active_fe_index lives on this
549  // face
550  face_dof_offsets[cell->face(face)->index()] =
551  next_free_face_slot;
552 
553  // Set the first and only slot for this face to
554  // active_fe_index of this face
555  face_dof_indices[next_free_face_slot] =
556  cell->active_fe_index();
557 
558  // The next dofs_per_face indices remain unset
559  // for the moment (i.e. at invalid_dof_index).
560  // Following this comes the stop index, which
561  // also is invalid_dof_index and therefore
562  // does not have to be explicitly set
563 
564  // Finally, move the current marker forward:
565  next_free_face_slot +=
566  1 // the active_fe_index
567  + dof_handler.get_fe(cell->active_fe_index())
568  .template n_dofs_per_object<
569  dim - 1>() // actual DoF indices
570  + 1; // the end marker
571  }
572  else
573  {
574  // There are two active_fe_indices that live on this
575  // face.
576  face_dof_offsets[cell->face(face)->index()] =
577  next_free_face_slot;
578 
579  // Store the two indices we will have to deal with.
580  // We sort these so that it does not matter which of
581  // the cells adjacent to this face we visit first. In
582  // sequential computations, this does not matter
583  // because the order in which we visit these cells is
584  // deterministic and always the same. But in parallel
585  // computations, we can get into trouble because two
586  // processors visit the cells in different order
587  // (because the mesh creation history on the two
588  // processors is different), and in that case it can
589  // happen that the order of active_fe_index values for
590  // a given face is different on the two processes,
591  // even though they agree on which two values need to
592  // be stored. Since the DoF unification on faces
593  // takes into account the order of the
594  // active_fe_indices, this leads to quite subtle bugs.
595  // We could fix this in the place where we do the DoF
596  // unification on cells, but it is better to just make
597  // sure that every process stores the exact same
598  // information (and in the same order) on each face.
599  unsigned int active_fe_indices[2] = {
600  cell->active_fe_index(),
601  cell->neighbor(face)->active_fe_index()};
602  if (active_fe_indices[1] < active_fe_indices[0])
603  std::swap(active_fe_indices[0],
604  active_fe_indices[1]);
605 
606  // Set first slot for this face to
607  // active_fe_index of this face
608  face_dof_indices[next_free_face_slot] =
609  active_fe_indices[0];
610 
611  // The next dofs_per_line/quad indices remain unset
612  // for the moment (i.e. at invalid_dof_index).
613  //
614  // Then comes the fe_index for the second slot:
615  face_dof_indices
616  [next_free_face_slot + 1 +
617  dof_handler.get_fe(active_fe_indices[0])
618  .template n_dofs_per_object<dim - 1>()] =
619  active_fe_indices[1];
620  // Then again a set of dofs that we need not
621  // set right now
622  //
623  // Following this comes the stop index, which
624  // also is invalid_dof_index and therefore
625  // does not have to be explicitly set
626 
627  // Finally, move the current marker forward:
628  next_free_face_slot +=
629  (1 // first active_fe_index
630  + dof_handler.get_fe(active_fe_indices[0])
631  .template n_dofs_per_object<
632  dim - 1>() // actual DoF indices
633  + 1 // second active_fe_index
634  + dof_handler.get_fe(active_fe_indices[1])
635  .template n_dofs_per_object<
636  dim - 1>() // actual DoF indices
637  + 1); // end marker
638  }
639 
640  // mark this face as visited
641  cell->face(face)->set_user_flag();
642  }
643 
644  // we should have moved the cursor for each level to the
645  // total number of dofs on that level. check that
646  Assert(next_free_face_slot == n_face_slots, ExcInternalError());
647 
648  // at the end, restore the user flags for the faces
649  switch (dim)
650  {
651  case 2:
652  {
653  const_cast<::Triangulation<dim, spacedim> &>(
654  *dof_handler.tria)
655  .load_user_flags_line(saved_face_user_flags);
656 
657  break;
658  }
659 
660  case 3:
661  {
662  const_cast<::Triangulation<dim, spacedim> &>(
663  *dof_handler.tria)
664  .load_user_flags_quad(saved_face_user_flags);
665 
666  break;
667  }
668 
669  default:
670  Assert(false, ExcNotImplemented());
671  }
672  }
673  }
674 
675 
676 
683  template <int spacedim>
684  static void reserve_space(DoFHandler<1, spacedim> &dof_handler)
685  {
686  Assert(dof_handler.fe_collection.size() > 0,
688  Assert(dof_handler.tria->n_levels() > 0,
689  ExcMessage("The current Triangulation must not be empty."));
690  Assert(dof_handler.tria->n_levels() == dof_handler.levels.size(),
691  ExcInternalError());
692 
693  reserve_space_release_space(dof_handler);
694 
695  Threads::TaskGroup<> tasks;
696  tasks +=
697  Threads::new_task(&reserve_space_cells<1, spacedim>, dof_handler);
698  tasks += Threads::new_task(&reserve_space_vertices<1, spacedim>,
699  dof_handler);
700  tasks.join_all();
701  }
702 
703 
704 
705  template <int spacedim>
706  static void reserve_space(DoFHandler<2, spacedim> &dof_handler)
707  {
708  Assert(dof_handler.fe_collection.size() > 0,
710  Assert(dof_handler.tria->n_levels() > 0,
711  ExcMessage("The current Triangulation must not be empty."));
712  Assert(dof_handler.tria->n_levels() == dof_handler.levels.size(),
713  ExcInternalError());
714 
715  reserve_space_release_space(dof_handler);
716 
717  Threads::TaskGroup<> tasks;
718  tasks +=
719  Threads::new_task(&reserve_space_cells<2, spacedim>, dof_handler);
720  tasks +=
721  Threads::new_task(&reserve_space_faces<2, spacedim>, dof_handler);
722  tasks += Threads::new_task(&reserve_space_vertices<2, spacedim>,
723  dof_handler);
724  tasks.join_all();
725  }
726 
727 
728 
729  template <int spacedim>
730  static void reserve_space(DoFHandler<3, spacedim> &dof_handler)
731  {
732  const unsigned int dim = 3;
733 
734  Assert(dof_handler.fe_collection.size() > 0,
736  Assert(dof_handler.tria->n_levels() > 0,
737  ExcMessage("The current Triangulation must not be empty."));
738  Assert(dof_handler.tria->n_levels() == dof_handler.levels.size(),
739  ExcInternalError());
740 
741  reserve_space_release_space(dof_handler);
742 
743  Threads::TaskGroup<> tasks;
744  tasks +=
745  Threads::new_task(&reserve_space_cells<3, spacedim>, dof_handler);
746  tasks +=
747  Threads::new_task(&reserve_space_faces<3, spacedim>, dof_handler);
748  tasks += Threads::new_task(&reserve_space_vertices<3, spacedim>,
749  dof_handler);
750 
751  // While the tasks above are running, we can turn to line dofs
752 
753  // the situation here is pretty much like with vertices:
754  // there can be an arbitrary number of finite elements
755  // associated with each line.
756  //
757  // the algorithm we use is somewhat similar to what we do in
758  // reserve_space_vertices()
759  {
760  // what we do first is to set up an array in which we
761  // record whether a line is associated with any of the
762  // given fe's, by setting a bit. in a later step, we
763  // then actually allocate memory for the required dofs
764  std::vector<std::vector<bool>> line_fe_association(
765  dof_handler.fe_collection.size(),
766  std::vector<bool>(dof_handler.tria->n_raw_lines(), false));
767 
768  for (const auto &cell : dof_handler.active_cell_iterators())
769  if (!cell->is_artificial())
770  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
771  ++l)
772  line_fe_association[cell->active_fe_index()]
773  [cell->line_index(l)] = true;
774 
775  // first check which of the lines is used at all,
776  // i.e. is associated with a finite element. we do this
777  // since not all lines may actually be used, in which
778  // case we do not have to allocate any memory at all
779  std::vector<bool> line_is_used(dof_handler.tria->n_raw_lines(),
780  false);
781  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
782  ++line)
783  for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
784  ++fe)
785  if (line_fe_association[fe][line] == true)
786  {
787  line_is_used[line] = true;
788  break;
789  }
790 
791  // next count how much memory we actually need. for each
792  // line, we need one slot per fe to store the fe_index,
793  // plus dofs_per_line for this fe. in addition, we need
794  // one slot as the end marker for the fe_indices. at the
795  // same time already fill the line_dofs_offsets field
796  dof_handler.faces->lines.dof_offsets.resize(
797  dof_handler.tria->n_raw_lines(), numbers::invalid_unsigned_int);
798 
799  unsigned int line_slots_needed = 0;
800  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
801  ++line)
802  if (line_is_used[line] == true)
803  {
804  dof_handler.faces->lines.dof_offsets[line] =
805  line_slots_needed;
806 
807  for (unsigned int fe = 0;
808  fe < dof_handler.fe_collection.size();
809  ++fe)
810  if (line_fe_association[fe][line] == true)
811  line_slots_needed +=
812  dof_handler.get_fe(fe).dofs_per_line + 1;
813  ++line_slots_needed;
814  }
815 
816  // now allocate the space we have determined we need,
817  // and set up the linked lists for each of the lines
818  dof_handler.faces->lines.dofs.resize(line_slots_needed,
820  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
821  ++line)
822  if (line_is_used[line] == true)
823  {
824  unsigned int pointer =
825  dof_handler.faces->lines.dof_offsets[line];
826  for (unsigned int fe = 0;
827  fe < dof_handler.fe_collection.size();
828  ++fe)
829  if (line_fe_association[fe][line] == true)
830  {
831  // if this line uses this fe, then set the
832  // fe_index and move the pointer ahead
833  dof_handler.faces->lines.dofs[pointer] = fe;
834  pointer += dof_handler.get_fe(fe).dofs_per_line + 1;
835  }
836  // finally place the end marker
837  dof_handler.faces->lines.dofs[pointer] =
839  }
840  }
841 
842  // Ensure that everything is done at this point.
843  tasks.join_all();
844  }
845 
846 
847 
851  template <int spacedim>
852  static unsigned int
854  {
855  return std::min(static_cast<types::global_dof_index>(
856  3 *
857  dof_handler.fe_collection.max_dofs_per_vertex() +
858  2 * dof_handler.fe_collection.max_dofs_per_line()),
859  dof_handler.n_dofs());
860  }
861 
862 
863 
864  template <int spacedim>
865  static unsigned int
867  {
868  // get these numbers by drawing pictures and counting...
869  // example:
870  // | | |
871  // --x-----x--x--X--
872  // | | | |
873  // | x--x--x
874  // | | | |
875  // --x--x--*--x--x--
876  // | | | |
877  // x--x--x |
878  // | | | |
879  // --X--x--x-----x--
880  // | | |
881  // x = vertices connected with center vertex *;
882  // = total of 19
883  //
884  // (the X vertices are connected with * if the vertices
885  // adjacent to X are hanging nodes) count lines -> 28 (don't
886  // forget to count mother and children separately!)
887  types::global_dof_index max_couplings;
888  switch (dof_handler.tria->max_adjacent_cells())
889  {
890  case 4:
891  max_couplings =
892  19 * dof_handler.fe_collection.max_dofs_per_vertex() +
893  28 * dof_handler.fe_collection.max_dofs_per_line() +
894  8 * dof_handler.fe_collection.max_dofs_per_quad();
895  break;
896  case 5:
897  max_couplings =
898  21 * dof_handler.fe_collection.max_dofs_per_vertex() +
899  31 * dof_handler.fe_collection.max_dofs_per_line() +
900  9 * dof_handler.fe_collection.max_dofs_per_quad();
901  break;
902  case 6:
903  max_couplings =
904  28 * dof_handler.fe_collection.max_dofs_per_vertex() +
905  42 * dof_handler.fe_collection.max_dofs_per_line() +
906  12 * dof_handler.fe_collection.max_dofs_per_quad();
907  break;
908  case 7:
909  max_couplings =
910  30 * dof_handler.fe_collection.max_dofs_per_vertex() +
911  45 * dof_handler.fe_collection.max_dofs_per_line() +
912  13 * dof_handler.fe_collection.max_dofs_per_quad();
913  break;
914  case 8:
915  max_couplings =
916  37 * dof_handler.fe_collection.max_dofs_per_vertex() +
917  56 * dof_handler.fe_collection.max_dofs_per_line() +
918  16 * dof_handler.fe_collection.max_dofs_per_quad();
919  break;
920  default:
921  Assert(false, ExcNotImplemented());
922  max_couplings = 0;
923  }
924  return std::min(max_couplings, dof_handler.n_dofs());
925  }
926 
927 
928 
929  template <int spacedim>
930  static unsigned int
932  {
933  // TODO:[?] Invent significantly better estimates than the ones in
934  // this function
935  // doing the same thing here is a rather complicated thing,
936  // compared to the 2d case, since it is hard to draw
937  // pictures with several refined hexahedra :-) so I
938  // presently only give a coarse estimate for the case that
939  // at most 8 hexes meet at each vertex
940  //
941  // can anyone give better estimate here?
942  const unsigned int max_adjacent_cells =
943  dof_handler.tria->max_adjacent_cells();
944 
945  types::global_dof_index max_couplings;
946  if (max_adjacent_cells <= 8)
947  max_couplings =
948  7 * 7 * 7 * dof_handler.fe_collection.max_dofs_per_vertex() +
949  7 * 6 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_line() +
950  9 * 4 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_quad() +
951  27 * dof_handler.fe_collection.max_dofs_per_hex();
952  else
953  {
954  Assert(false, ExcNotImplemented());
955  max_couplings = 0;
956  }
957 
958  return std::min(max_couplings, dof_handler.n_dofs());
959  }
960 
961 
962 
974  template <int dim, int spacedim>
975  static void
977  ::hp::DoFHandler<dim, spacedim> &dof_handler)
978  {
979  if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
980  dynamic_cast<
981  const ::parallel::shared::Triangulation<dim, spacedim>
982  *>(&dof_handler.get_triangulation()))
983  {
984  // we have a shared triangulation. in this case, every processor
985  // knows about all cells, but every processor only has knowledge
986  // about the active_fe_index on the cells it owns.
987  //
988  // we can create a complete set of active_fe_indices by letting
989  // every processor create a vector of indices for all cells,
990  // filling only those on the cells it owns and setting the indices
991  // on the other cells to zero. then we add all of these vectors
992  // up, and because every vector entry has exactly one processor
993  // that owns it, the sum is correct
994  std::vector<unsigned int> active_fe_indices(tr->n_active_cells(),
995  0u);
996  for (const auto &cell : dof_handler.active_cell_iterators())
997  if (cell->is_locally_owned())
998  active_fe_indices[cell->active_cell_index()] =
999  cell->active_fe_index();
1000 
1001  Utilities::MPI::sum(active_fe_indices,
1002  tr->get_communicator(),
1003  active_fe_indices);
1004 
1005  // now go back and fill the active_fe_index on all other
1006  // cells. we would like to call cell->set_active_fe_index(),
1007  // but that function does not allow setting these indices on
1008  // non-locally_owned cells. so we have to work around the
1009  // issue a little bit by accessing the underlying data
1010  // structures directly
1011  for (const auto &cell : dof_handler.active_cell_iterators())
1012  if (!cell->is_locally_owned())
1013  dof_handler.levels[cell->level()]->set_active_fe_index(
1014  cell->index(),
1015  active_fe_indices[cell->active_cell_index()]);
1016  }
1017  else if (const ::parallel::
1018  DistributedTriangulationBase<dim, spacedim> *tr =
1019  dynamic_cast<
1020  const ::parallel::
1021  DistributedTriangulationBase<dim, spacedim> *>(
1022  &dof_handler.get_triangulation()))
1023  {
1024  // For completely distributed meshes, use the function that is
1025  // able to move data from locally owned cells on one processor to
1026  // the corresponding ghost cells on others. To this end, we need
1027  // to have functions that can pack and unpack the data we want to
1028  // transport -- namely, the single unsigned int active_fe_index
1029  // objects
1030  auto pack =
1031  [](const typename ::hp::DoFHandler<dim, spacedim>::
1032  active_cell_iterator &cell) -> unsigned int {
1033  return cell->active_fe_index();
1034  };
1035 
1036  auto unpack =
1037  [&dof_handler](
1038  const typename ::hp::DoFHandler<dim, spacedim>::
1039  active_cell_iterator &cell,
1040  const unsigned int active_fe_index) -> void {
1041  // we would like to say
1042  // cell->set_active_fe_index(active_fe_index);
1043  // but this is not allowed on cells that are not
1044  // locally owned, and we are on a ghost cell
1045  dof_handler.levels[cell->level()]->set_active_fe_index(
1046  cell->index(), active_fe_index);
1047  };
1048 
1050  unsigned int,
1051  ::hp::DoFHandler<dim, spacedim>>(dof_handler,
1052  pack,
1053  unpack);
1054  }
1055  else
1056  {
1057  // a sequential triangulation. there is nothing we need to do here
1058  Assert(
1059  (dynamic_cast<
1060  const ::parallel::TriangulationBase<dim, spacedim> *>(
1061  &dof_handler.get_triangulation()) == nullptr),
1062  ExcInternalError());
1063  }
1064  }
1065 
1066 
1067 
1088  template <int dim, int spacedim>
1089  static void
1091  ::hp::DoFHandler<dim, spacedim> &dof_handler)
1092  {
1093  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1094 
1095  for (const auto &cell : dof_handler.active_cell_iterators())
1096  if (cell->is_locally_owned())
1097  {
1098  if (cell->refine_flag_set())
1099  {
1100  // Store the active_fe_index of each cell that will be
1101  // refined to and distribute it later on its children.
1102  // Pick their future index if flagged for p-refinement.
1103  fe_transfer->refined_cells_fe_index.insert(
1104  {cell, cell->future_fe_index()});
1105  }
1106  else if (cell->coarsen_flag_set())
1107  {
1108  // From all cells that will be coarsened, determine their
1109  // parent and calculate its proper active_fe_index, so that
1110  // it can be set after refinement. But first, check if that
1111  // particular cell has a parent at all.
1112  Assert(cell->level() > 0, ExcInternalError());
1113  const auto &parent = cell->parent();
1114 
1115  // Check if the active_fe_index for the current cell has
1116  // been determined already.
1117  if (fe_transfer->coarsened_cells_fe_index.find(parent) ==
1118  fe_transfer->coarsened_cells_fe_index.end())
1119  {
1120  // Find a suitable active_fe_index for the parent cell
1121  // based on the 'least dominant finite element' of its
1122  // children. Consider the childrens' hypothetical future
1123  // index when they have been flagged for p-refinement.
1124  std::set<unsigned int> fe_indices_children;
1125  for (unsigned int child_index = 0;
1126  child_index < parent->n_children();
1127  ++child_index)
1128  {
1129  const auto &sibling = parent->child(child_index);
1130  Assert(sibling->is_active() &&
1131  sibling->coarsen_flag_set(),
1133  dim>::ExcInconsistentCoarseningFlags());
1134 
1135  fe_indices_children.insert(
1136  sibling->future_fe_index());
1137  }
1138  Assert(!fe_indices_children.empty(),
1139  ExcInternalError());
1140 
1141  const unsigned int fe_index =
1142  dof_handler.fe_collection.find_dominated_fe_extended(
1143  fe_indices_children, /*codim=*/0);
1144 
1146  typename ::hp::FECollection<dim>::
1147  ExcNoDominatedFiniteElementAmongstChildren());
1148 
1149  fe_transfer->coarsened_cells_fe_index.insert(
1150  {parent, fe_index});
1151  }
1152  }
1153  else
1154  {
1155  // No h-refinement is scheduled for this cell.
1156  // However, it may have p-refinement indicators, so we
1157  // choose a new active_fe_index based on its flags.
1158  if (cell->future_fe_index_set() == true)
1159  fe_transfer->persisting_cells_fe_index.insert(
1160  {cell, cell->future_fe_index()});
1161  }
1162  }
1163  }
1164 
1165 
1166 
1171  template <int dim, int spacedim>
1172  static void
1174  ::hp::DoFHandler<dim, spacedim> &dof_handler)
1175  {
1176  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1177 
1178  // Set active_fe_indices on persisting cells.
1179  for (const auto &persist : fe_transfer->persisting_cells_fe_index)
1180  {
1181  const auto &cell = persist.first;
1182 
1183  if (cell->is_locally_owned())
1184  {
1185  Assert(cell->is_active(), ExcInternalError());
1186  cell->set_active_fe_index(persist.second);
1187  }
1188  }
1189 
1190  // Distribute active_fe_indices from all refined cells on their
1191  // respective children.
1192  for (const auto &refine : fe_transfer->refined_cells_fe_index)
1193  {
1194  const auto &parent = refine.first;
1195 
1196  for (unsigned int child_index = 0;
1197  child_index < parent->n_children();
1198  ++child_index)
1199  {
1200  const auto &child = parent->child(child_index);
1201  Assert(child->is_locally_owned() && child->is_active(),
1202  ExcInternalError());
1203  child->set_active_fe_index(refine.second);
1204  }
1205  }
1206 
1207  // Set active_fe_indices on coarsened cells that have been determined
1208  // before the actual coarsening happened.
1209  for (const auto &coarsen : fe_transfer->coarsened_cells_fe_index)
1210  {
1211  const auto &cell = coarsen.first;
1212  Assert(cell->is_locally_owned() && cell->is_active(),
1213  ExcInternalError());
1214  cell->set_active_fe_index(coarsen.second);
1215  }
1216  }
1217 
1218 
1229  template <int dim, int spacedim>
1230  static unsigned int
1233  const std::vector<unsigned int> & children_fe_indices,
1234  ::hp::FECollection<dim, spacedim> &fe_collection)
1235  {
1236  Assert(!children_fe_indices.empty(), ExcInternalError());
1237 
1238  // convert vector to set
1239  const std::set<unsigned int> children_fe_indices_set(
1240  children_fe_indices.begin(), children_fe_indices.end());
1241 
1242  const unsigned int dominated_fe_index =
1243  fe_collection.find_dominated_fe_extended(children_fe_indices_set,
1244  /*codim=*/0);
1245 
1246  Assert(dominated_fe_index != numbers::invalid_unsigned_int,
1247  typename ::hp::FECollection<
1248  dim>::ExcNoDominatedFiniteElementAmongstChildren());
1249 
1250  return dominated_fe_index;
1251  }
1252  };
1253  } // namespace DoFHandlerImplementation
1254  } // namespace hp
1255 } // namespace internal
1256 
1257 
1258 
1259 namespace hp
1260 {
1261  template <int dim, int spacedim>
1262  const unsigned int DoFHandler<dim, spacedim>::dimension;
1263 
1264  template <int dim, int spacedim>
1266 
1267 
1268 
1269  template <int dim, int spacedim>
1271  : tria(nullptr, typeid(*this).name())
1272  , faces(nullptr)
1273  {}
1274 
1275 
1276 
1277  template <int dim, int spacedim>
1279  const Triangulation<dim, spacedim> &tria)
1280  : tria(&tria, typeid(*this).name())
1281  , faces(nullptr)
1282  {
1283  setup_policy_and_listeners();
1284 
1285  create_active_fe_table();
1286  }
1287 
1288 
1289 
1290  template <int dim, int spacedim>
1292  {
1293  // unsubscribe as a listener to refinement of the underlying
1294  // triangulation
1295  for (auto &connection : tria_listeners)
1296  connection.disconnect();
1297  tria_listeners.clear();
1298 
1299  // ...and release allocated memory
1300  // virtual functions called in constructors and destructors never use the
1301  // override in a derived class
1302  // for clarity be explicit on which function is called
1304  }
1305 
1306 
1307 
1308  /*------------------------ Cell iterator functions ------------------------*/
1309 
1310 
1311  template <int dim, int spacedim>
1313  DoFHandler<dim, spacedim>::begin(const unsigned int level) const
1314  {
1315  return cell_iterator(*this->get_triangulation().begin(level), this);
1316  }
1317 
1318 
1319 
1320  template <int dim, int spacedim>
1323  {
1324  // level is checked in begin
1325  cell_iterator i = begin(level);
1326  if (i.state() != IteratorState::valid)
1327  return i;
1328  while (i->has_children())
1329  if ((++i).state() != IteratorState::valid)
1330  return i;
1331  return i;
1332  }
1333 
1334 
1335 
1336  template <int dim, int spacedim>
1339  {
1340  return cell_iterator(&this->get_triangulation(), -1, -1, this);
1341  }
1342 
1343 
1344 
1345  template <int dim, int spacedim>
1347  DoFHandler<dim, spacedim>::end(const unsigned int level) const
1348  {
1349  return (level == this->get_triangulation().n_levels() - 1 ?
1350  end() :
1351  begin(level + 1));
1352  }
1353 
1354 
1355 
1356  template <int dim, int spacedim>
1359  {
1360  return (level == this->get_triangulation().n_levels() - 1 ?
1362  begin_active(level + 1));
1363  }
1364 
1365 
1366 
1367  template <int dim, int spacedim>
1370  {
1372  begin(), end());
1373  }
1374 
1375 
1376 
1377  template <int dim, int spacedim>
1380  {
1381  return IteratorRange<
1382  typename DoFHandler<dim, spacedim>::active_cell_iterator>(begin_active(),
1383  end());
1384  }
1385 
1386 
1387 
1388  template <int dim, int spacedim>
1391  const unsigned int level) const
1392  {
1394  begin(level), end(level));
1395  }
1396 
1397 
1398 
1399  template <int dim, int spacedim>
1402  const unsigned int level) const
1403  {
1404  return IteratorRange<
1406  begin_active(level), end_active(level));
1407  }
1408 
1409 
1410 
1411  //------------------------------------------------------------------
1412 
1413 
1414  template <int dim, int spacedim>
1417  {
1418  Assert(fe_collection.size() > 0, ExcNoFESelected());
1419 
1420  std::unordered_set<types::global_dof_index> boundary_dofs;
1421  std::vector<types::global_dof_index> dofs_on_face;
1422  dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
1423 
1424  const IndexSet &owned_dofs = locally_owned_dofs();
1425 
1426  // loop over all faces to check whether they are at a
1427  // boundary. note that we need not take special care of single
1428  // lines in 3d (using @p{cell->has_boundary_lines}), since we do
1429  // not support boundaries of dimension dim-2, and so every
1430  // boundary line is also part of a boundary face.
1431  typename HpDoFHandler<dim, spacedim>::active_cell_iterator
1432  cell = this->begin_active(),
1433  endc = this->end();
1434  for (; cell != endc; ++cell)
1435  if (cell->is_locally_owned() && cell->at_boundary())
1436  {
1437  for (auto f : GeometryInfo<dim>::face_indices())
1438  if (cell->at_boundary(f))
1439  {
1440  const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
1441  dofs_on_face.resize(dofs_per_face);
1442 
1443  cell->face(f)->get_dof_indices(dofs_on_face,
1444  cell->active_fe_index());
1445  for (unsigned int i = 0; i < dofs_per_face; ++i)
1446  {
1447  const unsigned int global_idof_index = dofs_on_face[i];
1448  if (owned_dofs.is_element(global_idof_index))
1449  {
1450  boundary_dofs.insert(global_idof_index);
1451  }
1452  }
1453  }
1454  }
1455  return boundary_dofs.size();
1456  }
1457 
1458 
1459 
1460  template <int dim, int spacedim>
1463  const std::set<types::boundary_id> &boundary_ids) const
1464  {
1465  Assert(fe_collection.size() > 0, ExcNoFESelected());
1466  Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
1467  boundary_ids.end(),
1469 
1470  // same as above, but with additional checks for set of boundary
1471  // indicators
1472  std::unordered_set<types::global_dof_index> boundary_dofs;
1473  std::vector<types::global_dof_index> dofs_on_face;
1474  dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
1475 
1476  const IndexSet &owned_dofs = locally_owned_dofs();
1477 
1478  typename HpDoFHandler<dim, spacedim>::active_cell_iterator
1479  cell = this->begin_active(),
1480  endc = this->end();
1481  for (; cell != endc; ++cell)
1482  if (cell->is_locally_owned() && cell->at_boundary())
1483  {
1484  for (auto f : GeometryInfo<dim>::face_indices())
1485  if (cell->at_boundary(f) &&
1486  (boundary_ids.find(cell->face(f)->boundary_id()) !=
1487  boundary_ids.end()))
1488  {
1489  const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
1490  dofs_on_face.resize(dofs_per_face);
1491 
1492  cell->face(f)->get_dof_indices(dofs_on_face,
1493  cell->active_fe_index());
1494  for (unsigned int i = 0; i < dofs_per_face; ++i)
1495  {
1496  const unsigned int global_idof_index = dofs_on_face[i];
1497  if (owned_dofs.is_element(global_idof_index))
1498  {
1499  boundary_dofs.insert(global_idof_index);
1500  }
1501  }
1502  }
1503  }
1504  return boundary_dofs.size();
1505  }
1506 
1507 
1508 
1509  template <int dim, int spacedim>
1510  std::size_t
1512  {
1513  std::size_t mem =
1515  MemoryConsumption::memory_consumption(fe_collection) +
1521  MemoryConsumption::memory_consumption(vertex_dof_offsets));
1522  for (const auto &level : levels)
1525 
1526  return mem;
1527  }
1528 
1529 
1530 
1531  template <int dim, int spacedim>
1532  void
1534  const std::vector<unsigned int> &active_fe_indices)
1535  {
1536  Assert(active_fe_indices.size() == get_triangulation().n_active_cells(),
1537  ExcDimensionMismatch(active_fe_indices.size(),
1538  get_triangulation().n_active_cells()));
1539 
1540  create_active_fe_table();
1541  // we could set the values directly, since they are stored as
1542  // protected data of this object, but for simplicity we use the
1543  // cell-wise access. this way we also have to pass some debug-mode
1544  // tests which we would have to duplicate ourselves otherwise
1545  for (const auto &cell : active_cell_iterators())
1546  if (cell->is_locally_owned())
1547  cell->set_active_fe_index(active_fe_indices[cell->active_cell_index()]);
1548  }
1549 
1550 
1551 
1552  template <int dim, int spacedim>
1553  void
1555  std::vector<unsigned int> &active_fe_indices) const
1556  {
1557  active_fe_indices.resize(get_triangulation().n_active_cells());
1558 
1559  // we could try to extract the values directly, since they are
1560  // stored as protected data of this object, but for simplicity we
1561  // use the cell-wise access.
1562  for (const auto &cell : active_cell_iterators())
1563  if (!cell->is_artificial())
1564  active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
1565  }
1566 
1567 
1568 
1569  template <int dim, int spacedim>
1570  void
1572  const Triangulation<dim, spacedim> & tria,
1574  {
1575  clear();
1576 
1577  if (this->tria != &tria)
1578  {
1579  for (auto &connection : tria_listeners)
1580  connection.disconnect();
1581  tria_listeners.clear();
1582 
1583  this->tria = &tria;
1584 
1585  setup_policy_and_listeners();
1586  }
1587 
1588  create_active_fe_table();
1589 
1590  distribute_dofs(fe);
1591  }
1592 
1593 
1594 
1595  template <int dim, int spacedim>
1596  void
1598  {
1599  Assert(
1600  tria != nullptr,
1601  ExcMessage(
1602  "You need to set the Triangulation in the DoFHandler using initialize() or "
1603  "in the constructor before you can distribute DoFs."));
1604  Assert(tria->n_levels() > 0,
1605  ExcMessage("The Triangulation you are using is empty!"));
1606  Assert(ff.size() > 0, ExcMessage("The hp::FECollection given is empty!"));
1607 
1608  // don't create a new object if the one we have is already appropriate
1609  if (fe_collection != ff)
1610  fe_collection = hp::FECollection<dim, spacedim>(ff);
1611 
1612  // ensure that the active_fe_indices vectors are initialized correctly
1613  create_active_fe_table();
1614 
1615  // make sure every processor knows the active_fe_indices
1616  // on both its own cells and all ghost cells
1619 
1620  // make sure that the fe collection is large enough to
1621  // cover all fe indices presently in use on the mesh
1622  for (const auto &cell : active_cell_iterators())
1623  if (!cell->is_artificial())
1624  Assert(cell->active_fe_index() < fe_collection.size(),
1625  ExcInvalidFEIndex(cell->active_fe_index(),
1626  fe_collection.size()));
1627  }
1628 
1629 
1630 
1631  template <int dim, int spacedim>
1632  void
1635  {
1636  // assign the fe_collection and initialize all active_fe_indices
1637  set_fe(ff);
1638 
1639  // If an underlying shared::Tria allows artificial cells,
1640  // then save the current set of subdomain ids, and set
1641  // subdomain ids to the "true" owner of each cell. we later
1642  // restore these flags
1643  std::vector<types::subdomain_id> saved_subdomain_ids;
1645  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
1646  &get_triangulation()));
1647  if (shared_tria != nullptr && shared_tria->with_artificial_cells())
1648  {
1649  saved_subdomain_ids.resize(shared_tria->n_active_cells());
1650 
1651  const std::vector<types::subdomain_id> &true_subdomain_ids =
1652  shared_tria->get_true_subdomain_ids_of_cells();
1653 
1654  for (const auto &cell : shared_tria->active_cell_iterators())
1655  {
1656  const unsigned int index = cell->active_cell_index();
1657  saved_subdomain_ids[index] = cell->subdomain_id();
1658  cell->set_subdomain_id(true_subdomain_ids[index]);
1659  }
1660  }
1661 
1662  // then allocate space for all the other tables
1664  reserve_space(*this);
1665 
1666  // now undo the subdomain modification
1667  if (shared_tria != nullptr && shared_tria->with_artificial_cells())
1668  for (const auto &cell : shared_tria->active_cell_iterators())
1669  cell->set_subdomain_id(saved_subdomain_ids[cell->active_cell_index()]);
1670 
1671 
1672  // Clear user flags because we will need them. But first we save
1673  // them and make sure that we restore them later such that at the
1674  // end of this function the Triangulation will be in the same
1675  // state as it was at the beginning of this function.
1676  std::vector<bool> user_flags;
1677  tria->save_user_flags(user_flags);
1678  const_cast<Triangulation<dim, spacedim> &>(*tria).clear_user_flags();
1679 
1680 
1682 
1683  // Now for the real work:
1684  number_cache = policy->distribute_dofs();
1685 
1687 
1688  // do some housekeeping: compress indices
1689  {
1691  for (int level = levels.size() - 1; level >= 0; --level)
1692  tg += Threads::new_task(
1693  &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
1694  *levels[level],
1695  fe_collection);
1696  tg.join_all();
1697  }
1698 
1699  // finally restore the user flags
1700  const_cast<Triangulation<dim, spacedim> &>(*tria).load_user_flags(
1701  user_flags);
1702  }
1703 
1704 
1705 
1706  template <int dim, int spacedim>
1707  void
1709  {
1710  // connect functions to signals of the underlying triangulation
1711  tria_listeners.push_back(this->tria->signals.pre_refinement.connect(
1712  [this]() { this->pre_refinement_action(); }));
1713  tria_listeners.push_back(this->tria->signals.post_refinement.connect(
1714  [this]() { this->post_refinement_action(); }));
1715  tria_listeners.push_back(this->tria->signals.create.connect(
1716  [this]() { this->post_refinement_action(); }));
1717 
1718  // decide whether we need a sequential or a parallel shared/distributed
1719  // policy and attach corresponding callback functions dealing with the
1720  // transfer of active_fe_indices
1722  *>(&this->get_triangulation()))
1723  {
1724  policy = std_cxx14::make_unique<
1726  DoFHandler<dim, spacedim>>>(*this);
1727 
1728  // repartitioning signals
1729  tria_listeners.push_back(
1730  this->tria->signals.pre_distributed_repartition.connect([this]() {
1731  internal::hp::DoFHandlerImplementation::Implementation::
1732  ensure_absence_of_future_fe_indices<dim, spacedim>(*this);
1733  }));
1734  tria_listeners.push_back(
1735  this->tria->signals.pre_distributed_repartition.connect(
1736  [this]() { this->pre_distributed_active_fe_index_transfer(); }));
1737  tria_listeners.push_back(
1738  this->tria->signals.post_distributed_repartition.connect(
1739  [this] { this->post_distributed_active_fe_index_transfer(); }));
1740 
1741  // refinement signals
1742  tria_listeners.push_back(
1743  this->tria->signals.pre_distributed_refinement.connect(
1744  [this]() { this->pre_distributed_active_fe_index_transfer(); }));
1745  tria_listeners.push_back(
1746  this->tria->signals.post_distributed_refinement.connect(
1747  [this]() { this->post_distributed_active_fe_index_transfer(); }));
1748 
1749  // serialization signals
1750  tria_listeners.push_back(
1751  this->tria->signals.post_distributed_save.connect([this]() {
1752  this->post_distributed_serialization_of_active_fe_indices();
1753  }));
1754  }
1755  else if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim>
1756  *>(&this->get_triangulation()) != nullptr)
1757  {
1758  policy =
1759  std_cxx14::make_unique<internal::DoFHandlerImplementation::Policy::
1760  ParallelShared<DoFHandler<dim, spacedim>>>(
1761  *this);
1762 
1763  // partitioning signals
1764  tria_listeners.push_back(
1765  this->tria->signals.pre_partition.connect([this]() {
1766  internal::hp::DoFHandlerImplementation::Implementation::
1767  ensure_absence_of_future_fe_indices(*this);
1768  }));
1769 
1770  // refinement signals
1771  tria_listeners.push_back(this->tria->signals.pre_refinement.connect(
1772  [this] { this->pre_active_fe_index_transfer(); }));
1773  tria_listeners.push_back(this->tria->signals.post_refinement.connect(
1774  [this] { this->post_active_fe_index_transfer(); }));
1775  }
1776  else
1777  {
1778  policy =
1779  std_cxx14::make_unique<internal::DoFHandlerImplementation::Policy::
1780  Sequential<DoFHandler<dim, spacedim>>>(
1781  *this);
1782 
1783  // refinement signals
1784  tria_listeners.push_back(this->tria->signals.pre_refinement.connect(
1785  [this] { this->pre_active_fe_index_transfer(); }));
1786  tria_listeners.push_back(this->tria->signals.post_refinement.connect(
1787  [this] { this->post_active_fe_index_transfer(); }));
1788  }
1789  }
1790 
1791 
1792 
1793  template <int dim, int spacedim>
1794  void
1796  {
1797  // release memory
1798  clear_space();
1799  }
1800 
1801 
1802 
1803  template <int dim, int spacedim>
1804  void
1806  const std::vector<types::global_dof_index> &new_numbers)
1807  {
1808  Assert(levels.size() > 0,
1809  ExcMessage(
1810  "You need to distribute DoFs before you can renumber them."));
1811 
1812  AssertDimension(new_numbers.size(), n_locally_owned_dofs());
1813 
1814 #ifdef DEBUG
1815  // assert that the new indices are
1816  // consecutively numbered if we are
1817  // working on a single
1818  // processor. this doesn't need to
1819  // hold in the case of a parallel
1820  // mesh since we map the interval
1821  // [0...n_dofs()) into itself but
1822  // only globally, not on each
1823  // processor
1824  if (n_locally_owned_dofs() == n_dofs())
1825  {
1826  std::vector<types::global_dof_index> tmp(new_numbers);
1827  std::sort(tmp.begin(), tmp.end());
1828  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
1830  for (; p != tmp.end(); ++p, ++i)
1831  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
1832  }
1833  else
1834  for (const auto new_number : new_numbers)
1835  Assert(new_number < n_dofs(),
1836  ExcMessage(
1837  "New DoF index is not less than the total number of dofs."));
1838 #endif
1839 
1840  // uncompress the internal storage scheme of dofs on cells so that
1841  // we can access dofs in turns. uncompress in parallel, starting
1842  // with the most expensive levels (the highest ones)
1843  {
1845  for (int level = levels.size() - 1; level >= 0; --level)
1846  tg += Threads::new_task(
1847  &::internal::hp::DoFLevel::uncompress_data<dim, spacedim>,
1848  *levels[level],
1849  fe_collection);
1850  tg.join_all();
1851  }
1852 
1853  // do the renumbering
1854  number_cache = policy->renumber_dofs(new_numbers);
1855 
1856  // now re-compress the dof indices
1857  {
1859  for (int level = levels.size() - 1; level >= 0; --level)
1860  tg += Threads::new_task(
1861  &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
1862  *levels[level],
1863  fe_collection);
1864  tg.join_all();
1865  }
1866  }
1867 
1868 
1869 
1870  template <int dim, int spacedim>
1871  unsigned int
1873  {
1874  Assert(fe_collection.size() > 0, ExcNoFESelected());
1875  return ::internal::hp::DoFHandlerImplementation::Implementation::
1876  max_couplings_between_dofs(*this);
1877  }
1878 
1879 
1880 
1881  template <int dim, int spacedim>
1882  unsigned int
1884  {
1885  Assert(fe_collection.size() > 0, ExcNoFESelected());
1886 
1887  switch (dim)
1888  {
1889  case 1:
1890  return fe_collection.max_dofs_per_vertex();
1891  case 2:
1892  return (3 * fe_collection.max_dofs_per_vertex() +
1893  2 * fe_collection.max_dofs_per_line());
1894  case 3:
1895  // we need to take refinement of one boundary face into
1896  // consideration here; in fact, this function returns what
1897  // #max_coupling_between_dofs<2> returns
1898  //
1899  // we assume here, that only four faces meet at the boundary;
1900  // this assumption is not justified and needs to be fixed some
1901  // time. fortunately, omitting it for now does no harm since
1902  // the matrix will cry foul if its requirements are not
1903  // satisfied
1904  return (19 * fe_collection.max_dofs_per_vertex() +
1905  28 * fe_collection.max_dofs_per_line() +
1906  8 * fe_collection.max_dofs_per_quad());
1907  default:
1908  Assert(false, ExcNotImplemented());
1909  return 0;
1910  }
1911  }
1912 
1913 
1914 
1915  template <int dim, int spacedim>
1916  void
1918  {
1919  // Create sufficiently many hp::DoFLevels.
1920  while (levels.size() < tria->n_levels())
1921  levels.emplace_back(new ::internal::hp::DoFLevel);
1922 
1923  // then make sure that on each level we have the appropriate size
1924  // of active_fe_indices; preset them to zero, i.e. the default FE
1925  for (unsigned int level = 0; level < levels.size(); ++level)
1926  {
1927  if (levels[level]->active_fe_indices.size() == 0 &&
1928  levels[level]->future_fe_indices.size() == 0)
1929  {
1930  levels[level]->active_fe_indices.resize(tria->n_raw_cells(level),
1931  0);
1932  levels[level]->future_fe_indices.resize(
1933  tria->n_raw_cells(level),
1935  }
1936  else
1937  {
1938  // Either the active_fe_indices have size zero because
1939  // they were just created, or the correct size. Other
1940  // sizes indicate that something went wrong.
1941  Assert(levels[level]->active_fe_indices.size() ==
1942  tria->n_raw_cells(level) &&
1943  levels[level]->future_fe_indices.size() ==
1944  tria->n_raw_cells(level),
1945  ExcInternalError());
1946  }
1947 
1948  // it may be that the previous table was compressed; in that
1949  // case, restore the correct active_fe_index. the fact that
1950  // this no longer matches the indices in the table is of no
1951  // importance because the current function is called at a
1952  // point where we have to recreate the dof_indices tables in
1953  // the levels anyway
1954  levels[level]->normalize_active_fe_indices();
1955  }
1956  }
1957 
1958 
1959 
1960  template <int dim, int spacedim>
1961  void
1963  {
1964  create_active_fe_table();
1965  }
1966 
1967 
1968 
1969  template <int dim, int spacedim>
1970  void
1972  {
1973  // Normally only one level is added, but if this Triangulation
1974  // is created by copy_triangulation, it can be more than one level.
1975  while (levels.size() < tria->n_levels())
1976  levels.emplace_back(new ::internal::hp::DoFLevel);
1977 
1978  // Coarsening can lead to the loss of levels. Hence remove them.
1979  while (levels.size() > tria->n_levels())
1980  {
1981  // drop the last element. that also releases the memory pointed to
1982  levels.pop_back();
1983  }
1984 
1985  Assert(levels.size() == tria->n_levels(), ExcInternalError());
1986  for (unsigned int i = 0; i < levels.size(); ++i)
1987  {
1988  // Resize active_fe_indices vectors. Use zero indicator to extend.
1989  levels[i]->active_fe_indices.resize(tria->n_raw_cells(i), 0);
1990 
1991  // Resize future_fe_indices vectors. Make sure that all
1992  // future_fe_indices have been cleared after refinement happened.
1993  //
1994  // We have used future_fe_indices to update all active_fe_indices
1995  // before refinement happened, thus we are safe to clear them now.
1996  levels[i]->future_fe_indices.assign(
1997  tria->n_raw_cells(i),
1999  }
2000  }
2001 
2002 
2003 
2004  template <int dim, int spacedim>
2005  void
2007  {
2008  // Finite elements need to be assigned to each cell by calling
2009  // distribute_dofs() first to make this functionality available.
2010  if (fe_collection.size() > 0)
2011  {
2012  Assert(active_fe_index_transfer == nullptr, ExcInternalError());
2013 
2014  active_fe_index_transfer =
2015  std_cxx14::make_unique<ActiveFEIndexTransfer>();
2016 
2019  }
2020  }
2021 
2022 
2023 
2024  template <int dim, int spacedim>
2025  void
2027  {
2028 #ifndef DEAL_II_WITH_P4EST
2029  Assert(false,
2030  ExcMessage(
2031  "You are attempting to use a functionality that is only available "
2032  "if deal.II was configured to use p4est, but cmake did not find a "
2033  "valid p4est library."));
2034 #else
2035  // the implementation below requires a p:d:T currently
2036  Assert((dynamic_cast<
2038  &this->get_triangulation()) != nullptr),
2039  ExcNotImplemented());
2040 
2041  // Finite elements need to be assigned to each cell by calling
2042  // distribute_dofs() first to make this functionality available.
2043  if (fe_collection.size() > 0)
2044  {
2045  Assert(active_fe_index_transfer == nullptr, ExcInternalError());
2046 
2047  active_fe_index_transfer =
2048  std_cxx14::make_unique<ActiveFEIndexTransfer>();
2049 
2050  // If we work on a p::d::Triangulation, we have to transfer all
2051  // active_fe_indices since ownership of cells may change. We will
2052  // use our p::d::CellDataTransfer member to achieve this. Further,
2053  // we prepare the values in such a way that they will correspond to
2054  // the active_fe_indices on the new mesh.
2055 
2056  // Gather all current future_fe_indices.
2057  active_fe_index_transfer->active_fe_indices.resize(
2058  get_triangulation().n_active_cells(), numbers::invalid_unsigned_int);
2059 
2060  for (const auto &cell : active_cell_iterators())
2061  if (cell->is_locally_owned())
2062  active_fe_index_transfer
2063  ->active_fe_indices[cell->active_cell_index()] =
2064  cell->future_fe_index();
2065 
2066  // Create transfer object and attach to it.
2067  const auto *distributed_tria = dynamic_cast<
2069  &this->get_triangulation());
2070 
2071  active_fe_index_transfer->cell_data_transfer = std_cxx14::make_unique<
2072  parallel::distributed::
2073  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
2074  *distributed_tria,
2075  /*transfer_variable_size_data=*/false,
2076  /*refinement_strategy=*/
2077  &::AdaptationStrategies::Refinement::
2078  preserve<dim, spacedim, unsigned int>,
2079  /*coarsening_strategy=*/
2080  [this](
2081  const typename Triangulation<dim, spacedim>::cell_iterator &parent,
2082  const std::vector<unsigned int> &children_fe_indices)
2083  -> unsigned int {
2084  return ::internal::hp::DoFHandlerImplementation::
2085  Implementation::determine_fe_from_children<dim, spacedim>(
2086  parent, children_fe_indices, fe_collection);
2087  });
2088 
2089  active_fe_index_transfer->cell_data_transfer
2090  ->prepare_for_coarsening_and_refinement(
2091  active_fe_index_transfer->active_fe_indices);
2092  }
2093 #endif
2094  }
2095 
2096 
2097 
2098  template <int dim, int spacedim>
2099  void
2101  {
2102  // Finite elements need to be assigned to each cell by calling
2103  // distribute_dofs() first to make this functionality available.
2104  if (fe_collection.size() > 0)
2105  {
2106  Assert(active_fe_index_transfer != nullptr, ExcInternalError());
2107 
2110 
2111  // We have to distribute the information about active_fe_indices
2112  // of all cells (including the artificial ones) on all processors,
2113  // if a parallel::shared::Triangulation has been used.
2116 
2117  // Free memory.
2118  active_fe_index_transfer.reset();
2119  }
2120  }
2121 
2122 
2123 
2124  template <int dim, int spacedim>
2125  void
2127  {
2128 #ifndef DEAL_II_WITH_P4EST
2129  Assert(false,
2130  ExcMessage(
2131  "You are attempting to use a functionality that is only available "
2132  "if deal.II was configured to use p4est, but cmake did not find a "
2133  "valid p4est library."));
2134 #else
2135  // the implementation below requires a p:d:T currently
2136  Assert((dynamic_cast<
2138  &this->get_triangulation()) != nullptr),
2139  ExcNotImplemented());
2140 
2141  // Finite elements need to be assigned to each cell by calling
2142  // distribute_dofs() first to make this functionality available.
2143  if (fe_collection.size() > 0)
2144  {
2145  Assert(active_fe_index_transfer != nullptr, ExcInternalError());
2146 
2147  // Unpack active_fe_indices.
2148  active_fe_index_transfer->active_fe_indices.resize(
2149  get_triangulation().n_active_cells(), numbers::invalid_unsigned_int);
2150  active_fe_index_transfer->cell_data_transfer->unpack(
2151  active_fe_index_transfer->active_fe_indices);
2152 
2153  // Update all locally owned active_fe_indices.
2154  set_active_fe_indices(active_fe_index_transfer->active_fe_indices);
2155 
2156  // Update active_fe_indices on ghost cells.
2159 
2160  // Free memory.
2161  active_fe_index_transfer.reset();
2162  }
2163 #endif
2164  }
2165 
2166 
2167 
2168  template <int dim, int spacedim>
2169  void
2171  {
2172 #ifndef DEAL_II_WITH_P4EST
2173  Assert(false,
2174  ExcMessage(
2175  "You are attempting to use a functionality that is only available "
2176  "if deal.II was configured to use p4est, but cmake did not find a "
2177  "valid p4est library."));
2178 #else
2179  // the implementation below requires a p:d:T currently
2180  Assert((dynamic_cast<
2182  &this->get_triangulation()) != nullptr),
2183  ExcNotImplemented());
2184 
2185  // Finite elements need to be assigned to each cell by calling
2186  // distribute_dofs() first to make this functionality available.
2187  if (fe_collection.size() > 0)
2188  {
2189  Assert(active_fe_index_transfer == nullptr, ExcInternalError());
2190 
2191  active_fe_index_transfer =
2192  std_cxx14::make_unique<ActiveFEIndexTransfer>();
2193 
2194  // Create transfer object and attach to it.
2195  const auto *distributed_tria = dynamic_cast<
2197  &this->get_triangulation());
2198 
2199  active_fe_index_transfer->cell_data_transfer = std_cxx14::make_unique<
2200  parallel::distributed::
2201  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
2202  *distributed_tria,
2203  /*transfer_variable_size_data=*/false,
2204  /*refinement_strategy=*/
2205  &::AdaptationStrategies::Refinement::
2206  preserve<dim, spacedim, unsigned int>,
2207  /*coarsening_strategy=*/
2208  [this](
2209  const typename Triangulation<dim, spacedim>::cell_iterator &parent,
2210  const std::vector<unsigned int> &children_fe_indices)
2211  -> unsigned int {
2212  return ::internal::hp::DoFHandlerImplementation::
2213  Implementation::determine_fe_from_children<dim, spacedim>(
2214  parent, children_fe_indices, fe_collection);
2215  });
2216 
2217  // If we work on a p::d::Triangulation, we have to transfer all
2218  // active fe indices since ownership of cells may change.
2219 
2220  // Gather all current active_fe_indices
2221  get_active_fe_indices(active_fe_index_transfer->active_fe_indices);
2222 
2223  // Attach to transfer object
2224  active_fe_index_transfer->cell_data_transfer->prepare_for_serialization(
2225  active_fe_index_transfer->active_fe_indices);
2226  }
2227 #endif
2228  }
2229 
2230 
2231  template <int dim, int spacedim>
2232  void
2233  DoFHandler<dim,
2234  spacedim>::post_distributed_serialization_of_active_fe_indices()
2235  {
2236 #ifndef DEAL_II_WITH_P4EST
2237  Assert(false,
2238  ExcMessage(
2239  "You are attempting to use a functionality that is only available "
2240  "if deal.II was configured to use p4est, but cmake did not find a "
2241  "valid p4est library."));
2242 #else
2243  // the implementation below requires a p:d:T currently
2244  Assert((dynamic_cast<
2246  &this->get_triangulation()) != nullptr),
2247  ExcNotImplemented());
2248 
2249  if (fe_collection.size() > 0)
2250  {
2251  Assert(active_fe_index_transfer != nullptr, ExcInternalError());
2252 
2253  // Free memory.
2254  active_fe_index_transfer.reset();
2255  }
2256 #endif
2257  }
2258 
2259 
2260 
2261  template <int dim, int spacedim>
2262  void
2264  {
2265 #ifndef DEAL_II_WITH_P4EST
2266  Assert(false,
2267  ExcMessage(
2268  "You are attempting to use a functionality that is only available "
2269  "if deal.II was configured to use p4est, but cmake did not find a "
2270  "valid p4est library."));
2271 #else
2272  // the implementation below requires a p:d:T currently
2273  Assert((dynamic_cast<
2275  &this->get_triangulation()) != nullptr),
2276  ExcNotImplemented());
2277 
2278  // Finite elements need to be assigned to each cell by calling
2279  // distribute_dofs() first to make this functionality available.
2280  if (fe_collection.size() > 0)
2281  {
2282  Assert(active_fe_index_transfer == nullptr, ExcInternalError());
2283 
2284  active_fe_index_transfer =
2285  std_cxx14::make_unique<ActiveFEIndexTransfer>();
2286 
2287  // Create transfer object and attach to it.
2288  const auto *distributed_tria = dynamic_cast<
2290  &this->get_triangulation());
2291 
2292  active_fe_index_transfer->cell_data_transfer = std_cxx14::make_unique<
2293  parallel::distributed::
2294  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
2295  *distributed_tria,
2296  /*transfer_variable_size_data=*/false,
2297  /*refinement_strategy=*/
2298  &::AdaptationStrategies::Refinement::
2299  preserve<dim, spacedim, unsigned int>,
2300  /*coarsening_strategy=*/
2301  [this](
2302  const typename Triangulation<dim, spacedim>::cell_iterator &parent,
2303  const std::vector<unsigned int> &children_fe_indices)
2304  -> unsigned int {
2305  return ::internal::hp::DoFHandlerImplementation::
2306  Implementation::determine_fe_from_children<dim, spacedim>(
2307  parent, children_fe_indices, fe_collection);
2308  });
2309 
2310  // Unpack active_fe_indices.
2311  active_fe_index_transfer->active_fe_indices.resize(
2312  get_triangulation().n_active_cells(), numbers::invalid_unsigned_int);
2313  active_fe_index_transfer->cell_data_transfer->deserialize(
2314  active_fe_index_transfer->active_fe_indices);
2315 
2316  // Update all locally owned active_fe_indices.
2317  set_active_fe_indices(active_fe_index_transfer->active_fe_indices);
2318 
2319  // Update active_fe_indices on ghost cells.
2322 
2323  // Free memory.
2324  active_fe_index_transfer.reset();
2325  }
2326 #endif
2327  }
2328 
2329 
2330 
2331  template <int dim, int spacedim>
2332  template <int structdim>
2335  const unsigned int,
2336  const unsigned int,
2337  const unsigned int) const
2338  {
2339  Assert(false, ExcNotImplemented());
2341  }
2342 
2343 
2344 
2345  template <int dim, int spacedim>
2346  template <int structdim>
2347  void
2349  const unsigned int,
2350  const unsigned int,
2351  const unsigned int,
2352  const types::global_dof_index) const
2353  {
2354  Assert(false, ExcNotImplemented());
2355  }
2356 
2357 
2358 
2359  template <int dim, int spacedim>
2360  void
2362  {
2363  levels.clear();
2364  faces.reset();
2365 
2366  vertex_dofs = std::vector<types::global_dof_index>();
2367  vertex_dof_offsets = std::vector<unsigned int>();
2368  }
2369 } // namespace hp
2370 
2371 
2372 
2373 /*-------------- Explicit Instantiations -------------------------------*/
2374 #include "dof_handler.inst"
2375 
2376 
internal::hp::DoFHandlerImplementation::Implementation::ensure_absence_of_future_fe_indices
static void ensure_absence_of_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:94
IndexSet
Definition: index_set.h:74
internal::hp::DoFHandlerImplementation::Implementation::distribute_fe_indices_on_refined_cells
static void distribute_fe_indices_on_refined_cells(::hp::DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:1173
numbers::invalid_dof_index
const types::global_dof_index invalid_dof_index
Definition: types.h:206
parallel::shared::Triangulation
Definition: shared_tria.h:103
DoFHandler::faces
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > faces
Definition: dof_handler.h:1349
DoFHandler::active_cell_iterators_on_level
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
Definition: dof_handler.cc:1062
parallel::DistributedTriangulationBase
Definition: tria_base.h:352
tria_accessor.h
dof_faces.h
DoFHandler::levels
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > levels
Definition: dof_handler.h:1337
internal::hp::DoFHandlerImplementation::Implementation::max_couplings_between_dofs
static unsigned int max_couplings_between_dofs(const DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:866
internal::hp::DoFHandlerImplementation::Implementation
Definition: dof_handler.cc:85
internal::hp::DoFHandlerImplementation::Implementation::reserve_space
static void reserve_space(DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:706
hp::FECollection::size
unsigned int size() const
Definition: fe_collection.h:853
Triangulation::load_user_flags
void load_user_flags(std::istream &in)
Definition: tria.cc:11234
tria_levels.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Triangulation
Definition: tria.h:1109
tria.h
DoFHandler::end_active
active_cell_iterator end_active(const unsigned int level) const
Definition: dof_handler.cc:971
hp::DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:312
memory_consumption.h
tria_iterator.h
internal::DoFHandlerImplementation::Policy::ParallelDistributed
Definition: dof_handler_policy.h:221
hp::DoFHandler::fe_collection
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1051
internal::hp::DoFHandlerImplementation::Implementation::reserve_space_vertices
static void reserve_space_vertices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:160
GeometryInfo
Definition: geometry_info.h:1224
hp::DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
internal::hp::DoFLevel::invalid_active_fe_index
static const active_fe_index_type invalid_active_fe_index
Definition: dof_level.h:132
IteratorState::valid
@ valid
Iterator points to a valid object.
Definition: tria_iterator_base.h:38
Utilities::pack
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1209
Threads::new_task
Task< RT > new_task(const std::function< RT()> &function)
Definition: thread_management.h:1647
DoFHandler::renumber_dofs
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
Definition: dof_handler.cc:1363
hp::DoFHandler
Definition: dof_handler.h:203
DoFHandler::DoFHandler
DoFHandler()
Definition: dof_handler.cc:864
HpDoFHandler
#define HpDoFHandler
Definition: dof_handler.cc:57
DoFHandler::clear_space
void clear_space()
Definition: dof_handler.cc:1507
DoFHandler
Definition: dof_handler.h:205
DoFHandler::max_couplings_between_boundary_dofs
unsigned int max_couplings_between_boundary_dofs() const
Definition: dof_handler.cc:1470
GridRefinement::coarsen
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
Definition: grid_refinement.cc:88
DoFHandler::distribute_dofs
virtual void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
Definition: dof_handler.cc:1247
internal::hp::DoFLevel
Definition: dof_level.h:109
internal::TriangulationImplementation::n_active_cells
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12625
thread_management.h
level
unsigned int level
Definition: grid_out.cc:4355
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: dof_handler.cc:1030
DoFHandler::begin
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:922
DoFTools::ExcNoFESelected
static ::ExceptionBase & ExcNoFESelected()
DoFHandler::~DoFHandler
virtual ~DoFHandler() override
Definition: dof_handler.cc:870
DoFHandler::get_dof_index
types::global_dof_index get_dof_index(const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const
Definition: dof_handler.cc:1523
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:51
geometry_info.h
IteratorRange
Definition: iterator_range.h:129
internal::hp::DoFHandlerImplementation::Implementation::determine_fe_from_children
static unsigned int determine_fe_from_children(const typename Triangulation< dim, spacedim >::cell_iterator &, const std::vector< unsigned int > &children_fe_indices, ::hp::FECollection< dim, spacedim > &fe_collection)
Definition: dof_handler.cc:1231
grid_tools.h
Triangulation::clear_user_flags
void clear_user_flags()
Definition: tria.cc:11174
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
Utilities::unpack
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1353
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:268
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
dof_level.h
DoFHandler::vertex_dofs
std::vector< types::global_dof_index > vertex_dofs
Definition: dof_handler.h:1323
fe.h
parallel::distributed::Triangulation
Definition: tria.h:242
internal::hp::DoFHandlerImplementation::Implementation::max_couplings_between_dofs
static unsigned int max_couplings_between_dofs(const DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:931
shared_tria.h
IndexSet::is_element
bool is_element(const size_type index) const
Definition: index_set.h:1766
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
internal::hp::DoFHandlerImplementation::Implementation::max_couplings_between_dofs
static unsigned int max_couplings_between_dofs(const DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:853
MemorySpace::swap
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
Definition: memory_space.h:103
internal::hp::DoFHandlerImplementation::Implementation::reserve_space
static void reserve_space(DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:730
numbers::internal_face_boundary_id
const types::boundary_id internal_face_boundary_id
Definition: types.h:250
DoFTools::get_active_fe_indices
void get_active_fe_indices(const DoFHandlerType &dof_handler, std::vector< unsigned int > &active_fe_indices)
Definition: dof_tools.cc:1360
DoFHandler::cell_iterators
IteratorRange< cell_iterator > cell_iterators() const
Definition: dof_handler.cc:1021
hp::DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:325
hp::DoFHandler::levels
std::vector< std::unique_ptr<::internal::hp::DoFLevel > > levels
Definition: dof_handler.h:1180
Threads::TaskGroup::join_all
void join_all() const
Definition: thread_management.h:1833
DoFHandler::end
cell_iterator end() const
Definition: dof_handler.cc:951
DoFHandler::max_couplings_between_dofs
unsigned int max_couplings_between_dofs() const
Definition: dof_handler.cc:1460
parallel::Triangulation
TriangulationBase< dim, spacedim > Triangulation
Definition: tria_base.h:302
DoFHandler::cell_iterators_on_level
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
Definition: dof_handler.cc:1051
internal::hp::DoFHandlerImplementation::Implementation::reserve_space_cells
static void reserve_space_cells(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:271
unsigned int
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
tria.h
DoFHandler::max_dofs_per_face
unsigned int max_dofs_per_face(const DoFHandler< dim, spacedim > &dh)
DoFTools::ExcInvalidBoundaryIndicator
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
internal::hp::DoFIndicesOnFaces< 3 >::quads
internal::hp::DoFIndicesOnFacesOrEdges< 2 > quads
Definition: dof_faces.h:320
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
hp::DoFHandler::active_fe_index_transfer
std::unique_ptr< ActiveFEIndexTransfer > active_fe_index_transfer
Definition: dof_handler.h:1283
DoFHandler::memory_consumption
virtual std::size_t memory_consumption() const
Definition: dof_handler.cc:1203
internal::hp::DoFIndicesOnFacesOrEdges::dof_offsets
std::vector< unsigned int > dof_offsets
Definition: dof_faces.h:105
internal::hp::DoFIndicesOnFaces< 2 >
Definition: dof_faces.h:277
internal::hp::DoFHandlerImplementation::Implementation::reserve_space_faces
static void reserve_space_faces(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:365
DoFHandler::clear
virtual void clear()
Definition: dof_handler.cc:1352
GridTools::exchange_cell_data_to_ghosts
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack)
hp::FECollection::find_dominated_fe_extended
unsigned int find_dominated_fe_extended(const std::set< unsigned int > &fes, const unsigned int codim=0) const
Definition: fe_collection.cc:230
dof_accessor.h
Threads::TaskGroup
Definition: thread_management.h:1798
hp
Definition: hp.h:117
internal::hp::DoFHandlerImplementation::Implementation::reserve_space_release_space
static void reserve_space_release_space(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:114
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
Utilities::MPI::min
T min(const T &t, const MPI_Comm &mpi_communicator)
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
internal::hp::DoFIndicesOnFaces< 3 >
Definition: dof_faces.h:309
Triangulation::n_active_cells
unsigned int n_active_cells() const
Definition: tria.cc:12675
DoFHandler::set_dof_index
void set_dof_index(const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const
Definition: dof_handler.cc:1543
dof_handler_policy.h
parallel::shared::Triangulation::get_true_subdomain_ids_of_cells
const std::vector< types::subdomain_id > & get_true_subdomain_ids_of_cells() const
Definition: shared_tria.cc:333
internal
Definition: aligned_vector.h:369
memory.h
Triangulation::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12185
DoFHandler::n_boundary_dofs
types::global_dof_index n_boundary_dofs() const
Definition: dof_handler.cc:1089
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
internal::hp::DoFIndicesOnFacesOrEdges::dofs
std::vector< types::global_dof_index > dofs
Definition: dof_faces.h:111
DoFHandler::tria
SmartPointer< const Triangulation< dim, spacedim >, DoFHandler< dim, spacedim > > tria
Definition: dof_handler.h:1183
internal::hp::DoFHandlerImplementation::Implementation::reserve_space
static void reserve_space(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:684
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:935
GridRefinement::refine
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
Definition: grid_refinement.cc:41
hp::DoFHandler::DoFHandler
DoFHandler()
Definition: dof_handler.cc:1270
DoFHandler::fe_collection
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1190
TriaIterator
Definition: tria_iterator.h:578
internal::hp::DoFHandlerImplementation::Implementation::communicate_active_fe_indices
static void communicate_active_fe_indices(::hp::DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:976
DoFHandler::initialize
void initialize(const Triangulation< dim, spacedim > &tria, const FiniteElement< dim, spacedim > &fe)
Definition: dof_handler.cc:887
internal::hp::DoFHandlerImplementation::Implementation::collect_fe_indices_on_cells_to_be_refined
static void collect_fe_indices_on_cells_to_be_refined(::hp::DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:1090
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
DoFHandler::set_fe
virtual void set_fe(const FiniteElement< dim, spacedim > &fe)
Definition: dof_handler.cc:1235
parallel
Definition: distributed.h:416
parallel::shared::Triangulation::with_artificial_cells
bool with_artificial_cells() const
Definition: shared_tria.cc:324
hp::FECollection
Definition: fe_collection.h:54
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
int
dof_handler.h