Reference documentation for deal.II version 9.0.0
dof_handler.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/memory_consumption.h>
17 #include <deal.II/base/geometry_info.h>
18 #include <deal.II/base/std_cxx14/memory.h>
19 #include <deal.II/base/thread_management.h>
20 #include <deal.II/hp/dof_handler.h>
21 #include <deal.II/hp/dof_level.h>
22 #include <deal.II/hp/dof_faces.h>
23 #include <deal.II/dofs/dof_accessor.h>
24 #include <deal.II/dofs/dof_handler_policy.h>
25 #include <deal.II/grid/grid_tools.h>
26 #include <deal.II/grid/tria_accessor.h>
27 #include <deal.II/grid/tria_iterator.h>
28 #include <deal.II/grid/tria_levels.h>
29 #include <deal.II/grid/tria.h>
30 #include <deal.II/fe/fe.h>
31 #include <deal.II/distributed/shared_tria.h>
32 #include <deal.II/distributed/tria.h>
33 
34 #include <boost/serialization/array.hpp>
35 
36 #include <set>
37 #include <algorithm>
38 #include <functional>
39 
40 DEAL_II_NAMESPACE_OPEN
41 
42 // The following is necessary for compilation under Visual Studio which is unable to correctly
43 // distinguish between ::DoFHandler and ::hp::DoFHandler.
44 // Plus it makes code in dof_handler.cc easier to read.
45 #if defined(_MSC_VER) && (_MSC_VER >= 1800)
46 template <int dim, int spacedim> using HpDoFHandler = ::hp::DoFHandler<dim, spacedim>;
47 #else
48 // When using older Visual Studio or a different compiler just fall back.
49 #define HpDoFHandler DoFHandler
50 #endif
51 
52 namespace parallel
53 {
54  namespace distributed
55  {
56  template <int, int> class Triangulation;
57  }
58 }
59 
60 
61 
62 
63 namespace internal
64 {
65  namespace hp
66  {
67  namespace DoFHandlerImplementation
68  {
69  // access class ::hp::DoFHandler instead of namespace
70  // internal::hp::DoFHandler, etc
71  using ::hp::DoFHandler;
72 
78  {
83  template <int dim, int spacedim>
84  static
85  void
87  {
88  // Release all space except the active_fe_indices field
89  // which we have to back up before
90  {
91  std::vector<std::vector<DoFLevel::active_fe_index_type> >
92  active_fe_backup(dof_handler.levels.size ());
93  for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
94  active_fe_backup[level] = std::move(dof_handler.levels[level]->active_fe_indices);
95 
96  // delete all levels and set them up newly, since vectors
97  // are troublesome if you want to change their size
98  dof_handler.clear_space ();
99 
100  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
101  {
102  dof_handler.levels.emplace_back (new internal::hp::DoFLevel);
103  dof_handler.levels[level]->active_fe_indices = std::move(active_fe_backup[level]);
104  }
105 
106  if (dim > 1)
107  dof_handler.faces = std_cxx14::make_unique<internal::hp::DoFIndicesOnFaces<dim>> ();
108  }
109  }
110 
111 
112 
117  template <int dim, int spacedim>
118  static
119  void
121  {
122  // The final step in all of the reserve_space() functions is to set
123  // up vertex dof information. since vertices are sequentially
124  // numbered, what we do first is to set up an array in which
125  // we record whether a vertex is associated with any of the
126  // given fe's, by setting a bit. in a later step, we then
127  // actually allocate memory for the required dofs
128  //
129  // in the following, we only need to consider vertices that are
130  // adjacent to either a locally owned or a ghost cell; we never
131  // store anything on vertices that are only surrounded by
132  // artificial cells. so figure out that subset of vertices
133  // first
134  std::vector<bool> locally_used_vertices (dof_handler.tria->n_vertices(),
135  false);
136  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
137  cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
138  if (! cell->is_artificial())
139  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
140  locally_used_vertices[cell->vertex_index(v)] = true;
141 
142  std::vector<std::vector<bool> >
143  vertex_fe_association (dof_handler.fe_collection.size(),
144  std::vector<bool> (dof_handler.tria->n_vertices(), false));
145 
146  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
147  cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
148  if (! cell->is_artificial())
149  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
150  vertex_fe_association[cell->active_fe_index()][cell->vertex_index(v)]
151  = true;
152 
153  // in debug mode, make sure that each vertex is associated
154  // with at least one fe (note that except for unused
155  // vertices, all vertices are actually active). this is of
156  // course only true for vertices that are part of either
157  // ghost or locally owned cells
158 #ifdef DEBUG
159  for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
160  if (locally_used_vertices[v] == true)
161  if (dof_handler.tria->vertex_used(v) == true)
162  {
163  unsigned int fe=0;
164  for (; fe<dof_handler.fe_collection.size(); ++fe)
165  if (vertex_fe_association[fe][v] == true)
166  break;
167  Assert (fe != dof_handler.fe_collection.size(), ExcInternalError());
168  }
169 #endif
170 
171  // next count how much memory we actually need. for each
172  // vertex, we need one slot per fe to store the fe_index,
173  // plus dofs_per_vertex for this fe. in addition, we need
174  // one slot as the end marker for the fe_indices. at the
175  // same time already fill the vertex_dof_offsets field
176  dof_handler.vertex_dof_offsets.resize (dof_handler.tria->n_vertices(),
178 
179  unsigned int vertex_slots_needed = 0;
180  for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
181  if (dof_handler.tria->vertex_used(v) == true)
182  if (locally_used_vertices[v] == true)
183  {
184  dof_handler.vertex_dof_offsets[v] = vertex_slots_needed;
185 
186  for (unsigned int fe=0; fe<dof_handler.fe_collection.size(); ++fe)
187  if (vertex_fe_association[fe][v] == true)
188  vertex_slots_needed += dof_handler.get_fe(fe).dofs_per_vertex + 1;
189 
190  // don't forget the end_marker:
191  ++vertex_slots_needed;
192  }
193 
194  // now allocate the space we have determined we need, and
195  // set up the linked lists for each of the vertices
196  dof_handler.vertex_dofs.resize (vertex_slots_needed,
198  for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
199  if (dof_handler.tria->vertex_used(v) == true)
200  if (locally_used_vertices[v] == true)
201  {
202  unsigned int current_index = dof_handler.vertex_dof_offsets[v];
203  for (unsigned int fe=0; fe<dof_handler.fe_collection.size(); ++fe)
204  if (vertex_fe_association[fe][v] == true)
205  {
206  // if this vertex uses this fe, then set the
207  // fe_index and move the pointer ahead
208  dof_handler.vertex_dofs[current_index] = fe;
209  current_index += dof_handler.get_fe(fe).dofs_per_vertex + 1;
210  }
211  // finally place the end marker
212  dof_handler.vertex_dofs[current_index] = numbers::invalid_dof_index;
213  }
214  }
215 
216 
221  template <int dim, int spacedim>
222  static
223  void
225  {
226  // count how much space we need on each level for the cell
227  // dofs and set the dof_*_offsets data. initially set the
228  // latter to an invalid index, and only later set it to
229  // something reasonable for active dof_handler.cells
230  //
231  // note that for dof_handler.cells, the situation is simpler
232  // than for other (lower dimensional) objects since exactly
233  // one finite element is used for it
234  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
235  {
236  dof_handler.levels[level]->dof_offsets
237  = std::vector<DoFLevel::offset_type> (dof_handler.tria->n_raw_cells(level),
238  (DoFLevel::offset_type)(-1));
239  dof_handler.levels[level]->cell_cache_offsets
240  = std::vector<DoFLevel::offset_type> (dof_handler.tria->n_raw_cells(level),
241  (DoFLevel::offset_type)(-1));
242 
243  types::global_dof_index next_free_dof = 0;
244  types::global_dof_index cache_size = 0;
245  typename HpDoFHandler<dim,spacedim>::active_cell_iterator
246  cell=dof_handler.begin_active(level),
247  endc=dof_handler.end_active(level);
248  for (; cell!=endc; ++cell)
249  if (!cell->has_children()
250  &&
251  !cell->is_artificial())
252  {
253  dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
254  next_free_dof += cell->get_fe().template n_dofs_per_object<dim>();
255 
256  dof_handler.levels[level]->cell_cache_offsets[cell->index()] = cache_size;
257  cache_size += cell->get_fe().dofs_per_cell;
258  }
259 
260  dof_handler.levels[level]->dof_indices
261  = std::vector<types::global_dof_index> (next_free_dof,
263  dof_handler.levels[level]->cell_dof_indices_cache
264  = std::vector<types::global_dof_index> (cache_size,
266  }
267 
268  // safety check: make sure that the number of DoFs we
269  // allocated is actually correct (above we have also set the
270  // dof_*_offsets field, so we couldn't use this simpler
271  // algorithm)
272 #ifdef DEBUG
273  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
274  {
275  types::global_dof_index counter = 0;
276  typename HpDoFHandler<dim,spacedim>::active_cell_iterator
277  cell=dof_handler.begin_active(level),
278  endc=dof_handler.end_active(level);
279  for (; cell!=endc; ++cell)
280  if (!cell->has_children()
281  &&
282  !cell->is_artificial())
283  counter += cell->get_fe().template n_dofs_per_object<dim>();
284 
285  Assert (dof_handler.levels[level]->dof_indices.size() == counter,
286  ExcInternalError());
287 
288  // also check that the number of unassigned slots in the dof_offsets
289  // equals the number of cells on that level minus the number of
290  // active, non-artificial cells (because these are exactly the cells
291  // on which we do something)
292  unsigned int n_active_non_artificial_cells = 0;
293  for (cell=dof_handler.begin_active(level); cell!=endc; ++cell)
294  if (!cell->has_children()
295  &&
296  !cell->is_artificial())
297  ++n_active_non_artificial_cells;
298 
299  Assert (static_cast<unsigned int>
300  (std::count (dof_handler.levels[level]->dof_offsets.begin(),
301  dof_handler.levels[level]->dof_offsets.end(),
302  (DoFLevel::offset_type)(-1)))
303  ==
304  dof_handler.tria->n_raw_cells(level) - n_active_non_artificial_cells,
305  ExcInternalError());
306  }
307 #endif
308  }
309 
310 
311 
312 
317  template <int dim, int spacedim>
318  static
319  void
321  {
322  // make the code generic between lines and quads
323  std::vector<unsigned int> &face_dof_offsets
324  = (dim == 2
325  ?
326  reinterpret_cast<::internal::hp::DoFIndicesOnFaces<2>&>(*dof_handler.faces).lines.dof_offsets
327  :
328  reinterpret_cast<::internal::hp::DoFIndicesOnFaces<3>&>(*dof_handler.faces).quads.dof_offsets);
329 
330  std::vector<types::global_dof_index> &face_dof_indices
331  = (dim == 2
332  ?
333  reinterpret_cast<::internal::hp::DoFIndicesOnFaces<2>&>(*dof_handler.faces).lines.dofs
334  :
335  reinterpret_cast<::internal::hp::DoFIndicesOnFaces<3>&>(*dof_handler.faces).quads.dofs);
336 
337  // FACE DOFS
338  //
339  // count face dofs, then allocate as much space
340  // as we need and prime the linked list for faces (see the
341  // description in hp::DoFLevel) with the indices we will
342  // need. note that our task is more complicated than for the
343  // cell case above since two adjacent cells may have different
344  // active_fe_indices, in which case we need to allocate
345  // *two* sets of face dofs for the same face
346  //
347  // the way we do things is that we loop over all active
348  // cells (these are the ones that have DoFs only
349  // anyway) and all their faces. We note in the
350  // user flags whether we have previously visited a face and
351  // if so skip it (consequently, we have to save and later
352  // restore the face flags)
353  {
354  std::vector<bool> saved_face_user_flags;
355  switch (dim)
356  {
357  case 2:
358  {
359  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
360  .save_user_flags_line (saved_face_user_flags);
361  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
362  .clear_user_flags_line ();
363 
364  break;
365  }
366 
367  case 3:
368  {
369  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
370  .save_user_flags_quad (saved_face_user_flags);
371  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
372  .clear_user_flags_quad ();
373 
374  break;
375  }
376 
377  default:
378  Assert (false, ExcNotImplemented());
379  }
380 
381  // an array to hold how many slots (see the hp::DoFLevel
382  // class) we will have to store on each level
383  unsigned int n_face_slots = 0;
384 
385  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
386  cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
387  if (! cell->is_artificial())
388  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
389  if (! cell->face(face)->user_flag_set())
390  {
391  // ok, face has not been visited. so we need to
392  // allocate space for it. let's see how much we
393  // need: we need one set if a) there is no
394  // neighbor behind this face, or b) the neighbor
395  // is either coarser or finer than we are, or c)
396  // the neighbor is artificial, or d) the neighbor
397  // is neither coarser nor finer, nor is artificial,
398  // and just so happens to have the same active_fe_index :
399  if (cell->at_boundary(face)
400  ||
401  cell->face(face)->has_children()
402  ||
403  cell->neighbor_is_coarser(face)
404  ||
405  (!cell->at_boundary(face)
406  &&
407  cell->neighbor(face)->is_artificial())
408  ||
409  (!cell->at_boundary(face)
410  &&
411  !cell->neighbor(face)->is_artificial()
412  &&
413  (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
414  // ok, one set of dofs. that makes one active_fe_index, 1
415  // times dofs_per_face dofs, and one stop index
416  n_face_slots
417  += 1 + dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>() + 1;
418 
419  // otherwise we do indeed need two sets, i.e. two
420  // active_fe_indices, two sets of dofs, and one stop index:
421  else
422  n_face_slots
423  += (2
424  +
425  dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>()
426  +
427  dof_handler.get_fe(cell->neighbor(face)->active_fe_index())
428  .template n_dofs_per_object<dim-1>()
429  +
430  1);
431 
432  // mark this face as visited
433  cell->face(face)->set_user_flag ();
434  }
435 
436  // now that we know how many face dofs we will have to
437  // have on each level, allocate the memory. note that we
438  // allocate offsets for all faces, though only the active
439  // ones will have a non-invalid value later on
440  face_dof_offsets
441  = std::vector<unsigned int> (dof_handler.tria->n_raw_faces(),
442  (unsigned int)(-1));
443  face_dof_indices
444  = std::vector<types::global_dof_index> (n_face_slots,
446 
447  // with the memory now allocated, loop over the
448  // dof_handler.cells again and prime the _offset values as
449  // well as the fe_index fields
450  switch (dim)
451  {
452  case 2:
453  {
454  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
455  .clear_user_flags_line ();
456 
457  break;
458  }
459 
460  case 3:
461  {
462  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
463  .clear_user_flags_quad ();
464 
465  break;
466  }
467 
468  default:
469  Assert (false, ExcNotImplemented());
470  }
471 
472  unsigned int next_free_face_slot = 0;
473 
474  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
475  cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
476  if (! cell->is_artificial())
477  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
478  if (! cell->face(face)->user_flag_set())
479  {
480  // same decision tree as before
481  if (cell->at_boundary(face)
482  ||
483  cell->face(face)->has_children()
484  ||
485  cell->neighbor_is_coarser(face)
486  ||
487  (!cell->at_boundary(face)
488  &&
489  cell->neighbor(face)->is_artificial())
490  ||
491  (!cell->at_boundary(face)
492  &&
493  !cell->neighbor(face)->is_artificial()
494  &&
495  (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
496  {
497  face_dof_offsets[cell->face(face)->index()]
498  = next_free_face_slot;
499 
500  // set first slot for this face to
501  // active_fe_index of this face
502  face_dof_indices[next_free_face_slot]
503  = cell->active_fe_index();
504 
505  // the next dofs_per_face indices remain unset
506  // for the moment (i.e. at invalid_dof_index).
507  // following this comes the stop index, which
508  // also is invalid_dof_index and therefore
509  // does not have to be explicitly set
510 
511  // finally, mark those slots as used
512  next_free_face_slot
513  += dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>() + 2;
514  }
515  else
516  {
517  face_dof_offsets[cell->face(face)->index()]
518  = next_free_face_slot;
519 
520  // set first slot for this face to
521  // active_fe_index of this face
522  face_dof_indices[next_free_face_slot]
523  = cell->active_fe_index();
524 
525  // the next dofs_per_line/quad indices remain unset
526  // for the moment (i.e. at invalid_dof_index).
527  //
528  // then comes the fe_index for the neighboring
529  // cell:
530  face_dof_indices[next_free_face_slot
531  +
532  dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>()
533  +
534  1]
535  = cell->neighbor(face)->active_fe_index();
536  // then again a set of dofs that we need not
537  // set right now
538  //
539  // following this comes the stop index, which
540  // also is invalid_dof_index and therefore
541  // does not have to be explicitly set
542 
543  // finally, mark those slots as used
544  next_free_face_slot
545  += (dof_handler.get_fe(cell->active_fe_index()).template n_dofs_per_object<dim-1>()
546  +
547  dof_handler.get_fe(cell->neighbor(face)->active_fe_index())
548  .template n_dofs_per_object<dim-1>()
549  +
550  3);
551  }
552 
553  // mark this face as visited
554  cell->face(face)->set_user_flag ();
555  }
556 
557  // we should have moved the cursor for each level to the
558  // total number of dofs on that level. check that
559  Assert (next_free_face_slot == n_face_slots,
560  ExcInternalError());
561 
562  // at the end, restore the user flags for the faces
563  switch (dim)
564  {
565  case 2:
566  {
567  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
568  .load_user_flags_line (saved_face_user_flags);
569 
570  break;
571  }
572 
573  case 3:
574  {
575  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
576  .load_user_flags_quad (saved_face_user_flags);
577 
578  break;
579  }
580 
581  default:
582  Assert (false, ExcNotImplemented());
583  }
584 
585  }
586  }
587 
588 
595  template <int spacedim>
596  static
597  void
599  {
600  Assert (dof_handler.fe_collection.size() > 0,
602  Assert (dof_handler.tria->n_levels() > 0,
603  ExcMessage("The current Triangulation must not be empty."));
604  Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
605  ExcInternalError ());
606 
607  reserve_space_release_space (dof_handler);
608 
609  reserve_space_cells (dof_handler);
610  reserve_space_vertices (dof_handler);
611  }
612 
613 
614 
615  template <int spacedim>
616  static
617  void
619  {
620  Assert (dof_handler.fe_collection.size() > 0,
622  Assert (dof_handler.tria->n_levels() > 0,
623  ExcMessage("The current Triangulation must not be empty."));
624  Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
625  ExcInternalError ());
626 
627  reserve_space_release_space (dof_handler);
628 
629  // FIRST CELLS AND FACES
630  reserve_space_cells (dof_handler);
631  reserve_space_faces (dof_handler);
632 
633  // VERTEX DOFS
634  reserve_space_vertices (dof_handler);
635  }
636 
637 
638 
639  template <int spacedim>
640  static
641  void
643  {
644  const unsigned int dim = 3;
645 
646  Assert (dof_handler.fe_collection.size() > 0,
648  Assert (dof_handler.tria->n_levels() > 0,
649  ExcMessage("The current Triangulation must not be empty."));
650  Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
651  ExcInternalError ());
652 
653  reserve_space_release_space (dof_handler);
654 
655  // FIRST CELLS AND FACES
656  reserve_space_cells (dof_handler);
657  reserve_space_faces (dof_handler);
658 
659 
660  // LINE DOFS
661 
662  // the situation here is pretty much like with vertices:
663  // there can be an arbitrary number of finite elements
664  // associated with each line.
665  //
666  // the algorithm we use is somewhat similar to what we do in
667  // reserve_space_vertices()
668  if (true)
669  {
670  // what we do first is to set up an array in which we
671  // record whether a line is associated with any of the
672  // given fe's, by setting a bit. in a later step, we
673  // then actually allocate memory for the required dofs
674  std::vector<std::vector<bool> >
675  line_fe_association (dof_handler.fe_collection.size(),
676  std::vector<bool> (dof_handler.tria->n_raw_lines(),
677  false));
678 
679  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
680  cell=dof_handler.begin_active();
681  cell!=dof_handler.end(); ++cell)
682  for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
683  line_fe_association[cell->active_fe_index()][cell->line_index(l)]
684  = true;
685 
686  // first check which of the lines is used at all,
687  // i.e. is associated with a finite element. we do this
688  // since not all lines may actually be used, in which
689  // case we do not have to allocate any memory at all
690  std::vector<bool> line_is_used (dof_handler.tria->n_raw_lines(), false);
691  for (unsigned int line=0; line<dof_handler.tria->n_raw_lines(); ++line)
692  for (unsigned int fe=0; fe<dof_handler.fe_collection.size(); ++fe)
693  if (line_fe_association[fe][line] == true)
694  {
695  line_is_used[line] = true;
696  break;
697  }
698 
699  // next count how much memory we actually need. for each
700  // line, we need one slot per fe to store the fe_index,
701  // plus dofs_per_line for this fe. in addition, we need
702  // one slot as the end marker for the fe_indices. at the
703  // same time already fill the line_dofs_offsets field
704  dof_handler.faces->lines.dof_offsets
705  .resize (dof_handler.tria->n_raw_lines(),
707 
708  unsigned int line_slots_needed = 0;
709  for (unsigned int line=0; line<dof_handler.tria->n_raw_lines(); ++line)
710  if (line_is_used[line] == true)
711  {
712  dof_handler.faces->lines.dof_offsets[line] = line_slots_needed;
713 
714  for (unsigned int fe=0; fe<dof_handler.fe_collection.size(); ++fe)
715  if (line_fe_association[fe][line] == true)
716  line_slots_needed += dof_handler.get_fe(fe).dofs_per_line + 1;
717  ++line_slots_needed;
718  }
719 
720  // now allocate the space we have determined we need,
721  // and set up the linked lists for each of the lines
722  dof_handler.faces->lines.dofs.resize (line_slots_needed,
724  for (unsigned int line=0; line<dof_handler.tria->n_raw_lines(); ++line)
725  if (line_is_used[line] == true)
726  {
727  unsigned int pointer = dof_handler.faces->lines.dof_offsets[line];
728  for (unsigned int fe=0; fe<dof_handler.fe_collection.size(); ++fe)
729  if (line_fe_association[fe][line] == true)
730  {
731  // if this line uses this fe, then set the
732  // fe_index and move the pointer ahead
733  dof_handler.faces->lines.dofs[pointer] = fe;
734  pointer += dof_handler.get_fe(fe).dofs_per_line + 1;
735  }
736  // finally place the end marker
737  dof_handler.faces->lines.dofs[pointer] = numbers::invalid_dof_index;
738  }
739  }
740 
741 
742  // VERTEX DOFS
743  reserve_space_vertices (dof_handler);
744  }
745 
746 
750  template <int spacedim>
751  static
752  unsigned int
754  {
755  return std::min(static_cast<types::global_dof_index> (3*
756  dof_handler.fe_collection.max_dofs_per_vertex() +
757  2*dof_handler.fe_collection.max_dofs_per_line()),
758  dof_handler.n_dofs());
759  }
760 
761 
762 
763  template <int spacedim>
764  static
765  unsigned int
767  {
768  // get these numbers by drawing pictures and counting...
769  // example:
770  // | | |
771  // --x-----x--x--X--
772  // | | | |
773  // | x--x--x
774  // | | | |
775  // --x--x--*--x--x--
776  // | | | |
777  // x--x--x |
778  // | | | |
779  // --X--x--x-----x--
780  // | | |
781  // x = vertices connected with center vertex *;
782  // = total of 19
783  //
784  // (the X vertices are connected with * if the vertices
785  // adjacent to X are hanging nodes) count lines -> 28 (don't
786  // forget to count mother and children separately!)
787  types::global_dof_index max_couplings;
788  switch (dof_handler.tria->max_adjacent_cells())
789  {
790  case 4:
791  max_couplings=19*dof_handler.fe_collection.max_dofs_per_vertex() +
792  28*dof_handler.fe_collection.max_dofs_per_line() +
793  8*dof_handler.fe_collection.max_dofs_per_quad();
794  break;
795  case 5:
796  max_couplings=21*dof_handler.fe_collection.max_dofs_per_vertex() +
797  31*dof_handler.fe_collection.max_dofs_per_line() +
798  9*dof_handler.fe_collection.max_dofs_per_quad();
799  break;
800  case 6:
801  max_couplings=28*dof_handler.fe_collection.max_dofs_per_vertex() +
802  42*dof_handler.fe_collection.max_dofs_per_line() +
803  12*dof_handler.fe_collection.max_dofs_per_quad();
804  break;
805  case 7:
806  max_couplings=30*dof_handler.fe_collection.max_dofs_per_vertex() +
807  45*dof_handler.fe_collection.max_dofs_per_line() +
808  13*dof_handler.fe_collection.max_dofs_per_quad();
809  break;
810  case 8:
811  max_couplings=37*dof_handler.fe_collection.max_dofs_per_vertex() +
812  56*dof_handler.fe_collection.max_dofs_per_line() +
813  16*dof_handler.fe_collection.max_dofs_per_quad();
814  break;
815  default:
816  Assert (false, ExcNotImplemented());
817  max_couplings=0;
818  };
819  return std::min(max_couplings,dof_handler.n_dofs());
820  }
821 
822 
823  template <int spacedim>
824  static
825  unsigned int
827  {
828 //TODO:[?] Invent significantly better estimates than the ones in this function
829  // doing the same thing here is a rather complicated thing,
830  // compared to the 2d case, since it is hard to draw
831  // pictures with several refined hexahedra :-) so I
832  // presently only give a coarse estimate for the case that
833  // at most 8 hexes meet at each vertex
834  //
835  // can anyone give better estimate here?
836  const unsigned int max_adjacent_cells = dof_handler.tria->max_adjacent_cells();
837 
838  types::global_dof_index max_couplings;
839  if (max_adjacent_cells <= 8)
840  max_couplings=7*7*7*dof_handler.fe_collection.max_dofs_per_vertex() +
841  7*6*7*3*dof_handler.fe_collection.max_dofs_per_line() +
842  9*4*7*3*dof_handler.fe_collection.max_dofs_per_quad() +
843  27*dof_handler.fe_collection.max_dofs_per_hex();
844  else
845  {
846  Assert (false, ExcNotImplemented());
847  max_couplings=0;
848  }
849 
850  return std::min(max_couplings,dof_handler.n_dofs());
851  }
852 
853 
859  template <int dim, int spacedim>
860  static
862  {
864  dynamic_cast<const parallel::shared::Triangulation<dim, spacedim>*> (&dof_handler.get_triangulation()))
865  {
866  // we have a shared triangulation. in this case, every processor
867  // knows about all cells, but every processor only has knowledge
868  // about the active_fe_index on the cells it owns.
869  //
870  // we can create a complete set of active_fe_indices by letting
871  // every processor create a vector of indices for all cells,
872  // filling only those on the cells it owns and setting the indices
873  // on the other cells to zero. then we add all of these vectors
874  // up, and because every vector entry has exactly one processor
875  // that owns it, the sum is correct
876  std::vector<unsigned int> active_fe_indices (tr->n_active_cells(),
877  (unsigned int)0);
878  for (const auto &cell : dof_handler.active_cell_iterators())
879  if (cell->is_locally_owned())
880  active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
881 
882  Utilities::MPI::sum (active_fe_indices,
883  tr->get_communicator(),
884  active_fe_indices);
885 
886  // now go back and fill the active_fe_index on ghost
887  // cells. we would like to call cell->set_active_fe_index(),
888  // but that function does not allow setting these indices on
889  // non-locally_owned cells. so we have to work around the
890  // issue a little bit by accessing the underlying data
891  // structures directly
892  for (auto cell : dof_handler.active_cell_iterators())
893  if (cell->is_ghost())
894  dof_handler.levels[cell->level()]->
895  set_active_fe_index(cell->index(),
896  active_fe_indices[cell->active_cell_index()]);
897  }
899  dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim>*> (&dof_handler.get_triangulation()))
900  {
901  // For completely distributed meshes, use the function that is able to move
902  // data from locally owned cells on one processor to the corresponding
903  // ghost cells on others. To this end, we need to have functions that can pack
904  // and unpack the data we want to transport -- namely, the single unsigned int
905  // active_fe_index objects
906  auto pack
907  = [] (const typename ::hp::DoFHandler<dim,spacedim>::active_cell_iterator &cell) -> unsigned int
908  {
909  return cell->active_fe_index();
910  };
911 
912  auto unpack
913  = [&dof_handler] (const typename ::hp::DoFHandler<dim,spacedim>::active_cell_iterator &cell,
914  const unsigned int &active_fe_index) -> void
915  {
916  // we would like to say
917  // cell->set_active_fe_index(active_fe_index);
918  // but this is not allowed on cells that are not
919  // locally owned, and we are on a ghost cell
920  dof_handler.levels[cell->level()]->
921  set_active_fe_index(cell->index(),
922  active_fe_index);
923  };
924 
925  GridTools::exchange_cell_data_to_ghosts<unsigned int, ::hp::DoFHandler<dim,spacedim>>
926  (dof_handler, pack, unpack);
927  }
928  else
929  {
930  // a sequential triangulation. there is nothing we need to do here
931  Assert ((dynamic_cast<const parallel::Triangulation<dim, spacedim>*> (&dof_handler.get_triangulation())
932  == nullptr),
933  ExcInternalError());
934  }
935  }
936  };
937  }
938  }
939 }
940 
941 
942 namespace hp
943 {
944  template <int dim, int spacedim>
945  const unsigned int DoFHandler<dim,spacedim>::dimension;
946 
947  template <int dim, int spacedim>
949 
950  template <int dim, int spacedim>
952 
953 
954 
955  template <int dim, int spacedim>
957  :
958  tria(nullptr, typeid(*this).name()),
959  faces (nullptr)
960  {}
961 
962 
963  template <int dim, int spacedim>
965  :
966  tria(&tria, typeid(*this).name()),
967  faces (nullptr)
968  {
969  setup_policy_and_listeners ();
970 
971  create_active_fe_table ();
972  }
973 
974 
975  template <int dim, int spacedim>
977  {
978  // unsubscribe as a listener to refinement of the underlying
979  // triangulation
980  for (unsigned int i=0; i<tria_listeners.size(); ++i)
981  tria_listeners[i].disconnect ();
982  tria_listeners.clear ();
983 
984  // ...and release allocated memory
985  // virtual functions called in constructors and destructors never use the
986  // override in a derived class
987  // for clarity be explicit on which function is called
989  }
990 
991 
992  /*------------------------ Cell iterator functions ------------------------*/
993 
994  template <int dim, int spacedim>
996  DoFHandler<dim, spacedim>::begin(const unsigned int level) const
997  {
998  return cell_iterator (*this->get_triangulation().begin(level),
999  this);
1000  }
1001 
1002 
1003 
1004  template <int dim, int spacedim>
1006  DoFHandler<dim,spacedim>::begin_active (const unsigned int level) const
1007  {
1008  // level is checked in begin
1009  cell_iterator i = begin (level);
1010  if (i.state() != IteratorState::valid)
1011  return i;
1012  while (i->has_children())
1013  if ((++i).state() != IteratorState::valid)
1014  return i;
1015  return i;
1016  }
1017 
1018 
1019 
1020  template <int dim, int spacedim>
1023  {
1024  return cell_iterator (&this->get_triangulation(),
1025  -1,
1026  -1,
1027  this);
1028  }
1029 
1030 
1031  template <int dim, int spacedim>
1033  DoFHandler<dim,spacedim>::end (const unsigned int level) const
1034  {
1035  return (level == this->get_triangulation().n_levels()-1 ?
1036  end() :
1037  begin (level+1));
1038  }
1039 
1040 
1041  template <int dim, int spacedim>
1043  DoFHandler<dim, spacedim>::end_active (const unsigned int level) const
1044  {
1045  return (level == this->get_triangulation().n_levels()-1 ?
1046  active_cell_iterator(end()) :
1047  begin_active (level+1));
1048  }
1049 
1050 
1051 
1052  template <int dim, int spacedim>
1055  {
1056  return
1058  (begin(), end());
1059  }
1060 
1061 
1062  template <int dim, int spacedim>
1065  {
1066  return
1068  (begin_active(), end());
1069  }
1070 
1071 
1072 
1073  template <int dim, int spacedim>
1076  {
1077  return
1079  (begin(level), end(level));
1080  }
1081 
1082 
1083 
1084  template <int dim, int spacedim>
1087  {
1088  return
1090  (begin_active(level), end_active(level));
1091  }
1092 
1093 
1094 
1095 
1096 //------------------------------------------------------------------
1097 
1098 
1099  template <int dim, int spacedim>
1101  {
1102  Assert(fe_collection.size() > 0, ExcNoFESelected());
1103 
1104  std::set<types::global_dof_index> boundary_dofs;
1105  std::vector<types::global_dof_index> dofs_on_face;
1106  dofs_on_face.reserve (this->get_fe_collection().max_dofs_per_face());
1107 
1108  // loop over all faces to check whether they are at a
1109  // boundary. note that we need not take special care of single
1110  // lines in 3d (using @p{cell->has_boundary_lines}), since we do
1111  // not support boundaries of dimension dim-2, and so every
1112  // boundary line is also part of a boundary face.
1113  typename HpDoFHandler<dim,spacedim>::active_cell_iterator cell = this->begin_active (),
1114  endc = this->end();
1115  for (; cell!=endc; ++cell)
1116  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
1117  if (cell->at_boundary(f))
1118  {
1119  const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
1120  dofs_on_face.resize (dofs_per_face);
1121 
1122  cell->face(f)->get_dof_indices (dofs_on_face,
1123  cell->active_fe_index());
1124  for (unsigned int i=0; i<dofs_per_face; ++i)
1125  boundary_dofs.insert(dofs_on_face[i]);
1126  }
1127  return boundary_dofs.size();
1128  }
1129 
1130 
1131 
1132  template <int dim, int spacedim>
1134  DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_ids) const
1135  {
1136  Assert(fe_collection.size() > 0, ExcNoFESelected());
1137  Assert (boundary_ids.find (numbers::internal_face_boundary_id) == boundary_ids.end(),
1139 
1140  // same as above, but with additional checks for set of boundary
1141  // indicators
1142  std::set<types::global_dof_index> boundary_dofs;
1143  std::vector<types::global_dof_index> dofs_on_face;
1144  dofs_on_face.reserve (this->get_fe_collection().max_dofs_per_face());
1145 
1146  typename HpDoFHandler<dim,spacedim>::active_cell_iterator cell = this->begin_active (),
1147  endc = this->end();
1148  for (; cell!=endc; ++cell)
1149  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
1150  if (cell->at_boundary(f) &&
1151  (boundary_ids.find(cell->face(f)->boundary_id()) !=
1152  boundary_ids.end()))
1153  {
1154  const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
1155  dofs_on_face.resize (dofs_per_face);
1156 
1157  cell->face(f)->get_dof_indices (dofs_on_face,
1158  cell->active_fe_index());
1159  for (unsigned int i=0; i<dofs_per_face; ++i)
1160  boundary_dofs.insert(dofs_on_face[i]);
1161  }
1162  return boundary_dofs.size();
1163  }
1164 
1165 
1166 
1167  template <>
1169  {
1170  Assert(false,ExcNotImplemented());
1171  return 0;
1172  }
1173 
1174 
1175 
1176  template <>
1177  template <typename number>
1179  {
1180  Assert(false,ExcNotImplemented());
1181  return 0;
1182  }
1183 
1184 
1185 
1186  template <>
1187  types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
1188  {
1189  Assert(false,ExcNotImplemented());
1190  return 0;
1191  }
1192 
1193 
1194 
1195  template <int dim, int spacedim>
1196  std::size_t
1198  {
1199  std::size_t mem = (MemoryConsumption::memory_consumption (tria) +
1200  MemoryConsumption::memory_consumption (fe_collection) +
1204  MemoryConsumption::memory_consumption (number_cache) +
1206  MemoryConsumption::memory_consumption (vertex_dof_offsets) +
1207  MemoryConsumption::memory_consumption (has_children));
1208  for (unsigned int i=0; i<levels.size(); ++i)
1209  mem += MemoryConsumption::memory_consumption (*levels[i]);
1210  mem += MemoryConsumption::memory_consumption (*faces);
1211 
1212  return mem;
1213  }
1214 
1215 
1216 
1217 
1218 
1219 
1220  template <int dim, int spacedim>
1221  void DoFHandler<dim,spacedim>::set_active_fe_indices (const std::vector<unsigned int> &active_fe_indices)
1222  {
1223  Assert(active_fe_indices.size()==get_triangulation().n_active_cells(),
1224  ExcDimensionMismatch(active_fe_indices.size(), get_triangulation().n_active_cells()));
1225 
1226  create_active_fe_table ();
1227  // we could set the values directly, since they are stored as
1228  // protected data of this object, but for simplicity we use the
1229  // cell-wise access. this way we also have to pass some debug-mode
1230  // tests which we would have to duplicate ourselves otherwise
1231  active_cell_iterator cell=begin_active(),
1232  endc=end();
1233  for (unsigned int i=0; cell!=endc; ++cell, ++i)
1234  if (cell->is_locally_owned())
1235  cell->set_active_fe_index(active_fe_indices[i]);
1236  }
1237 
1238 
1239 
1240  template <int dim, int spacedim>
1241  void DoFHandler<dim,spacedim>::get_active_fe_indices (std::vector<unsigned int> &active_fe_indices) const
1242  {
1243  active_fe_indices.resize(get_triangulation().n_active_cells());
1244 
1245  // we could try to extract the values directly, since they are
1246  // stored as protected data of this object, but for simplicity we
1247  // use the cell-wise access.
1248  active_cell_iterator cell=begin_active(),
1249  endc=end();
1250  for (unsigned int i=0; cell!=endc; ++cell, ++i)
1251  active_fe_indices[i] = cell->active_fe_index();
1252  }
1253 
1254  template <int dim, int spacedim>
1257  {
1258  if (this->tria != &tria)
1259  {
1260  for (unsigned int i=0; i<tria_listeners.size(); ++i)
1261  tria_listeners[i].disconnect ();
1262  tria_listeners.clear ();
1263 
1264  this->tria = &tria;
1265 
1266  setup_policy_and_listeners ();
1267  }
1268 
1269  create_active_fe_table ();
1270 
1271  distribute_dofs (fe);
1272  }
1273 
1274 
1275  template <int dim, int spacedim>
1277  {
1278  Assert(tria!=nullptr,
1279  ExcMessage("You need to set the Triangulation in the DoFHandler using initialize() or "
1280  "in the constructor before you can distribute DoFs."));
1281  Assert (tria->n_levels() > 0,
1282  ExcMessage("The Triangulation you are using is empty!"));
1283  Assert (ff.size() > 0,
1284  ExcMessage("The hp::FECollection given is empty!"));
1285 
1286  // don't create a new object if the one we have is already appropriate
1287  if (fe_collection != ff)
1288  fe_collection = hp::FECollection<dim,spacedim>(ff);
1289 
1290  // at the beginning, make sure every processor knows the
1291  // active_fe_indices on both its own cells and all ghost cells
1293 
1294  // If an underlying shared::Tria allows artificial cells,
1295  // then save the current set of subdomain ids, and set
1296  // subdomain ids to the "true" owner of each cell. we later
1297  // restore these flags
1298  std::vector<types::subdomain_id> saved_subdomain_ids;
1299  if (const parallel::shared::Triangulation<dim, spacedim> *shared_tria =
1300  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim>*> (&get_triangulation())))
1301  if (shared_tria->with_artificial_cells())
1302  {
1303  saved_subdomain_ids.resize (shared_tria->n_active_cells());
1304 
1305  typename parallel::shared::Triangulation<dim,spacedim>::active_cell_iterator
1306  cell = get_triangulation().begin_active(),
1307  endc = get_triangulation().end();
1308 
1309  const std::vector<types::subdomain_id> &true_subdomain_ids
1310  = shared_tria->get_true_subdomain_ids_of_cells();
1311 
1312  for (unsigned int index=0; cell != endc; ++cell, ++index)
1313  {
1314  saved_subdomain_ids[index] = cell->subdomain_id();
1315  cell->set_subdomain_id(true_subdomain_ids[index]);
1316  }
1317  }
1318 
1319  // ensure that the active_fe_indices vectors are
1320  // initialized correctly.
1321  create_active_fe_table ();
1322 
1323  // up front make sure that the fe collection is large enough to
1324  // cover all fe indices presently in use on the mesh
1325  for (active_cell_iterator cell = begin_active(); cell != end(); ++cell)
1326  if (cell->is_locally_owned())
1327  Assert (cell->active_fe_index() < fe_collection.size(),
1328  ExcInvalidFEIndex (cell->active_fe_index(),
1329  fe_collection.size()));
1330 
1331 
1332  // then allocate space for all the other tables
1334 
1335 
1336  // now undo the subdomain modification
1337  if (const parallel::shared::Triangulation<dim, spacedim> *shared_tria =
1338  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim>*> (&get_triangulation())))
1339  if (shared_tria->with_artificial_cells())
1340  {
1341  typename parallel::shared::Triangulation<dim,spacedim>::active_cell_iterator
1342  cell = get_triangulation().begin_active(),
1343  endc = get_triangulation().end();
1344 
1345  for (unsigned int index=0; cell != endc; ++cell, ++index)
1346  cell->set_subdomain_id(saved_subdomain_ids[index]);
1347  }
1348 
1349 
1350  // Clear user flags because we will need them. But first we save
1351  // them and make sure that we restore them later such that at the
1352  // end of this function the Triangulation will be in the same
1353  // state as it was at the beginning of this function.
1354  std::vector<bool> user_flags;
1355  tria->save_user_flags(user_flags);
1356  const_cast<Triangulation<dim,spacedim> &>(*tria).clear_user_flags ();
1357 
1358 
1360 
1361  // Now for the real work:
1362  number_cache = policy->distribute_dofs ();
1363 
1365 
1366  // do some housekeeping: compress indices
1367  {
1369  for (int level=levels.size()-1; level>=0; --level)
1370  tg += Threads::new_task (&::internal::hp::DoFLevel::compress_data<dim,spacedim>,
1371  *levels[level], fe_collection);
1372  tg.join_all ();
1373  }
1374 
1375  // finally restore the user flags
1376  const_cast<Triangulation<dim,spacedim> &>(*tria).load_user_flags(user_flags);
1377  }
1378 
1379 
1380  template <int dim, int spacedim>
1382  {
1383  // decide whether we need a sequential or a parallel shared/distributed policy
1384  if (dynamic_cast<const parallel::shared::Triangulation< dim, spacedim>*> (&*this->tria) != nullptr)
1385  policy = std_cxx14::make_unique<internal::DoFHandlerImplementation::Policy::ParallelShared<DoFHandler<dim,spacedim> >> (*this);
1386  else if (dynamic_cast<const parallel::distributed::Triangulation< dim, spacedim >*> (&*this->tria) != nullptr)
1387  policy = std_cxx14::make_unique<internal::DoFHandlerImplementation::Policy::ParallelDistributed<DoFHandler<dim,spacedim> >> (*this);
1388  else
1389  policy = std_cxx14::make_unique<internal::DoFHandlerImplementation::Policy::Sequential<DoFHandler<dim,spacedim> >> (*this);
1390 
1391  tria_listeners.push_back
1392  (this->tria->signals.pre_refinement
1394  std::ref(*this))));
1395  tria_listeners.push_back
1396  (this->tria->signals.post_refinement
1398  std::ref(*this))));
1399  tria_listeners.push_back
1400  (this->tria->signals.create
1402  std::ref(*this))));
1403  }
1404 
1405 
1406  template <int dim, int spacedim>
1408  {
1409  // release memory
1410  clear_space ();
1411  }
1412 
1413 
1414 
1415  template <int dim, int spacedim>
1416  void DoFHandler<dim,spacedim>::renumber_dofs (const std::vector<types::global_dof_index> &new_numbers)
1417  {
1418  Assert(levels.size()>0, ExcMessage("You need to distribute DoFs before you can renumber them."));
1419 
1420  AssertDimension (new_numbers.size(), n_locally_owned_dofs());
1421 
1422 #ifdef DEBUG
1423  // assert that the new indices are
1424  // consecutively numbered if we are
1425  // working on a single
1426  // processor. this doesn't need to
1427  // hold in the case of a parallel
1428  // mesh since we map the interval
1429  // [0...n_dofs()) into itself but
1430  // only globally, not on each
1431  // processor
1432  if (n_locally_owned_dofs() == n_dofs())
1433  {
1434  std::vector<types::global_dof_index> tmp(new_numbers);
1435  std::sort (tmp.begin(), tmp.end());
1436  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
1438  for (; p!=tmp.end(); ++p, ++i)
1439  Assert (*p == i, ExcNewNumbersNotConsecutive(i));
1440  }
1441  else
1442  for (types::global_dof_index i=0; i<new_numbers.size(); ++i)
1443  Assert (new_numbers[i] < n_dofs(),
1444  ExcMessage ("New DoF index is not less than the total number of dofs."));
1445 #endif
1446 
1447  // uncompress the internal storage scheme of dofs on cells so that
1448  // we can access dofs in turns. uncompress in parallel, starting
1449  // with the most expensive levels (the highest ones)
1450  {
1452  for (int level=levels.size()-1; level>=0; --level)
1453  tg += Threads::new_task (&::internal::hp::DoFLevel::uncompress_data<dim,spacedim>,
1454  *levels[level], fe_collection);
1455  tg.join_all ();
1456  }
1457 
1458  // do the renumbering
1459  number_cache = policy->renumber_dofs(new_numbers);
1460 
1461  // now re-compress the dof indices
1462  {
1464  for (int level=levels.size()-1; level>=0; --level)
1465  tg += Threads::new_task (&::internal::hp::DoFLevel::compress_data<dim,spacedim>,
1466  *levels[level], fe_collection);
1467  tg.join_all ();
1468  }
1469  }
1470 
1471 
1472 
1473  template <int dim, int spacedim>
1474  unsigned int
1476  {
1477  Assert(fe_collection.size() > 0, ExcNoFESelected());
1478  return ::internal::hp::DoFHandlerImplementation::Implementation::max_couplings_between_dofs (*this);
1479  }
1480 
1481 
1482 
1483  template <int dim, int spacedim>
1484  unsigned int
1486  {
1487  Assert (fe_collection.size() > 0, ExcNoFESelected());
1488 
1489  switch (dim)
1490  {
1491  case 1:
1492  return fe_collection.max_dofs_per_vertex();
1493  case 2:
1494  return (3*fe_collection.max_dofs_per_vertex()
1495  +
1496  2*fe_collection.max_dofs_per_line());
1497  case 3:
1498  // we need to take refinement of one boundary face into
1499  // consideration here; in fact, this function returns what
1500  // #max_coupling_between_dofs<2> returns
1501  //
1502  // we assume here, that only four faces meet at the boundary;
1503  // this assumption is not justified and needs to be fixed some
1504  // time. fortunately, omitting it for now does no harm since
1505  // the matrix will cry foul if its requirements are not
1506  // satisfied
1507  return (19*fe_collection.max_dofs_per_vertex() +
1508  28*fe_collection.max_dofs_per_line() +
1509  8*fe_collection.max_dofs_per_quad());
1510  default:
1511  Assert (false, ExcNotImplemented());
1512  return 0;
1513  }
1514  }
1515 
1516 
1517 
1518  template <int dim, int spacedim>
1520  {
1521  // Create sufficiently many hp::DoFLevels.
1522  while (levels.size () < tria->n_levels ())
1523  levels.emplace_back (new ::internal::hp::DoFLevel);
1524 
1525  // then make sure that on each level we have the appropriate size
1526  // of active_fe_indices; preset them to zero, i.e. the default FE
1527  for (unsigned int level=0; level<levels.size(); ++level)
1528  {
1529  if (levels[level]->active_fe_indices.size () == 0)
1530  levels[level]->active_fe_indices.resize (tria->n_raw_cells(level),
1531  0);
1532  else
1533  {
1534  // Either the active_fe_indices have size zero because
1535  // they were just created, or the correct size. Other
1536  // sizes indicate that something went wrong.
1537  Assert (levels[level]->active_fe_indices.size () ==
1538  tria->n_raw_cells(level),
1539  ExcInternalError ());
1540  }
1541 
1542  // it may be that the previous table was compressed; in that
1543  // case, restore the correct active_fe_index. the fact that
1544  // this no longer matches the indices in the table is of no
1545  // importance because the current function is called at a
1546  // point where we have to recreate the dof_indices tables in
1547  // the levels anyway
1548  levels[level]->normalize_active_fe_indices ();
1549  }
1550  }
1551 
1552 
1553  template <int dim, int spacedim>
1555  {
1556  create_active_fe_table ();
1557 
1558  // Remember if the cells already have children. That will make the
1559  // transfer of the active_fe_index to the finer levels easier.
1560  Assert (has_children.size () == 0, ExcInternalError ());
1561  for (unsigned int i=0; i<levels.size(); ++i)
1562  {
1563  const unsigned int cells_on_level = tria->n_raw_cells(i);
1564  std::unique_ptr<std::vector<bool> > has_children_level(new std::vector<bool> (cells_on_level));
1565 
1566  // Check for each cell, if it has children. in 1d, we don't
1567  // store refinement cases, so use the 'children' vector
1568  // instead
1569  if (dim == 1)
1570  std::transform (tria->levels[i]->cells.children.begin (),
1571  tria->levels[i]->cells.children.end (),
1572  has_children_level->begin (),
1573  std::bind(std::not_equal_to<int>(),
1574  std::placeholders::_1,
1575  -1));
1576  else
1577  std::transform (tria->levels[i]->cells.refinement_cases.begin (),
1578  tria->levels[i]->cells.refinement_cases.end (),
1579  has_children_level->begin (),
1580  std::bind (std::not_equal_to<unsigned char>(),
1581  std::placeholders::_1,
1582  static_cast<unsigned char>(RefinementCase<dim>::no_refinement)));
1583 
1584  has_children.emplace_back (std::move(has_children_level));
1585  }
1586  }
1587 
1588 
1589 
1590  template <int dim, int spacedim>
1591  void
1593  {
1594  Assert (has_children.size () == levels.size (), ExcInternalError ());
1595 
1596  // Normally only one level is added, but if this Triangulation
1597  // is created by copy_triangulation, it can be more than one level.
1598  while (levels.size () < tria->n_levels ())
1599  levels.emplace_back (new ::internal::hp::DoFLevel);
1600 
1601  // Coarsening can lead to the loss of levels. Hence remove them.
1602  while (levels.size () > tria->n_levels ())
1603  {
1604  // drop the last element. that also releases the memory pointed to
1605  levels.pop_back ();
1606  }
1607 
1608  Assert(levels.size () == tria->n_levels (), ExcInternalError());
1609 
1610  // Resize active_fe_indices vectors. use zero indicator to extend
1611  for (unsigned int i=0; i<levels.size(); ++i)
1612  levels[i]->active_fe_indices.resize (tria->n_raw_cells(i), 0);
1613 
1614  // if a finite element collection has already been set, then
1615  // actually try to set active_fe_indices for child cells of
1616  // refined cells to the active_fe_index of the mother cell. if no
1617  // finite element collection has been assigned yet, then all
1618  // indicators are zero anyway, and there is no point trying to set
1619  // anything (besides, we would trip over an assertion in
1620  // set_active_fe_index)
1621  if (fe_collection.size() > 0)
1622  {
1623  cell_iterator cell = begin(),
1624  endc = end ();
1625  for (; cell != endc; ++cell)
1626  {
1627  // Look if the cell got children during refinement by
1628  // checking whether it has children now but didn't have
1629  // children before refinement (the has_children array is
1630  // set in pre-refinement action)
1631  //
1632  // Note: Although one level is added to the DoFHandler
1633  // levels, when the triangulation got one, for the buffer
1634  // has_children this new level is not required, because
1635  // the cells on the finest level never have
1636  // children. Hence cell->has_children () will always
1637  // return false on that level, which would cause shortcut
1638  // evaluation of the following expression. Thus an index
1639  // error in has_children should never occur.
1640  if (cell->has_children () &&
1641  !(*has_children [cell->level ()])[cell->index ()])
1642  {
1643  // Set active_fe_index in children to the same value
1644  // as in the parent cell. we can't access the
1645  // active_fe_index in the parent cell any more through
1646  // cell->active_fe_index() since that function is not
1647  // allowed for inactive cells, but we can access this
1648  // information from the DoFLevels directly
1649  //
1650  // we don't have to set the active_fe_index for ghost
1651  // cells -- these will be exchanged automatically
1652  // upon distribute_dofs()
1653  for (unsigned int i = 0; i < cell->n_children(); ++i)
1654  if (cell->child (i)->is_locally_owned())
1655  cell->child (i)->set_active_fe_index
1656  (levels[cell->level()]->active_fe_index (cell->index()));
1657  }
1658  }
1659  }
1660 
1661  // Free buffer objects
1662  has_children.clear ();
1663  }
1664 
1665 
1666  template <int dim, int spacedim>
1667  template <int structdim>
1669  DoFHandler<dim,spacedim>::get_dof_index (const unsigned int,
1670  const unsigned int,
1671  const unsigned int,
1672  const unsigned int) const
1673  {
1674  Assert (false, ExcNotImplemented());
1676  }
1677 
1678 
1679  template <int dim, int spacedim>
1680  template <int structdim>
1681  void
1682  DoFHandler<dim,spacedim>::set_dof_index (const unsigned int,
1683  const unsigned int,
1684  const unsigned int,
1685  const unsigned int,
1686  const types::global_dof_index) const
1687  {
1688  Assert (false, ExcNotImplemented());
1689  }
1690 
1691 
1692  template <int dim, int spacedim>
1694  {
1695  levels.clear ();
1696  faces.reset ();
1697 
1698  vertex_dofs = std::move(std::vector<types::global_dof_index>());
1699  vertex_dof_offsets = std::move (std::vector<unsigned int>());
1700  }
1701 }
1702 
1703 
1704 
1705 /*-------------- Explicit Instantiations -------------------------------*/
1706 #include "dof_handler.inst"
1707 
1708 
1709 DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
Definition: dof_handler.cc:864
unsigned int max_couplings_between_dofs() const
static const unsigned int invalid_unsigned_int
Definition: types.h:173
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1248
void clear_user_flags()
Definition: tria.cc:9729
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:723
unsigned int offset_type
Definition: dof_level.h:109
Task< RT > new_task(const std::function< RT()> &function)
std::vector< std::unique_ptr<::internal::hp::DoFLevel > > levels
Definition: dof_handler.h:886
void load_user_flags(std::istream &in)
Definition: tria.cc:9786
unsigned int boundary_id
Definition: types.h:110
cell_iterator end() const
Definition: dof_handler.cc:751
virtual void clear()
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > levels
Definition: dof_handler.h:1173
unsigned int max_dofs_per_face(const DoFHandler< dim, spacedim > &dh)
SmartPointer< const Triangulation< dim, spacedim >, DoFHandler< dim, spacedim > > tria
Definition: dof_handler.h:1032
unsigned int size() const
virtual std::size_t memory_consumption() const
Definition: dof_handler.cc:964
static void reserve_space_cells(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:224
void initialize(const Triangulation< dim, spacedim > &tria, const FiniteElement< dim, spacedim > &fe)
Definition: dof_handler.cc:699
ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:186
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:735
static ::ExceptionBase & ExcMessage(std::string arg1)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:88
#define Assert(cond, exc)
Definition: exceptions.h:1142
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:110
active_cell_iterator end_active(const unsigned int level) const
Definition: dof_handler.cc:773
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
types::global_dof_index n_boundary_dofs() const
Definition: dof_handler.cc:889
static void reserve_space_release_space(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:86
types::global_dof_index n_dofs() const
ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:198
static unsigned int max_couplings_between_dofs(const DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:753
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > faces
Definition: dof_handler.h:1182
std::vector< unsigned int > dof_offsets
Definition: dof_faces.h:104
Definition: hp.h:102
virtual ~DoFHandler()
Definition: dof_handler.cc:681
internal::hp::DoFIndicesOnFacesOrEdges< 2 > quads
Definition: dof_faces.h:313
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
const Triangulation< dim, spacedim > & get_triangulation() const
void clear_space()
virtual void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
Definition: dof_handler.cc:993
static void reserve_space_vertices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:120
void save_user_flags(std::ostream &out) const
Definition: tria.cc:9739
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
Definition: dof_handler.cc:853
static ::ExceptionBase & ExcNotImplemented()
Iterator points to a valid object.
const types::boundary_id internal_face_boundary_id
Definition: types.h:219
unsigned int max_couplings_between_boundary_dofs() const
static void communicate_active_fe_indices(::hp::DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:861
static void reserve_space(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:598
const types::global_dof_index invalid_dof_index
Definition: types.h:187
static ::ExceptionBase & ExcNoFESelected()
IteratorRange< cell_iterator > cell_iterators() const
Definition: dof_handler.cc:820
static void reserve_space_faces(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:320
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:1161
static ::ExceptionBase & ExcInternalError()
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1039