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