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