Reference documentation for deal.II version 8.5.1
dof_handler.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2017 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/thread_management.h>
19 #include <deal.II/base/std_cxx11/bind.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/grid/tria_accessor.h>
25 #include <deal.II/grid/tria_iterator.h>
26 #include <deal.II/grid/tria_levels.h>
27 #include <deal.II/grid/tria.h>
28 #include <deal.II/fe/fe.h>
29 #include <deal.II/distributed/shared_tria.h>
30 #include <deal.II/distributed/tria.h>
31 
32 #include <set>
33 #include <algorithm>
34 #include <functional>
35 
36 DEAL_II_NAMESPACE_OPEN
37 
38 // The following is necessary for compilation under Visual Studio which is unable to correctly
39 // distinguish between ::DoFHandler and ::hp::DoFHandler.
40 // Plus it makes code in dof_handler.cc easier to read.
41 // Requires C++11 support which is in Visual Studio 2013 and newer.
42 #if defined(_MSC_VER) && (_MSC_VER >= 1800)
43 template <int dim, int spacedim> using HpDoFHandler = ::hp::DoFHandler<dim, spacedim>;
44 #else
45 // When using older Visual Studio or a different compiler just fall back.
46 #define HpDoFHandler DoFHandler
47 #endif
48 
49 namespace parallel
50 {
51  namespace distributed
52  {
53  template <int, int> class Triangulation;
54  }
55 }
56 
57 
58 namespace internal
59 {
60  namespace hp
61  {
62  typedef
63  std::vector<std::pair<unsigned int, unsigned int> > DoFIdentities;
64 
65 
81  template <int structdim, int dim, int spacedim>
82  void
83  ensure_existence_of_dof_identities (const FiniteElement<dim,spacedim> &fe1,
84  const FiniteElement<dim,spacedim> &fe2,
85  std_cxx11::shared_ptr<DoFIdentities> &identities)
86  {
87  // see if we need to fill this
88  // entry, or whether it already
89  // exists
90  if (identities.get() == 0)
91  {
92  switch (structdim)
93  {
94  case 0:
95  {
96  identities =
97  std_cxx11::shared_ptr<DoFIdentities>
98  (new DoFIdentities(fe1.hp_vertex_dof_identities(fe2)));
99  break;
100  }
101 
102  case 1:
103  {
104  identities =
105  std_cxx11::shared_ptr<DoFIdentities>
106  (new DoFIdentities(fe1.hp_line_dof_identities(fe2)));
107  break;
108  }
109 
110  case 2:
111  {
112  identities =
113  std_cxx11::shared_ptr<DoFIdentities>
114  (new DoFIdentities(fe1.hp_quad_dof_identities(fe2)));
115  break;
116  }
117 
118  default:
119  Assert (false, ExcNotImplemented());
120  }
121 
122  // double check whether the
123  // newly created entries
124  // make any sense at all
125  for (unsigned int i=0; i<identities->size(); ++i)
126  {
127  Assert ((*identities)[i].first < fe1.template n_dofs_per_object<structdim>(),
128  ExcInternalError());
129  Assert ((*identities)[i].second < fe2.template n_dofs_per_object<structdim>(),
130  ExcInternalError());
131  }
132  }
133  }
134 
135 
136 
146  template <int dim, int spacedim, typename iterator>
147  unsigned int
148  get_most_dominating_fe_index (const iterator &object)
149  {
150  unsigned int dominating_fe_index = 0;
151  for (; dominating_fe_index<object->n_active_fe_indices();
152  ++dominating_fe_index)
153  {
154  const FiniteElement<dim, spacedim> &this_fe
155  = object->get_fe (object->nth_active_fe_index(dominating_fe_index));
156 
157  FiniteElementDomination::Domination
158  domination = FiniteElementDomination::either_element_can_dominate;
159  for (unsigned int other_fe_index=0;
160  other_fe_index<object->n_active_fe_indices();
161  ++other_fe_index)
162  if (other_fe_index != dominating_fe_index)
163  {
164  const FiniteElement<dim, spacedim>
165  &that_fe
166  = object->get_fe (object->nth_active_fe_index(other_fe_index));
167 
168  domination = domination &
169  this_fe.compare_for_face_domination(that_fe);
170  }
171 
172  // see if this element is
173  // able to dominate all the
174  // other ones, and if so
175  // take it
176  if ((domination == FiniteElementDomination::this_element_dominates)
177  ||
178  (domination == FiniteElementDomination::either_element_can_dominate)
179  ||
180  (domination == FiniteElementDomination::no_requirements))
181  break;
182  }
183 
184  // check that we have
185  // found one such fe
186  if (dominating_fe_index != object->n_active_fe_indices())
187  {
188  // return the finite element
189  // index used on it. note
190  // that only a single fe can
191  // be active on such subfaces
192  return object->nth_active_fe_index(dominating_fe_index);
193  }
194  else
195  {
196  // if we couldn't find the most dominating object
197  return numbers::invalid_unsigned_int;
198  }
199  }
200  }
201 }
202 
203 
204 
205 namespace internal
206 {
207  namespace hp
208  {
209  namespace DoFHandler
210  {
211  // access class
212  // ::hp::DoFHandler instead of
213  // namespace internal::hp::DoFHandler, etc
214  using ::hp::DoFHandler;
215 
220  struct Implementation
221  {
229  template<int dim, int spacedim>
230  static
231  void
232  reserve_space_vertices (DoFHandler<dim,spacedim> &dof_handler)
233  {
234  // The final step is allocating
235  // memory is to set up vertex dof
236  // information. since vertices
237  // are sequentially numbered,
238  // what we do first is to set up
239  // an array in which we record
240  // whether a vertex is associated
241  // with any of the given fe's, by
242  // setting a bit. in a later
243  // step, we then actually
244  // allocate memory for the
245  // required dofs
246  std::vector<std::vector<bool> >
247  vertex_fe_association (dof_handler.finite_elements->size(),
248  std::vector<bool> (dof_handler.tria->n_vertices(), false));
249 
250  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
251  cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
252  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
253  vertex_fe_association[cell->active_fe_index()][cell->vertex_index(v)]
254  = true;
255 
256  // in debug mode, make sure
257  // that each vertex is
258  // associated with at least one
259  // fe (note that except for
260  // unused vertices, all
261  // vertices are actually
262  // active)
263 #ifdef DEBUG
264  for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
265  if (dof_handler.tria->vertex_used(v) == true)
266  {
267  unsigned int fe=0;
268  for (; fe<dof_handler.finite_elements->size(); ++fe)
269  if (vertex_fe_association[fe][v] == true)
270  break;
271  Assert (fe != dof_handler.finite_elements->size(), ExcInternalError());
272  }
273 #endif
274 
275  // next count how much memory
276  // we actually need. for each
277  // vertex, we need one slot per
278  // fe to store the fe_index,
279  // plus dofs_per_vertex for
280  // this fe. in addition, we
281  // need one slot as the end
282  // marker for the
283  // fe_indices. at the same time
284  // already fill the
285  // vertex_dofs_offsets field
286  dof_handler.vertex_dofs_offsets.resize (dof_handler.tria->n_vertices(),
287  numbers::invalid_dof_index);
288 
289  unsigned int vertex_slots_needed = 0;
290  for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
291  if (dof_handler.tria->vertex_used(v) == true)
292  {
293  dof_handler.vertex_dofs_offsets[v] = vertex_slots_needed;
294 
295  for (unsigned int fe=0; fe<dof_handler.finite_elements->size(); ++fe)
296  if (vertex_fe_association[fe][v] == true)
297  vertex_slots_needed += (*dof_handler.finite_elements)[fe].dofs_per_vertex + 1;
298  ++vertex_slots_needed;
299  }
300 
301  // now allocate the space we
302  // have determined we need, and
303  // set up the linked lists for
304  // each of the vertices
305  dof_handler.vertex_dofs.resize (vertex_slots_needed,
306  DoFHandler<dim,spacedim>::invalid_dof_index);
307  for (unsigned int v=0; v<dof_handler.tria->n_vertices(); ++v)
308  if (dof_handler.tria->vertex_used(v) == true)
309  {
310  types::global_dof_index pointer = dof_handler.vertex_dofs_offsets[v];
311  for (unsigned int fe=0; fe<dof_handler.finite_elements->size(); ++fe)
312  if (vertex_fe_association[fe][v] == true)
313  {
314  // if this vertex
315  // uses this fe,
316  // then set the
317  // fe_index and
318  // move the pointer
319  // ahead
320  dof_handler.vertex_dofs[pointer] = fe;
321  pointer += (*dof_handler.finite_elements)[fe].dofs_per_vertex + 1;
322  }
323  // finally place the end
324  // marker
325  dof_handler.vertex_dofs[pointer] = numbers::invalid_dof_index;
326  }
327  }
328 
329 
330 
345  template <int spacedim>
346  static
347  types::global_dof_index
348  distribute_dofs_on_cell (const typename ::hp::DoFHandler<1,spacedim>::active_cell_iterator &cell,
349  types::global_dof_index next_free_dof)
350  {
351  const unsigned int dim = 1;
352 
353  const FiniteElement<dim,spacedim> &fe = cell->get_fe();
354  const unsigned int fe_index = cell->active_fe_index ();
355 
356  // number dofs on vertices. to do
357  // so, check whether dofs for
358  // this vertex have been
359  // distributed and for the
360  // present fe (only check the
361  // first dof), and if this isn't
362  // the case distribute new ones
363  // there
364  if (fe.dofs_per_vertex > 0)
365  for (unsigned int vertex=0; vertex<GeometryInfo<1>::vertices_per_cell; ++vertex)
366  if (cell->vertex_dof_index(vertex, 0, fe_index) ==
367  DoFHandler<dim,spacedim>::invalid_dof_index)
368  for (unsigned int d=0; d<fe.dofs_per_vertex; ++d, ++next_free_dof)
369  cell->set_vertex_dof_index (vertex, d, next_free_dof, fe_index);
370 
371  // finally for the line. this one
372  // shouldn't be numbered yet
373  if (fe.dofs_per_line > 0)
374  {
375  Assert ((cell->dof_index(0, fe_index) ==
376  DoFHandler<dim,spacedim>::invalid_dof_index),
377  ExcInternalError());
378 
379  for (unsigned int d=0; d<fe.dofs_per_line; ++d, ++next_free_dof)
380  cell->set_dof_index (d, next_free_dof, fe_index);
381  }
382 
383  // note that this cell has been processed
384  cell->set_user_flag ();
385 
386  return next_free_dof;
387  }
388 
389 
390  template <int spacedim>
391  static
392  types::global_dof_index
393  distribute_dofs_on_cell (const typename ::hp::DoFHandler<2,spacedim>::active_cell_iterator &cell,
394  types::global_dof_index next_free_dof)
395  {
396  const unsigned int dim = 2;
397 
398  const FiniteElement<dim,spacedim> &fe = cell->get_fe();
399  const unsigned int fe_index = cell->active_fe_index ();
400 
401  // number dofs on vertices. to do
402  // so, check whether dofs for
403  // this vertex have been
404  // distributed and for the
405  // present fe (only check the
406  // first dof), and if this isn't
407  // the case distribute new ones
408  // there
409  if (fe.dofs_per_vertex > 0)
410  for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_cell; ++vertex)
411  if (cell->vertex_dof_index(vertex, 0, fe_index) ==
412  DoFHandler<dim,spacedim>::invalid_dof_index)
413  for (unsigned int d=0; d<fe.dofs_per_vertex; ++d, ++next_free_dof)
414  cell->set_vertex_dof_index (vertex, d, next_free_dof, fe_index);
415 
416  // next the sides. do the
417  // same as above: check whether
418  // the line is already numbered
419  // for the present fe_index, and
420  // if not do it
421  if (fe.dofs_per_line > 0)
422  for (unsigned int l=0; l<GeometryInfo<2>::lines_per_cell; ++l)
423  {
424  typename HpDoFHandler<dim,spacedim>::line_iterator
425  line = cell->line(l);
426 
427  if (line->dof_index(0,fe_index) ==
428  DoFHandler<dim,spacedim>::invalid_dof_index)
429  for (unsigned int d=0; d<fe.dofs_per_line; ++d, ++next_free_dof)
430  line->set_dof_index (d, next_free_dof, fe_index);
431  }
432 
433 
434  // finally for the quad. this one
435  // shouldn't be numbered yet
436  if (fe.dofs_per_quad > 0)
437  {
438  Assert ((cell->dof_index(0, fe_index) ==
439  DoFHandler<dim,spacedim>::invalid_dof_index),
440  ExcInternalError());
441 
442  for (unsigned int d=0; d<fe.dofs_per_quad; ++d, ++next_free_dof)
443  cell->set_dof_index (d, next_free_dof, fe_index);
444  }
445 
446  // note that this cell has been processed
447  cell->set_user_flag ();
448 
449  return next_free_dof;
450  }
451 
452 
453  template <int spacedim>
454  static
455  types::global_dof_index
456  distribute_dofs_on_cell (const typename ::hp::DoFHandler<3,spacedim>::active_cell_iterator &cell,
457  types::global_dof_index next_free_dof)
458  {
459  const unsigned int dim = 3;
460 
461  const FiniteElement<dim,spacedim> &fe = cell->get_fe();
462  const unsigned int fe_index = cell->active_fe_index ();
463 
464  // number dofs on vertices. to do
465  // so, check whether dofs for
466  // this vertex have been
467  // distributed and for the
468  // present fe (only check the
469  // first dof), and if this isn't
470  // the case distribute new ones
471  // there
472  if (fe.dofs_per_vertex > 0)
473  for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
474  if (cell->vertex_dof_index(vertex, 0, fe_index) ==
475  DoFHandler<dim,spacedim>::invalid_dof_index)
476  for (unsigned int d=0; d<fe.dofs_per_vertex; ++d, ++next_free_dof)
477  cell->set_vertex_dof_index (vertex, d, next_free_dof, fe_index);
478 
479  // next the four lines. do the
480  // same as above: check whether
481  // the line is already numbered
482  // for the present fe_index, and
483  // if not do it
484  if (fe.dofs_per_line > 0)
485  for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
486  {
487  typename HpDoFHandler<dim,spacedim>::line_iterator
488  line = cell->line(l);
489 
490  if (line->dof_index(0,fe_index) ==
491  DoFHandler<dim,spacedim>::invalid_dof_index)
492  for (unsigned int d=0; d<fe.dofs_per_line; ++d, ++next_free_dof)
493  line->set_dof_index (d, next_free_dof, fe_index);
494  }
495 
496  // same for quads
497  if (fe.dofs_per_quad > 0)
498  for (unsigned int q=0; q<GeometryInfo<3>::quads_per_cell; ++q)
499  {
500  typename HpDoFHandler<dim,spacedim>::quad_iterator
501  quad = cell->quad(q);
502 
503  if (quad->dof_index(0,fe_index) ==
504  DoFHandler<dim,spacedim>::invalid_dof_index)
505  for (unsigned int d=0; d<fe.dofs_per_quad; ++d, ++next_free_dof)
506  quad->set_dof_index (d, next_free_dof, fe_index);
507  }
508 
509 
510  // finally for the hex. this one
511  // shouldn't be numbered yet
512  if (fe.dofs_per_hex > 0)
513  {
514  Assert ((cell->dof_index(0, fe_index) ==
515  DoFHandler<dim,spacedim>::invalid_dof_index),
516  ExcInternalError());
517 
518  for (unsigned int d=0; d<fe.dofs_per_hex; ++d, ++next_free_dof)
519  cell->set_dof_index (d, next_free_dof, fe_index);
520  }
521 
522  // note that this cell has been processed
523  cell->set_user_flag ();
524 
525  return next_free_dof;
526  }
527 
528 
538  template <int spacedim>
539  static
540  void
541  reserve_space (DoFHandler<1,spacedim> &dof_handler)
542  {
543  const unsigned int dim = 1;
544 
545  typedef DoFHandler<dim,spacedim> BaseClass;
546 
547  Assert (dof_handler.finite_elements != 0,
548  typename BaseClass::ExcNoFESelected());
549  Assert (dof_handler.finite_elements->size() > 0,
550  typename BaseClass::ExcNoFESelected());
551  Assert (dof_handler.tria->n_levels() > 0,
552  typename
553  BaseClass::ExcInvalidTriangulation());
554  Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
555  ExcInternalError ());
556 
557  // Release all space except the
558  // active_fe_indices field which
559  // we have to backup before
560  {
561  std::vector<std::vector<DoFLevel::active_fe_index_type> >
562  active_fe_backup(dof_handler.levels.size ());
563  for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
564  std::swap (dof_handler.levels[level]->active_fe_indices,
565  active_fe_backup[level]);
566 
567  // delete all levels and set them up
568  // newly, since vectors are
569  // troublesome if you want to change
570  // their size
571  dof_handler.clear_space ();
572 
573  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
574  {
575  dof_handler.levels.push_back (new internal::hp::DoFLevel);
576  std::swap (active_fe_backup[level],
577  dof_handler.levels[level]->active_fe_indices);
578  }
579  }
580 
581  // LINE (CELL) DOFs
582 
583  // count how much space we need
584  // on each level for the cell
585  // dofs and set the
586  // dof_*_offsets
587  // data. initially set the latter
588  // to an invalid index, and only
589  // later set it to something
590  // reasonable for active dof_handler.cells
591  //
592  // note that for dof_handler.cells, the
593  // situation is simpler than for
594  // other (lower dimensional)
595  // objects since exactly one
596  // finite element is used for it
597  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
598  {
599  dof_handler.levels[level]->dof_offsets
600  = std::vector<DoFLevel::offset_type> (
601  dof_handler.tria->n_raw_lines(level),
602  (DoFLevel::offset_type)(-1));
603  dof_handler.levels[level]->cell_cache_offsets
604  = std::vector<DoFLevel::offset_type> (
605  dof_handler.tria->n_raw_lines(level),
606  (DoFLevel::offset_type)(-1));
607 
608  types::global_dof_index next_free_dof = 0;
609  types::global_dof_index cache_size = 0;
610  typename HpDoFHandler<dim,spacedim>::active_cell_iterator
611  cell=dof_handler.begin_active(level),
612  endc=dof_handler.end_active(level);
613  for (; cell!=endc; ++cell)
614  if (!cell->has_children())
615  {
616  dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
617  next_free_dof += cell->get_fe().dofs_per_line;
618 
619  dof_handler.levels[level]->cell_cache_offsets[cell->index()] = cache_size;
620  cache_size += cell->get_fe().dofs_per_cell;
621  }
622 
623  dof_handler.levels[level]->dof_indices
624  = std::vector<types::global_dof_index> (next_free_dof,
625  DoFHandler<dim,spacedim>::invalid_dof_index);
626  dof_handler.levels[level]->cell_dof_indices_cache
627  = std::vector<types::global_dof_index> (cache_size,
628  DoFHandler<dim,spacedim>::invalid_dof_index);
629  }
630 
631  // safety check: make sure that
632  // the number of DoFs we
633  // allocated is actually correct
634  // (above we have also set the
635  // dof_*_offsets field, so
636  // we couldn't use this simpler
637  // algorithm)
638 #ifdef DEBUG
639  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
640  {
641  types::global_dof_index counter = 0;
642  typename HpDoFHandler<dim,spacedim>::active_cell_iterator
643  cell=dof_handler.begin_active(level),
644  endc=dof_handler.end_active(level);
645  for (; cell!=endc; ++cell)
646  if (!cell->has_children())
647  counter += cell->get_fe().dofs_per_line;
648 
649  Assert (dof_handler.levels[level]->dof_indices.size() == counter,
650  ExcInternalError());
651  Assert (static_cast<unsigned int>
652  (std::count (dof_handler.levels[level]->dof_offsets.begin(),
653  dof_handler.levels[level]->dof_offsets.end(),
654  (DoFLevel::offset_type)(-1)))
655  ==
656  dof_handler.tria->n_raw_lines(level) - dof_handler.tria->n_active_lines(level),
657  ExcInternalError());
658  }
659 #endif
660 
661 
662  // VERTEX DOFS
663  reserve_space_vertices (dof_handler);
664  }
665 
666 
667  template <int spacedim>
668  static
669  void
670  reserve_space (DoFHandler<2,spacedim> &dof_handler)
671  {
672  const unsigned int dim = 2;
673 
674  typedef DoFHandler<dim,spacedim> BaseClass;
675 
676  Assert (dof_handler.finite_elements != 0,
677  typename BaseClass::ExcNoFESelected());
678  Assert (dof_handler.finite_elements->size() > 0,
679  typename BaseClass::ExcNoFESelected());
680  Assert (dof_handler.tria->n_levels() > 0,
681  typename BaseClass::ExcInvalidTriangulation());
682  Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
683  ExcInternalError ());
684 
685  // Release all space except the
686  // active_fe_indices field which
687  // we have to backup before
688  {
689  std::vector<std::vector<DoFLevel::active_fe_index_type> >
690  active_fe_backup(dof_handler.levels.size ());
691  for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
692  std::swap (dof_handler.levels[level]->active_fe_indices,
693  active_fe_backup[level]);
694 
695  // delete all levels and set them up
696  // newly, since vectors are
697  // troublesome if you want to change
698  // their size
699  dof_handler.clear_space ();
700 
701  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
702  {
703  dof_handler.levels.push_back (new internal::hp::DoFLevel);
704  std::swap (active_fe_backup[level],
705  dof_handler.levels[level]->active_fe_indices);
706  }
707  dof_handler.faces = new internal::hp::DoFIndicesOnFaces<2>;
708  }
709 
710  // QUAD (CELL) DOFs
711 
712  // count how much space we need
713  // on each level for the cell
714  // dofs and set the
715  // dof_*_offsets
716  // data. initially set the latter
717  // to an invalid index, and only
718  // later set it to something
719  // reasonable for active dof_handler.cells
720  //
721  // note that for dof_handler.cells, the
722  // situation is simpler than for
723  // other (lower dimensional)
724  // objects since exactly one
725  // finite element is used for it
726  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
727  {
728  dof_handler.levels[level]->dof_offsets
729  = std::vector<DoFLevel::offset_type> (
730  dof_handler.tria->n_raw_quads(level),
731  (DoFLevel::offset_type)(-1));
732  dof_handler.levels[level]->cell_cache_offsets
733  = std::vector<DoFLevel::offset_type> (
734  dof_handler.tria->n_raw_quads(level),
735  (DoFLevel::offset_type)(-1));
736 
737  types::global_dof_index next_free_dof = 0;
738  types::global_dof_index cache_size = 0;
739  typename HpDoFHandler<dim,spacedim>::active_cell_iterator
740  cell=dof_handler.begin_active(level),
741  endc=dof_handler.end_active(level);
742  for (; cell!=endc; ++cell)
743  if (!cell->has_children())
744  {
745  dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
746  next_free_dof += cell->get_fe().dofs_per_quad;
747 
748  dof_handler.levels[level]->cell_cache_offsets[cell->index()] = cache_size;
749  cache_size += cell->get_fe().dofs_per_cell;
750  }
751 
752  dof_handler.levels[level]->dof_indices
753  = std::vector<types::global_dof_index> (next_free_dof,
754  DoFHandler<dim,spacedim>::invalid_dof_index);
755  dof_handler.levels[level]->cell_dof_indices_cache
756  = std::vector<types::global_dof_index> (cache_size,
757  DoFHandler<dim,spacedim>::invalid_dof_index);
758  }
759 
760  // safety check: make sure that
761  // the number of DoFs we
762  // allocated is actually correct
763  // (above we have also set the
764  // dof_*_offsets field, so
765  // we couldn't use this simpler
766  // algorithm)
767 #ifdef DEBUG
768  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
769  {
770  types::global_dof_index counter = 0;
771  typename HpDoFHandler<dim,spacedim>::active_cell_iterator
772  cell=dof_handler.begin_active(level),
773  endc=dof_handler.end_active(level);
774  for (; cell!=endc; ++cell)
775  if (!cell->has_children())
776  counter += cell->get_fe().dofs_per_quad;
777 
778  Assert (dof_handler.levels[level]->dof_indices.size() == counter,
779  ExcInternalError());
780  Assert (static_cast<unsigned int>
781  (std::count (dof_handler.levels[level]->dof_offsets.begin(),
782  dof_handler.levels[level]->dof_offsets.end(),
783  (DoFLevel::offset_type)(-1)))
784  ==
785  dof_handler.tria->n_raw_quads(level) - dof_handler.tria->n_active_quads(level),
786  ExcInternalError());
787  }
788 #endif
789 
790 
791  // LINE DOFS
792  //
793  // same here: count line dofs,
794  // then allocate as much space as
795  // we need and prime the linked
796  // list for lines (see the
797  // description in hp::DoFLevel)
798  // with the indices we will
799  // need. note that our task is
800  // more complicated since two
801  // adjacent dof_handler.cells may have
802  // different active_fe_indices,
803  // in which case we need to
804  // allocate *two* sets of line
805  // dofs for the same line
806  //
807  // the way we do things is that
808  // we loop over all active dof_handler.cells
809  // (these are the ones that have
810  // DoFs only anyway) and all
811  // their dof_handler.faces. We note in the
812  // user flags whether we have
813  // previously visited a face and
814  // if so skip it (consequently,
815  // we have to save and later
816  // restore the line flags)
817  {
818  std::vector<bool> saved_line_user_flags;
819  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
820  .save_user_flags_line (saved_line_user_flags);
821  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
822  .clear_user_flags_line ();
823 
824  // an array to hold how many
825  // slots (see the hp::DoFLevel
826  // class) we will have to store
827  // on each level
828  unsigned int n_line_slots = 0;
829 
830  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
831  cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
832  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
833  if (! cell->face(face)->user_flag_set())
834  {
835  // ok, face has not been
836  // visited. so we need to
837  // allocate space for it. let's
838  // see how much we need: we need
839  // one set if a) there is no
840  // neighbor behind this face, or
841  // b) the neighbor is either
842  // coarser or finer than we are,
843  // or c) the neighbor is neither
844  // coarser nor finer, but has
845  // happens to have the same
846  // active_fe_index:
847  if (cell->at_boundary(face)
848  ||
849  cell->face(face)->has_children()
850  ||
851  cell->neighbor_is_coarser(face)
852  ||
853  (!cell->at_boundary(face)
854  &&
855  (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
856  // ok, one set of
857  // dofs. that makes
858  // one index, 1 times
859  // dofs_per_line
860  // dofs, and one stop
861  // index
862  n_line_slots
863  += (*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_line + 2;
864 
865  // otherwise we do
866  // indeed need two
867  // sets, i.e. two
868  // indices, two sets of
869  // dofs, and one stop
870  // index:
871  else
872  n_line_slots
873  += ((*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_line
874  +
875  (*dof_handler.finite_elements)[cell->neighbor(face)->active_fe_index()]
876  .dofs_per_line
877  +
878  3);
879 
880  // mark this face as
881  // visited
882  cell->face(face)->set_user_flag ();
883  }
884 
885  // now that we know how many
886  // line dofs we will have to
887  // have on each level, allocate
888  // the memory. note that we
889  // allocate offsets for all
890  // lines, though only the
891  // active ones will have a
892  // non-invalid value later on
893  dof_handler.faces->lines.dof_offsets
894  = std::vector<unsigned int> (dof_handler.tria->n_raw_lines(),
895  (unsigned int)(-1));
896  dof_handler.faces->lines.dofs
897  = std::vector<types::global_dof_index> (n_line_slots,
898  DoFHandler<dim,spacedim>::invalid_dof_index);
899 
900  // with the memory now
901  // allocated, loop over the
902  // dof_handler.cells again and prime the
903  // _offset values as well as
904  // the fe_index fields
905  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
906  .clear_user_flags_line ();
907 
908  unsigned int next_free_line_slot = 0;
909 
910  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
911  cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
912  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
913  if (! cell->face(face)->user_flag_set())
914  {
915  // same decision tree
916  // as before
917  if (cell->at_boundary(face)
918  ||
919  cell->face(face)->has_children()
920  ||
921  cell->neighbor_is_coarser(face)
922  ||
923  (!cell->at_boundary(face)
924  &&
925  (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
926  {
927  dof_handler.faces
928  ->lines.dof_offsets[cell->face(face)->index()]
929  = next_free_line_slot;
930 
931  // set first slot
932  // for this line to
933  // active_fe_index
934  // of this face
935  dof_handler.faces
936  ->lines.dofs[next_free_line_slot]
937  = cell->active_fe_index();
938 
939  // the next
940  // dofs_per_line
941  // indices remain
942  // unset for the
943  // moment (i.e. at
944  // invalid_dof_index).
945  // following this
946  // comes the stop
947  // index, which
948  // also is
949  // invalid_dof_index
950  // and therefore
951  // does not have to
952  // be explicitly
953  // set
954 
955  // finally, mark
956  // those slots as
957  // used
958  next_free_line_slot
959  += (*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_line + 2;
960  }
961  else
962  {
963  dof_handler.faces
964  ->lines.dof_offsets[cell->face(face)->index()]
965  = next_free_line_slot;
966 
967  // set first slot
968  // for this line to
969  // active_fe_index
970  // of this face
971  dof_handler.faces
972  ->lines.dofs[next_free_line_slot]
973  = cell->active_fe_index();
974 
975  // the next
976  // dofs_per_line
977  // indices remain
978  // unset for the
979  // moment (i.e. at
980  // invalid_dof_index).
981  //
982  // then comes the
983  // fe_index for the
984  // neighboring
985  // cell:
986  dof_handler.faces
987  ->lines.dofs[next_free_line_slot
988  +
989  (*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_line
990  +
991  1]
992  = cell->neighbor(face)->active_fe_index();
993  // then again a set
994  // of dofs that we
995  // need not set
996  // right now
997  //
998  // following this
999  // comes the stop
1000  // index, which
1001  // also is
1002  // invalid_dof_index
1003  // and therefore
1004  // does not have to
1005  // be explicitly
1006  // set
1007 
1008  // finally, mark
1009  // those slots as
1010  // used
1011  next_free_line_slot
1012  += ((*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_line
1013  +
1014  (*dof_handler.finite_elements)[cell->neighbor(face)->active_fe_index()]
1015  .dofs_per_line
1016  +
1017  3);
1018  }
1019 
1020  // mark this face as
1021  // visited
1022  cell->face(face)->set_user_flag ();
1023  }
1024 
1025  // we should have moved the
1026  // cursor for each level to the
1027  // total number of dofs on that
1028  // level. check that
1029  Assert (next_free_line_slot == n_line_slots,
1030  ExcInternalError());
1031 
1032  // at the end, restore the user
1033  // flags for the lines
1034  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
1035  .load_user_flags_line (saved_line_user_flags);
1036  }
1037 
1038 
1039  // VERTEX DOFS
1040  reserve_space_vertices (dof_handler);
1041  }
1042 
1043 
1044  template <int spacedim>
1045  static
1046  void
1047  reserve_space (DoFHandler<3,spacedim> &dof_handler)
1048  {
1049  const unsigned int dim = 3;
1050 
1051  typedef DoFHandler<dim,spacedim> BaseClass;
1052 
1053  Assert (dof_handler.finite_elements != 0,
1054  typename BaseClass::ExcNoFESelected());
1055  Assert (dof_handler.finite_elements->size() > 0,
1056  typename BaseClass::ExcNoFESelected());
1057  Assert (dof_handler.tria->n_levels() > 0,
1058  typename BaseClass::ExcInvalidTriangulation());
1059  Assert (dof_handler.tria->n_levels() == dof_handler.levels.size (),
1060  ExcInternalError ());
1061 
1062  // Release all space except the
1063  // active_fe_indices field which
1064  // we have to backup before
1065  {
1066  std::vector<std::vector<DoFLevel::active_fe_index_type> >
1067  active_fe_backup(dof_handler.levels.size ());
1068  for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
1069  std::swap (dof_handler.levels[level]->active_fe_indices,
1070  active_fe_backup[level]);
1071 
1072  // delete all levels and set them up
1073  // newly, since vectors are
1074  // troublesome if you want to change
1075  // their size
1076  dof_handler.clear_space ();
1077 
1078  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
1079  {
1080  dof_handler.levels.push_back (new internal::hp::DoFLevel);
1081  std::swap (active_fe_backup[level],
1082  dof_handler.levels[level]->active_fe_indices);
1083  }
1084  dof_handler.faces = new internal::hp::DoFIndicesOnFaces<3>;
1085  }
1086 
1087  // HEX (CELL) DOFs
1088 
1089  // count how much space we need
1090  // on each level for the cell
1091  // dofs and set the
1092  // dof_*_offsets
1093  // data. initially set the latter
1094  // to an invalid index, and only
1095  // later set it to something
1096  // reasonable for active dof_handler.cells
1097  //
1098  // note that for dof_handler.cells, the
1099  // situation is simpler than for
1100  // other (lower dimensional)
1101  // objects since exactly one
1102  // finite element is used for it
1103  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
1104  {
1105  dof_handler.levels[level]->dof_offsets
1106  = std::vector<DoFLevel::offset_type> (
1107  dof_handler.tria->n_raw_hexs(level),
1108  (DoFLevel::offset_type)(-1));
1109  dof_handler.levels[level]->cell_cache_offsets
1110  = std::vector<DoFLevel::offset_type> (
1111  dof_handler.tria->n_raw_hexs(level),
1112  (DoFLevel::offset_type)(-1));
1113 
1114  types::global_dof_index next_free_dof = 0;
1115  types::global_dof_index cache_size = 0;
1116  typename HpDoFHandler<dim,spacedim>::active_cell_iterator
1117  cell=dof_handler.begin_active(level),
1118  endc=dof_handler.end_active(level);
1119  for (; cell!=endc; ++cell)
1120  if (!cell->has_children())
1121  {
1122  dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
1123  next_free_dof += cell->get_fe().dofs_per_hex;
1124 
1125  dof_handler.levels[level]->cell_cache_offsets[cell->index()] = cache_size;
1126  cache_size += cell->get_fe().dofs_per_cell;
1127  }
1128 
1129  dof_handler.levels[level]->dof_indices
1130  = std::vector<types::global_dof_index> (next_free_dof,
1131  DoFHandler<dim,spacedim>::invalid_dof_index);
1132  dof_handler.levels[level]->cell_dof_indices_cache
1133  = std::vector<types::global_dof_index> (cache_size,
1134  DoFHandler<dim,spacedim>::invalid_dof_index);
1135  }
1136 
1137  // safety check: make sure that
1138  // the number of DoFs we
1139  // allocated is actually correct
1140  // (above we have also set the
1141  // dof_*_offsets field, so
1142  // we couldn't use this simpler
1143  // algorithm)
1144 #ifdef DEBUG
1145  for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
1146  {
1147  types::global_dof_index counter = 0;
1148  typename HpDoFHandler<dim,spacedim>::active_cell_iterator
1149  cell=dof_handler.begin_active(level),
1150  endc=dof_handler.end_active(level);
1151  for (; cell!=endc; ++cell)
1152  if (!cell->has_children())
1153  counter += cell->get_fe().dofs_per_hex;
1154 
1155  Assert (dof_handler.levels[level]->dof_indices.size() == counter,
1156  ExcInternalError());
1157  Assert (static_cast<unsigned int>
1158  (std::count (dof_handler.levels[level]->dof_offsets.begin(),
1159  dof_handler.levels[level]->dof_offsets.end(),
1160  (DoFLevel::offset_type)(-1)))
1161  ==
1162  dof_handler.tria->n_raw_hexs(level) - dof_handler.tria->n_active_hexs(level),
1163  ExcInternalError());
1164  }
1165 #endif
1166 
1167 
1168  // QUAD DOFS
1169  //
1170  // same here: count quad dofs,
1171  // then allocate as much space as
1172  // we need and prime the linked
1173  // list for quad (see the
1174  // description in hp::DoFLevel)
1175  // with the indices we will
1176  // need. note that our task is
1177  // more complicated since two
1178  // adjacent dof_handler.cells may have
1179  // different active_fe_indices,
1180  // in which case we need to
1181  // allocate *two* sets of line
1182  // dofs for the same line
1183  //
1184  // the way we do things is that
1185  // we loop over all active dof_handler.cells
1186  // (these are the ones that have
1187  // DoFs only anyway) and all
1188  // their dof_handler.faces. We note in the
1189  // user flags whether we have
1190  // previously visited a face and
1191  // if so skip it (consequently,
1192  // we have to save and later
1193  // restore the line flags)
1194  {
1195  std::vector<bool> saved_quad_user_flags;
1196  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
1197  .save_user_flags_quad (saved_quad_user_flags);
1198  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
1199  .clear_user_flags_quad ();
1200 
1201  // examine, how how many
1202  // slots (see the hp::DoFLevel
1203  // class) we will have to store
1204  unsigned int n_quad_slots = 0;
1205 
1206  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
1207  cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
1208  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1209  if (! cell->face(face)->user_flag_set())
1210  {
1211  // ok, face has not been
1212  // visited. so we need to
1213  // allocate space for
1214  // it. let's see how much
1215  // we need: we need one
1216  // set if a) there is no
1217  // neighbor behind this
1218  // face, or b) the
1219  // neighbor is not on the
1220  // same level or further
1221  // refined, or c) the
1222  // neighbor is on the
1223  // same level, but
1224  // happens to have the
1225  // same active_fe_index:
1226  if (cell->at_boundary(face)
1227  ||
1228  cell->face(face)->has_children()
1229  ||
1230  cell->neighbor_is_coarser(face)
1231  ||
1232  (!cell->at_boundary(face)
1233  &&
1234  (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
1235  // ok, one set of
1236  // dofs. that makes
1237  // one index, 1 times
1238  // dofs_per_quad
1239  // dofs, and one stop
1240  // index
1241  n_quad_slots
1242  += (*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_quad + 2;
1243 
1244  // otherwise we do
1245  // indeed need two
1246  // sets, i.e. two
1247  // indices, two sets of
1248  // dofs, and one stop
1249  // index:
1250  else
1251  n_quad_slots
1252  += ((*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_quad
1253  +
1254  (*dof_handler.finite_elements)[cell->neighbor(face)->active_fe_index()]
1255  .dofs_per_quad
1256  +
1257  3);
1258 
1259  // mark this face as
1260  // visited
1261  cell->face(face)->set_user_flag ();
1262  }
1263 
1264  // now that we know how many
1265  // quad dofs we will have to
1266  // have, allocate
1267  // the memory. note that we
1268  // allocate offsets for all
1269  // quads, though only the
1270  // active ones will have a
1271  // non-invalid value later on
1272  if (true)
1273  {
1274  dof_handler.faces->quads.dof_offsets
1275  = std::vector<unsigned int>
1276  (dof_handler.tria->n_raw_quads(),
1277  (unsigned int)(-1));
1278  dof_handler.faces->quads.dofs
1279  = std::vector<types::global_dof_index> (n_quad_slots,
1280  DoFHandler<dim,spacedim>::invalid_dof_index);
1281  }
1282 
1283  // with the memory now
1284  // allocated, loop over the
1285  // dof_handler.cells again and prime the
1286  // _offset values as well as
1287  // the fe_index fields
1288  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
1289  .clear_user_flags_quad ();
1290 
1291  unsigned int next_free_quad_slot = 0;
1292 
1293  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
1294  cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
1295  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1296  if (! cell->face(face)->user_flag_set())
1297  {
1298  // same decision tree
1299  // as before
1300  if (cell->at_boundary(face)
1301  ||
1302  cell->face(face)->has_children()
1303  ||
1304  cell->neighbor_is_coarser(face)
1305  ||
1306  (!cell->at_boundary(face)
1307  &&
1308  (cell->active_fe_index() == cell->neighbor(face)->active_fe_index())))
1309  {
1310  dof_handler.faces
1311  ->quads.dof_offsets[cell->face(face)->index()]
1312  = next_free_quad_slot;
1313 
1314  // set first slot
1315  // for this quad to
1316  // active_fe_index
1317  // of this face
1318  dof_handler.faces
1319  ->quads.dofs[next_free_quad_slot]
1320  = cell->active_fe_index();
1321 
1322  // the next
1323  // dofs_per_quad
1324  // indices remain
1325  // unset for the
1326  // moment (i.e. at
1327  // invalid_dof_index).
1328  // following this
1329  // comes the stop
1330  // index, which
1331  // also is
1332  // invalid_dof_index
1333  // and therefore
1334  // does not have to
1335  // be explicitly
1336  // set
1337 
1338  // finally, mark
1339  // those slots as
1340  // used
1341  next_free_quad_slot
1342  += (*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_quad + 2;
1343  }
1344  else
1345  {
1346  dof_handler.faces
1347  ->quads.dof_offsets[cell->face(face)->index()]
1348  = next_free_quad_slot;
1349 
1350  // set first slot
1351  // for this quad to
1352  // active_fe_index
1353  // of this face
1354  dof_handler.faces
1355  ->quads.dofs[next_free_quad_slot]
1356  = cell->active_fe_index();
1357 
1358  // the next
1359  // dofs_per_quad
1360  // indices remain
1361  // unset for the
1362  // moment (i.e. at
1363  // invalid_dof_index).
1364  //
1365  // then comes the
1366  // fe_index for the
1367  // neighboring
1368  // cell:
1369  dof_handler.faces
1370  ->quads.dofs[next_free_quad_slot
1371  +
1372  (*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_quad
1373  +
1374  1]
1375  = cell->neighbor(face)->active_fe_index();
1376  // then again a set
1377  // of dofs that we
1378  // need not set
1379  // right now
1380  //
1381  // following this
1382  // comes the stop
1383  // index, which
1384  // also is
1385  // invalid_dof_index
1386  // and therefore
1387  // does not have to
1388  // be explicitly
1389  // set
1390 
1391  // finally, mark
1392  // those slots as
1393  // used
1394  next_free_quad_slot
1395  += ((*dof_handler.finite_elements)[cell->active_fe_index()].dofs_per_quad
1396  +
1397  (*dof_handler.finite_elements)[cell->neighbor(face)->active_fe_index()]
1398  .dofs_per_quad
1399  +
1400  3);
1401  }
1402 
1403  // mark this face as
1404  // visited
1405  cell->face(face)->set_user_flag ();
1406  }
1407 
1408  // we should have moved the
1409  // cursor to the total number
1410  // of dofs. check that
1411  Assert (next_free_quad_slot == n_quad_slots,
1412  ExcInternalError());
1413 
1414  // at the end, restore the user
1415  // flags for the quads
1416  const_cast<::Triangulation<dim,spacedim>&>(*dof_handler.tria)
1417  .load_user_flags_quad (saved_quad_user_flags);
1418  }
1419 
1420 
1421  // LINE DOFS
1422 
1423  // the situation here is pretty
1424  // much like with vertices: there
1425  // can be an arbitrary number of
1426  // finite elements associated
1427  // with each line.
1428  //
1429  // the algorithm we use is
1430  // somewhat similar to what we do
1431  // in reserve_space_vertices()
1432  if (true)
1433  {
1434  // what we do first is to set up
1435  // an array in which we record
1436  // whether a line is associated
1437  // with any of the given fe's, by
1438  // setting a bit. in a later
1439  // step, we then actually
1440  // allocate memory for the
1441  // required dofs
1442  std::vector<std::vector<bool> >
1443  line_fe_association (dof_handler.finite_elements->size(),
1444  std::vector<bool> (dof_handler.tria->n_raw_lines(),
1445  false));
1446 
1447  for (typename HpDoFHandler<dim,spacedim>::active_cell_iterator
1448  cell=dof_handler.begin_active();
1449  cell!=dof_handler.end(); ++cell)
1450  for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
1451  line_fe_association[cell->active_fe_index()][cell->line_index(l)]
1452  = true;
1453 
1454  // first check which of the
1455  // lines is used at all,
1456  // i.e. is associated with a
1457  // finite element. we do this
1458  // since not all lines may
1459  // actually be used, in which
1460  // case we do not have to
1461  // allocate any memory at
1462  // all
1463  std::vector<bool> line_is_used (dof_handler.tria->n_raw_lines(), false);
1464  for (unsigned int line=0; line<dof_handler.tria->n_raw_lines(); ++line)
1465  for (unsigned int fe=0; fe<dof_handler.finite_elements->size(); ++fe)
1466  if (line_fe_association[fe][line] == true)
1467  {
1468  line_is_used[line] = true;
1469  break;
1470  }
1471 
1472  // next count how much memory
1473  // we actually need. for each
1474  // line, we need one slot per
1475  // fe to store the fe_index,
1476  // plus dofs_per_line for
1477  // this fe. in addition, we
1478  // need one slot as the end
1479  // marker for the
1480  // fe_indices. at the same
1481  // time already fill the
1482  // line_dofs_offsets field
1483  dof_handler.faces->lines.dof_offsets
1484  .resize (dof_handler.tria->n_raw_lines(),
1485  numbers::invalid_unsigned_int);
1486 
1487  unsigned int line_slots_needed = 0;
1488  for (unsigned int line=0; line<dof_handler.tria->n_raw_lines(); ++line)
1489  if (line_is_used[line] == true)
1490  {
1491  dof_handler.faces->lines.dof_offsets[line] = line_slots_needed;
1492 
1493  for (unsigned int fe=0; fe<dof_handler.finite_elements->size(); ++fe)
1494  if (line_fe_association[fe][line] == true)
1495  line_slots_needed += (*dof_handler.finite_elements)[fe].dofs_per_line + 1;
1496  ++line_slots_needed;
1497  }
1498 
1499  // now allocate the space we
1500  // have determined we need, and
1501  // set up the linked lists for
1502  // each of the lines
1503  dof_handler.faces->lines.dofs.resize (line_slots_needed,
1504  DoFHandler<dim,spacedim>::invalid_dof_index);
1505  for (unsigned int line=0; line<dof_handler.tria->n_raw_lines(); ++line)
1506  if (line_is_used[line] == true)
1507  {
1508  unsigned int pointer = dof_handler.faces->lines.dof_offsets[line];
1509  for (unsigned int fe=0; fe<dof_handler.finite_elements->size(); ++fe)
1510  if (line_fe_association[fe][line] == true)
1511  {
1512  // if this line
1513  // uses this fe,
1514  // then set the
1515  // fe_index and
1516  // move the
1517  // pointer ahead
1518  dof_handler.faces->lines.dofs[pointer] = fe;
1519  pointer += (*dof_handler.finite_elements)[fe].dofs_per_line + 1;
1520  }
1521  // finally place the end
1522  // marker
1523  dof_handler.faces->lines.dofs[pointer] = numbers::invalid_dof_index;
1524  }
1525  }
1526 
1527 
1528 
1529  // VERTEX DOFS
1530  reserve_space_vertices (dof_handler);
1531  }
1532 
1533 
1538  template <int spacedim>
1539  static
1540  unsigned int
1541  max_couplings_between_dofs (const DoFHandler<1,spacedim> &dof_handler)
1542  {
1543  return std::min(static_cast<types::global_dof_index> (3*
1544  dof_handler.finite_elements->max_dofs_per_vertex() +
1545  2*dof_handler.finite_elements->max_dofs_per_line()),
1546  dof_handler.n_dofs());
1547  }
1548 
1549 
1550 
1551  template <int spacedim>
1552  static
1553  unsigned int
1554  max_couplings_between_dofs (const DoFHandler<2,spacedim> &dof_handler)
1555  {
1556  // get these numbers by drawing pictures
1557  // and counting...
1558  // example:
1559  // | | |
1560  // --x-----x--x--X--
1561  // | | | |
1562  // | x--x--x
1563  // | | | |
1564  // --x--x--*--x--x--
1565  // | | | |
1566  // x--x--x |
1567  // | | | |
1568  // --X--x--x-----x--
1569  // | | |
1570  // x = vertices connected with center vertex *;
1571  // = total of 19
1572  // (the X vertices are connected with * if
1573  // the vertices adjacent to X are hanging
1574  // nodes)
1575  // count lines -> 28 (don't forget to count
1576  // mother and children separately!)
1577  types::global_dof_index max_couplings;
1578  switch (dof_handler.tria->max_adjacent_cells())
1579  {
1580  case 4:
1581  max_couplings=19*dof_handler.finite_elements->max_dofs_per_vertex() +
1582  28*dof_handler.finite_elements->max_dofs_per_line() +
1583  8*dof_handler.finite_elements->max_dofs_per_quad();
1584  break;
1585  case 5:
1586  max_couplings=21*dof_handler.finite_elements->max_dofs_per_vertex() +
1587  31*dof_handler.finite_elements->max_dofs_per_line() +
1588  9*dof_handler.finite_elements->max_dofs_per_quad();
1589  break;
1590  case 6:
1591  max_couplings=28*dof_handler.finite_elements->max_dofs_per_vertex() +
1592  42*dof_handler.finite_elements->max_dofs_per_line() +
1593  12*dof_handler.finite_elements->max_dofs_per_quad();
1594  break;
1595  case 7:
1596  max_couplings=30*dof_handler.finite_elements->max_dofs_per_vertex() +
1597  45*dof_handler.finite_elements->max_dofs_per_line() +
1598  13*dof_handler.finite_elements->max_dofs_per_quad();
1599  break;
1600  case 8:
1601  max_couplings=37*dof_handler.finite_elements->max_dofs_per_vertex() +
1602  56*dof_handler.finite_elements->max_dofs_per_line() +
1603  16*dof_handler.finite_elements->max_dofs_per_quad();
1604  break;
1605  default:
1606  Assert (false, ExcNotImplemented());
1607  max_couplings=0;
1608  };
1609  return std::min(max_couplings,dof_handler.n_dofs());
1610  }
1611 
1612 
1613  template <int spacedim>
1614  static
1615  unsigned int
1616  max_couplings_between_dofs (const DoFHandler<3,spacedim> &dof_handler)
1617  {
1618 //TODO:[?] Invent significantly better estimates than the ones in this function
1619  // doing the same thing here is a rather
1620  // complicated thing, compared to the 2d
1621  // case, since it is hard to draw pictures
1622  // with several refined hexahedra :-) so I
1623  // presently only give a coarse estimate
1624  // for the case that at most 8 hexes meet
1625  // at each vertex
1626  //
1627  // can anyone give better estimate here?
1628  const unsigned int max_adjacent_cells = dof_handler.tria->max_adjacent_cells();
1629 
1630  types::global_dof_index max_couplings;
1631  if (max_adjacent_cells <= 8)
1632  max_couplings=7*7*7*dof_handler.finite_elements->max_dofs_per_vertex() +
1633  7*6*7*3*dof_handler.finite_elements->max_dofs_per_line() +
1634  9*4*7*3*dof_handler.finite_elements->max_dofs_per_quad() +
1635  27*dof_handler.finite_elements->max_dofs_per_hex();
1636  else
1637  {
1638  Assert (false, ExcNotImplemented());
1639  max_couplings=0;
1640  }
1641 
1642  return std::min(max_couplings,dof_handler.n_dofs());
1643  }
1644  };
1645  }
1646  }
1647 }
1648 
1649 
1650 namespace hp
1651 {
1652  template<int dim, int spacedim>
1653  const unsigned int DoFHandler<dim,spacedim>::dimension;
1654 
1655  template<int dim, int spacedim>
1656  const types::global_dof_index DoFHandler<dim,spacedim>::invalid_dof_index;
1657 
1658  template<int dim, int spacedim>
1659  const unsigned int DoFHandler<dim,spacedim>::default_fe_index;
1660 
1661 
1662 
1663  template<int dim, int spacedim>
1664  DoFHandler<dim,spacedim>::DoFHandler (const Triangulation<dim,spacedim> &tria)
1665  :
1666  tria(&tria, typeid(*this).name()),
1667  faces (NULL)
1668  {
1669  Assert ((dynamic_cast<const parallel::distributed::Triangulation< dim, spacedim >*>
1670  (&tria)
1671  == 0),
1672  ExcMessage ("The given triangulation is parallel distributed but "
1673  "this class does not currently support this."));
1674 
1675  create_active_fe_table ();
1676 
1677  tria_listeners.push_back
1678  (tria.signals.pre_refinement
1679  .connect (std_cxx11::bind (&DoFHandler<dim,spacedim>::pre_refinement_action,
1680  std_cxx11::ref(*this))));
1681  tria_listeners.push_back
1682  (tria.signals.post_refinement
1683  .connect (std_cxx11::bind (&DoFHandler<dim,spacedim>::post_refinement_action,
1684  std_cxx11::ref(*this))));
1685  tria_listeners.push_back
1686  (tria.signals.create
1687  .connect (std_cxx11::bind (&DoFHandler<dim,spacedim>::post_refinement_action,
1688  std_cxx11::ref(*this))));
1689  }
1690 
1691 
1692  template<int dim, int spacedim>
1693  DoFHandler<dim,spacedim>::~DoFHandler ()
1694  {
1695  // unsubscribe as a listener to refinement
1696  // of the underlying triangulation
1697  for (unsigned int i=0; i<tria_listeners.size(); ++i)
1698  tria_listeners[i].disconnect ();
1699  tria_listeners.clear ();
1700 
1701  // ...and release allocated memory
1702  clear ();
1703  }
1704 
1705 
1706  /*------------------------ Cell iterator functions ------------------------*/
1707 
1708  template <int dim, int spacedim>
1709  typename DoFHandler<dim,spacedim>::cell_iterator
1710  DoFHandler<dim, spacedim>::begin(const unsigned int level) const
1711  {
1712  return cell_iterator (*this->get_triangulation().begin(level),
1713  this);
1714  }
1715 
1716 
1717 
1718  template <int dim, int spacedim>
1719  typename DoFHandler<dim,spacedim>::active_cell_iterator
1720  DoFHandler<dim,spacedim>::begin_active (const unsigned int level) const
1721  {
1722  // level is checked in begin
1723  cell_iterator i = begin (level);
1724  if (i.state() != IteratorState::valid)
1725  return i;
1726  while (i->has_children())
1727  if ((++i).state() != IteratorState::valid)
1728  return i;
1729  return i;
1730  }
1731 
1732 
1733 
1734  template <int dim, int spacedim>
1735  typename DoFHandler<dim,spacedim>::cell_iterator
1736  DoFHandler<dim,spacedim>::end () const
1737  {
1738  return cell_iterator (&this->get_triangulation(),
1739  -1,
1740  -1,
1741  this);
1742  }
1743 
1744 
1745  template <int dim, int spacedim>
1746  typename DoFHandler<dim,spacedim>::cell_iterator
1747  DoFHandler<dim,spacedim>::end (const unsigned int level) const
1748  {
1749  return (level == this->get_triangulation().n_levels()-1 ?
1750  end() :
1751  begin (level+1));
1752  }
1753 
1754 
1755  template <int dim, int spacedim>
1756  typename DoFHandler<dim, spacedim>::active_cell_iterator
1757  DoFHandler<dim, spacedim>::end_active (const unsigned int level) const
1758  {
1759  return (level == this->get_triangulation().n_levels()-1 ?
1760  active_cell_iterator(end()) :
1761  begin_active (level+1));
1762  }
1763 
1764 
1765 
1766  template <int dim, int spacedim>
1767  IteratorRange<typename DoFHandler<dim, spacedim>::cell_iterator>
1768  DoFHandler<dim, spacedim>::cell_iterators () const
1769  {
1770  return
1771  IteratorRange<typename DoFHandler<dim, spacedim>::cell_iterator>
1772  (begin(), end());
1773  }
1774 
1775 
1776  template <int dim, int spacedim>
1777  IteratorRange<typename DoFHandler<dim, spacedim>::active_cell_iterator>
1778  DoFHandler<dim, spacedim>::active_cell_iterators () const
1779  {
1780  return
1781  IteratorRange<typename DoFHandler<dim, spacedim>::active_cell_iterator>
1782  (begin_active(), end());
1783  }
1784 
1785 
1786 
1787  template <int dim, int spacedim>
1788  IteratorRange<typename DoFHandler<dim, spacedim>::cell_iterator>
1789  DoFHandler<dim, spacedim>::cell_iterators_on_level (const unsigned int level) const
1790  {
1791  return
1792  IteratorRange<typename DoFHandler<dim, spacedim>::cell_iterator>
1793  (begin(level), end(level));
1794  }
1795 
1796 
1797 
1798  template <int dim, int spacedim>
1799  IteratorRange<typename DoFHandler<dim, spacedim>::active_cell_iterator>
1800  DoFHandler<dim, spacedim>::active_cell_iterators_on_level (const unsigned int level) const
1801  {
1802  return
1803  IteratorRange<typename DoFHandler<dim, spacedim>::active_cell_iterator>
1804  (begin_active(level), end_active(level));
1805  }
1806 
1807 
1808 
1809 
1810 //------------------------------------------------------------------
1811 
1812 
1813  template <>
1814  types::global_dof_index DoFHandler<1>::n_boundary_dofs () const
1815  {
1816  Assert (finite_elements != 0, ExcNoFESelected());
1817 
1818  DoFHandler<1,1>::cell_iterator cell;
1819  types::global_dof_index n = 0;
1820 
1821  // search left-most cell
1822  cell = this->begin_active();
1823  while (!cell->at_boundary(0))
1824  cell = cell->neighbor(0);
1825  n += cell->get_fe().dofs_per_vertex;
1826 
1827  // same with right-most cell
1828  cell = this->begin_active();
1829  while (!cell->at_boundary(1))
1830  cell = cell->neighbor(1);
1831  n += cell->get_fe().dofs_per_vertex;
1832 
1833  return n;
1834  }
1835 
1836 
1837 
1838  template <>
1839  template <typename number>
1840  types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::map<types::boundary_id, const Function<1,number>*> &boundary_ids) const
1841  {
1842  Assert (finite_elements != 0, ExcNoFESelected());
1843 
1844  // check that only boundary
1845  // indicators 0 and 1 are allowed
1846  // in 1d
1847  for (typename std::map<types::boundary_id, const Function<1,number>*>::const_iterator i=boundary_ids.begin();
1848  i!=boundary_ids.end(); ++i)
1849  Assert ((i->first == 0) || (i->first == 1),
1850  ExcInvalidBoundaryIndicator());
1851 
1852  DoFHandler<1,1>::active_cell_iterator cell;
1853  types::global_dof_index n = 0;
1854 
1855  // search left-most cell
1856  if (boundary_ids.find (0) != boundary_ids.end())
1857  {
1858  cell = this->begin_active();
1859  while (!cell->at_boundary(0))
1860  cell = cell->neighbor(0);
1861  n += cell->get_fe().dofs_per_vertex;
1862  }
1863 
1864  // same with right-most cell
1865  if (boundary_ids.find (1) != boundary_ids.end())
1866  {
1867  cell = this->begin_active();
1868  while (!cell->at_boundary(1))
1869  cell = cell->neighbor(1);
1870  n += cell->get_fe().dofs_per_vertex;
1871  }
1872 
1873  return n;
1874  }
1875 
1876 
1877 
1878  template <>
1879  types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_ids) const
1880  {
1881  Assert (finite_elements != 0, ExcNoFESelected());
1882 
1883  // check that only boundary
1884  // indicators 0 and 1 are allowed
1885  // in 1d
1886  for (std::set<types::boundary_id>::const_iterator i=boundary_ids.begin();
1887  i!=boundary_ids.end(); ++i)
1888  Assert ((*i == 0) || (*i == 1),
1889  ExcInvalidBoundaryIndicator());
1890 
1891  DoFHandler<1,1>::active_cell_iterator cell;
1892  types::global_dof_index n = 0;
1893 
1894  // search left-most cell
1895  if (boundary_ids.find (0) != boundary_ids.end())
1896  {
1897  cell = this->begin_active();
1898  while (!cell->at_boundary(0))
1899  cell = cell->neighbor(0);
1900  n += cell->get_fe().dofs_per_vertex;
1901  }
1902 
1903  // same with right-most cell
1904  if (boundary_ids.find (1) != boundary_ids.end())
1905  {
1906  cell = this->begin_active();
1907  while (!cell->at_boundary(1))
1908  cell = cell->neighbor(1);
1909  n += cell->get_fe().dofs_per_vertex;
1910  }
1911 
1912  return n;
1913  }
1914 
1915 
1916  template <>
1917  types::global_dof_index DoFHandler<1,2>::n_boundary_dofs () const
1918  {
1919  Assert(false,ExcNotImplemented());
1920  return 0;
1921  }
1922 
1923  template <>
1924  template <typename number>
1925  types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::map<types::boundary_id, const Function<2,number>*> &) const
1926  {
1927  Assert(false,ExcNotImplemented());
1928  return 0;
1929  }
1930 
1931  template <>
1932  types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &) const
1933  {
1934  Assert(false,ExcNotImplemented());
1935  return 0;
1936  }
1937 
1938 
1939 
1940  template <>
1941  types::global_dof_index DoFHandler<1,3>::n_boundary_dofs () const
1942  {
1943  Assert(false,ExcNotImplemented());
1944  return 0;
1945  }
1946 
1947  template <>
1948  template <typename number>
1949  types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const std::map<types::boundary_id, const Function<3,number>*> &) const
1950  {
1951  Assert(false,ExcNotImplemented());
1952  return 0;
1953  }
1954 
1955  template <>
1956  types::global_dof_index DoFHandler<1,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
1957  {
1958  Assert(false,ExcNotImplemented());
1959  return 0;
1960  }
1961 
1962 
1963  template<int dim, int spacedim>
1964  types::global_dof_index DoFHandler<dim,spacedim>::n_boundary_dofs () const
1965  {
1966  Assert (finite_elements != 0, ExcNoFESelected());
1967 
1968  std::set<types::global_dof_index> boundary_dofs;
1969  std::vector<types::global_dof_index> dofs_on_face;
1970  dofs_on_face.reserve (this->get_fe ().max_dofs_per_face());
1971 
1972  // loop over all faces to check
1973  // whether they are at a
1974  // boundary. note that we need not
1975  // take special care of single
1976  // lines in 3d (using
1977  // @p{cell->has_boundary_lines}),
1978  // since we do not support
1979  // boundaries of dimension dim-2,
1980  // and so every boundary line is
1981  // also part of a boundary face.
1982  typename HpDoFHandler<dim,spacedim>::active_cell_iterator cell = this->begin_active (),
1983  endc = this->end();
1984  for (; cell!=endc; ++cell)
1985  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
1986  if (cell->at_boundary(f))
1987  {
1988  const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
1989  dofs_on_face.resize (dofs_per_face);
1990 
1991  cell->face(f)->get_dof_indices (dofs_on_face,
1992  cell->active_fe_index());
1993  for (unsigned int i=0; i<dofs_per_face; ++i)
1994  boundary_dofs.insert(dofs_on_face[i]);
1995  };
1996  return boundary_dofs.size();
1997  }
1998 
1999 
2000 
2001  template<int dim, int spacedim>
2002  types::global_dof_index
2003  DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_ids) const
2004  {
2005  Assert (finite_elements != 0, ExcNoFESelected());
2006  Assert (boundary_ids.find (numbers::internal_face_boundary_id) == boundary_ids.end(),
2007  ExcInvalidBoundaryIndicator());
2008 
2009  // same as above, but with
2010  // additional checks for set of
2011  // boundary indicators
2012  std::set<types::global_dof_index> boundary_dofs;
2013  std::vector<types::global_dof_index> dofs_on_face;
2014  dofs_on_face.reserve (this->get_fe ().max_dofs_per_face());
2015 
2016  typename HpDoFHandler<dim,spacedim>::active_cell_iterator cell = this->begin_active (),
2017  endc = this->end();
2018  for (; cell!=endc; ++cell)
2019  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
2020  if (cell->at_boundary(f) &&
2021  (boundary_ids.find(cell->face(f)->boundary_id()) !=
2022  boundary_ids.end()))
2023  {
2024  const unsigned int dofs_per_face = cell->get_fe().dofs_per_face;
2025  dofs_on_face.resize (dofs_per_face);
2026 
2027  cell->face(f)->get_dof_indices (dofs_on_face,
2028  cell->active_fe_index());
2029  for (unsigned int i=0; i<dofs_per_face; ++i)
2030  boundary_dofs.insert(dofs_on_face[i]);
2031  };
2032  return boundary_dofs.size();
2033  }
2034 
2035 
2036 
2037  template <>
2038  types::global_dof_index DoFHandler<2,3>::n_boundary_dofs () const
2039  {
2040  Assert(false,ExcNotImplemented());
2041  return 0;
2042  }
2043 
2044 
2045 
2046  template <>
2047  template <typename number>
2048  types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::map<types::boundary_id, const Function<3,number>*> &) const
2049  {
2050  Assert(false,ExcNotImplemented());
2051  return 0;
2052  }
2053 
2054 
2055 
2056  template <>
2057  types::global_dof_index DoFHandler<2,3>::n_boundary_dofs (const std::set<types::boundary_id> &) const
2058  {
2059  Assert(false,ExcNotImplemented());
2060  return 0;
2061  }
2062 
2063 
2064 
2065  template<int dim, int spacedim>
2066  std::size_t
2067  DoFHandler<dim,spacedim>::memory_consumption () const
2068  {
2069  std::size_t mem = (MemoryConsumption::memory_consumption (tria) +
2070  MemoryConsumption::memory_consumption (finite_elements) +
2071  MemoryConsumption::memory_consumption (tria) +
2072  MemoryConsumption::memory_consumption (levels) +
2073  MemoryConsumption::memory_consumption (*faces) +
2074  MemoryConsumption::memory_consumption (number_cache) +
2075  MemoryConsumption::memory_consumption (vertex_dofs) +
2076  MemoryConsumption::memory_consumption (vertex_dofs_offsets) +
2077  MemoryConsumption::memory_consumption (has_children));
2078  for (unsigned int i=0; i<levels.size(); ++i)
2079  mem += MemoryConsumption::memory_consumption (*levels[i]);
2080  mem += MemoryConsumption::memory_consumption (*faces);
2081 
2082  return mem;
2083  }
2084 
2085 
2086 
2087  template<int dim, int spacedim>
2088  void
2089  DoFHandler<dim,spacedim>::
2090  compute_vertex_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const
2091  {
2092  // Note: we may wish to have
2093  // something here similar to what
2094  // we do for lines and quads,
2095  // namely that we only identify
2096  // dofs for any fe towards the
2097  // most dominating one. however,
2098  // it is not clear whether this
2099  // is actually necessary for
2100  // vertices at all, I can't think
2101  // of a finite element that would
2102  // make that necessary...
2103  Table<2,std_cxx11::shared_ptr<::internal::hp::DoFIdentities> >
2104  vertex_dof_identities (get_fe().size(),
2105  get_fe().size());
2106 
2107  // loop over all vertices and
2108  // see which one we need to
2109  // work on
2110  for (unsigned int vertex_index=0; vertex_index<get_tria().n_vertices();
2111  ++vertex_index)
2112  {
2113  const unsigned int n_active_fe_indices
2114  = ::internal::DoFAccessor::Implementation::
2115  n_active_vertex_fe_indices (*this, vertex_index);
2116  if (n_active_fe_indices > 1)
2117  {
2118  const unsigned int
2119  first_fe_index
2120  = ::internal::DoFAccessor::Implementation::
2121  nth_active_vertex_fe_index (*this, vertex_index, 0);
2122 
2123  // loop over all the
2124  // other FEs with which
2125  // we want to identify
2126  // the DoF indices of
2127  // the first FE of
2128  for (unsigned int f=1; f<n_active_fe_indices; ++f)
2129  {
2130  const unsigned int
2131  other_fe_index
2132  = ::internal::DoFAccessor::Implementation::
2133  nth_active_vertex_fe_index (*this, vertex_index, f);
2134 
2135  // make sure the
2136  // entry in the
2137  // equivalence
2138  // table exists
2139  ::internal::hp::ensure_existence_of_dof_identities<0>
2140  (get_fe()[first_fe_index],
2141  get_fe()[other_fe_index],
2142  vertex_dof_identities[first_fe_index][other_fe_index]);
2143 
2144  // then loop
2145  // through the
2146  // identities we
2147  // have. first get
2148  // the global
2149  // numbers of the
2150  // dofs we want to
2151  // identify and
2152  // make sure they
2153  // are not yet
2154  // constrained to
2155  // anything else,
2156  // except for to
2157  // each other. use
2158  // the rule that we
2159  // will always
2160  // constrain the
2161  // dof with the
2162  // higher fe
2163  // index to the
2164  // one with the
2165  // lower, to avoid
2166  // circular
2167  // reasoning.
2168  ::internal::hp::DoFIdentities &identities
2169  = *vertex_dof_identities[first_fe_index][other_fe_index];
2170  for (unsigned int i=0; i<identities.size(); ++i)
2171  {
2172  const types::global_dof_index lower_dof_index
2173  = ::internal::DoFAccessor::Implementation::
2174  get_vertex_dof_index (*this,
2175  vertex_index,
2176  first_fe_index,
2177  identities[i].first);
2178  const types::global_dof_index higher_dof_index
2179  = ::internal::DoFAccessor::Implementation::
2180  get_vertex_dof_index (*this,
2181  vertex_index,
2182  other_fe_index,
2183  identities[i].second);
2184 
2185  Assert ((new_dof_indices[higher_dof_index] ==
2186  numbers::invalid_dof_index)
2187  ||
2188  (new_dof_indices[higher_dof_index] ==
2189  lower_dof_index),
2190  ExcInternalError());
2191 
2192  new_dof_indices[higher_dof_index] = lower_dof_index;
2193  }
2194  }
2195  }
2196  }
2197  }
2198 
2199 
2200  template <>
2201  void
2202  DoFHandler<1,1>::
2203  compute_line_dof_identities (std::vector<types::global_dof_index> &) const
2204  {}
2205 
2206 
2207 
2208  template <>
2209  void
2210  DoFHandler<1,2>::
2211  compute_line_dof_identities (std::vector<types::global_dof_index> &) const
2212  {}
2213 
2214  template <>
2215  void
2216  DoFHandler<1,3>::
2217  compute_line_dof_identities (std::vector<types::global_dof_index> &) const
2218  {}
2219 
2220 
2221  template<int dim, int spacedim>
2222  void
2223  DoFHandler<dim,spacedim>::
2224  compute_line_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const
2225  {
2226  // we will mark lines that we have already treated, so first save and clear
2227  // the user flags on lines and later restore them
2228  std::vector<bool> user_flags;
2229  this->get_triangulation().save_user_flags_line(user_flags);
2230  const_cast<Triangulation<dim,spacedim> &>(this->get_triangulation()).clear_user_flags_line ();
2231 
2232  // An implementation of the algorithm described in the hp paper, including
2233  // the modification mentioned later in the "complications in 3-d" subsections
2234  //
2235  // as explained there, we do something only if there are exactly 2 finite
2236  // elements associated with an object. if there is only one, then there is
2237  // nothing to do anyway, and if there are 3 or more, then we can get into
2238  // trouble. note that this only happens for lines in 3d and higher, and for
2239  // quads only in 4d and higher, so this isn't a particularly frequent case
2240  //
2241  // there is one case, however, that we would like to handle (see, for
2242  // example, the hp/crash_15 testcase): if we have FESystem(FE_Q(2),FE_DGQ(i))
2243  // elements for a bunch of values 'i', then we should be able to handle this
2244  // because we can simply unify *all* dofs, not only a some. so what we do
2245  // is to first treat all pairs of finite elements that have *identical* dofs,
2246  // and then only deal with those that are not identical of which we can
2247  // handle at most 2
2248  Table<2,std_cxx11::shared_ptr<internal::hp::DoFIdentities> >
2249  line_dof_identities (finite_elements->size(),
2250  finite_elements->size());
2251 
2252  for (active_cell_iterator cell=begin_active(); cell!=end(); ++cell)
2253  for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
2254  if (cell->line(l)->user_flag_set() == false)
2255  {
2256  const line_iterator line = cell->line(l);
2257  line->set_user_flag ();
2258 
2259  unsigned int unique_sets_of_dofs
2260  = line->n_active_fe_indices();
2261 
2262  // do a first loop over all sets of dofs and do identity
2263  // uniquification
2264  for (unsigned int f=0; f<line->n_active_fe_indices(); ++f)
2265  for (unsigned int g=f+1; g<line->n_active_fe_indices(); ++g)
2266  {
2267  const unsigned int fe_index_1 = line->nth_active_fe_index (f),
2268  fe_index_2 = line->nth_active_fe_index (g);
2269 
2270  if (((*finite_elements)[fe_index_1].dofs_per_line
2271  ==
2272  (*finite_elements)[fe_index_2].dofs_per_line)
2273  &&
2274  ((*finite_elements)[fe_index_1].dofs_per_line > 0))
2275  {
2276  internal::hp::ensure_existence_of_dof_identities<1>
2277  ((*finite_elements)[fe_index_1],
2278  (*finite_elements)[fe_index_2],
2279  line_dof_identities[fe_index_1][fe_index_2]);
2280  // see if these sets of dofs are identical. the first
2281  // condition for this is that indeed there are n identities
2282  if (line_dof_identities[fe_index_1][fe_index_2]->size()
2283  ==
2284  (*finite_elements)[fe_index_1].dofs_per_line)
2285  {
2286  unsigned int i=0;
2287  for (; i<(*finite_elements)[fe_index_1].dofs_per_line; ++i)
2288  if (((*(line_dof_identities[fe_index_1][fe_index_2]))[i].first != i)
2289  &&
2290  ((*(line_dof_identities[fe_index_1][fe_index_2]))[i].second != i))
2291  // not an identity
2292  break;
2293 
2294  if (i == (*finite_elements)[fe_index_1].dofs_per_line)
2295  {
2296  // The line dofs (i.e., the ones interior to a line) of these two finite elements are identical.
2297  // Note that there could be situations when one element still dominates another, e.g.:
2298  // FE_Q(2) x FE_Nothing(dominate) vs
2299  // FE_Q(2) x FE_Q(1)
2300 
2301  --unique_sets_of_dofs;
2302 
2303  for (unsigned int j=0; j<(*finite_elements)[fe_index_1].dofs_per_line; ++j)
2304  {
2305  const types::global_dof_index master_dof_index
2306  = line->dof_index (j, fe_index_1);
2307  const types::global_dof_index slave_dof_index
2308  = line->dof_index (j, fe_index_2);
2309 
2310  // if master dof was already constrained,
2311  // constrain to that one, otherwise constrain
2312  // slave to master
2313  if (new_dof_indices[master_dof_index] !=
2314  numbers::invalid_dof_index)
2315  {
2316  Assert (new_dof_indices[new_dof_indices[master_dof_index]] ==
2317  numbers::invalid_dof_index,
2318  ExcInternalError());
2319 
2320  new_dof_indices[slave_dof_index]
2321  = new_dof_indices[master_dof_index];
2322  }
2323  else
2324  {
2325  Assert ((new_dof_indices[master_dof_index] ==
2326  numbers::invalid_dof_index)
2327  ||
2328  (new_dof_indices[slave_dof_index] ==
2329  master_dof_index),
2330  ExcInternalError());
2331 
2332  new_dof_indices[slave_dof_index] = master_dof_index;
2333  }
2334  }
2335  }
2336  }
2337  }
2338  }
2339 
2340  // if at this point, there is only one unique set of dofs left, then
2341  // we have taken care of everything above. if there are two, then we
2342  // need to deal with them here. if there are more, then we punt, as
2343  // described in the paper (and mentioned above)
2344 //TODO: The check for 'dim==2' was inserted by intuition. It fixes
2345 // the previous problems with @ref step_27 "step-27" in 3D. But an explanation
2346 // for this is still required, and what we do here is not what we
2347 // describe in the paper!.
2348  if ((unique_sets_of_dofs == 2) && (dim == 2))
2349  {
2350  // find out which is the most dominating finite element of the
2351  // ones that are used on this line
2352  const unsigned int most_dominating_fe_index
2353  = internal::hp::get_most_dominating_fe_index<dim,spacedim> (line);
2354 
2355  // if we found the most dominating element, then use this to eliminate some of
2356  // the degrees of freedom by identification. otherwise, the code that computes
2357  // hanging node constraints will have to deal with it by computing
2358  // appropriate constraints along this face/edge
2359  if (most_dominating_fe_index != numbers::invalid_unsigned_int)
2360  {
2361  const unsigned int n_active_fe_indices
2362  = line->n_active_fe_indices ();
2363 
2364  // loop over the indices of all the finite elements that are not
2365  // dominating, and identify their dofs to the most dominating
2366  // one
2367  for (unsigned int f=0; f<n_active_fe_indices; ++f)
2368  if (line->nth_active_fe_index (f) !=
2369  most_dominating_fe_index)
2370  {
2371  const unsigned int
2372  other_fe_index = line->nth_active_fe_index (f);
2373 
2374  internal::hp::ensure_existence_of_dof_identities<1>
2375  ((*finite_elements)[most_dominating_fe_index],
2376  (*finite_elements)[other_fe_index],
2377  line_dof_identities[most_dominating_fe_index][other_fe_index]);
2378 
2379  internal::hp::DoFIdentities &identities
2380  = *line_dof_identities[most_dominating_fe_index][other_fe_index];
2381  for (unsigned int i=0; i<identities.size(); ++i)
2382  {
2383  const types::global_dof_index master_dof_index
2384  = line->dof_index (identities[i].first, most_dominating_fe_index);
2385  const types::global_dof_index slave_dof_index
2386  = line->dof_index (identities[i].second, other_fe_index);
2387 
2388  Assert ((new_dof_indices[master_dof_index] ==
2389  numbers::invalid_dof_index)
2390  ||
2391  (new_dof_indices[slave_dof_index] ==
2392  master_dof_index),
2393  ExcInternalError());
2394 
2395  new_dof_indices[slave_dof_index] = master_dof_index;
2396  }
2397  }
2398  }
2399  }
2400  }
2401 
2402  // finally restore the user flags
2403  const_cast<Triangulation<dim,spacedim> &>(this->get_triangulation())
2404  .load_user_flags_line(user_flags);
2405  }
2406 
2407 
2408 
2409  template <int dim, int spacedim>
2410  void
2411  DoFHandler<dim,spacedim>::
2412  compute_quad_dof_identities (std::vector<types::global_dof_index> &) const
2413  {
2414  // this function should only be called for dim<3 where there are
2415  // no quad dof identies. for dim>=3, the specialization below should
2416  // take care of it
2417  Assert (dim < 3, ExcInternalError());
2418  }
2419 
2420 
2421  template <>
2422  void
2423  DoFHandler<3,3>::
2424  compute_quad_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const
2425  {
2426  const int dim = 3;
2427  const int spacedim = 3;
2428 
2429  // we will mark quads that we
2430  // have already treated, so first
2431  // save and clear the user flags
2432  // on quads and later restore
2433  // them
2434  std::vector<bool> user_flags;
2435  this->get_triangulation().save_user_flags_quad(user_flags);
2436  const_cast<Triangulation<dim,spacedim> &>(this->get_triangulation()).clear_user_flags_quad ();
2437 
2438  // An implementation of the
2439  // algorithm described in the hp
2440  // paper, including the
2441  // modification mentioned later
2442  // in the "complications in 3-d"
2443  // subsections
2444  //
2445  // as explained there, we do
2446  // something only if there are
2447  // exactly 2 finite elements
2448  // associated with an object. if
2449  // there is only one, then there
2450  // is nothing to do anyway, and
2451  // if there are 3 or more, then
2452  // we can get into trouble. note
2453  // that this only happens for
2454  // lines in 3d and higher, and
2455  // for quads only in 4d and
2456  // higher, so this isn't a
2457  // particularly frequent case
2458  Table<2,std_cxx11::shared_ptr<internal::hp::DoFIdentities> >
2459  quad_dof_identities (finite_elements->size(),
2460  finite_elements->size());
2461 
2462  for (active_cell_iterator cell=begin_active(); cell!=end(); ++cell)
2463  for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
2464  if ((cell->quad(q)->user_flag_set() == false)
2465  &&
2466  (cell->quad(q)->n_active_fe_indices() == 2))
2467  {
2468  const quad_iterator quad = cell->quad(q);
2469  quad->set_user_flag ();
2470 
2471  // find out which is the
2472  // most dominating finite
2473  // element of the ones that
2474  // are used on this quad
2475  const unsigned int most_dominating_fe_index
2476  = internal::hp::get_most_dominating_fe_index<dim,spacedim> (quad);
2477 
2478  // if we found the most dominating element, then use this to eliminate some of
2479  // the degrees of freedom by identification. otherwise, the code that computes
2480  // hanging node constraints will have to deal with it by computing
2481  // appropriate constraints along this face/edge
2482  if (most_dominating_fe_index != numbers::invalid_unsigned_int)
2483  {
2484  const unsigned int n_active_fe_indices
2485  = quad->n_active_fe_indices ();
2486 
2487  // loop over the indices of
2488  // all the finite elements
2489  // that are not dominating,
2490  // and identify their dofs
2491  // to the most dominating
2492  // one
2493  for (unsigned int f=0; f<n_active_fe_indices; ++f)
2494  if (quad->nth_active_fe_index (f) !=
2495  most_dominating_fe_index)
2496  {
2497  const unsigned int
2498  other_fe_index = quad->nth_active_fe_index (f);
2499 
2500  internal::hp::ensure_existence_of_dof_identities<2>
2501  ((*finite_elements)[most_dominating_fe_index],
2502  (*finite_elements)[other_fe_index],
2503  quad_dof_identities[most_dominating_fe_index][other_fe_index]);
2504 
2505  internal::hp::DoFIdentities &identities
2506  = *quad_dof_identities[most_dominating_fe_index][other_fe_index];
2507  for (unsigned int i=0; i<identities.size(); ++i)
2508  {
2509  const types::global_dof_index master_dof_index
2510  = quad->dof_index (identities[i].first, most_dominating_fe_index);
2511  const types::global_dof_index slave_dof_index
2512  = quad->dof_index (identities[i].second, other_fe_index);
2513 
2514  Assert ((new_dof_indices[master_dof_index] ==
2515  numbers::invalid_dof_index)
2516  ||
2517  (new_dof_indices[slave_dof_index] ==
2518  master_dof_index),
2519  ExcInternalError());
2520 
2521  new_dof_indices[slave_dof_index] = master_dof_index;
2522  }
2523  }
2524  }
2525  }
2526 
2527  // finally restore the user flags
2528  const_cast<Triangulation<dim,spacedim> &>(this->get_triangulation())
2529  .load_user_flags_quad(user_flags);
2530  }
2531 
2532 
2533 
2534  template <int dim, int spacedim>
2535  void DoFHandler<dim,spacedim>::set_active_fe_indices (const std::vector<unsigned int> &active_fe_indices)
2536  {
2537  Assert(active_fe_indices.size()==get_tria().n_active_cells(),
2538  ExcDimensionMismatch(active_fe_indices.size(), get_tria().n_active_cells()));
2539 
2540  create_active_fe_table ();
2541  // we could set the values directly, since
2542  // they are stored as protected data of
2543  // this object, but for simplicity we use
2544  // the cell-wise access. this way we also
2545  // have to pass some debug-mode tests which
2546  // we would have to duplicate ourselves
2547  // otherwise
2548  active_cell_iterator cell=begin_active(),
2549  endc=end();
2550  for (unsigned int i=0; cell!=endc; ++cell, ++i)
2551  cell->set_active_fe_index(active_fe_indices[i]);
2552  }
2553 
2554 
2555 
2556  template <int dim, int spacedim>
2557  void DoFHandler<dim,spacedim>::get_active_fe_indices (std::vector<unsigned int> &active_fe_indices) const
2558  {
2559  active_fe_indices.resize(get_tria().n_active_cells());
2560 
2561  // we could try to extract the values directly, since
2562  // they are stored as protected data of
2563  // this object, but for simplicity we use
2564  // the cell-wise access.
2565  active_cell_iterator cell=begin_active(),
2566  endc=end();
2567  for (unsigned int i=0; cell!=endc; ++cell, ++i)
2568  active_fe_indices[i]=cell->active_fe_index();
2569  }
2570 
2571 
2572 
2573  template<int dim, int spacedim>
2574  void DoFHandler<dim,spacedim>::distribute_dofs (const hp::FECollection<dim,spacedim> &ff)
2575  {
2576  Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
2577 
2578  finite_elements = &ff;
2579 
2580  // This call ensures that the
2581  // active_fe_indices vectors are
2582  // initialized correctly.
2583  create_active_fe_table ();
2584 
2585  // up front make sure that the fe
2586  // collection is large enough to
2587  // cover all fe indices presently
2588  // in use on the mesh
2589  for (active_cell_iterator cell = begin_active(); cell != end(); ++cell)
2590  Assert (cell->active_fe_index() < finite_elements->size(),
2591  ExcInvalidFEIndex (cell->active_fe_index(),
2592  finite_elements->size()));
2593 
2594 
2595  // then allocate space for all
2596  // the other tables
2597  ::internal::hp::DoFHandler::Implementation::reserve_space (*this);
2598 
2599  // Clear user flags because we will
2600  // need them. But first we save
2601  // them and make sure that we
2602  // restore them later such that at
2603  // the end of this function the
2604  // Triangulation will be in the
2605  // same state as it was at the
2606  // beginning of this function.
2607  std::vector<bool> user_flags;
2608  tria->save_user_flags(user_flags);
2609  const_cast<Triangulation<dim,spacedim> &>(*tria).clear_user_flags ();
2610 
2611 
2613 
2614  // Step 1: distribute DoFs on all
2615  // active entities
2616  {
2617  types::global_dof_index next_free_dof = 0;
2618  active_cell_iterator cell = begin_active(),
2619  endc = end();
2620 
2621  for (; cell != endc; ++cell)
2622  next_free_dof
2623  = ::internal::hp::DoFHandler::Implementation::distribute_dofs_on_cell<spacedim> (cell,
2624  next_free_dof);
2625 
2626  number_cache.n_global_dofs = next_free_dof;
2627  }
2628 
2629 
2631 
2632  // Step 2: identify certain dofs
2633  // if the finite element tells us
2634  // that they should have the same
2635  // value. only pertinent for
2636  // faces and other
2637  // lower-dimensional objects
2638  // where elements come together
2639  std::vector<types::global_dof_index>
2640  constrained_indices (number_cache.n_global_dofs, numbers::invalid_dof_index);
2641  compute_vertex_dof_identities (constrained_indices);
2642  compute_line_dof_identities (constrained_indices);
2643  compute_quad_dof_identities (constrained_indices);
2644 
2645  // loop over all dofs and assign
2646  // new numbers to those which are
2647  // not constrained
2648  std::vector<types::global_dof_index>
2649  new_dof_indices (number_cache.n_global_dofs, numbers::invalid_dof_index);
2650  types::global_dof_index next_free_dof = 0;
2651  for (types::global_dof_index i=0; i<number_cache.n_global_dofs; ++i)
2652  if (constrained_indices[i] == numbers::invalid_dof_index)
2653  {
2654  new_dof_indices[i] = next_free_dof;
2655  ++next_free_dof;
2656  }
2657 
2658  // then loop over all those that
2659  // are constrained and record the
2660  // new dof number for those:
2661  for (types::global_dof_index i=0; i<number_cache.n_global_dofs; ++i)
2662  if (constrained_indices[i] != numbers::invalid_dof_index)
2663  {
2664  Assert (new_dof_indices[constrained_indices[i]] !=
2665  numbers::invalid_dof_index,
2666  ExcInternalError());
2667 
2668  new_dof_indices[i] = new_dof_indices[constrained_indices[i]];
2669  }
2670 
2671  for (types::global_dof_index i=0; i<number_cache.n_global_dofs; ++i)
2672  {
2673  Assert (new_dof_indices[i] != numbers::invalid_dof_index,
2674  ExcInternalError());
2675  Assert (new_dof_indices[i] < next_free_dof,
2676  ExcInternalError());
2677  }
2678 
2679  // finally, do the renumbering
2680  // and set the number of actually
2681  // used dof indices
2682  renumber_dofs_internal (new_dof_indices, ::internal::int2type<dim>());
2683 
2684  // now set the elements of the
2685  // number cache appropriately
2686  number_cache.n_global_dofs = next_free_dof;
2687  number_cache.n_locally_owned_dofs = number_cache.n_global_dofs;
2688 
2689  if (dynamic_cast<const parallel::shared::Triangulation< dim, spacedim >*>
2690  (&this->get_triangulation())
2691  == 0)
2692  {
2693  number_cache.locally_owned_dofs
2694  = IndexSet (number_cache.n_global_dofs);
2695  number_cache.locally_owned_dofs.add_range (0,
2696  number_cache.n_global_dofs);
2697  Assert (number_cache.n_global_dofs < std::numeric_limits<unsigned int>::max (),
2698  ExcMessage ("Global number of degrees of freedom is too large."));
2699  number_cache.n_locally_owned_dofs_per_processor
2700  = std::vector<types::global_dof_index> (1,
2701  (types::global_dof_index) number_cache.n_global_dofs);
2702  }
2703  else
2704  {
2705  AssertThrow(false, ExcNotImplemented() );
2706  //number_cache.locally_owned_dofs = ::DoFTools::locally_owned_dofs_with_subdomain(this,tria->locally_owned_subdomain() );
2707  //TODO: update n_locally_owned_dofs_per_processor as well
2708  }
2709 
2710  number_cache.locally_owned_dofs_per_processor
2711  = std::vector<IndexSet> (1,
2712  number_cache.locally_owned_dofs);
2713 
2714  // update the cache used for cell dof indices and compress the data on the levels. do
2715  // the latter on separate tasks to gain parallelism, starting with the highest
2716  // level (there is most to do there, so start it first)
2717  for (active_cell_iterator cell = begin_active();
2718  cell != end(); ++cell)
2719  cell->update_cell_dof_indices_cache ();
2720 
2721  {
2722  Threads::TaskGroup<> tg;
2723  for (int level=levels.size()-1; level>=0; --level)
2724  tg += Threads::new_task (&::internal::hp::DoFLevel::compress_data<dim,spacedim>,
2725  *levels[level], *finite_elements);
2726  tg.join_all ();
2727  }
2728 
2729  // finally restore the user flags
2730  const_cast<Triangulation<dim,spacedim> &>(*tria).load_user_flags(user_flags);
2731  }
2732 
2733 
2734 
2735  template<int dim, int spacedim>
2736  void DoFHandler<dim,spacedim>::clear ()
2737  {
2738  // release lock to old fe
2739  finite_elements = 0;
2740 
2741  // release memory
2742  clear_space ();
2743  }
2744 
2745 
2746 
2747  template<int dim, int spacedim>
2748  void DoFHandler<dim,spacedim>::renumber_dofs (const std::vector<types::global_dof_index> &new_numbers)
2749  {
2750  Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
2751 #ifdef DEBUG
2752  // assert that the new indices are
2753  // consecutively numbered
2754  if (true)
2755  {
2756  std::vector<types::global_dof_index> tmp(new_numbers);
2757  std::sort (tmp.begin(), tmp.end());
2758  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2759  types::global_dof_index i = 0;
2760  for (; p!=tmp.end(); ++p, ++i)
2761  Assert (*p == i, ExcNewNumbersNotConsecutive(i));
2762  }
2763 #endif
2764 
2765  // uncompress the internal storage scheme of dofs on cells
2766  // so that we can access dofs in turns. uncompress in parallel, starting
2767  // with the most expensive levels (the highest ones)
2768  {
2769  Threads::TaskGroup<> tg;
2770  for (int level=levels.size()-1; level>=0; --level)
2771  tg += Threads::new_task (&::internal::hp::DoFLevel::uncompress_data<dim,spacedim>,
2772  *levels[level], *finite_elements);
2773  tg.join_all ();
2774  }
2775 
2776  // do the renumbering
2777  renumber_dofs_internal (new_numbers, ::internal::int2type<dim>());
2778 
2779  // update the cache used for cell dof indices
2780  for (active_cell_iterator cell = begin_active();
2781  cell != end(); ++cell)
2782  cell->update_cell_dof_indices_cache ();
2783 
2784  // now re-compress the dof indices
2785  {
2786  Threads::TaskGroup<> tg;
2787  for (int level=levels.size()-1; level>=0; --level)
2788  tg += Threads::new_task (&::internal::hp::DoFLevel::compress_data<dim,spacedim>,
2789  *levels[level], *finite_elements);
2790  tg.join_all ();
2791  }
2792  }
2793 
2794 
2795 
2796  template<int dim, int spacedim>
2797  void
2798  DoFHandler<dim,spacedim>::
2799  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
2800  ::internal::int2type<0>)
2801  {
2802  Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
2803 
2804  for (unsigned int vertex_index=0; vertex_index<get_tria().n_vertices();
2805  ++vertex_index)
2806  {
2807  const unsigned int n_active_fe_indices
2808  = ::internal::DoFAccessor::Implementation::
2809  n_active_vertex_fe_indices (*this, vertex_index);
2810 
2811  for (unsigned int f=0; f<n_active_fe_indices; ++f)
2812  {
2813  const unsigned int fe_index
2814  = ::internal::DoFAccessor::Implementation::
2815  nth_active_vertex_fe_index (*this, vertex_index, f);
2816 
2817  for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_vertex; ++d)
2818  {
2819  const types::global_dof_index vertex_dof_index
2820  = ::internal::DoFAccessor::Implementation::
2821  get_vertex_dof_index(*this,
2822  vertex_index,
2823  fe_index,
2824  d);
2825  ::internal::DoFAccessor::Implementation::
2826  set_vertex_dof_index (*this,
2827  vertex_index,
2828  fe_index,
2829  d,
2830  new_numbers[vertex_dof_index]);
2831  }
2832  }
2833  }
2834  }
2835 
2836 
2837 
2838  template<int dim, int spacedim>
2839  void
2840  DoFHandler<dim,spacedim>::
2841  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
2842  ::internal::int2type<1>)
2843  {
2844  Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
2845 
2846  renumber_dofs_internal (new_numbers, internal::int2type<0>());
2847 
2848  // save user flags on lines so we
2849  // can use them to mark lines
2850  // we've already treated
2851  std::vector<bool> saved_line_user_flags;
2852  const_cast<::Triangulation<dim,spacedim>&>(*tria)
2853  .save_user_flags_line (saved_line_user_flags);
2854  const_cast<::Triangulation<dim,spacedim>&>(*tria)
2855  .clear_user_flags_line ();
2856 
2857  for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
2858  for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
2859  if (cell->line(l)->user_flag_set() == false)
2860  {
2861  const line_iterator line = cell->line(l);
2862  line->set_user_flag();
2863 
2864  const unsigned int n_active_fe_indices
2865  = line->n_active_fe_indices ();
2866 
2867  for (unsigned int f=0; f<n_active_fe_indices; ++f)
2868  {
2869  const unsigned int fe_index
2870  = line->nth_active_fe_index (f);
2871 
2872  for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_line; ++d)
2873  line->set_dof_index (d,
2874  new_numbers[line->dof_index(d,fe_index)],
2875  fe_index);
2876  }
2877  }
2878 
2879  // at the end, restore the user
2880  // flags for the lines
2881  const_cast<::Triangulation<dim,spacedim>&>(*tria)
2882  .load_user_flags_line (saved_line_user_flags);
2883  }
2884 
2885 
2886 
2887 //TODO: Merge the following three functions -- they are identical
2888  template<>
2889  void
2890  DoFHandler<2,2>::
2891  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
2892  ::internal::int2type<2>)
2893  {
2894  const unsigned int dim = 2;
2895  const unsigned int spacedim = 2;
2896 
2897  Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
2898 
2899  renumber_dofs_internal (new_numbers, internal::int2type<1>());
2900 
2901  // save user flags on quads so we
2902  // can use them to mark quads
2903  // we've already treated
2904  std::vector<bool> saved_quad_user_flags;
2905  const_cast<::Triangulation<dim,spacedim>&>(*tria)
2906  .save_user_flags_quad (saved_quad_user_flags);
2907  const_cast<::Triangulation<dim,spacedim>&>(*tria)
2908  .clear_user_flags_quad ();
2909 
2910  for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
2911  for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
2912  if (cell->quad(q)->user_flag_set() == false)
2913  {
2914  const quad_iterator quad = cell->quad(q);
2915  quad->set_user_flag();
2916 
2917  const unsigned int n_active_fe_indices
2918  = quad->n_active_fe_indices ();
2919 
2920  for (unsigned int f=0; f<n_active_fe_indices; ++f)
2921  {
2922  const unsigned int fe_index
2923  = quad->nth_active_fe_index (f);
2924 
2925  for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_quad; ++d)
2926  quad->set_dof_index (d,
2927  new_numbers[quad->dof_index(d,fe_index)],
2928  fe_index);
2929  }
2930  }
2931 
2932  // at the end, restore the user
2933  // flags for the quads
2934  const_cast<::Triangulation<dim,spacedim>&>(*tria)
2935  .load_user_flags_quad (saved_quad_user_flags);
2936  }
2937 
2938 
2939 
2940  template<>
2941  void
2942  DoFHandler<2,3>::
2943  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
2944  ::internal::int2type<2>)
2945  {
2946  const unsigned int dim = 2;
2947  const unsigned int spacedim = 3;
2948 
2949  Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
2950 
2951  renumber_dofs_internal (new_numbers, internal::int2type<1>());
2952 
2953  // save user flags on quads so we
2954  // can use them to mark quads
2955  // we've already treated
2956  std::vector<bool> saved_quad_user_flags;
2957  const_cast<::Triangulation<dim,spacedim>&>(*tria)
2958  .save_user_flags_quad (saved_quad_user_flags);
2959  const_cast<::Triangulation<dim,spacedim>&>(*tria)
2960  .clear_user_flags_quad ();
2961 
2962  for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
2963  for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
2964  if (cell->quad(q)->user_flag_set() == false)
2965  {
2966  const quad_iterator quad = cell->quad(q);
2967  quad->set_user_flag();
2968 
2969  const unsigned int n_active_fe_indices
2970  = quad->n_active_fe_indices ();
2971 
2972  for (unsigned int f=0; f<n_active_fe_indices; ++f)
2973  {
2974  const unsigned int fe_index
2975  = quad->nth_active_fe_index (f);
2976 
2977  for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_quad; ++d)
2978  quad->set_dof_index (d,
2979  new_numbers[quad->dof_index(d,fe_index)],
2980  fe_index);
2981  }
2982  }
2983 
2984  // at the end, restore the user
2985  // flags for the quads
2986  const_cast<::Triangulation<dim,spacedim>&>(*tria)
2987  .load_user_flags_quad (saved_quad_user_flags);
2988  }
2989 
2990 
2991  template<>
2992  void
2993  DoFHandler<3,3>::
2994  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
2995  ::internal::int2type<2>)
2996  {
2997  const unsigned int dim = 3;
2998  const unsigned int spacedim = 3;
2999 
3000  Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
3001 
3002  renumber_dofs_internal (new_numbers, internal::int2type<1>());
3003 
3004  // save user flags on quads so we
3005  // can use them to mark quads
3006  // we've already treated
3007  std::vector<bool> saved_quad_user_flags;
3008  const_cast<::Triangulation<dim,spacedim>&>(*tria)
3009  .save_user_flags_quad (saved_quad_user_flags);
3010  const_cast<::Triangulation<dim,spacedim>&>(*tria)
3011  .clear_user_flags_quad ();
3012 
3013  for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
3014  for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
3015  if (cell->quad(q)->user_flag_set() == false)
3016  {
3017  const quad_iterator quad = cell->quad(q);
3018  quad->set_user_flag();
3019 
3020  const unsigned int n_active_fe_indices
3021  = quad->n_active_fe_indices ();
3022 
3023  for (unsigned int f=0; f<n_active_fe_indices; ++f)
3024  {
3025  const unsigned int fe_index
3026  = quad->nth_active_fe_index (f);
3027 
3028  for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_quad; ++d)
3029  quad->set_dof_index (d,
3030  new_numbers[quad->dof_index(d,fe_index)],
3031  fe_index);
3032  }
3033  }
3034 
3035  // at the end, restore the user
3036  // flags for the quads
3037  const_cast<::Triangulation<dim,spacedim>&>(*tria)
3038  .load_user_flags_quad (saved_quad_user_flags);
3039  }
3040 
3041 
3042  template<>
3043  void
3044  DoFHandler<3,3>::
3045  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
3046  ::internal::int2type<3>)
3047  {
3048  const unsigned int dim = 3;
3049  const unsigned int spacedim = 3;
3050 
3051  Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
3052 
3053  renumber_dofs_internal (new_numbers, internal::int2type<2>());
3054 
3055  // save user flags on hexes so we
3056  // can use them to mark hexes
3057  // we've already treated
3058  std::vector<bool> saved_hex_user_flags;
3059  const_cast<::Triangulation<dim,spacedim>&>(*tria)
3060  .save_user_flags_hex (saved_hex_user_flags);
3061  const_cast<::Triangulation<dim,spacedim>&>(*tria)
3062  .clear_user_flags_hex ();
3063 
3064  // we're in 3d, so hexes are also
3065  // cells. stick with the same
3066  // kind of notation as in the
3067  // previous functions, though
3068  for (active_cell_iterator cell = begin_active(); cell!=end(); ++cell)
3069  if (cell->user_flag_set() == false)
3070  {
3071  const hex_iterator hex = cell;
3072  hex->set_user_flag();
3073 
3074  const unsigned int n_active_fe_indices
3075  = hex->n_active_fe_indices ();
3076 
3077  for (unsigned int f=0; f<n_active_fe_indices; ++f)
3078  {
3079  const unsigned int fe_index
3080  = hex->nth_active_fe_index (f);
3081 
3082  for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_hex; ++d)
3083  hex->set_dof_index (d,
3084  new_numbers[hex->dof_index(d,fe_index)],
3085  fe_index);
3086  }
3087  }
3088 
3089  // at the end, restore the user
3090  // flags for the hexs
3091  const_cast<::Triangulation<dim,spacedim>&>(*tria)
3092  .load_user_flags_hex (saved_hex_user_flags);
3093  }
3094 
3095 
3096 
3097  template <int dim, int spacedim>
3098  unsigned int
3099  DoFHandler<dim, spacedim>::max_couplings_between_dofs () const
3100  {
3101  Assert (finite_elements != 0, ExcNoFESelected());
3102  return ::internal::hp::DoFHandler::Implementation::max_couplings_between_dofs (*this);
3103  }
3104 
3105 
3106 
3107  template <int dim, int spacedim>
3108  unsigned int
3109  DoFHandler<dim,spacedim>::max_couplings_between_boundary_dofs () const
3110  {
3111  Assert (finite_elements != 0, ExcNoFESelected());
3112 
3113  switch (dim)
3114  {
3115  case 1:
3116  return finite_elements->max_dofs_per_vertex();
3117  case 2:
3118  return (3*finite_elements->max_dofs_per_vertex()
3119  +
3120  2*finite_elements->max_dofs_per_line());
3121  case 3:
3122  // we need to take refinement of
3123  // one boundary face into consideration
3124  // here; in fact, this function returns
3125  // what #max_coupling_between_dofs<2>
3126  // returns
3127  //
3128  // we assume here, that only four faces
3129  // meet at the boundary; this assumption
3130  // is not justified and needs to be
3131  // fixed some time. fortunately, omitting
3132  // it for now does no harm since the
3133  // matrix will cry foul if its requirements
3134  // are not satisfied
3135  return (19*finite_elements->max_dofs_per_vertex() +
3136  28*finite_elements->max_dofs_per_line() +
3137  8*finite_elements->max_dofs_per_quad());
3138  default:
3139  Assert (false, ExcNotImplemented());
3140  return 0;
3141  }
3142  }
3143 
3144 
3145 
3146  template<int dim, int spacedim>
3147  void DoFHandler<dim,spacedim>::create_active_fe_table ()
3148  {
3149  // Create sufficiently many
3150  // hp::DoFLevels.
3151  while (levels.size () < tria->n_levels ())
3152  levels.push_back (new ::internal::hp::DoFLevel);
3153 
3154  // then make sure that on each
3155  // level we have the appropriate
3156  // size of active_fe_indices;
3157  // preset them to zero, i.e. the
3158  // default FE
3159  for (unsigned int level=0; level<levels.size(); ++level)
3160  {
3161  if (levels[level]->active_fe_indices.size () == 0)
3162  levels[level]->active_fe_indices.resize (tria->n_raw_cells(level),
3163  0);
3164  else
3165  {
3166  // Either the
3167  // active_fe_indices have
3168  // size zero because they
3169  // were just created, or
3170  // the correct
3171  // size. Other sizes
3172  // indicate that
3173  // something went wrong.
3174  Assert (levels[level]->active_fe_indices.size () ==
3175  tria->n_raw_cells(level),
3176  ExcInternalError ());
3177  }
3178 
3179  // it may be that the previous table was compressed; in that
3180  // case, restore the correct active_fe_index. the fact that
3181  // this no longer matches the indices in the table is of no
3182  // importance because the current function is called at a
3183  // point where we have to recreate the dof_indices tables in
3184  // the levels anyway
3185  levels[level]->normalize_active_fe_indices ();
3186  }
3187  }
3188 
3189 
3190  template <int dim, int spacedim>
3191  void DoFHandler<dim,spacedim>::pre_refinement_action ()
3192  {
3193  create_active_fe_table ();
3194 
3195  // Remember if the cells already have
3196  // children. That will make the transfer
3197  // of the active_fe_index to the finer
3198  // levels easier.
3199  Assert (has_children.size () == 0, ExcInternalError ());
3200  for (unsigned int i=0; i<levels.size(); ++i)
3201  {
3202  const unsigned int cells_on_level = tria->n_raw_cells(i);
3203  std::vector<bool> *has_children_level =
3204  new std::vector<bool> (cells_on_level);
3205 
3206  // Check for each cell, if it has children. in 1d,
3207  // we don't store refinement cases, so use the 'children'
3208  // vector instead
3209  if (dim == 1)
3210  std::transform (tria->levels[i]->cells.children.begin (),
3211  tria->levels[i]->cells.children.end (),
3212  has_children_level->begin (),
3213  std_cxx11::bind(std::not_equal_to<int>(),
3214  std_cxx11::_1,
3215  -1));
3216  else
3217  std::transform (tria->levels[i]->cells.refinement_cases.begin (),
3218  tria->levels[i]->cells.refinement_cases.end (),
3219  has_children_level->begin (),
3220  std_cxx11::bind (std::not_equal_to<unsigned char>(),
3221  std_cxx11::_1,
3222  static_cast<unsigned char>(RefinementCase<dim>::no_refinement)));
3223 
3224  has_children.push_back (has_children_level);
3225  }
3226  }
3227 
3228 
3229 
3230  template<int dim, int spacedim>
3231  void
3232  DoFHandler<dim,spacedim>::post_refinement_action ()
3233  {
3234  Assert (has_children.size () == levels.size (), ExcInternalError ());
3235 
3236  // Normally only one level is added, but if this Triangulation
3237  // is created by copy_triangulation, it can be more than one level.
3238  while (levels.size () < tria->n_levels ())
3239  levels.push_back (new ::internal::hp::DoFLevel);
3240 
3241  // Coarsening can lead to the loss
3242  // of levels. Hence remove them.
3243  while (levels.size () > tria->n_levels ())
3244  {
3245  delete levels[levels.size ()-1];
3246  levels.pop_back ();
3247  }
3248 
3249  Assert(levels.size () == tria->n_levels (), ExcInternalError());
3250 
3251  // Resize active_fe_indices
3252  // vectors. use zero indicator to
3253  // extend
3254  for (unsigned int i=0; i<levels.size(); ++i)
3255  levels[i]->active_fe_indices.resize (tria->n_raw_cells(i), 0);
3256 
3257  // if a finite element collection
3258  // has already been set, then
3259  // actually try to set
3260  // active_fe_indices for child
3261  // cells of refined cells to the
3262  // active_fe_index of the mother
3263  // cell. if no finite element
3264  // collection has been assigned
3265  // yet, then all indicators are
3266  // zero anyway, and there is no
3267  // point trying to set anything
3268  // (besides, we would trip over
3269  // an assertion in
3270  // set_active_fe_index)
3271  if (finite_elements != 0)
3272  {
3273  cell_iterator cell = begin(),
3274  endc = end ();
3275  for (; cell != endc; ++cell)
3276  {
3277  // Look if the cell got children during refinement by
3278  // checking whether it has children now but didn't have
3279  // children before refinement (the has_children array is
3280  // set in pre-refinement action)
3281  //
3282  // Note: Although one level is added to
3283  // the DoFHandler levels, when the
3284  // triangulation got one, for the buffer
3285  // has_children this new level is not
3286  // required, because the cells on the
3287  // finest level never have children. Hence
3288  // cell->has_children () will always return
3289  // false on that level, which would cause
3290  // shortcut evaluation of the following
3291  // expression. Thus an index error in
3292  // has_children should never occur.
3293  if (cell->has_children () &&
3294  !(*has_children [cell->level ()])[cell->index ()])
3295  {
3296  // Set active_fe_index in children to the same value
3297  // as in the parent cell. we can't access the
3298  // active_fe_index in the parent cell any more through
3299  // cell->active_fe_index() since that function is not
3300  // allowed for inactive cells, but we can access this
3301  // information from the DoFLevels directly
3302  for (unsigned int i = 0; i < cell->n_children(); ++i)
3303  cell->child (i)->set_active_fe_index
3304  (levels[cell->level()]->active_fe_index (cell->index()));
3305  }
3306  }
3307  }
3308 
3309  // Free buffer objects
3310  std::vector<std::vector<bool> *>::iterator children_level;
3311  for (children_level = has_children.begin ();
3312  children_level != has_children.end ();
3313  ++children_level)
3314  delete (*children_level);
3315  has_children.clear ();
3316  }
3317 
3318 
3319  template <int dim, int spacedim>
3320  template <int structdim>
3321  types::global_dof_index
3322  DoFHandler<dim,spacedim>::get_dof_index (const unsigned int,
3323  const unsigned int,
3324  const unsigned int,
3325  const unsigned int) const
3326  {
3327  Assert (false, ExcNotImplemented());
3328  return numbers::invalid_dof_index;
3329  }
3330 
3331 
3332  template <int dim, int spacedim>
3333  template <int structdim>
3334  void
3335  DoFHandler<dim,spacedim>::set_dof_index (const unsigned int,
3336  const unsigned int,
3337  const unsigned int,
3338  const unsigned int,
3339  const types::global_dof_index) const
3340  {
3341  Assert (false, ExcNotImplemented());
3342  }
3343 
3344 
3345  template<int dim, int spacedim>
3346  void DoFHandler<dim,spacedim>::clear_space ()
3347  {
3348  for (unsigned int i=0; i<levels.size(); ++i)
3349  delete levels[i];
3350  levels.resize (0);
3351  delete faces;
3352  faces = NULL;
3353 
3354  {
3355  std::vector<types::global_dof_index> tmp;
3356  std::swap (vertex_dofs, tmp);
3357  }
3358 
3359  {
3360  std::vector<types::global_dof_index> tmp;
3361  std::swap (vertex_dofs_offsets, tmp);
3362  }
3363  }
3364 }
3365 
3366 
3367 
3368 /*-------------- Explicit Instantiations -------------------------------*/
3369 #include "dof_handler.inst"
3370 
3371 
3372 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:869
DoFHandler::max_couplings_between_dofs
unsigned int max_couplings_between_dofs() const
Definition: dof_handler.cc:1431
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:170
Triangulation::load_user_flags_line
void load_user_flags_line(std::istream &in)
Definition: tria.cc:10024
hp::DoFHandler::DoFHandler
DoFHandler(const Triangulation< dim, spacedim > &tria)
Definition: dof_handler.cc:1664
Triangulation::clear_user_flags
void clear_user_flags()
Definition: tria.cc:9889
DoFHandler::begin
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:728
internal::hp::DoFLevel::offset_type
unsigned int offset_type
Definition: dof_level.h:109
Triangulation::save_user_flags_quad
void save_user_flags_quad(std::ostream &out) const
Definition: tria.cc:10123
Triangulation::load_user_flags
void load_user_flags(std::istream &in)
Definition: tria.cc:9946
DoFHandler::end
cell_iterator end() const
Definition: dof_handler.cc:756
DoFHandler::clear
virtual void clear()
Definition: dof_handler.cc:1180
FiniteElementData::dofs_per_quad
const unsigned int dofs_per_quad
Definition: fe_base.h:250
FiniteElement::hp_quad_dof_identities
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:923
Triangulation::clear_user_flags_quad
void clear_user_flags_quad()
Definition: tria.cc:9849
DoFHandler::max_dofs_per_face
unsigned int max_dofs_per_face(const DoFHandler< dim, spacedim > &dh)
DoFHandler::levels
std::vector<::internal::DoFHandler::DoFLevel< dim > * > levels
Definition: dof_handler.h:1073
DoFHandler::active_cell_iterator
ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:228
DoFHandler::tria
SmartPointer< const Triangulation< dim, spacedim >, DoFHandler< dim, spacedim > > tria
Definition: dof_handler.h:919
DoFHandler::max_dofs_per_vertex
unsigned int max_dofs_per_vertex(const DoFHandler< dim, spacedim > &dh)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:369
DoFHandler::memory_consumption
virtual std::size_t memory_consumption() const
Definition: dof_handler.cc:1050
internal::hp::DoFHandler::Implementation::max_couplings_between_dofs
static unsigned int max_couplings_between_dofs(const DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:1541
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe() const
FiniteElementData::dofs_per_line
const unsigned int dofs_per_line
Definition: fe_base.h:244
Triangulation::load_user_flags_quad
void load_user_flags_quad(std::istream &in)
Definition: tria.cc:10134
hp::DoFHandler::active_cell_iterator
ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:163
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:740
hp::DoFHandler::get_fe
const hp::FECollection< dim, spacedim > & get_fe() const
Triangulation::load_user_flags_hex
void load_user_flags_hex(std::istream &in)
Definition: tria.cc:10197
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
Triangulation::save_user_flags_line
void save_user_flags_line(std::ostream &out) const
Definition: tria.cc:10013
Triangulation::clear_user_flags_hex
void clear_user_flags_hex()
Definition: tria.cc:9881
types::global_dof_index
unsigned int global_dof_index
Definition: types.h:88
internal::hp::DoFHandler::Implementation::distribute_dofs_on_cell
static types::global_dof_index distribute_dofs_on_cell(const typename ::hp::DoFHandler< 1, spacedim >::active_cell_iterator &cell, types::global_dof_index next_free_dof)
Definition: dof_handler.cc:348
FiniteElementData::dofs_per_hex
const unsigned int dofs_per_hex
Definition: fe_base.h:256
Assert
#define Assert(cond, exc)
Definition: exceptions.h:313
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: dof_handler.cc:835
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
DoFHandler::end_active
active_cell_iterator end_active(const unsigned int level) const
Definition: dof_handler.cc:778
DoFHandler::n_boundary_dofs
types::global_dof_index n_boundary_dofs() const
Definition: dof_handler.cc:975
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
internal::hp::DoFHandler::Implementation::reserve_space_vertices
static void reserve_space_vertices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:232
DoFRenumbering::ExcRenumberingIncomplete
static ::ExceptionBase & ExcRenumberingIncomplete()
DoFHandler::faces
::internal::DoFHandler::DoFFaces< dim > * faces
Definition: dof_handler.h:1082
hp::DoFHandler::cell_iterator
ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:197
FiniteElement::hp_vertex_dof_identities
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:901
internal::hp::DoFHandler::Implementation::reserve_space
static void reserve_space(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:541
Triangulation::save_user_flags_hex
void save_user_flags_hex(std::ostream &out) const
Definition: tria.cc:10186
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
hp
Definition: hp.h:102
DoFHandler::~DoFHandler
virtual ~DoFHandler()
Definition: dof_handler.cc:688
DoFHandler::renumber_dofs
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
Definition: dof_handler.cc:1194
MemoryConsumption::memory_consumption
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:274
FiniteElement::hp_line_dof_identities
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:912
DoFHandler::clear_space
void clear_space()
Definition: dof_handler.cc:1478
DoFHandler::distribute_dofs
virtual void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
Definition: dof_handler.cc:1079
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:858
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
IteratorState::valid
Iterator points to a valid object.
Definition: tria_iterator_base.h:39
FiniteElementData::dofs_per_vertex
const unsigned int dofs_per_vertex
Definition: fe_base.h:238
Table
Definition: table.h:33
types::boundary_id
unsigned char boundary_id
Definition: types.h:110
numbers::internal_face_boundary_id
const types::boundary_id internal_face_boundary_id
Definition: types.h:216
DoFHandler::max_couplings_between_boundary_dofs
unsigned int max_couplings_between_boundary_dofs() const
Definition: dof_handler.cc:1440
Triangulation::clear_user_flags_line
void clear_user_flags_line()
Definition: tria.cc:9817
numbers::invalid_dof_index
const types::global_dof_index invalid_dof_index
Definition: types.h:184
DoFTools::ExcNoFESelected
static ::ExceptionBase & ExcNoFESelected()
DoFHandler::cell_iterators
IteratorRange< cell_iterator > cell_iterators() const
Definition: dof_handler.cc:825
FiniteElement::compare_for_face_domination
virtual FiniteElementDomination::Domination compare_for_face_domination(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:934
DoFHandler::vertex_dofs
std::vector< types::global_dof_index > vertex_dofs
Definition: dof_handler.h:1061
Threads::new_task
Task< RT > new_task(const std_cxx11::function< RT()> &function)
Definition: thread_management.h:3211
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Generated by   doxygen 1.8.14