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