Reference documentation for deal.II version 8.5.0
dof_handler.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2016 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/dofs/dof_handler.h>
18 #include <deal.II/dofs/dof_handler_policy.h>
19 #include <deal.II/dofs/dof_levels.h>
20 #include <deal.II/dofs/dof_faces.h>
21 #include <deal.II/dofs/dof_accessor.h>
22 #include <deal.II/grid/tria_accessor.h>
23 #include <deal.II/grid/tria_iterator.h>
24 #include <deal.II/grid/tria_levels.h>
25 #include <deal.II/grid/tria.h>
26 #include <deal.II/base/geometry_info.h>
27 #include <deal.II/fe/fe.h>
28 #include <deal.II/distributed/shared_tria.h>
29 #include <deal.II/distributed/tria.h>
30 
31 #include <set>
32 #include <algorithm>
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 
37 //TODO[WB]: do not use a plain pointer for DoFHandler::faces, but rather an
38 //unique_ptr or some such thing. alternatively, why not use the DoFFaces object
39 //right away?
40 
41 template<int dim, int spacedim>
42 const unsigned int DoFHandler<dim,spacedim>::dimension;
43 
44 template<int dim, int spacedim>
46 
47 template <int dim, int spacedim>
49 
50 template <int dim, int spacedim>
52 
53 
54 // reference the invalid_dof_index variable explicitly to work around
55 // a bug in the icc8 compiler
56 namespace internal
57 {
58  template <int dim, int spacedim>
59  const types::global_dof_index *dummy ()
60  {
62  }
63 }
64 
65 
66 
67 namespace internal
68 {
69  template<int dim, int spacedim>
70  std::string policy_to_string(const ::internal::DoFHandler::Policy::PolicyBase<dim,spacedim> &policy)
71  {
72  std::string policy_name;
73  if (dynamic_cast<const typename ::internal::DoFHandler::Policy::Sequential<dim,spacedim>*>(&policy))
74  policy_name = "Policy::Sequential<";
75  else if (dynamic_cast<const typename ::internal::DoFHandler::Policy::ParallelDistributed<dim,spacedim>*>(&policy))
76  policy_name = "Policy::ParallelDistributed<";
77  else if (dynamic_cast<const typename ::internal::DoFHandler::Policy::ParallelShared<dim,spacedim>*>(&policy))
78  policy_name = "Policy::ParallelShared<";
79  else
81  policy_name += Utilities::int_to_string(dim)+
82  ","+Utilities::int_to_string(spacedim)+">";
83  return policy_name;
84  }
85 
86 
87  namespace DoFHandler
88  {
89  // access class
90  // ::DoFHandler instead of
91  // namespace internal::DoFHandler
92  using ::DoFHandler;
93 
94 
100  {
105  template <int spacedim>
106  static
107  unsigned int
109  {
110  return std::min(static_cast<types::global_dof_index>(3*dof_handler.selected_fe->dofs_per_vertex +
111  2*dof_handler.selected_fe->dofs_per_line),
112  dof_handler.n_dofs());
113  }
114 
115 
116 
117  template <int spacedim>
118  static
119  unsigned int
121  {
122 
123  // get these numbers by drawing pictures
124  // and counting...
125  // example:
126  // | | |
127  // --x-----x--x--X--
128  // | | | |
129  // | x--x--x
130  // | | | |
131  // --x--x--*--x--x--
132  // | | | |
133  // x--x--x |
134  // | | | |
135  // --X--x--x-----x--
136  // | | |
137  // x = vertices connected with center vertex *;
138  // = total of 19
139  // (the X vertices are connected with * if
140  // the vertices adjacent to X are hanging
141  // nodes)
142  // count lines -> 28 (don't forget to count
143  // mother and children separately!)
144  types::global_dof_index max_couplings;
145  switch (dof_handler.tria->max_adjacent_cells())
146  {
147  case 4:
148  max_couplings=19*dof_handler.selected_fe->dofs_per_vertex +
149  28*dof_handler.selected_fe->dofs_per_line +
150  8*dof_handler.selected_fe->dofs_per_quad;
151  break;
152  case 5:
153  max_couplings=21*dof_handler.selected_fe->dofs_per_vertex +
154  31*dof_handler.selected_fe->dofs_per_line +
155  9*dof_handler.selected_fe->dofs_per_quad;
156  break;
157  case 6:
158  max_couplings=28*dof_handler.selected_fe->dofs_per_vertex +
159  42*dof_handler.selected_fe->dofs_per_line +
160  12*dof_handler.selected_fe->dofs_per_quad;
161  break;
162  case 7:
163  max_couplings=30*dof_handler.selected_fe->dofs_per_vertex +
164  45*dof_handler.selected_fe->dofs_per_line +
165  13*dof_handler.selected_fe->dofs_per_quad;
166  break;
167  case 8:
168  max_couplings=37*dof_handler.selected_fe->dofs_per_vertex +
169  56*dof_handler.selected_fe->dofs_per_line +
170  16*dof_handler.selected_fe->dofs_per_quad;
171  break;
172 
173  // the following numbers are not based on actual counting but by
174  // extrapolating the number sequences from the previous ones (for
175  // example, for dofs_per_vertex, the sequence above is 19, 21, 28,
176  // 30, 37, and is continued as follows):
177  case 9:
178  max_couplings=39*dof_handler.selected_fe->dofs_per_vertex +
179  59*dof_handler.selected_fe->dofs_per_line +
180  17*dof_handler.selected_fe->dofs_per_quad;
181  break;
182  case 10:
183  max_couplings=46*dof_handler.selected_fe->dofs_per_vertex +
184  70*dof_handler.selected_fe->dofs_per_line +
185  20*dof_handler.selected_fe->dofs_per_quad;
186  break;
187  case 11:
188  max_couplings=48*dof_handler.selected_fe->dofs_per_vertex +
189  73*dof_handler.selected_fe->dofs_per_line +
190  21*dof_handler.selected_fe->dofs_per_quad;
191  break;
192  case 12:
193  max_couplings=55*dof_handler.selected_fe->dofs_per_vertex +
194  84*dof_handler.selected_fe->dofs_per_line +
195  24*dof_handler.selected_fe->dofs_per_quad;
196  break;
197  case 13:
198  max_couplings=57*dof_handler.selected_fe->dofs_per_vertex +
199  87*dof_handler.selected_fe->dofs_per_line +
200  25*dof_handler.selected_fe->dofs_per_quad;
201  break;
202  case 14:
203  max_couplings=63*dof_handler.selected_fe->dofs_per_vertex +
204  98*dof_handler.selected_fe->dofs_per_line +
205  28*dof_handler.selected_fe->dofs_per_quad;
206  break;
207  case 15:
208  max_couplings=65*dof_handler.selected_fe->dofs_per_vertex +
209  103*dof_handler.selected_fe->dofs_per_line +
210  29*dof_handler.selected_fe->dofs_per_quad;
211  break;
212  case 16:
213  max_couplings=72*dof_handler.selected_fe->dofs_per_vertex +
214  114*dof_handler.selected_fe->dofs_per_line +
215  32*dof_handler.selected_fe->dofs_per_quad;
216  break;
217 
218  default:
219  Assert (false, ExcNotImplemented());
220  max_couplings=0;
221  }
222  return std::min(max_couplings,dof_handler.n_dofs());
223  }
224 
225 
226  template <int spacedim>
227  static
228  unsigned int
230  {
231 //TODO:[?] Invent significantly better estimates than the ones in this function
232 
233  // doing the same thing here is a
234  // rather complicated thing, compared
235  // to the 2d case, since it is hard
236  // to draw pictures with several
237  // refined hexahedra :-) so I
238  // presently only give a coarse
239  // estimate for the case that at most
240  // 8 hexes meet at each vertex
241  //
242  // can anyone give better estimate
243  // here?
244  const unsigned int max_adjacent_cells
245  = dof_handler.tria->max_adjacent_cells();
246 
247  types::global_dof_index max_couplings;
248  if (max_adjacent_cells <= 8)
249  max_couplings=7*7*7*dof_handler.selected_fe->dofs_per_vertex +
250  7*6*7*3*dof_handler.selected_fe->dofs_per_line +
251  9*4*7*3*dof_handler.selected_fe->dofs_per_quad +
252  27*dof_handler.selected_fe->dofs_per_hex;
253  else
254  {
255  Assert (false, ExcNotImplemented());
256  max_couplings=0;
257  }
258 
259  return std::min(max_couplings,dof_handler.n_dofs());
260  }
261 
262 
272  template <int spacedim>
273  static
275  {
276  dof_handler.vertex_dofs
277  .resize(dof_handler.tria->n_vertices() *
278  dof_handler.selected_fe->dofs_per_vertex,
280 
281  for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
282  {
283  dof_handler.levels
284  .push_back (new internal::DoFHandler::DoFLevel<1>);
285 
286  dof_handler.levels.back()->dof_object.dofs
287  .resize (dof_handler.tria->n_raw_cells(i) *
288  dof_handler.selected_fe->dofs_per_line,
290 
291  dof_handler.levels.back()->cell_dof_indices_cache
292  .resize (dof_handler.tria->n_raw_cells(i) *
293  dof_handler.selected_fe->dofs_per_cell,
295  }
296  }
297 
298 
299  template <int spacedim>
300  static
301  void reserve_space (DoFHandler<2,spacedim> &dof_handler)
302  {
303  dof_handler.vertex_dofs
304  .resize(dof_handler.tria->n_vertices() *
305  dof_handler.selected_fe->dofs_per_vertex,
307 
308  for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
309  {
310  dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<2>);
311 
312  dof_handler.levels.back()->dof_object.dofs
313  .resize (dof_handler.tria->n_raw_cells(i) *
314  dof_handler.selected_fe->dofs_per_quad,
316 
317  dof_handler.levels.back()->cell_dof_indices_cache
318  .resize (dof_handler.tria->n_raw_cells(i) *
319  dof_handler.selected_fe->dofs_per_cell,
321  }
322 
323  dof_handler.faces = new internal::DoFHandler::DoFFaces<2>;
324  // avoid access to n_raw_lines when there are no cells
325  if (dof_handler.tria->n_cells() > 0)
326  {
327  dof_handler.faces->lines.dofs
328  .resize (dof_handler.tria->n_raw_lines() *
329  dof_handler.selected_fe->dofs_per_line,
331  }
332  }
333 
334 
335  template <int spacedim>
336  static
337  void reserve_space (DoFHandler<3,spacedim> &dof_handler)
338  {
339  dof_handler.vertex_dofs
340  .resize(dof_handler.tria->n_vertices() *
341  dof_handler.selected_fe->dofs_per_vertex,
343 
344  for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
345  {
346  dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<3>);
347 
348  dof_handler.levels.back()->dof_object.dofs
349  .resize (dof_handler.tria->n_raw_cells(i) *
350  dof_handler.selected_fe->dofs_per_hex,
352 
353  dof_handler.levels.back()->cell_dof_indices_cache
354  .resize (dof_handler.tria->n_raw_cells(i) *
355  dof_handler.selected_fe->dofs_per_cell,
357  }
358  dof_handler.faces = new internal::DoFHandler::DoFFaces<3>;
359 
360  // avoid access to n_raw_lines when there are no cells
361  if (dof_handler.tria->n_cells() > 0)
362  {
363  dof_handler.faces->lines.dofs
364  .resize (dof_handler.tria->n_raw_lines() *
365  dof_handler.selected_fe->dofs_per_line,
367  dof_handler.faces->quads.dofs
368  .resize (dof_handler.tria->n_raw_quads() *
369  dof_handler.selected_fe->dofs_per_quad,
371  }
372  }
373 
374  template<int spacedim>
375  static
376  void reserve_space_mg (DoFHandler<1, spacedim> &dof_handler)
377  {
378  Assert (dof_handler.get_triangulation().n_levels () > 0, ExcMessage ("Invalid triangulation"));
379  dof_handler.clear_mg_space ();
380 
381  const ::Triangulation<1, spacedim> &tria = dof_handler.get_triangulation();
382  const unsigned int &dofs_per_line = dof_handler.get_fe ().dofs_per_line;
383  const unsigned int &n_levels = tria.n_levels ();
384 
385  for (unsigned int i = 0; i < n_levels; ++i)
386  {
387  dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<1>);
388  dof_handler.mg_levels.back ()->dof_object.dofs = std::vector<types::global_dof_index> (tria.n_raw_lines (i) * dofs_per_line, DoFHandler<1>::invalid_dof_index);
389  }
390 
391  const unsigned int &n_vertices = tria.n_vertices ();
392 
393  dof_handler.mg_vertex_dofs.resize (n_vertices);
394 
395  std::vector<unsigned int> max_level (n_vertices, 0);
396  std::vector<unsigned int> min_level (n_vertices, n_levels);
397 
398  for (typename ::Triangulation<1, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell)
399  {
400  const unsigned int level = cell->level ();
401 
402  for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex)
403  {
404  const unsigned int vertex_index = cell->vertex_index (vertex);
405 
406  if (min_level[vertex_index] > level)
407  min_level[vertex_index] = level;
408 
409  if (max_level[vertex_index] < level)
410  max_level[vertex_index] = level;
411  }
412  }
413 
414  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
415  if (tria.vertex_used (vertex))
416  {
417  Assert (min_level[vertex] < n_levels, ExcInternalError ());
418  Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ());
419  dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], dof_handler.get_fe ().dofs_per_vertex);
420  }
421 
422  else
423  {
424  Assert (min_level[vertex] == n_levels, ExcInternalError ());
425  Assert (max_level[vertex] == 0, ExcInternalError ());
426  dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0);
427  }
428  }
429 
430  template<int spacedim>
431  static
432  void reserve_space_mg (DoFHandler<2, spacedim> &dof_handler)
433  {
434  Assert (dof_handler.get_triangulation().n_levels () > 0, ExcMessage ("Invalid triangulation"));
435  dof_handler.clear_mg_space ();
436 
437  const ::FiniteElement<2, spacedim> &fe = dof_handler.get_fe ();
438  const ::Triangulation<2, spacedim> &tria = dof_handler.get_triangulation();
439  const unsigned int &n_levels = tria.n_levels ();
440 
441  for (unsigned int i = 0; i < n_levels; ++i)
442  {
443  dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<2>);
444  dof_handler.mg_levels.back ()->dof_object.dofs = std::vector<types::global_dof_index> (tria.n_raw_quads (i) * fe.dofs_per_quad, DoFHandler<2>::invalid_dof_index);
445  }
446 
447  dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<2>;
448  dof_handler.mg_faces->lines.dofs = std::vector<types::global_dof_index> (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<2>::invalid_dof_index);
449 
450  const unsigned int &n_vertices = tria.n_vertices ();
451 
452  dof_handler.mg_vertex_dofs.resize (n_vertices);
453 
454  std::vector<unsigned int> max_level (n_vertices, 0);
455  std::vector<unsigned int> min_level (n_vertices, n_levels);
456 
457  for (typename ::Triangulation<2, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell)
458  {
459  const unsigned int level = cell->level ();
460 
461  for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex)
462  {
463  const unsigned int vertex_index = cell->vertex_index (vertex);
464 
465  if (min_level[vertex_index] > level)
466  min_level[vertex_index] = level;
467 
468  if (max_level[vertex_index] < level)
469  max_level[vertex_index] = level;
470  }
471  }
472 
473  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
474  if (tria.vertex_used (vertex))
475  {
476  Assert (min_level[vertex] < n_levels, ExcInternalError ());
477  Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ());
478  dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex);
479  }
480 
481  else
482  {
483  Assert (min_level[vertex] == n_levels, ExcInternalError ());
484  Assert (max_level[vertex] == 0, ExcInternalError ());
485  dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0);
486  }
487  }
488 
489  template<int spacedim>
490  static
491  void reserve_space_mg (DoFHandler<3, spacedim> &dof_handler)
492  {
493  Assert (dof_handler.get_triangulation().n_levels () > 0, ExcMessage ("Invalid triangulation"));
494  dof_handler.clear_mg_space ();
495 
496  const ::FiniteElement<3, spacedim> &fe = dof_handler.get_fe ();
497  const ::Triangulation<3, spacedim> &tria = dof_handler.get_triangulation();
498  const unsigned int &n_levels = tria.n_levels ();
499 
500  for (unsigned int i = 0; i < n_levels; ++i)
501  {
502  dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<3>);
503  dof_handler.mg_levels.back ()->dof_object.dofs = std::vector<types::global_dof_index> (tria.n_raw_hexs (i) * fe.dofs_per_hex, DoFHandler<3>::invalid_dof_index);
504  }
505 
506  dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<3>;
507  dof_handler.mg_faces->lines.dofs = std::vector<types::global_dof_index> (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<3>::invalid_dof_index);
508  dof_handler.mg_faces->quads.dofs = std::vector<types::global_dof_index> (tria.n_raw_quads () * fe.dofs_per_quad, DoFHandler<3>::invalid_dof_index);
509 
510  const unsigned int &n_vertices = tria.n_vertices ();
511 
512  dof_handler.mg_vertex_dofs.resize (n_vertices);
513 
514  std::vector<unsigned int> max_level (n_vertices, 0);
515  std::vector<unsigned int> min_level (n_vertices, n_levels);
516 
517  for (typename ::Triangulation<3, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell)
518  {
519  const unsigned int level = cell->level ();
520 
521  for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex)
522  {
523  const unsigned int vertex_index = cell->vertex_index (vertex);
524 
525  if (min_level[vertex_index] > level)
526  min_level[vertex_index] = level;
527 
528  if (max_level[vertex_index] < level)
529  max_level[vertex_index] = level;
530  }
531  }
532 
533  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
534  if (tria.vertex_used (vertex))
535  {
536  Assert (min_level[vertex] < n_levels, ExcInternalError ());
537  Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ());
538  dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex);
539  }
540 
541  else
542  {
543  Assert (min_level[vertex] == n_levels, ExcInternalError ());
544  Assert (max_level[vertex] == 0, ExcInternalError ());
545  dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0);
546  }
547  }
548 
549 
550 
551  template<int spacedim>
552  static
554  get_dof_index (
555  const DoFHandler<1, spacedim> &dof_handler,
558  const unsigned int obj_index,
559  const unsigned int fe_index,
560  const unsigned int local_index,
561  const int2type<1>)
562  {
563  return mg_level.dof_object.get_dof_index (dof_handler, obj_index, fe_index, local_index);
564  }
565 
566  template<int spacedim>
567  static
569  get_dof_index (const DoFHandler<2, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<2> &, internal::DoFHandler::DoFFaces<2> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>)
570  {
571  return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index);
572  }
573 
574  template<int spacedim>
575  static
577  get_dof_index (const DoFHandler<2, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<2> &mg_level, internal::DoFHandler::DoFFaces<2> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>)
578  {
579  return mg_level.dof_object.get_dof_index (dof_handler, obj_index, fe_index, local_index);
580  }
581 
582  template<int spacedim>
583  static
585  get_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &, internal::DoFHandler::DoFFaces<3> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>)
586  {
587  return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index);
588  }
589 
590  template<int spacedim>
591  static
593  get_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &, internal::DoFHandler::DoFFaces<3> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>)
594  {
595  return mg_faces.quads.get_dof_index (dof_handler, obj_index, fe_index, local_index);
596  }
597 
598  template<int spacedim>
599  static
601  get_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &mg_level, internal::DoFHandler::DoFFaces<3> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<3>)
602  {
603  return mg_level.dof_object.get_dof_index (dof_handler, obj_index, fe_index, local_index);
604  }
605 
606  template<int spacedim>
607  static
608  void set_dof_index (const DoFHandler<1, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<1> &mg_level, internal::DoFHandler::DoFFaces<1> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>)
609  {
610  mg_level.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
611  }
612 
613  template<int spacedim>
614  static
615  void set_dof_index (const DoFHandler<2, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<2> &, internal::DoFHandler::DoFFaces<2> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>)
616  {
617  mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
618  }
619 
620  template<int spacedim>
621  static
622  void set_dof_index (const DoFHandler<2, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<2> &mg_level, internal::DoFHandler::DoFFaces<2> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>)
623  {
624  mg_level.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
625  }
626 
627  template<int spacedim>
628  static
629  void set_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &, internal::DoFHandler::DoFFaces<3> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>)
630  {
631  mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
632  }
633 
634  template<int spacedim>
635  static
636  void set_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &, internal::DoFHandler::DoFFaces<3> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>)
637  {
638  mg_faces.quads.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
639  }
640 
641  template<int spacedim>
642  static
643  void set_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &mg_level, internal::DoFHandler::DoFFaces<3> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<3>)
644  {
645  mg_level.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
646  }
647  };
648  }
649 }
650 
651 
652 
653 template<int dim, int spacedim>
655  :
656  tria(&tria, typeid(*this).name()),
657  selected_fe(0, typeid(*this).name()),
658  faces(NULL),
659  mg_faces (NULL)
660 {
661  // decide whether we need a
662  // sequential or a parallel
663  // distributed policy
664  if (dynamic_cast<const parallel::shared::Triangulation< dim, spacedim>*>
665  (&tria)
666  != 0)
668  else if (dynamic_cast<const parallel::distributed::Triangulation< dim, spacedim >*>
669  (&tria)
670  == 0)
672  else
674 }
675 
676 
677 template<int dim, int spacedim>
679  :
680  tria(0, typeid(*this).name()),
681  selected_fe(0, typeid(*this).name()),
682  faces(NULL),
683  mg_faces (NULL)
684 {}
685 
686 
687 template <int dim, int spacedim>
689 {
690  // release allocated memory
691  clear ();
692 }
693 
694 
695 template<int dim, int spacedim>
696 void
699  const FiniteElement<dim,spacedim> &fe)
700 {
701  tria = &t;
702  faces = 0;
703  number_cache.n_global_dofs = 0;
704 
705  // decide whether we need a
706  // sequential or a parallel
707  // distributed policy
708  if (dynamic_cast<const parallel::shared::Triangulation< dim, spacedim>*>
709  (&t)
710  != 0)
712  else if (dynamic_cast<const parallel::distributed::Triangulation< dim, spacedim >*>
713  (&t)
714  == 0)
716  else
718 
719  distribute_dofs(fe);
720 }
721 
722 
723 
724 /*------------------------ Cell iterator functions ------------------------*/
725 
726 template <int dim, int spacedim>
728 DoFHandler<dim,spacedim>::begin (const unsigned int level) const
729 {
730  typename Triangulation<dim,spacedim>::cell_iterator cell = this->get_triangulation().begin(level);
731  if (cell == this->get_triangulation().end(level))
732  return end(level);
733  return cell_iterator (*cell, this);
734 }
735 
736 
737 
738 template <int dim, int spacedim>
740 DoFHandler<dim,spacedim>::begin_active (const unsigned int level) const
741 {
742  // level is checked in begin
743  cell_iterator i = begin (level);
744  if (i.state() != IteratorState::valid)
745  return i;
746  while (i->has_children())
747  if ((++i).state() != IteratorState::valid)
748  return i;
749  return i;
750 }
751 
752 
753 
754 template <int dim, int spacedim>
757 {
758  return cell_iterator (&this->get_triangulation(),
759  -1,
760  -1,
761  this);
762 }
763 
764 
765 template <int dim, int spacedim>
767 DoFHandler<dim,spacedim>::end (const unsigned int level) const
768 {
769  typename Triangulation<dim,spacedim>::cell_iterator cell = this->get_triangulation().end(level);
770  if (cell.state() != IteratorState::valid)
771  return end();
772  return cell_iterator (*cell, this);
773 }
774 
775 
776 template <int dim, int spacedim>
778 DoFHandler<dim, spacedim>::end_active (const unsigned int level) const
779 {
780  typename Triangulation<dim,spacedim>::cell_iterator cell = this->get_triangulation().end_active(level);
781  if (cell.state() != IteratorState::valid)
782  return active_cell_iterator(end());
783  return active_cell_iterator (*cell, this);
784 }
785 
786 
787 
788 template <int dim, int spacedim>
789 typename DoFHandler<dim, spacedim>::level_cell_iterator
790 DoFHandler<dim, spacedim>::begin_mg (const unsigned int level) const
791 {
792  // Assert(this->has_level_dofs(), ExcMessage("You can only iterate over mg "
793  // "levels if mg dofs got distributed."));
794  typename Triangulation<dim,spacedim>::cell_iterator cell = this->get_triangulation().begin(level);
795  if (cell == this->get_triangulation().end(level))
796  return end_mg(level);
797  return level_cell_iterator (*cell, this);
798 }
799 
800 
801 template <int dim, int spacedim>
802 typename DoFHandler<dim, spacedim>::level_cell_iterator
803 DoFHandler<dim, spacedim>::end_mg (const unsigned int level) const
804 {
805  // Assert(this->has_level_dofs(), ExcMessage("You can only iterate over mg "
806  // "levels if mg dofs got distributed."));
807  typename Triangulation<dim,spacedim>::cell_iterator cell = this->get_triangulation().end(level);
808  if (cell.state() != IteratorState::valid)
809  return end();
810  return level_cell_iterator (*cell, this);
811 }
812 
813 
814 template <int dim, int spacedim>
815 typename DoFHandler<dim, spacedim>::level_cell_iterator
817 {
818  return level_cell_iterator (&this->get_triangulation(), -1, -1, this);
819 }
820 
821 
822 
823 template <int dim, int spacedim>
826 {
827  return
829  (begin(), end());
830 }
831 
832 
833 template <int dim, int spacedim>
836 {
837  return
839  (begin_active(), end());
840 }
841 
842 
843 
844 template <int dim, int spacedim>
847 {
848  return
850  (begin_mg(), end_mg());
851 }
852 
853 
854 
855 
856 template <int dim, int spacedim>
859 {
860  return
862  (begin(level), end(level));
863 }
864 
865 
866 
867 template <int dim, int spacedim>
870 {
871  return
873  (begin_active(level), end_active(level));
874 }
875 
876 
877 
878 template <int dim, int spacedim>
881 {
882  return
884  (begin_mg(level), end_mg(level));
885 }
886 
887 
888 
889 //---------------------------------------------------------------------------
890 
891 
892 
893 template <>
895 {
896  return 2*get_fe().dofs_per_vertex;
897 }
898 
899 
900 
901 template <>
902 template <typename number>
904 {
905  // check that only boundary
906  // indicators 0 and 1 are allowed
907  // in 1d
908  for (typename std::map<types::boundary_id, const Function<1,number>*>::const_iterator i=boundary_ids.begin();
909  i!=boundary_ids.end(); ++i)
910  Assert ((i->first == 0) || (i->first == 1),
912 
913  return boundary_ids.size()*get_fe().dofs_per_vertex;
914 }
915 
916 
917 
918 template <>
919 types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_ids) const
920 {
921  // check that only boundary
922  // indicators 0 and 1 are allowed
923  // in 1d
924  for (std::set<types::boundary_id>::const_iterator i=boundary_ids.begin();
925  i!=boundary_ids.end(); ++i)
926  Assert ((*i == 0) || (*i == 1),
928 
929  return boundary_ids.size()*get_fe().dofs_per_vertex;
930 }
931 
932 
933 template <>
935 {
936  return 2*get_fe().dofs_per_vertex;
937 }
938 
939 
940 
941 template <>
942 template <typename number>
944 {
945  // check that only boundary
946  // indicators 0 and 1 are allowed
947  // in 1d
948  for (typename std::map<types::boundary_id, const Function<2,number>*>::const_iterator i=boundary_ids.begin();
949  i!=boundary_ids.end(); ++i)
950  Assert ((i->first == 0) || (i->first == 1),
952 
953  return boundary_ids.size()*get_fe().dofs_per_vertex;
954 }
955 
956 
957 
958 template <>
959 types::global_dof_index DoFHandler<1,2>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_ids) const
960 {
961  // check that only boundary
962  // indicators 0 and 1 are allowed
963  // in 1d
964  for (std::set<types::boundary_id>::const_iterator i=boundary_ids.begin();
965  i!=boundary_ids.end(); ++i)
966  Assert ((*i == 0) || (*i == 1),
968 
969  return boundary_ids.size()*get_fe().dofs_per_vertex;
970 }
971 
972 
973 
974 template<int dim, int spacedim>
976 {
977  std::set<int> boundary_dofs;
978 
979  const unsigned int dofs_per_face = get_fe().dofs_per_face;
980  std::vector<types::global_dof_index> dofs_on_face(dofs_per_face);
981 
982  // loop over all faces of all cells
983  // and see whether they are at a
984  // boundary. note (i) that we visit
985  // interior faces twice (which we
986  // don't care about) but exterior
987  // faces only once as is
988  // appropriate, and (ii) that we
989  // need not take special care of
990  // single lines (using
991  // @p{cell->has_boundary_lines}),
992  // since we do not support
993  // boundaries of dimension dim-2,
994  // and so every boundary line is
995  // also part of a boundary face.
996  active_cell_iterator cell = begin_active (),
997  endc = end();
998  for (; cell!=endc; ++cell)
999  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
1000  if (cell->at_boundary(f))
1001  {
1002  cell->face(f)->get_dof_indices (dofs_on_face);
1003  for (unsigned int i=0; i<dofs_per_face; ++i)
1004  boundary_dofs.insert(dofs_on_face[i]);
1005  }
1006 
1007  return boundary_dofs.size();
1008 }
1009 
1010 
1011 
1012 template <int dim, int spacedim>
1014 DoFHandler<dim,spacedim>::n_boundary_dofs (const std::set<types::boundary_id> &boundary_ids) const
1015 {
1016  Assert (boundary_ids.find (numbers::internal_face_boundary_id) == boundary_ids.end(),
1018 
1019  std::set<types::global_dof_index> boundary_dofs;
1020 
1021  const unsigned int dofs_per_face = get_fe().dofs_per_face;
1022  std::vector<types::global_dof_index> dofs_on_face(dofs_per_face);
1023 
1024  // same as in the previous
1025  // function, but with a different
1026  // check for the boundary indicator
1027  active_cell_iterator cell = begin_active (),
1028  endc = end();
1029  for (; cell!=endc; ++cell)
1030  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
1031  if (cell->at_boundary(f)
1032  &&
1033  (std::find (boundary_ids.begin(),
1034  boundary_ids.end(),
1035  cell->face(f)->boundary_id()) !=
1036  boundary_ids.end()))
1037  {
1038  cell->face(f)->get_dof_indices (dofs_on_face);
1039  for (unsigned int i=0; i<dofs_per_face; ++i)
1040  boundary_dofs.insert(dofs_on_face[i]);
1041  }
1042 
1043  return boundary_dofs.size();
1044 }
1045 
1046 
1047 
1048 template<int dim, int spacedim>
1049 std::size_t
1051 {
1052  std::size_t mem = (MemoryConsumption::memory_consumption (tria) +
1054  MemoryConsumption::memory_consumption (block_info_object) +
1058  sizeof (number_cache) +
1059  MemoryConsumption::memory_consumption (mg_number_cache) +
1061  for (unsigned int i=0; i<levels.size(); ++i)
1062  mem += MemoryConsumption::memory_consumption (*levels[i]);
1063 
1064  for (unsigned int level = 0; level < mg_levels.size (); ++level)
1065  mem += mg_levels[level]->memory_consumption ();
1066 
1067  if (mg_faces != 0)
1068  mem += MemoryConsumption::memory_consumption (*mg_faces);
1069 
1070  for (unsigned int i = 0; i < mg_vertex_dofs.size (); ++i)
1071  mem += sizeof (MGVertexDoFs) + (1 + mg_vertex_dofs[i].get_finest_level () - mg_vertex_dofs[i].get_coarsest_level ()) * sizeof (types::global_dof_index);
1072 
1073  return mem;
1074 }
1075 
1076 
1077 
1078 template<int dim, int spacedim>
1080 {
1081  selected_fe = &ff;
1082 
1083  // delete all levels and set them
1084  // up newly. note that we still
1085  // have to allocate space for all
1086  // degrees of freedom on this mesh
1087  // (including ghost and cells that
1088  // are entirely stored on different
1089  // processors), though we may not
1090  // assign numbers to some of them
1091  // (i.e. they will remain at
1092  // invalid_dof_index). We need to
1093  // allocate the space because we
1094  // will want to be able to query
1095  // the dof_indices on each cell,
1096  // and simply be told that we don't
1097  // know them on some cell (i.e. get
1098  // back invalid_dof_index)
1099  clear_space ();
1101 
1102  // hand things off to the policy
1103  policy->distribute_dofs (*this,number_cache);
1104 
1105  // initialize the block info object
1106  // only if this is a sequential
1107  // triangulation. it doesn't work
1108  // correctly yet if it is parallel
1109  if (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&*tria) == 0)
1110  block_info_object.initialize(*this, false, true);
1111 }
1112 
1113 
1114 template<int dim, int spacedim>
1116 {
1117  (void)fe;
1118  Assert(levels.size()>0, ExcMessage("Distribute active DoFs using distribute_dofs() before calling distribute_mg_dofs()."));
1119 
1120  const FiniteElement<dim, spacedim> *old_fe = selected_fe;
1121  (void)old_fe;
1122  Assert(old_fe == &fe, ExcMessage("You are required to use the same FE for level and active DoFs!") );
1123 
1126  ExcMessage("The mesh smoothing requirement 'limit_level_difference_at_vertices' has to be set for using multigrid!"));
1127 
1128  clear_mg_space();
1129 
1130  internal::DoFHandler::Implementation::reserve_space_mg (*this);
1132  if (!dist_tr)
1133  mg_number_cache.resize((*tria).n_levels());
1134  else
1135  mg_number_cache.resize(dist_tr->n_global_levels());
1136 
1137  policy->distribute_mg_dofs (*this, mg_number_cache);
1138 
1139  // initialize the block info object
1140  // only if this is a sequential
1141  // triangulation. it doesn't work
1142  // correctly yet if it is parallel
1143  if (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&*tria) == 0)
1144  block_info_object.initialize (*this, true, false);
1145 }
1146 
1147 template<int dim, int spacedim>
1149 {
1150  //TODO: move this to distribute_mg_dofs and remove function
1151 }
1152 
1153 template<int dim, int spacedim>
1155 {
1156  for (unsigned int i = 0; i < mg_levels.size (); ++i)
1157  delete mg_levels[i];
1158 
1159  mg_levels.clear ();
1160  delete mg_faces;
1161  mg_faces = NULL;
1162 
1163  std::vector<MGVertexDoFs> tmp;
1164 
1165  std::swap (mg_vertex_dofs, tmp);
1166 
1167  mg_number_cache.clear();
1168 }
1169 
1170 
1171 template<int dim, int spacedim>
1173 {
1174  block_info_object.initialize_local(*this);
1175 }
1176 
1177 
1178 
1179 template<int dim, int spacedim>
1181 {
1182  // release lock to old fe
1183  selected_fe = 0;
1184 
1185  // release memory
1186  clear_space ();
1187  clear_mg_space ();
1188 }
1189 
1190 
1191 
1192 template <int dim, int spacedim>
1193 void
1194 DoFHandler<dim,spacedim>::renumber_dofs (const std::vector<types::global_dof_index> &new_numbers)
1195 {
1196  Assert(levels.size()>0, ExcMessage("You need to distribute DoFs before you can renumber them."));
1197 
1198  Assert (new_numbers.size() == n_locally_owned_dofs(),
1200 
1201 #ifdef DEBUG
1202  // assert that the new indices are
1203  // consecutively numbered if we are
1204  // working on a single
1205  // processor. this doesn't need to
1206  // hold in the case of a parallel
1207  // mesh since we map the interval
1208  // [0...n_dofs()) into itself but
1209  // only globally, not on each
1210  // processor
1211  if (n_locally_owned_dofs() == n_dofs())
1212  {
1213  std::vector<types::global_dof_index> tmp(new_numbers);
1214  std::sort (tmp.begin(), tmp.end());
1215  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
1217  for (; p!=tmp.end(); ++p, ++i)
1218  Assert (*p == i, ExcNewNumbersNotConsecutive(i));
1219  }
1220  else
1221  for (types::global_dof_index i=0; i<new_numbers.size(); ++i)
1222  Assert (new_numbers[i] < n_dofs(),
1223  ExcMessage ("New DoF index is not less than the total number of dofs."));
1224 #endif
1225 
1226  policy->renumber_dofs (new_numbers, *this,number_cache);
1227 }
1228 
1229 
1230 template <int dim, int spacedim>
1231 void
1233  const std::vector<types::global_dof_index> &)
1234 {
1235  Assert(false, ExcNotImplemented());
1236 }
1237 
1238 
1239 template<>
1240 void DoFHandler<1>::renumber_dofs (const unsigned int level,
1241  const std::vector<types::global_dof_index> &new_numbers)
1242 {
1243  Assert(mg_levels.size()>0 && levels.size()>0,
1244  ExcMessage("You need to distribute active and level DoFs before you can renumber level DoFs."));
1245  Assert (new_numbers.size() == n_dofs(level), DoFHandler<1>::ExcRenumberingIncomplete());
1246 
1247  // note that we can not use cell iterators
1248  // in this function since then we would
1249  // renumber the dofs on the interface of
1250  // two cells more than once. Anyway, this
1251  // ways it's not only more correct but also
1252  // faster
1253  for (std::vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
1254  i!=mg_vertex_dofs.end(); ++i)
1255  // if the present vertex lives on
1256  // the present level
1257  if ((i->get_coarsest_level() <= level) &&
1258  (i->get_finest_level() >= level))
1259  for (unsigned int d=0; d<this->get_fe().dofs_per_vertex; ++d)
1260  i->set_index (level, d, new_numbers[i->get_index (level, d)]);
1261 
1262  for (std::vector<types::global_dof_index>::iterator i=mg_levels[level]->dof_object.dofs.begin();
1263  i!=mg_levels[level]->dof_object.dofs.end(); ++i)
1264  {
1266  {
1267  Assert(*i<new_numbers.size(), ExcInternalError());
1268  *i = new_numbers[*i];
1269  }
1270  }
1271 }
1272 
1273 
1274 
1275 template<>
1276 void DoFHandler<2>::renumber_dofs (const unsigned int level,
1277  const std::vector<types::global_dof_index> &new_numbers)
1278 {
1279  Assert(mg_levels.size()>0 && levels.size()>0,
1280  ExcMessage("You need to distribute active and level DoFs before you can renumber level DoFs."));
1281  Assert (new_numbers.size() == n_dofs(level),
1283 
1284  for (std::vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
1285  i!=mg_vertex_dofs.end(); ++i)
1286  // if the present vertex lives on
1287  // the present level
1288  if ((i->get_coarsest_level() <= level) &&
1289  (i->get_finest_level() >= level))
1290  for (unsigned int d=0; d<this->get_fe().dofs_per_vertex; ++d)
1291  i->set_index (level, d, new_numbers[i->get_index (level, d)]);
1292 
1293  if (this->get_fe().dofs_per_line > 0)
1294  {
1295  // save user flags as they will be modified
1296  std::vector<bool> user_flags;
1297  this->get_triangulation().save_user_flags(user_flags);
1298  const_cast<Triangulation<2> &>(this->get_triangulation()).clear_user_flags ();
1299 
1300  // flag all lines adjacent to cells of the current
1301  // level, as those lines logically belong to the same
1302  // level as the cell, at least for for isotropic
1303  // refinement
1304  level_cell_iterator cell, endc = end(level);
1305  for (cell = begin(level); cell != endc; ++cell)
1306  for (unsigned int line=0; line < GeometryInfo<2>::faces_per_cell; ++line)
1307  cell->face(line)->set_user_flag();
1308 
1309  for (cell_iterator cell = begin(); cell != end(); ++cell)
1310  for (unsigned int l=0; l<GeometryInfo<2>::lines_per_cell; ++l)
1311  if (cell->line(l)->user_flag_set())
1312  {
1313  for (unsigned int d=0; d<this->get_fe().dofs_per_line; ++d)
1314  cell->line(l)->set_mg_dof_index (level, d,
1315  new_numbers[cell->line(l)->mg_dof_index(level, d)]);
1316  cell->line(l)->clear_user_flag();
1317  }
1318  // finally, restore user flags
1319  const_cast<Triangulation<2> &>(this->get_triangulation()).load_user_flags (user_flags);
1320  }
1321 
1322  for (std::vector<types::global_dof_index>::iterator i=mg_levels[level]->dof_object.dofs.begin();
1323  i!=mg_levels[level]->dof_object.dofs.end(); ++i)
1324  {
1326  {
1327  Assert(*i<new_numbers.size(), ExcInternalError());
1328  *i = new_numbers[*i];
1329  }
1330  }
1331 }
1332 
1333 
1334 
1335 template<>
1336 void DoFHandler<3>::renumber_dofs (const unsigned int level,
1337  const std::vector<types::global_dof_index> &new_numbers)
1338 {
1339  Assert(mg_levels.size()>0 && levels.size()>0,
1340  ExcMessage("You need to distribute active and level DoFs before you can renumber level DoFs."));
1341  Assert (new_numbers.size() == n_dofs(level),
1343 
1344  for (std::vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
1345  i!=mg_vertex_dofs.end(); ++i)
1346  // if the present vertex lives on
1347  // the present level
1348  if ((i->get_coarsest_level() <= level) &&
1349  (i->get_finest_level() >= level))
1350  for (unsigned int d=0; d<this->get_fe().dofs_per_vertex; ++d)
1351  i->set_index (level, d, new_numbers[i->get_index (level, d)]);
1352 
1353  // LINE DoFs
1354  if (this->get_fe().dofs_per_line > 0)
1355  {
1356  // save user flags as they will be modified
1357  std::vector<bool> user_flags;
1358  this->get_triangulation().save_user_flags(user_flags);
1359  const_cast<Triangulation<3> &>(this->get_triangulation()).clear_user_flags ();
1360 
1361  // flag all lines adjacent to cells of the current
1362  // level, as those lines logically belong to the same
1363  // level as the cell, at least for for isotropic
1364  // refinement
1365  level_cell_iterator cell, endc = end(level);
1366  for (cell = begin(level) ; cell != endc ; ++cell)
1367  for (unsigned int line=0; line < GeometryInfo<3>::lines_per_cell; ++line)
1368  cell->line(line)->set_user_flag();
1369 
1370 
1371  for (cell = begin(level); cell != endc; ++cell)
1372  for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
1373  if (cell->line(l)->user_flag_set())
1374  {
1375  for (unsigned int d=0; d<this->get_fe().dofs_per_line; ++d)
1376  cell->line(l)->set_mg_dof_index (level, d,
1377  new_numbers[cell->line(l)->mg_dof_index(level, d)]);
1378  cell->line(l)->clear_user_flag();
1379  }
1380  // finally, restore user flags
1381  const_cast<Triangulation<3> &>(this->get_triangulation()).load_user_flags (user_flags);
1382  }
1383 
1384  // QUAD DoFs
1385  if (this->get_fe().dofs_per_quad > 0)
1386  {
1387  // save user flags as they will be modified
1388  std::vector<bool> user_flags;
1389  this->get_triangulation().save_user_flags(user_flags);
1390  const_cast<Triangulation<3> &>(this->get_triangulation()).clear_user_flags ();
1391 
1392  // flag all quads adjacent to cells of the current
1393  // level, as those lines logically belong to the same
1394  // level as the cell, at least for for isotropic
1395  // refinement
1396  level_cell_iterator cell, endc = end(level);
1397  for (cell = begin(level) ; cell != endc; ++cell)
1398  for (unsigned int quad=0; quad < GeometryInfo<3>::faces_per_cell; ++quad)
1399  cell->face(quad)->set_user_flag();
1400 
1401  for (cell = begin(level); cell != endc; ++cell)
1402  for (unsigned int q=0; q<GeometryInfo<3>::quads_per_cell; ++q)
1403  if (cell->quad(q)->user_flag_set())
1404  {
1405  for (unsigned int d=0; d<this->get_fe().dofs_per_quad; ++d)
1406  cell->quad(q)->set_mg_dof_index (level, d,
1407  new_numbers[cell->quad(q)->mg_dof_index(level, d)]);
1408  cell->quad(q)->clear_user_flag();
1409  }
1410  // finally, restore user flags
1411  const_cast<Triangulation<3> &>(this->get_triangulation()).load_user_flags (user_flags);
1412  }
1413 
1414  //HEX DoFs
1415  for (std::vector<types::global_dof_index>::iterator i=mg_levels[level]->dof_object.dofs.begin();
1416  i!=mg_levels[level]->dof_object.dofs.end(); ++i)
1417  {
1419  {
1420  Assert(*i<new_numbers.size(), ExcInternalError());
1421  *i = new_numbers[*i];
1422  }
1423  }
1424 }
1425 
1426 
1427 
1428 
1429 template <int dim, int spacedim>
1430 unsigned int
1432 {
1434 }
1435 
1436 
1437 
1438 template <int dim, int spacedim>
1439 unsigned int
1441 {
1442  switch (dim)
1443  {
1444  case 1:
1445  return get_fe().dofs_per_vertex;
1446  case 2:
1447  return (3*get_fe().dofs_per_vertex +
1448  2*get_fe().dofs_per_line);
1449  case 3:
1450  // we need to take refinement of
1451  // one boundary face into
1452  // consideration here; in fact,
1453  // this function returns what
1454  // #max_coupling_between_dofs<2>
1455  // returns
1456  //
1457  // we assume here, that only four
1458  // faces meet at the boundary;
1459  // this assumption is not
1460  // justified and needs to be
1461  // fixed some time. fortunately,
1462  // omitting it for now does no
1463  // harm since the matrix will cry
1464  // foul if its requirements are
1465  // not satisfied
1466  return (19*get_fe().dofs_per_vertex +
1467  28*get_fe().dofs_per_line +
1468  8*get_fe().dofs_per_quad);
1469  default:
1470  Assert (false, ExcNotImplemented());
1472  }
1473 }
1474 
1475 
1476 
1477 template<int dim, int spacedim>
1479 {
1480  for (unsigned int i=0; i<levels.size(); ++i)
1481  delete levels[i];
1482  levels.resize (0);
1483 
1484  delete faces;
1485  faces = 0;
1486 
1487  std::vector<types::global_dof_index> tmp;
1488  std::swap (vertex_dofs, tmp);
1489 
1490  number_cache.clear ();
1491 }
1492 
1493 template<int dim, int spacedim>
1494 template<int structdim>
1497  const unsigned int obj_level,
1498  const unsigned int obj_index,
1499  const unsigned int fe_index,
1500  const unsigned int local_index) const
1501 {
1502  return internal::DoFHandler::Implementation::get_dof_index (*this, *this->mg_levels[obj_level],
1503  *this->mg_faces, obj_index,
1504  fe_index, local_index,
1506 }
1507 
1508 template<int dim, int spacedim>
1509 template<int structdim>
1510 void DoFHandler<dim, spacedim>::set_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const
1511 {
1512  internal::DoFHandler::Implementation::set_dof_index (*this, *this->mg_levels[obj_level], *this->mg_faces, obj_index, fe_index, local_index, global_index, internal::int2type<structdim> ());
1513 }
1514 
1515 
1516 template<int dim, int spacedim>
1517 DoFHandler<dim, spacedim>::MGVertexDoFs::MGVertexDoFs (): coarsest_level (numbers::invalid_unsigned_int), finest_level (0), indices (0), indices_offset (0)
1518 {
1519 }
1520 
1521 
1522 template<int dim, int spacedim>
1524 {
1525  delete[] indices;
1526  delete[] indices_offset;
1527 }
1528 
1529 template<int dim, int spacedim>
1530 void DoFHandler<dim, spacedim>::MGVertexDoFs::init (const unsigned int cl, const unsigned int fl, const unsigned int dofs_per_vertex)
1531 {
1532  if (indices != 0)
1533  {
1534  delete[] indices;
1535  indices = 0;
1536  }
1537 
1538  if (indices_offset != 0)
1539  {
1540  delete[] indices_offset;
1541  indices_offset = 0;
1542  }
1543 
1544  coarsest_level = cl;
1545  finest_level = fl;
1546 
1547  if (cl > fl)
1548  return;
1549 
1550  const unsigned int n_levels = finest_level - coarsest_level + 1;
1551  const unsigned int n_indices = n_levels * dofs_per_vertex;
1552 
1553  indices = new types::global_dof_index[n_indices];
1554  Assert (indices != 0, ExcNoMemory ());
1555 
1556  for (unsigned int i = 0; i < n_indices; ++i)
1558 
1559  indices_offset = new types::global_dof_index[n_levels];
1560  Assert (indices != 0, ExcNoMemory ());
1561 
1562  for (unsigned int i = 0; i < n_levels; ++i)
1563  indices_offset[i] = i * dofs_per_vertex;
1564 }
1565 
1566 template<int dim, int spacedim>
1568 {
1569  return coarsest_level;
1570 }
1571 
1572 template<int dim, int spacedim>
1574 {
1575  return finest_level;
1576 }
1577 
1578 
1579 /*-------------- Explicit Instantiations -------------------------------*/
1580 #include "dof_handler.inst"
1581 
1582 
1583 DEAL_II_NAMESPACE_CLOSE
virtual unsigned int n_global_levels() const
Definition: tria_base.cc:125
unsigned int max_couplings_between_boundary_dofs() const
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1067
active_cell_iterator end_active(const unsigned int level) const
Definition: dof_handler.cc:778
static const unsigned int invalid_unsigned_int
Definition: types.h:170
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:740
Iterator points to a valid object.
static::ExceptionBase & ExcRenumberingIncomplete()
level_cell_iterator begin_mg(const unsigned int level=0) const
Definition: dof_handler.cc:790
virtual void clear()
const Triangulation< dim, spacedim > & get_triangulation() const
SmartPointer< const FiniteElement< dim, spacedim >, DoFHandler< dim, spacedim > > selected_fe
Definition: dof_handler.h:930
std::vector<::internal::DoFHandler::DoFLevel< dim > * > levels
Definition: dof_handler.h:1073
ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:228
SmartPointer< const Triangulation< dim, spacedim >, DoFHandler< dim, spacedim > > tria
Definition: dof_handler.h:919
#define AssertThrow(cond, exc)
Definition: exceptions.h:369
unsigned int max_couplings_between_dofs() const
void initialize(const Triangulation< dim, spacedim > &tria, const FiniteElement< dim, spacedim > &fe)
Definition: dof_handler.cc:697
virtual void distribute_mg_dofs(const FiniteElement< dim, spacedim > &fe)
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
Definition: dof_handler.cc:869
static::ExceptionBase & ExcMessage(std::string arg1)
level_cell_iterator end_mg() const
Definition: dof_handler.cc:816
ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:256
unsigned int global_dof_index
Definition: types.h:88
#define Assert(cond, exc)
Definition: exceptions.h:313
types::global_dof_index n_dofs() const
void set_dof_index(const ::DoFHandler< dh_dim, spacedim > &dof_handler, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index)
Definition: dof_objects.cc:41
unsigned int get_finest_level() const
types::global_dof_index get_dof_index(const ::DoFHandler< dh_dim, spacedim > &dof_handler, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const
Definition: dof_objects.h:189
unsigned int get_coarsest_level() const
static::ExceptionBase & ExcRenumberingIncomplete()
::internal::DoFHandler::DoFFaces< dim > * faces
Definition: dof_handler.h:1082
IteratorState::IteratorStates state() const
internal::DoFHandler::DoFObjects< 1 > lines
Definition: dof_faces.h:116
DoFObjects< dim > dof_object
Definition: dof_levels.h:82
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:85
cell_iterator end() const
Definition: dof_handler.cc:756
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
virtual ~DoFHandler()
Definition: dof_handler.cc:688
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
Definition: dof_handler.cc:858
IteratorRange< cell_iterator > cell_iterators() const
Definition: dof_handler.cc:825
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
void initialize_local_block_info()
static unsigned int max_couplings_between_dofs(const DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:108
internal::DoFHandler::DoFObjects< 1 > lines
Definition: dof_faces.h:146
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:728
IteratorRange< level_cell_iterator > mg_cell_iterators() const
Definition: dof_handler.cc:846
types::global_dof_index n_boundary_dofs() const
Definition: dof_handler.cc:975
void clear_space()
internal::DoFHandler::DoFObjects< 2 > quads
Definition: dof_faces.h:151
IteratorRange< level_cell_iterator > mg_cell_iterators_on_level(const unsigned int level) const
Definition: dof_handler.cc:880
void init(const unsigned int coarsest_level, const unsigned int finest_level, const unsigned int dofs_per_vertex)
virtual void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
static::ExceptionBase & ExcInvalidBoundaryIndicator()
static::ExceptionBase & ExcNotImplemented()
unsigned char boundary_id
Definition: types.h:110
const types::boundary_id internal_face_boundary_id
Definition: types.h:216
virtual std::size_t memory_consumption() const
static void reserve_space(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:274
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: dof_handler.cc:835
std::vector< types::global_dof_index > vertex_dofs
Definition: dof_handler.h:1061
const FiniteElement< dim, spacedim > & get_fe() const
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static::ExceptionBase & ExcInternalError()