Reference documentation for deal.II version 8.5.1
dof_handler_policy.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 
17 //TODO [TH]: renumber DoFs for multigrid is not done yet
18 
19 #include <deal.II/base/geometry_info.h>
20 #include <deal.II/base/utilities.h>
21 #include <deal.II/base/memory_consumption.h>
22 #include <deal.II/grid/tria.h>
23 #include <deal.II/grid/tria_iterator.h>
24 #include <deal.II/dofs/dof_handler.h>
25 #include <deal.II/dofs/dof_accessor.h>
26 #include <deal.II/dofs/dof_handler_policy.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 #include <numeric>
34 
35 DEAL_II_NAMESPACE_OPEN
36 
37 
38 namespace internal
39 {
40  namespace DoFHandler
41  {
42  namespace Policy
43  {
44  // use class ::DoFHandler instead
45  // of namespace internal::DoFHandler in
46  // the following
47  using ::DoFHandler;
48 
49  struct Implementation
50  {
51 
52  /* -------------- distribute_dofs functionality ------------- */
53 
65  template <int spacedim>
66  static
67  types::global_dof_index
68  distribute_dofs_on_cell (const DoFHandler<1,spacedim> &dof_handler,
69  const typename DoFHandler<1,spacedim>::active_cell_iterator &cell,
70  types::global_dof_index next_free_dof)
71  {
72 
73  // distribute dofs of vertices
74  if (dof_handler.get_fe().dofs_per_vertex > 0)
75  for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
76  {
77  if (cell->vertex_dof_index (v,0) ==
78  DoFHandler<1,spacedim>::invalid_dof_index)
79  for (unsigned int d=0;
80  d<dof_handler.get_fe().dofs_per_vertex; ++d)
81  {
82  Assert ((cell->vertex_dof_index (v,d) ==
83  DoFHandler<1,spacedim>::invalid_dof_index),
84  ExcInternalError());
85  cell->set_vertex_dof_index (v, d, next_free_dof++);
86  }
87  else
88  for (unsigned int d=0;
89  d<dof_handler.get_fe().dofs_per_vertex; ++d)
90  Assert ((cell->vertex_dof_index (v,d) !=
91  DoFHandler<1,spacedim>::invalid_dof_index),
92  ExcInternalError());
93  }
94 
95  // dofs of line
96  for (unsigned int d=0;
97  d<dof_handler.get_fe().dofs_per_line; ++d)
98  cell->set_dof_index (d, next_free_dof++);
99 
100  return next_free_dof;
101  }
102 
103 
104 
105  template <int spacedim>
106  static
107  types::global_dof_index
108  distribute_dofs_on_cell (const DoFHandler<2,spacedim> &dof_handler,
109  const typename DoFHandler<2,spacedim>::active_cell_iterator &cell,
110  types::global_dof_index next_free_dof)
111  {
112  if (dof_handler.get_fe().dofs_per_vertex > 0)
113  // number dofs on vertices
114  for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_cell; ++vertex)
115  // check whether dofs for this
116  // vertex have been distributed
117  // (only check the first dof)
118  if (cell->vertex_dof_index(vertex, 0) == DoFHandler<2,spacedim>::invalid_dof_index)
119  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_vertex; ++d)
120  cell->set_vertex_dof_index (vertex, d, next_free_dof++);
121 
122  // for the four sides
123  if (dof_handler.get_fe().dofs_per_line > 0)
124  for (unsigned int side=0; side<GeometryInfo<2>::faces_per_cell; ++side)
125  {
126  const typename DoFHandler<2,spacedim>::line_iterator
127  line = cell->line(side);
128 
129  // distribute dofs if necessary:
130  // check whether line dof is already
131  // numbered (check only first dof)
132  if (line->dof_index(0) == DoFHandler<2,spacedim>::invalid_dof_index)
133  // if not: distribute dofs
134  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_line; ++d)
135  line->set_dof_index (d, next_free_dof++);
136  }
137 
138 
139  // dofs of quad
140  if (dof_handler.get_fe().dofs_per_quad > 0)
141  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_quad; ++d)
142  cell->set_dof_index (d, next_free_dof++);
143 
144  return next_free_dof;
145  }
146 
147 
148  template <int spacedim>
149  static
150  types::global_dof_index
151  distribute_dofs_on_cell (const DoFHandler<3,spacedim> &dof_handler,
152  const typename DoFHandler<3,spacedim>::active_cell_iterator &cell,
153  types::global_dof_index next_free_dof)
154  {
155  if (dof_handler.get_fe().dofs_per_vertex > 0)
156  // number dofs on vertices
157  for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
158  // check whether dofs for this
159  // vertex have been distributed
160  // (only check the first dof)
161  if (cell->vertex_dof_index(vertex, 0) == DoFHandler<3,spacedim>::invalid_dof_index)
162  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_vertex; ++d)
163  cell->set_vertex_dof_index (vertex, d, next_free_dof++);
164 
165  // for the lines
166  if (dof_handler.get_fe().dofs_per_line > 0)
167  for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
168  {
169  const typename DoFHandler<3,spacedim>::line_iterator
170  line = cell->line(l);
171 
172  // distribute dofs if necessary:
173  // check whether line dof is already
174  // numbered (check only first dof)
175  if (line->dof_index(0) == DoFHandler<3,spacedim>::invalid_dof_index)
176  // if not: distribute dofs
177  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_line; ++d)
178  line->set_dof_index (d, next_free_dof++);
179  }
180 
181  // for the quads
182  if (dof_handler.get_fe().dofs_per_quad > 0)
183  for (unsigned int q=0; q<GeometryInfo<3>::quads_per_cell; ++q)
184  {
185  const typename DoFHandler<3,spacedim>::quad_iterator
186  quad = cell->quad(q);
187 
188  // distribute dofs if necessary:
189  // check whether quad dof is already
190  // numbered (check only first dof)
191  if (quad->dof_index(0) == DoFHandler<3,spacedim>::invalid_dof_index)
192  // if not: distribute dofs
193  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_quad; ++d)
194  quad->set_dof_index (d, next_free_dof++);
195  }
196 
197 
198  // dofs of hex
199  if (dof_handler.get_fe().dofs_per_hex > 0)
200  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_hex; ++d)
201  cell->set_dof_index (d, next_free_dof++);
202 
203  return next_free_dof;
204  }
205 
206 
213  template <int dim, int spacedim>
214  static
215  types::global_dof_index
216  distribute_dofs (const types::global_dof_index offset,
217  const types::subdomain_id subdomain_id,
218  DoFHandler<dim,spacedim> &dof_handler)
219  {
220  Assert (dof_handler.get_triangulation().n_levels() > 0,
221  ExcMessage("Empty triangulation"));
222 
223  types::global_dof_index next_free_dof = offset;
224  typename DoFHandler<dim,spacedim>::active_cell_iterator
225  cell = dof_handler.begin_active(),
226  endc = dof_handler.end();
227 
228  for (; cell != endc; ++cell)
229  if ((subdomain_id == numbers::invalid_subdomain_id)
230  ||
231  (cell->subdomain_id() == subdomain_id))
232  next_free_dof
233  = Implementation::distribute_dofs_on_cell (dof_handler,
234  cell,
235  next_free_dof);
236 
237  // update the cache used for cell dof indices
238  for (cell = dof_handler.begin_active(); cell != endc; ++cell)
239  if (!cell->is_artificial())
240  cell->update_cell_dof_indices_cache ();
241 
242  return next_free_dof;
243  }
244 
245 
264  // These three function
265  // have an unused
266  // DoFHandler object as
267  // their first
268  // argument. Without it,
269  // the file was not
270  // compileable under gcc
271  // 4.4.5 (Debian).
272  template <int spacedim>
273  static
274  unsigned int
275  distribute_mg_dofs_on_cell (const DoFHandler<1,spacedim> &,
276  typename DoFHandler<1,spacedim>::level_cell_iterator &cell,
277  unsigned int next_free_dof)
278  {
279  const unsigned int dim = 1;
280 
281  // distribute dofs of vertices
282  if (cell->get_fe().dofs_per_vertex > 0)
283  for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
284  {
285  typename DoFHandler<dim,spacedim>::level_cell_iterator neighbor = cell->neighbor(v);
286 
287  if (neighbor.state() == IteratorState::valid)
288  {
289  // has neighbor already been processed?
290  if (neighbor->user_flag_set() &&
291  (neighbor->level() == cell->level()))
292  // copy dofs if the neighbor is on
293  // the same level (only then are
294  // mg dofs the same)
295  {
296  if (v==0)
297  for (unsigned int d=0; d<cell->get_fe().dofs_per_vertex; ++d)
298  cell->set_mg_vertex_dof_index (cell->level(), 0, d,
299  neighbor->mg_vertex_dof_index (cell->level(), 1, d));
300  else
301  for (unsigned int d=0; d<cell->get_fe().dofs_per_vertex; ++d)
302  cell->set_mg_vertex_dof_index (cell->level(), 1, d,
303  neighbor->mg_vertex_dof_index (cell->level(), 0, d));
304 
305  // next neighbor
306  continue;
307  };
308  };
309 
310  // otherwise: create dofs newly
311  for (unsigned int d=0; d<cell->get_fe().dofs_per_vertex; ++d)
312  cell->set_mg_vertex_dof_index (cell->level(), v, d, next_free_dof++);
313  };
314 
315  // dofs of line
316  if (cell->get_fe().dofs_per_line > 0)
317  for (unsigned int d=0; d<cell->get_fe().dofs_per_line; ++d)
318  cell->set_mg_dof_index (cell->level(), d, next_free_dof++);
319 
320  // note that this cell has been processed
321  cell->set_user_flag ();
322 
323  return next_free_dof;
324  }
325 
326 
327  template <int spacedim>
328  static
329  unsigned int
330  distribute_mg_dofs_on_cell (const DoFHandler<2,spacedim> &,
331  typename DoFHandler<2,spacedim>::level_cell_iterator &cell,
332  unsigned int next_free_dof)
333  {
334  const unsigned int dim = 2;
335  if (cell->get_fe().dofs_per_vertex > 0)
336  // number dofs on vertices
337  for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_cell; ++vertex)
338  // check whether dofs for this
339  // vertex have been distributed
340  // (only check the first dof)
341  if (cell->mg_vertex_dof_index(cell->level(), vertex, 0) == DoFHandler<2>::invalid_dof_index)
342  for (unsigned int d=0; d<cell->get_fe().dofs_per_vertex; ++d)
343  cell->set_mg_vertex_dof_index (cell->level(), vertex, d, next_free_dof++);
344 
345  // for the four sides
346  if (cell->get_fe().dofs_per_line > 0)
347  for (unsigned int side=0; side<GeometryInfo<2>::faces_per_cell; ++side)
348  {
349  typename DoFHandler<dim,spacedim>::line_iterator line = cell->line(side);
350 
351  // distribute dofs if necessary:
352  // check whether line dof is already
353  // numbered (check only first dof)
354  if (line->mg_dof_index(cell->level(), 0) == DoFHandler<2>::invalid_dof_index)
355  // if not: distribute dofs
356  for (unsigned int d=0; d<cell->get_fe().dofs_per_line; ++d)
357  line->set_mg_dof_index (cell->level(), d, next_free_dof++);
358  };
359 
360 
361  // dofs of quad
362  if (cell->get_fe().dofs_per_quad > 0)
363  for (unsigned int d=0; d<cell->get_fe().dofs_per_quad; ++d)
364  cell->set_mg_dof_index (cell->level(), d, next_free_dof++);
365 
366 
367  // note that this cell has been processed
368  cell->set_user_flag ();
369 
370  return next_free_dof;
371  }
372 
373 
374  template <int spacedim>
375  static
376  unsigned int
377  distribute_mg_dofs_on_cell (const DoFHandler<3,spacedim> &,
378  typename DoFHandler<3,spacedim>::level_cell_iterator &cell,
379  unsigned int next_free_dof)
380  {
381  const unsigned int dim = 3;
382  if (cell->get_fe().dofs_per_vertex > 0)
383  // number dofs on vertices
384  for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
385  // check whether dofs for this
386  // vertex have been distributed
387  // (only check the first dof)
388  if (cell->mg_vertex_dof_index(cell->level(), vertex, 0) == DoFHandler<3>::invalid_dof_index)
389  for (unsigned int d=0; d<cell->get_fe().dofs_per_vertex; ++d)
390  cell->set_mg_vertex_dof_index (cell->level(), vertex, d, next_free_dof++);
391 
392  // for the lines
393  if (cell->get_fe().dofs_per_line > 0)
394  for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
395  {
396  typename DoFHandler<dim,spacedim>::line_iterator line = cell->line(l);
397 
398  // distribute dofs if necessary:
399  // check whether line dof is already
400  // numbered (check only first dof)
401  if (line->mg_dof_index(cell->level(), 0) == DoFHandler<3>::invalid_dof_index)
402  // if not: distribute dofs
403  for (unsigned int d=0; d<cell->get_fe().dofs_per_line; ++d)
404  line->set_mg_dof_index (cell->level(), d, next_free_dof++);
405  };
406 
407  // for the quads
408  if (cell->get_fe().dofs_per_quad > 0)
409  for (unsigned int q=0; q<GeometryInfo<3>::quads_per_cell; ++q)
410  {
411  typename DoFHandler<dim,spacedim>::quad_iterator quad = cell->quad(q);
412 
413  // distribute dofs if necessary:
414  // check whether line dof is already
415  // numbered (check only first dof)
416  if (quad->mg_dof_index(cell->level(), 0) == DoFHandler<3>::invalid_dof_index)
417  // if not: distribute dofs
418  for (unsigned int d=0; d<cell->get_fe().dofs_per_quad; ++d)
419  quad->set_mg_dof_index (cell->level(), d, next_free_dof++);
420  };
421 
422 
423  // dofs of cell
424  if (cell->get_fe().dofs_per_hex > 0)
425  for (unsigned int d=0; d<cell->get_fe().dofs_per_hex; ++d)
426  cell->set_mg_dof_index (cell->level(), d, next_free_dof++);
427 
428 
429  // note that this cell has
430  // been processed
431  cell->set_user_flag ();
432 
433  return next_free_dof;
434  }
435 
436 
437  template <int dim, int spacedim>
438  static
439  unsigned int
440  distribute_dofs_on_level (const unsigned int offset,
441  const types::subdomain_id level_subdomain_id,
442  DoFHandler<dim,spacedim> &dof_handler,
443  const unsigned int level)
444  {
445  const ::Triangulation<dim,spacedim> &tria
446  = dof_handler.get_triangulation();
447  Assert (tria.n_levels() > 0, ExcMessage("Empty triangulation"));
448  if (level>=tria.n_levels())
449  return 0; //this is allowed for multigrid
450 
451  // Clear user flags because we will
452  // need them. But first we save
453  // them and make sure that we
454  // restore them later such that at
455  // the end of this function the
456  // Triangulation will be in the
457  // same state as it was at the
458  // beginning of this function.
459  std::vector<bool> user_flags;
460  tria.save_user_flags(user_flags);
461  const_cast<::Triangulation<dim,spacedim> &>(tria).clear_user_flags ();
462 
463  unsigned int next_free_dof = offset;
464  typename DoFHandler<dim,spacedim>::level_cell_iterator
465  cell = dof_handler.begin(level),
466  endc = dof_handler.end(level);
467 
468  for (; cell != endc; ++cell)
469  if ((level_subdomain_id == numbers::invalid_subdomain_id)
470  ||
471  (cell->level_subdomain_id() == level_subdomain_id))
472  next_free_dof
473  = Implementation::distribute_mg_dofs_on_cell (dof_handler, cell, next_free_dof);
474 
475 // // update the cache used
476 // // for cell dof indices
477 // for (typename DoFHandler<dim,spacedim>::level_cell_iterator
478 // cell = dof_handler.begin(); cell != dof_handler.end(); ++cell)
479 // if (cell->subdomain_id() != numbers::artificial_subdomain_id)
480 // cell->update_cell_dof_indices_cache ();
481 
482  // finally restore the user flags
483  const_cast<::Triangulation<dim,spacedim> &>(tria).load_user_flags(user_flags);
484 
485  return next_free_dof;
486  }
487 
488 
489  /* --------------------- renumber_dofs functionality ---------------- */
490 
491 
508  template <int spacedim>
509  static
510  void
511  renumber_dofs (const std::vector<types::global_dof_index> &new_numbers,
512  const IndexSet &,
513  DoFHandler<1,spacedim> &dof_handler,
514  const bool check_validity)
515  {
516  // note that we can not use cell
517  // iterators in this function since
518  // then we would renumber the dofs on
519  // the interface of two cells more
520  // than once. Anyway, this way it's
521  // not only more correct but also
522  // faster; note, however, that dof
523  // numbers may be invalid_dof_index,
524  // namely when the appropriate
525  // vertex/line/etc is unused
526  for (std::vector<types::global_dof_index>::iterator
527  i=dof_handler.vertex_dofs.begin();
528  i!=dof_handler.vertex_dofs.end(); ++i)
529  if (*i != DoFHandler<1,spacedim>::invalid_dof_index)
530  *i = new_numbers[*i];
531  else if (check_validity)
532  // if index is
533  // invalid_dof_index:
534  // check if this one
535  // really is unused
536  Assert (dof_handler.get_triangulation()
537  .vertex_used((i-dof_handler.vertex_dofs.begin()) /
538  dof_handler.selected_fe->dofs_per_vertex)
539  == false,
540  ExcInternalError ());
541 
542  for (unsigned int level=0; level<dof_handler.levels.size(); ++level)
543  for (std::vector<types::global_dof_index>::iterator
544  i=dof_handler.levels[level]->dof_object.dofs.begin();
545  i!=dof_handler.levels[level]->dof_object.dofs.end(); ++i)
546  if (*i != DoFHandler<1,spacedim>::invalid_dof_index)
547  *i = new_numbers[*i];
548 
549  // update the cache
550  // used for cell dof
551  // indices
552  for (typename DoFHandler<1,spacedim>::level_cell_iterator
553  cell = dof_handler.begin();
554  cell != dof_handler.end(); ++cell)
555  cell->update_cell_dof_indices_cache ();
556  }
557 
558  template <int spacedim>
559  static
560  void
561  renumber_mg_dofs (const std::vector<::types::global_dof_index> &new_numbers,
562  const IndexSet &indices,
563  DoFHandler<1,spacedim> &dof_handler,
564  const unsigned int level,
565  const bool check_validity)
566  {
567  for (typename std::vector<typename DoFHandler<1,spacedim>::MGVertexDoFs>::iterator
568  i=dof_handler.mg_vertex_dofs.begin();
569  i!=dof_handler.mg_vertex_dofs.end();
570  ++i)
571  // if the present vertex lives on
572  // the current level
573  if ((i->get_coarsest_level() <= level) &&
574  (i->get_finest_level() >= level))
575  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_vertex; ++d)
576  {
577  ::types::global_dof_index idx = i->get_index (level, d);
578  if (idx != DoFHandler<1>::invalid_dof_index)
579  i->set_index (level, d,
580  (indices.n_elements() == 0)?
581  (new_numbers[idx]) :
582  (new_numbers[indices.index_within_set(idx)]));
583 
584  if (check_validity)
585  Assert(idx != DoFHandler<1>::invalid_dof_index, ExcInternalError ());
586  }
587 
588 
589  for (std::vector<types::global_dof_index>::iterator
590  i=dof_handler.mg_levels[level]->dof_object.dofs.begin();
591  i!=dof_handler.mg_levels[level]->dof_object.dofs.end();
592  ++i)
593  {
594  if (*i != DoFHandler<1>::invalid_dof_index)
595  {
596  Assert(*i<new_numbers.size(), ExcInternalError());
597  *i = (indices.n_elements() == 0)?
598  (new_numbers[*i]) :
599  (new_numbers[indices.index_within_set(*i)]);
600  }
601  }
602  }
603 
604 
605  template <int spacedim>
606  static
607  void
608  renumber_dofs (const std::vector<types::global_dof_index> &new_numbers,
609  const IndexSet &indices,
610  DoFHandler<2,spacedim> &dof_handler,
611  const bool check_validity)
612  {
613  // note that we can not use cell
614  // iterators in this function since
615  // then we would renumber the dofs on
616  // the interface of two cells more
617  // than once. Anyway, this way it's
618  // not only more correct but also
619  // faster; note, however, that dof
620  // numbers may be invalid_dof_index,
621  // namely when the appropriate
622  // vertex/line/etc is unused
623  for (std::vector<types::global_dof_index>::iterator
624  i=dof_handler.vertex_dofs.begin();
625  i!=dof_handler.vertex_dofs.end(); ++i)
626  if (*i != DoFHandler<2,spacedim>::invalid_dof_index)
627  *i = (indices.n_elements() == 0)?
628  (new_numbers[*i]) :
629  (new_numbers[indices.index_within_set(*i)]);
630  else if (check_validity)
631  // if index is invalid_dof_index:
632  // check if this one really is
633  // unused
634  Assert (dof_handler.get_triangulation()
635  .vertex_used((i-dof_handler.vertex_dofs.begin()) /
636  dof_handler.selected_fe->dofs_per_vertex)
637  == false,
638  ExcInternalError ());
639 
640  for (std::vector<types::global_dof_index>::iterator
641  i=dof_handler.faces->lines.dofs.begin();
642  i!=dof_handler.faces->lines.dofs.end(); ++i)
643  if (*i != DoFHandler<2,spacedim>::invalid_dof_index)
644  *i = ((indices.n_elements() == 0) ?
645  new_numbers[*i] :
646  new_numbers[indices.index_within_set(*i)]);
647 
648  for (unsigned int level=0; level<dof_handler.levels.size(); ++level)
649  {
650  for (std::vector<types::global_dof_index>::iterator
651  i=dof_handler.levels[level]->dof_object.dofs.begin();
652  i!=dof_handler.levels[level]->dof_object.dofs.end(); ++i)
653  if (*i != DoFHandler<2,spacedim>::invalid_dof_index)
654  *i = ((indices.n_elements() == 0) ?
655  new_numbers[*i] :
656  new_numbers[indices.index_within_set(*i)]);
657  }
658 
659  // update the cache
660  // used for cell dof
661  // indices
662  for (typename DoFHandler<2,spacedim>::level_cell_iterator
663  cell = dof_handler.begin();
664  cell != dof_handler.end(); ++cell)
665  cell->update_cell_dof_indices_cache ();
666  }
667 
668  template <int spacedim>
669  static
670  void
671  renumber_mg_dofs (const std::vector<::types::global_dof_index> &new_numbers,
672  const IndexSet &indices,
673  DoFHandler<2,spacedim> &dof_handler,
674  const unsigned int level,
675  const bool check_validity)
676  {
677  if (level>=dof_handler.get_triangulation().n_levels())
678  return;
679  for (typename std::vector<typename DoFHandler<2,spacedim>::MGVertexDoFs>::iterator i=dof_handler.mg_vertex_dofs.begin();
680  i!=dof_handler.mg_vertex_dofs.end(); ++i)
681  // if the present vertex lives on
682  // the present level
683  if ((i->get_coarsest_level() <= level) &&
684  (i->get_finest_level() >= level))
685  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_vertex; ++d)
686  {
687  ::types::global_dof_index idx =i->get_index (level, d/*, dof_handler.get_fe().dofs_per_vertex*/);
688  if (idx != DoFHandler<1>::invalid_dof_index)
689  i->set_index (level, d/*, dof_handler.get_fe().dofs_per_vertex*/,
690  ((indices.n_elements() == 0) ?
691  new_numbers[idx] :
692  new_numbers[indices.index_within_set(idx)]));
693 
694  if (check_validity)
695  Assert(idx != DoFHandler<2>::invalid_dof_index, ExcInternalError ());
696  }
697 
698  if (dof_handler.get_fe().dofs_per_line > 0)
699  {
700  // save user flags as they will be modified
701  std::vector<bool> user_flags;
702  dof_handler.get_triangulation().save_user_flags(user_flags);
703  const_cast<::Triangulation<2,spacedim> &>(dof_handler.get_triangulation()).clear_user_flags ();
704 
705  // flag all lines adjacent to cells of the current
706  // level, as those lines logically belong to the same
707  // level as the cell, at least for for isotropic
708  // refinement
709  typename DoFHandler<2,spacedim>::level_cell_iterator cell,
710  endc = dof_handler.end(level);
711  for (cell = dof_handler.begin(level); cell != endc; ++cell)
712  for (unsigned int line=0; line < GeometryInfo<2>::faces_per_cell; ++line)
713  cell->face(line)->set_user_flag();
714 
715  for (typename DoFHandler<2,spacedim>::cell_iterator cell = dof_handler.begin();
716  cell != dof_handler.end(); ++cell)
717  for (unsigned int l=0; l<GeometryInfo<2>::lines_per_cell; ++l)
718  if (cell->line(l)->user_flag_set())
719  {
720  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_line; ++d)
721  {
722  ::types::global_dof_index idx = cell->line(l)->mg_dof_index(level, d);
723  if (idx != DoFHandler<1>::invalid_dof_index)
724  cell->line(l)->set_mg_dof_index (level, d, ((indices.n_elements() == 0) ?
725  new_numbers[idx] :
726  new_numbers[indices.index_within_set(idx)]));
727  if (check_validity)
728  Assert(idx != DoFHandler<2>::invalid_dof_index, ExcInternalError ());
729  }
730  cell->line(l)->clear_user_flag();
731  }
732  // finally, restore user flags
733  const_cast<::Triangulation<2,spacedim> &>(dof_handler.get_triangulation()).load_user_flags (user_flags);
734  }
735 
736  for (std::vector<types::global_dof_index>::iterator i=dof_handler.mg_levels[level]->dof_object.dofs.begin();
737  i!=dof_handler.mg_levels[level]->dof_object.dofs.end(); ++i)
738  {
739  if (*i != DoFHandler<2>::invalid_dof_index)
740  {
741  Assert(*i<new_numbers.size(), ExcInternalError());
742  *i = ((indices.n_elements() == 0) ?
743  new_numbers[*i] :
744  new_numbers[indices.index_within_set(*i)]);
745  }
746  }
747 
748  }
749 
750 
751  template <int spacedim>
752  static
753  void
754  renumber_dofs (const std::vector<types::global_dof_index> &new_numbers,
755  const IndexSet &indices,
756  DoFHandler<3,spacedim> &dof_handler,
757  const bool check_validity)
758  {
759  // note that we can not use cell
760  // iterators in this function since
761  // then we would renumber the dofs on
762  // the interface of two cells more
763  // than once. Anyway, this way it's
764  // not only more correct but also
765  // faster; note, however, that dof
766  // numbers may be invalid_dof_index,
767  // namely when the appropriate
768  // vertex/line/etc is unused
769  for (std::vector<types::global_dof_index>::iterator
770  i=dof_handler.vertex_dofs.begin();
771  i!=dof_handler.vertex_dofs.end(); ++i)
772  if (*i != DoFHandler<3,spacedim>::invalid_dof_index)
773  *i = ((indices.n_elements() == 0) ?
774  new_numbers[*i] :
775  new_numbers[indices.index_within_set(*i)]);
776  else if (check_validity)
777  // if index is invalid_dof_index:
778  // check if this one really is
779  // unused
780  Assert (dof_handler.get_triangulation()
781  .vertex_used((i-dof_handler.vertex_dofs.begin()) /
782  dof_handler.selected_fe->dofs_per_vertex)
783  == false,
784  ExcInternalError ());
785 
786  for (std::vector<types::global_dof_index>::iterator
787  i=dof_handler.faces->lines.dofs.begin();
788  i!=dof_handler.faces->lines.dofs.end(); ++i)
789  if (*i != DoFHandler<3,spacedim>::invalid_dof_index)
790  *i = ((indices.n_elements() == 0) ?
791  new_numbers[*i] :
792  new_numbers[indices.index_within_set(*i)]);
793  for (std::vector<types::global_dof_index>::iterator
794  i=dof_handler.faces->quads.dofs.begin();
795  i!=dof_handler.faces->quads.dofs.end(); ++i)
796  if (*i != DoFHandler<3,spacedim>::invalid_dof_index)
797  *i = ((indices.n_elements() == 0) ?
798  new_numbers[*i] :
799  new_numbers[indices.index_within_set(*i)]);
800 
801  for (unsigned int level=0; level<dof_handler.levels.size(); ++level)
802  {
803  for (std::vector<types::global_dof_index>::iterator
804  i=dof_handler.levels[level]->dof_object.dofs.begin();
805  i!=dof_handler.levels[level]->dof_object.dofs.end(); ++i)
806  if (*i != DoFHandler<3,spacedim>::invalid_dof_index)
807  *i = ((indices.n_elements() == 0) ?
808  new_numbers[*i] :
809  new_numbers[indices.index_within_set(*i)]);
810  }
811 
812  // update the cache
813  // used for cell dof
814  // indices
815  for (typename DoFHandler<3,spacedim>::level_cell_iterator
816  cell = dof_handler.begin();
817  cell != dof_handler.end(); ++cell)
818  cell->update_cell_dof_indices_cache ();
819  }
820 
821  template <int spacedim>
822  static
823  void
824  renumber_mg_dofs (const std::vector<::types::global_dof_index> &new_numbers,
825  const IndexSet &indices,
826  DoFHandler<3,spacedim> &dof_handler,
827  const unsigned int level,
828  const bool check_validity)
829  {
830  if (level>=dof_handler.get_triangulation().n_levels())
831  return;
832  for (typename std::vector<typename DoFHandler<3,spacedim>::MGVertexDoFs>::iterator i=dof_handler.mg_vertex_dofs.begin();
833  i!=dof_handler.mg_vertex_dofs.end(); ++i)
834  // if the present vertex lives on
835  // the present level
836  if ((i->get_coarsest_level() <= level) &&
837  (i->get_finest_level() >= level))
838  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_vertex; ++d)
839  {
840  ::types::global_dof_index idx =i->get_index (level, d);
841  if (idx != DoFHandler<3>::invalid_dof_index)
842  i->set_index (level, d,
843  ((indices.n_elements() == 0) ?
844  new_numbers[idx] :
845  new_numbers[indices.index_within_set(idx)]));
846 
847  if (check_validity)
848  Assert(idx != DoFHandler<3>::invalid_dof_index, ExcInternalError ());
849  }
850 
851  if (dof_handler.get_fe().dofs_per_line > 0 ||
852  dof_handler.get_fe().dofs_per_quad > 0)
853  {
854  // save user flags as they will be modified
855  std::vector<bool> user_flags;
856  dof_handler.get_triangulation().save_user_flags(user_flags);
857  const_cast<::Triangulation<3,spacedim> &>(dof_handler.get_triangulation()).clear_user_flags ();
858 
859  // flag all lines adjacent to cells of the current level, as
860  // those lines logically belong to the same level as the cell,
861  // at least for isotropic refinement
862  typename DoFHandler<3,spacedim>::level_cell_iterator cell,
863  endc = dof_handler.end(level);
864  for (cell = dof_handler.begin(level); cell != endc; ++cell)
865  for (unsigned int line=0; line < GeometryInfo<3>::lines_per_cell; ++line)
866  cell->line(line)->set_user_flag();
867 
868  for (typename DoFHandler<3,spacedim>::cell_iterator cell = dof_handler.begin();
869  cell != dof_handler.end(); ++cell)
870  for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
871  if (cell->line(l)->user_flag_set())
872  {
873  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_line; ++d)
874  {
875  ::types::global_dof_index idx = cell->line(l)->mg_dof_index(level, d);
876  if (idx != DoFHandler<3>::invalid_dof_index)
877  cell->line(l)->set_mg_dof_index (level, d, ((indices.n_elements() == 0) ?
878  new_numbers[idx] :
879  new_numbers[indices.index_within_set(idx)]));
880  if (check_validity)
881  Assert(idx != DoFHandler<3>::invalid_dof_index, ExcInternalError ());
882  }
883  cell->line(l)->clear_user_flag();
884  }
885 
886  // flag all quads adjacent to cells of the current level, as
887  // those quads logically belong to the same level as the cell,
888  // at least for isotropic refinement
889  for (cell = dof_handler.begin(level); cell != endc; ++cell)
890  for (unsigned int quad=0; quad < GeometryInfo<3>::quads_per_cell; ++quad)
891  cell->quad(quad)->set_user_flag();
892 
893  for (typename DoFHandler<3,spacedim>::cell_iterator cell = dof_handler.begin();
894  cell != dof_handler.end(); ++cell)
895  for (unsigned int l=0; l<GeometryInfo<3>::quads_per_cell; ++l)
896  if (cell->quad(l)->user_flag_set())
897  {
898  for (unsigned int d=0; d<dof_handler.get_fe().dofs_per_quad; ++d)
899  {
900  ::types::global_dof_index idx = cell->quad(l)->mg_dof_index(level, d);
901  if (idx != DoFHandler<1>::invalid_dof_index)
902  cell->quad(l)->set_mg_dof_index (level, d, ((indices.n_elements() == 0) ?
903  new_numbers[idx] :
904  new_numbers[indices.index_within_set(idx)]));
905  if (check_validity)
906  Assert(idx != DoFHandler<3>::invalid_dof_index, ExcInternalError ());
907  }
908  cell->quad(l)->clear_user_flag();
909  }
910 
911  // finally, restore user flags
912  const_cast<::Triangulation<3,spacedim> &>(dof_handler.get_triangulation()).load_user_flags (user_flags);
913  }
914 
915  for (std::vector<types::global_dof_index>::iterator i=dof_handler.mg_levels[level]->dof_object.dofs.begin();
916  i!=dof_handler.mg_levels[level]->dof_object.dofs.end(); ++i)
917  {
918  if (*i != DoFHandler<3>::invalid_dof_index)
919  {
920  Assert(*i<new_numbers.size(), ExcInternalError());
921  *i = ((indices.n_elements() == 0) ?
922  new_numbers[*i] :
923  new_numbers[indices.index_within_set(*i)]);
924  }
925  }
926  }
927 
928 
929  };
930 
931 
932 
933  /* --------------------- class PolicyBase ---------------- */
934 
935  template <int dim, int spacedim>
936  PolicyBase<dim,spacedim>::~PolicyBase ()
937  {}
938 
939 
940  /* --------------------- class Sequential ---------------- */
941 
942 
943  template <int dim, int spacedim>
944  void
945  Sequential<dim,spacedim>::
946  distribute_dofs (DoFHandler<dim,spacedim> &dof_handler,
947  NumberCache &number_cache_current ) const
948  {
949  const types::global_dof_index n_dofs =
950  Implementation::distribute_dofs (0,
951  numbers::invalid_subdomain_id,
952  dof_handler);
953 
954  // now set the elements of the
955  // number cache appropriately
956  NumberCache number_cache;
957  number_cache.n_global_dofs = n_dofs;
958  number_cache.n_locally_owned_dofs = number_cache.n_global_dofs;
959 
960  number_cache.locally_owned_dofs
961  = IndexSet (number_cache.n_global_dofs);
962  number_cache.locally_owned_dofs.add_range (0,
963  number_cache.n_global_dofs);
964  number_cache.locally_owned_dofs.compress();
965 
966  number_cache.n_locally_owned_dofs_per_processor
967  = std::vector<types::global_dof_index> (1,
968  number_cache.n_global_dofs);
969 
970  number_cache.locally_owned_dofs_per_processor
971  = std::vector<IndexSet> (1,
972  number_cache.locally_owned_dofs);
973  number_cache_current = number_cache;
974  }
975 
976 
977  template <int dim, int spacedim>
978  void
979  Sequential<dim,spacedim>::
980  distribute_mg_dofs (DoFHandler<dim,spacedim> &dof_handler,
981  std::vector<NumberCache> &number_caches) const
982  {
983  std::vector<bool> user_flags;
984 
985  dof_handler.get_triangulation().save_user_flags (user_flags);
986  const_cast<::Triangulation<dim, spacedim>&>(dof_handler.get_triangulation()).clear_user_flags ();
987 
988  for (unsigned int level = 0; level < dof_handler.get_triangulation().n_levels(); ++level)
989  {
990  types::global_dof_index next_free_dof = Implementation::distribute_dofs_on_level(0, numbers::invalid_subdomain_id, dof_handler, level);
991 
992  number_caches[level].n_global_dofs = next_free_dof;
993  number_caches[level].n_locally_owned_dofs = next_free_dof;
994  number_caches[level].locally_owned_dofs = complete_index_set(next_free_dof);
995  number_caches[level].locally_owned_dofs_per_processor.resize(1);
996  number_caches[level].locally_owned_dofs_per_processor[0] = complete_index_set(next_free_dof);
997  number_caches[level].n_locally_owned_dofs_per_processor.resize(1);
998  number_caches[level].n_locally_owned_dofs_per_processor[0] = next_free_dof;
999  }
1000  const_cast<::Triangulation<dim, spacedim>&>(dof_handler.get_triangulation()).load_user_flags (user_flags);
1001  }
1002 
1003  template <int dim, int spacedim>
1004  void
1005  Sequential<dim,spacedim>::
1006  renumber_dofs (const std::vector<types::global_dof_index> &new_numbers,
1007  ::DoFHandler<dim,spacedim> &dof_handler,
1008  NumberCache &number_cache_current) const
1009  {
1010  Implementation::renumber_dofs (new_numbers, IndexSet(0),
1011  dof_handler, true);
1012 
1013  // in the sequential case,
1014  // the number cache should
1015  // not have changed but we
1016  // have to set the elements
1017  // of the structure
1018  // appropriately anyway
1019  NumberCache number_cache;
1020  number_cache.n_global_dofs = dof_handler.n_dofs();
1021  number_cache.n_locally_owned_dofs = number_cache.n_global_dofs;
1022 
1023  number_cache.locally_owned_dofs
1024  = IndexSet (number_cache.n_global_dofs);
1025  number_cache.locally_owned_dofs.add_range (0,
1026  number_cache.n_global_dofs);
1027  number_cache.locally_owned_dofs.compress();
1028 
1029  number_cache.n_locally_owned_dofs_per_processor
1030  = std::vector<types::global_dof_index> (1,
1031  number_cache.n_global_dofs);
1032 
1033  number_cache.locally_owned_dofs_per_processor
1034  = std::vector<IndexSet> (1,
1035  number_cache.locally_owned_dofs);
1036  number_cache_current = number_cache;
1037  }
1038 
1039  /* --------------------- class ParallelShared ---------------- */
1040 
1041  template <int dim, int spacedim>
1042  void
1043  ParallelShared<dim,spacedim>::
1044  distribute_dofs (DoFHandler<dim,spacedim> &dof_handler,
1045  NumberCache &number_cache) const
1046  {
1047  // If the underlying shared::Tria allows artifical cells, we need to do
1048  // some tricks here to make Sequential algorithms play nicely.
1049  // Namely, we first restore original partition (without artificial cells)
1050  // and then turn artificial cells on at the end of this function.
1051  const parallel::shared::Triangulation<dim, spacedim> *tr =
1052  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim>*> (&dof_handler.get_triangulation()));
1053  Assert(tr != 0, ExcInternalError());
1054  typename parallel::shared::Triangulation<dim,spacedim>::active_cell_iterator
1055  cell = dof_handler.get_triangulation().begin_active(),
1056  endc = dof_handler.get_triangulation().end();
1057  std::vector<types::subdomain_id> current_subdomain_ids(tr->n_active_cells());
1058  const std::vector<types::subdomain_id> &true_subdomain_ids = tr->get_true_subdomain_ids_of_cells();
1059  if (tr->with_artificial_cells())
1060  for (unsigned int index=0; cell != endc; cell++, index++)
1061  {
1062  current_subdomain_ids[index] = cell->subdomain_id();
1063  cell->set_subdomain_id(true_subdomain_ids[index]);
1064  }
1065 
1066  Sequential<dim,spacedim>::distribute_dofs (dof_handler,number_cache);
1067  DoFRenumbering::subdomain_wise (dof_handler);
1068  // dofrenumbering will reset subdomains, this is ugly but we need to do it again:
1069  cell = tr->begin_active();
1070  if (tr->with_artificial_cells())
1071  for (unsigned int index=0; cell != endc; cell++, index++)
1072  cell->set_subdomain_id(true_subdomain_ids[index]);
1073 
1074  number_cache.locally_owned_dofs_per_processor = DoFTools::locally_owned_dofs_per_subdomain (dof_handler);
1075  number_cache.locally_owned_dofs = number_cache.locally_owned_dofs_per_processor[dof_handler.get_triangulation().locally_owned_subdomain()];
1076  number_cache.n_locally_owned_dofs_per_processor.resize (number_cache.locally_owned_dofs_per_processor.size());
1077  for (unsigned int i = 0; i < number_cache.n_locally_owned_dofs_per_processor.size(); i++)
1078  number_cache.n_locally_owned_dofs_per_processor[i] = number_cache.locally_owned_dofs_per_processor[i].n_elements();
1079  number_cache.n_locally_owned_dofs = number_cache.n_locally_owned_dofs_per_processor[dof_handler.get_triangulation().locally_owned_subdomain()];
1080 
1081  // restore current subdomain ids
1082  cell = tr->begin_active();
1083  if (tr->with_artificial_cells())
1084  for (unsigned int index=0; cell != endc; cell++, index++)
1085  cell->set_subdomain_id(current_subdomain_ids[index]);
1086  }
1087 
1088  template <int dim, int spacedim>
1089  void
1090  ParallelShared<dim,spacedim>::
1091  distribute_mg_dofs (DoFHandler<dim,spacedim> &dof_handler,
1092  std::vector<NumberCache> &number_caches) const
1093  {
1094  // first, call the sequential function to distribute dofs
1095  Sequential<dim,spacedim>::distribute_mg_dofs (dof_handler, number_caches);
1096  // now we need to update the number cache.
1097  // This part is not yet implemented.
1098  AssertThrow(false,ExcNotImplemented());
1099  }
1100 
1101  template <int dim, int spacedim>
1102  void
1103  ParallelShared<dim,spacedim>::
1104  renumber_dofs (const std::vector<types::global_dof_index> &new_numbers,
1105  ::DoFHandler<dim,spacedim> &dof_handler,
1106  NumberCache &number_cache) const
1107  {
1108 
1109 #ifndef DEAL_II_WITH_MPI
1110  (void)new_numbers;
1111  (void)dof_handler;
1112  (void)number_cache;
1113  Assert (false, ExcNotImplemented());
1114 #else
1115  // Similar to distribute_dofs() we need to have a special treatment in
1116  // case artificial cells are present.
1117  const parallel::shared::Triangulation<dim, spacedim> *tr =
1118  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim>*> (&dof_handler.get_triangulation()));
1119  Assert(tr != 0, ExcInternalError());
1120  typename parallel::shared::Triangulation<dim,spacedim>::active_cell_iterator
1121  cell = dof_handler.get_triangulation().begin_active(),
1122  endc = dof_handler.get_triangulation().end();
1123  std::vector<types::subdomain_id> current_subdomain_ids(tr->n_active_cells());
1124  const std::vector<types::subdomain_id> &true_subdomain_ids = tr->get_true_subdomain_ids_of_cells();
1125  if (tr->with_artificial_cells())
1126  for (unsigned int index=0; cell != endc; cell++, index++)
1127  {
1128  current_subdomain_ids[index] = cell->subdomain_id();
1129  cell->set_subdomain_id(true_subdomain_ids[index]);
1130  }
1131 
1132  std::vector<types::global_dof_index> global_gathered_numbers (dof_handler.n_dofs (), 0);
1133  // as we call DoFRenumbering::subdomain_wise (dof_handler) from distribute_dofs(),
1134  // we need to support sequential-like input.
1135  // Distributed-like input from, for example, component_wise renumbering is also supported.
1136  if (new_numbers.size () == dof_handler.n_dofs ())
1137  {
1138  global_gathered_numbers = new_numbers;
1139  }
1140  else
1141  {
1142  Assert(new_numbers.size() == dof_handler.locally_owned_dofs().n_elements(),
1143  ExcInternalError());
1144  const unsigned int n_cpu = Utilities::MPI::n_mpi_processes (tr->get_communicator ());
1145  std::vector<types::global_dof_index> gathered_new_numbers (dof_handler.n_dofs (), 0);
1146  Assert(Utilities::MPI::this_mpi_process (tr->get_communicator ()) ==
1147  dof_handler.get_triangulation().locally_owned_subdomain (),
1148  ExcInternalError())
1149 
1150  //gather new numbers among processors into one vector
1151  {
1152  std::vector<types::global_dof_index> new_numbers_copy (new_numbers);
1153  // displs:
1154  // Entry i specifies the displacement (relative to recvbuf )
1155  // at which to place the incoming data from process i
1156  // rcounts:
1157  // containing the number of elements that are to be received from each process
1158  std::vector<int> displs(n_cpu),
1159  rcounts(n_cpu);
1160  types::global_dof_index shift = 0;
1161  //set rcounts based on new_numbers:
1162  int cur_count = new_numbers_copy.size ();
1163  int ierr = MPI_Allgather (&cur_count, 1, MPI_INT,
1164  &rcounts[0], 1, MPI_INT,
1165  tr->get_communicator ());
1166  AssertThrowMPI(ierr);
1167 
1168  for (unsigned int i = 0; i < n_cpu; i++)
1169  {
1170  displs[i] = shift;
1171  shift += rcounts[i];
1172  }
1173  Assert(((int)new_numbers_copy.size()) ==
1174  rcounts[Utilities::MPI::this_mpi_process (tr->get_communicator ())],
1175  ExcInternalError());
1176  ierr = MPI_Allgatherv (&new_numbers_copy[0], new_numbers_copy.size (),
1177  DEAL_II_DOF_INDEX_MPI_TYPE,
1178  &gathered_new_numbers[0], &rcounts[0],
1179  &displs[0],
1180  DEAL_II_DOF_INDEX_MPI_TYPE,
1181  tr->get_communicator ());
1182  AssertThrowMPI(ierr);
1183  }
1184 
1185  // put new numbers according to the current locally_owned_dofs_per_processor IndexSets
1186  types::global_dof_index shift = 0;
1187  // flag_1 and flag_2 are
1188  // used to control that there is a
1189  // one-to-one relation between old and new DoFs.
1190  std::vector<unsigned int> flag_1 (dof_handler.n_dofs (), 0),
1191  flag_2 (dof_handler.n_dofs (), 0);
1192  for (unsigned int i = 0; i < n_cpu; i++)
1193  {
1194  const IndexSet &iset =
1195  number_cache.locally_owned_dofs_per_processor[i];
1196  for (types::global_dof_index ind = 0;
1197  ind < iset.n_elements (); ind++)
1198  {
1199  const types::global_dof_index target = iset.nth_index_in_set (ind);
1200  const types::global_dof_index value = gathered_new_numbers[shift + ind];
1201  Assert(target < dof_handler.n_dofs(), ExcInternalError());
1202  Assert(value < dof_handler.n_dofs(), ExcInternalError());
1203  global_gathered_numbers[target] = value;
1204  flag_1[target]++;
1205  flag_2[value]++;
1206  }
1207  shift += iset.n_elements ();
1208  }
1209 
1210  Assert(*std::max_element(flag_1.begin(), flag_1.end()) == 1,
1211  ExcInternalError());
1212  Assert(*std::min_element(flag_1.begin(), flag_1.end()) == 1,
1213  ExcInternalError());
1214  Assert((*std::max_element(flag_2.begin(), flag_2.end())) == 1,
1215  ExcInternalError());
1216  Assert((*std::min_element(flag_2.begin(), flag_2.end())) == 1,
1217  ExcInternalError());
1218  }
1219  Sequential<dim, spacedim>::renumber_dofs (global_gathered_numbers, dof_handler, number_cache);
1220  // correct number_cache:
1221  number_cache.locally_owned_dofs_per_processor =
1222  DoFTools::locally_owned_dofs_per_subdomain (dof_handler);
1223  number_cache.locally_owned_dofs =
1224  number_cache.locally_owned_dofs_per_processor[dof_handler.get_triangulation().locally_owned_subdomain ()];
1225  // sequential renumbering returns a vector of size 1 here,
1226  // correct this:
1227  number_cache.n_locally_owned_dofs_per_processor.resize(number_cache.locally_owned_dofs_per_processor.size());
1228  for (unsigned int i = 0;
1229  i < number_cache.n_locally_owned_dofs_per_processor.size (); i++)
1230  number_cache.n_locally_owned_dofs_per_processor[i] = number_cache.locally_owned_dofs_per_processor[i].n_elements ();
1231 
1232  number_cache.n_locally_owned_dofs =
1233  number_cache.n_locally_owned_dofs_per_processor[dof_handler.get_triangulation().locally_owned_subdomain ()];
1234 
1235  // restore artificial cells
1236  cell = tr->begin_active();
1237  if (tr->with_artificial_cells())
1238  for (unsigned int index=0; cell != endc; cell++, index++)
1239  cell->set_subdomain_id(current_subdomain_ids[index]);
1240 #endif
1241  }
1242 
1243  /* --------------------- class ParallelDistributed ---------------- */
1244 
1245 #ifdef DEAL_II_WITH_P4EST
1246 
1247  namespace
1248  {
1249  template <int dim>
1250  struct types
1251  {
1252 
1264  struct cellinfo
1265  {
1266  std::vector<unsigned int> tree_index;
1267  std::vector<typename ::internal::p4est::types<dim>::quadrant> quadrants;
1268  std::vector<::types::global_dof_index> dofs;
1269 
1270  unsigned int bytes_for_buffer () const
1271  {
1272  return (sizeof(unsigned int)*2 +
1273  tree_index.size() * sizeof(unsigned int) +
1274  quadrants.size() * sizeof(typename ::internal::p4est
1275  ::types<dim>::quadrant) +
1276  dofs.size() * sizeof(::types::global_dof_index));
1277  }
1278 
1279  void pack_data (std::vector<char> &buffer) const
1280  {
1281  buffer.resize(bytes_for_buffer());
1282 
1283  char *ptr = &buffer[0];
1284 
1285  const unsigned int num_cells = tree_index.size();
1286  const unsigned int num_dofs = dofs.size();
1287  std::memcpy(ptr, &num_cells, sizeof(unsigned int));
1288  ptr += sizeof(unsigned int);
1289  std::memcpy(ptr, &num_dofs, sizeof(unsigned int));
1290  ptr += sizeof(unsigned int);
1291 
1292  std::memcpy(ptr,
1293  &tree_index[0],
1294  num_cells*sizeof(unsigned int));
1295  ptr += num_cells*sizeof(unsigned int);
1296 
1297  std::memcpy(ptr,
1298  &quadrants[0],
1299  num_cells * sizeof(typename ::internal::p4est::
1300  types<dim>::quadrant));
1301  ptr += num_cells*sizeof(typename ::internal::p4est::types<dim>::
1302  quadrant);
1303 
1304  if (num_dofs>0)
1305  std::memcpy(ptr,
1306  &dofs[0],
1307  dofs.size() * sizeof(::types::global_dof_index));
1308  ptr += dofs.size() * sizeof(::types::global_dof_index);
1309 
1310  Assert (ptr == &buffer[0]+buffer.size(),
1311  ExcInternalError());
1312  }
1313 
1314  void unpack_data(const std::vector<char> &buffer)
1315  {
1316  const char *ptr = &buffer[0];
1317  unsigned int num_cells;
1318  unsigned int num_dofs;
1319  memcpy(&num_cells, ptr, sizeof(unsigned int));
1320  ptr+=sizeof(unsigned int);
1321  memcpy(&num_dofs, ptr, sizeof(unsigned int));
1322  ptr+=sizeof(unsigned int);
1323 
1324  tree_index.resize(num_cells);
1325  std::memcpy(&tree_index[0],
1326  ptr,
1327  num_cells*sizeof(unsigned int));
1328  ptr += num_cells*sizeof(unsigned int);
1329 
1330  quadrants.resize(num_cells);
1331  std::memcpy(&quadrants[0],
1332  ptr,
1333  num_cells * sizeof(typename ::internal::p4est::
1334  types<dim>::quadrant));
1335  ptr += num_cells*sizeof(typename ::internal::p4est::types<dim>::
1336  quadrant);
1337 
1338  dofs.resize(num_dofs);
1339  if (num_dofs>0)
1340  std::memcpy(&dofs[0],
1341  ptr,
1342  dofs.size() * sizeof(::types::global_dof_index));
1343  ptr += dofs.size() * sizeof(::types::global_dof_index);
1344 
1345  Assert (ptr == &buffer[0]+buffer.size(),
1346  ExcInternalError());
1347  }
1348 
1349  };
1350 
1351  };
1352 
1353 
1354 
1355  template <int dim, int spacedim>
1356  void
1357  fill_dofindices_recursively (const typename parallel::distributed::Triangulation<dim,spacedim> &tria,
1358  const unsigned int tree_index,
1359  const typename DoFHandler<dim,spacedim>::level_cell_iterator &dealii_cell,
1360  const typename ::internal::p4est::types<dim>::quadrant &p4est_cell,
1361  const std::map<unsigned int, std::set<::types::subdomain_id> > &vertices_with_ghost_neighbors,
1362  std::map<::types::subdomain_id, typename types<dim>::cellinfo> &needs_to_get_cell)
1363  {
1364  // see if we have to
1365  // recurse...
1366  if (dealii_cell->has_children())
1367  {
1368  typename ::internal::p4est::types<dim>::quadrant
1369  p4est_child[GeometryInfo<dim>::max_children_per_cell];
1370  internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
1371 
1372 
1373  for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
1374  fill_dofindices_recursively<dim,spacedim>(tria,
1375  tree_index,
1376  dealii_cell->child(c),
1377  p4est_child[c],
1378  vertices_with_ghost_neighbors,
1379  needs_to_get_cell);
1380  return;
1381  }
1382 
1383  // we're at a leaf cell. see if
1384  // the cell is flagged as
1385  // interesting. note that we
1386  // have only flagged our own
1387  // cells before
1388  if (dealii_cell->user_flag_set() && !dealii_cell->is_ghost())
1389  {
1390  Assert (!dealii_cell->is_artificial(), ExcInternalError());
1391 
1392  // check each vertex if
1393  // it is interesting and
1394  // push dofindices if yes
1395  std::set<::types::subdomain_id> send_to;
1396  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
1397  {
1398  const std::map<unsigned int, std::set<::types::subdomain_id> >::const_iterator
1399  neighbor_subdomains_of_vertex
1400  = vertices_with_ghost_neighbors.find (dealii_cell->vertex_index(v));
1401 
1402  if (neighbor_subdomains_of_vertex ==
1403  vertices_with_ghost_neighbors.end())
1404  continue;
1405 
1406  Assert(neighbor_subdomains_of_vertex->second.size()!=0,
1407  ExcInternalError());
1408 
1409  send_to.insert(neighbor_subdomains_of_vertex->second.begin(),
1410  neighbor_subdomains_of_vertex->second.end());
1411  }
1412 
1413  if (send_to.size() > 0)
1414  {
1415  // this cell's dof_indices
1416  // need to be sent to
1417  // someone
1418  std::vector<::types::global_dof_index>
1419  local_dof_indices (dealii_cell->get_fe().dofs_per_cell);
1420  dealii_cell->get_dof_indices (local_dof_indices);
1421 
1422  for (std::set<::types::subdomain_id>::iterator it=send_to.begin();
1423  it!=send_to.end(); ++it)
1424  {
1425  const ::types::subdomain_id subdomain = *it;
1426 
1427  // get an iterator
1428  // to what needs to
1429  // be sent to that
1430  // subdomain (if
1431  // already exists),
1432  // or create such
1433  // an object
1434  typename std::map<::types::subdomain_id, typename types<dim>::cellinfo>::iterator
1435  p
1436  = needs_to_get_cell.insert (std::make_pair(subdomain,
1437  typename types<dim>::cellinfo()))
1438  .first;
1439 
1440  p->second.tree_index.push_back(tree_index);
1441  p->second.quadrants.push_back(p4est_cell);
1442 
1443  p->second.dofs.push_back(dealii_cell->get_fe().dofs_per_cell);
1444  p->second.dofs.insert(p->second.dofs.end(),
1445  local_dof_indices.begin(),
1446  local_dof_indices.end());
1447 
1448  }
1449  }
1450  }
1451  }
1452 
1453  template <int dim, int spacedim>
1454  void
1455  get_mg_dofindices_recursively (
1456  const parallel::distributed::Triangulation<dim,spacedim> &tria,
1457  const typename ::internal::p4est::types<dim>::quadrant &p4est_cell,
1458  const typename DoFHandler<dim,spacedim>::level_cell_iterator &dealii_cell,
1459  const typename ::internal::p4est::types<dim>::quadrant &quadrant,
1460  typename types<dim>::cellinfo &cell_info)
1461  {
1462  if (internal::p4est::quadrant_is_equal<dim>(p4est_cell, quadrant))
1463  {
1464  // why would somebody request a cell that is not ours?
1465  Assert(dealii_cell->level_subdomain_id()==tria.locally_owned_subdomain(), ExcInternalError());
1466 
1467 
1468  std::vector<::types::global_dof_index>
1469  local_dof_indices (dealii_cell->get_fe().dofs_per_cell);
1470  dealii_cell->get_mg_dof_indices (local_dof_indices);
1471 
1472  cell_info.dofs.push_back(dealii_cell->get_fe().dofs_per_cell);
1473  cell_info.dofs.insert(cell_info.dofs.end(),
1474  local_dof_indices.begin(),
1475  local_dof_indices.end());
1476  return; // we are done
1477  }
1478 
1479  if (! dealii_cell->has_children())
1480  return;
1481 
1482  if (! internal::p4est::quadrant_is_ancestor<dim> (p4est_cell, quadrant))
1483  return;
1484 
1485  typename ::internal::p4est::types<dim>::quadrant
1486  p4est_child[GeometryInfo<dim>::max_children_per_cell];
1487  internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
1488 
1489  for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
1490  get_mg_dofindices_recursively<dim,spacedim> (tria, p4est_child[c],
1491  dealii_cell->child(c),
1492  quadrant, cell_info);
1493  }
1494 
1495 
1496  template <int dim, int spacedim>
1497  void
1498  find_marked_mg_ghost_cells_recursively(const typename parallel::distributed::Triangulation<dim,spacedim> &tria,
1499  const unsigned int tree_index,
1500  const typename DoFHandler<dim,spacedim>::level_cell_iterator &dealii_cell,
1501  const typename ::internal::p4est::types<dim>::quadrant &p4est_cell,
1502  std::map<::types::subdomain_id, typename types<dim>::cellinfo> &neighbor_cell_list)
1503  {
1504  // recurse...
1505  if (dealii_cell->has_children())
1506  {
1507  typename ::internal::p4est::types<dim>::quadrant
1508  p4est_child[GeometryInfo<dim>::max_children_per_cell];
1509  internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
1510 
1511 
1512  for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
1513  find_marked_mg_ghost_cells_recursively<dim,spacedim>(tria,
1514  tree_index,
1515  dealii_cell->child(c),
1516  p4est_child[c],
1517  neighbor_cell_list);
1518  }
1519 
1520  if (dealii_cell->user_flag_set() && dealii_cell->level_subdomain_id() != tria.locally_owned_subdomain())
1521  {
1522  typename types<dim>::cellinfo &cell_info = neighbor_cell_list[dealii_cell->level_subdomain_id()];
1523 
1524  cell_info.tree_index.push_back(tree_index);
1525  cell_info.quadrants.push_back(p4est_cell);
1526  }
1527  }
1528 
1529 
1530  template <int dim, int spacedim>
1531  void
1532  set_mg_dofindices_recursively (
1533  const parallel::distributed::Triangulation<dim,spacedim> &tria,
1534  const typename ::internal::p4est::types<dim>::quadrant &p4est_cell,
1535  const typename DoFHandler<dim,spacedim>::level_cell_iterator &dealii_cell,
1536  const typename ::internal::p4est::types<dim>::quadrant &quadrant,
1537  ::types::global_dof_index *dofs)
1538  {
1539  if (internal::p4est::quadrant_is_equal<dim>(p4est_cell, quadrant))
1540  {
1541  Assert(dealii_cell->level_subdomain_id()!=::numbers::artificial_subdomain_id, ExcInternalError());
1542 
1543  // update dof indices of cell
1544  std::vector<::types::global_dof_index>
1545  dof_indices (dealii_cell->get_fe().dofs_per_cell);
1546  dealii_cell->get_mg_dof_indices(dof_indices);
1547 
1548  bool complete = true;
1549  for (unsigned int i=0; i<dof_indices.size(); ++i)
1550  if (dofs[i] != DoFHandler<dim,spacedim>::invalid_dof_index)
1551  {
1552  Assert((dof_indices[i] ==
1553  (DoFHandler<dim,spacedim>::invalid_dof_index))
1554  ||
1555  (dof_indices[i]==dofs[i]),
1556  ExcInternalError());
1557  dof_indices[i]=dofs[i];
1558  }
1559  else
1560  complete=false;
1561 
1562  if (!complete)
1563  const_cast
1564  <typename DoFHandler<dim,spacedim>::level_cell_iterator &>
1565  (dealii_cell)->set_user_flag();
1566  else
1567  const_cast
1568  <typename DoFHandler<dim,spacedim>::level_cell_iterator &>
1569  (dealii_cell)->clear_user_flag();
1570 
1571  const_cast
1572  <typename DoFHandler<dim,spacedim>::level_cell_iterator &>
1573  (dealii_cell)->set_mg_dof_indices(dof_indices);
1574  return;
1575  }
1576 
1577  if (! dealii_cell->has_children())
1578  return;
1579 
1580  if (! internal::p4est::quadrant_is_ancestor<dim> (p4est_cell, quadrant))
1581  return;
1582 
1583  typename ::internal::p4est::types<dim>::quadrant
1584  p4est_child[GeometryInfo<dim>::max_children_per_cell];
1585  internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
1586 
1587  for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
1588  set_mg_dofindices_recursively<dim,spacedim> (tria, p4est_child[c],
1589  dealii_cell->child(c),
1590  quadrant, dofs);
1591 
1592  }
1593 
1594 
1595 
1596  template <int dim, int spacedim>
1597  void
1598  communicate_mg_ghost_cells(const typename parallel::distributed::Triangulation<dim,spacedim> &tria,
1599  DoFHandler<dim,spacedim> &dof_handler,
1600  const std::vector<::types::global_dof_index> &coarse_cell_to_p4est_tree_permutation,
1601  const std::vector<::types::global_dof_index> &p4est_tree_to_coarse_cell_permutation)
1602  {
1603  // build list of cells to request for each neighbor
1604  std::set<::types::subdomain_id> level_ghost_owners = tria.level_ghost_owners();
1605  typedef std::map<::types::subdomain_id, typename types<dim>::cellinfo> cellmap_t;
1606  cellmap_t neighbor_cell_list;
1607  for (std::set<::types::subdomain_id>::iterator it = level_ghost_owners.begin();
1608  it != level_ghost_owners.end();
1609  ++it)
1610  neighbor_cell_list.insert(std::make_pair(*it, typename types<dim>::cellinfo()));
1611 
1612  for (typename DoFHandler<dim,spacedim>::level_cell_iterator
1613  cell = dof_handler.begin(0);
1614  cell != dof_handler.end(0);
1615  ++cell)
1616  {
1617  typename ::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
1618  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
1619 
1620  find_marked_mg_ghost_cells_recursively<dim,spacedim>
1621  (tria,
1622  coarse_cell_to_p4est_tree_permutation[cell->index()],
1623  cell,
1624  p4est_coarse_cell,
1625  neighbor_cell_list);
1626  }
1627  Assert(level_ghost_owners.size() == neighbor_cell_list.size(), ExcInternalError());
1628 
1629  //* send our requests:
1630  std::vector<std::vector<char> > sendbuffers (level_ghost_owners.size());
1631  std::vector<MPI_Request> requests (level_ghost_owners.size());
1632 
1633  unsigned int idx=0;
1634  for (typename cellmap_t::iterator it = neighbor_cell_list.begin();
1635  it!=neighbor_cell_list.end();
1636  ++it, ++idx)
1637  {
1638  // pack all the data into
1639  // the buffer for this
1640  // recipient and send
1641  // it. keep data around
1642  // till we can make sure
1643  // that the packet has been
1644  // received
1645  it->second.pack_data (sendbuffers[idx]);
1646  const int ierr = MPI_Isend(sendbuffers[idx].data(), sendbuffers[idx].size(),
1647  MPI_BYTE, it->first,
1648  1100101, tria.get_communicator(), &requests[idx]);
1649  AssertThrowMPI(ierr);
1650  }
1651 
1652  //* receive requests and reply
1653  std::vector<std::vector<char> > reply_buffers (level_ghost_owners.size());
1654  std::vector<MPI_Request> reply_requests (level_ghost_owners.size());
1655 
1656  for (unsigned int idx=0; idx<level_ghost_owners.size(); ++idx)
1657  {
1658  std::vector<char> receive;
1659  typename types<dim>::cellinfo cellinfo;
1660 
1661  MPI_Status status;
1662  int len;
1663  int ierr = MPI_Probe(MPI_ANY_SOURCE, 1100101, tria.get_communicator(), &status);
1664  AssertThrowMPI(ierr);
1665  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
1666  AssertThrowMPI(ierr);
1667  receive.resize(len);
1668 
1669  char *ptr = &receive[0];
1670  ierr = MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG,
1671  tria.get_communicator(), &status);
1672  AssertThrowMPI(ierr);
1673 
1674  cellinfo.unpack_data(receive);
1675 
1676  // store the dof indices for each cell
1677  for (unsigned int c=0; c<cellinfo.tree_index.size(); ++c)
1678  {
1679  typename DoFHandler<dim,spacedim>::level_cell_iterator
1680  cell (&dof_handler.get_triangulation(),
1681  0,
1682  p4est_tree_to_coarse_cell_permutation[cellinfo.tree_index[c]],
1683  &dof_handler);
1684 
1685  typename ::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
1686  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
1687 
1688  get_mg_dofindices_recursively<dim,spacedim> (tria,
1689  p4est_coarse_cell,
1690  cell,
1691  cellinfo.quadrants[c],
1692  cellinfo);
1693  }
1694 
1695  //send reply
1696  cellinfo.pack_data(reply_buffers[idx]);
1697  ierr = MPI_Isend(&(reply_buffers[idx])[0], reply_buffers[idx].size(),
1698  MPI_BYTE, status.MPI_SOURCE,
1699  1100102, tria.get_communicator(), &reply_requests[idx]);
1700  AssertThrowMPI(ierr);
1701  }
1702 
1703  // * finally receive the replies
1704  for (unsigned int idx=0; idx<level_ghost_owners.size(); ++idx)
1705  {
1706  std::vector<char> receive;
1707  typename types<dim>::cellinfo cellinfo;
1708 
1709  MPI_Status status;
1710  int len;
1711  int ierr = MPI_Probe(MPI_ANY_SOURCE, 1100102, tria.get_communicator(), &status);
1712  AssertThrowMPI(ierr);
1713  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
1714  AssertThrowMPI(ierr);
1715  receive.resize(len);
1716 
1717  char *ptr = &receive[0];
1718  ierr = MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG,
1719  tria.get_communicator(), &status);
1720  AssertThrowMPI(ierr);
1721 
1722  cellinfo.unpack_data(receive);
1723  if (cellinfo.tree_index.size()==0)
1724  continue;
1725 
1726  // set the dof indices for each cell
1727  ::types::global_dof_index *dofs = cellinfo.dofs.data();
1728  for (unsigned int c=0; c<cellinfo.tree_index.size(); ++c, dofs+=1+dofs[0])
1729  {
1730  typename DoFHandler<dim,spacedim>::level_cell_iterator
1731  cell (&tria,
1732  0,
1733  p4est_tree_to_coarse_cell_permutation[cellinfo.tree_index[c]],
1734  &dof_handler);
1735 
1736  typename ::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
1737  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
1738 
1739  Assert(cell->get_fe().dofs_per_cell==dofs[0], ExcInternalError());
1740 
1741  set_mg_dofindices_recursively<dim,spacedim> (tria,
1742  p4est_coarse_cell,
1743  cell,
1744  cellinfo.quadrants[c],
1745  dofs+1);
1746  }
1747  }
1748 
1749  // complete all sends, so that we can
1750  // safely destroy the buffers.
1751  if (requests.size() > 0)
1752  {
1753  const int ierr = MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE);
1754  AssertThrowMPI(ierr);
1755  }
1756  if (reply_requests.size() > 0)
1757  {
1758  const int ierr = MPI_Waitall(reply_requests.size(), &reply_requests[0], MPI_STATUSES_IGNORE);
1759  AssertThrowMPI(ierr);
1760  }
1761 
1762  }
1763 
1764 
1765  template <int spacedim>
1766  void
1767  communicate_mg_ghost_cells(const typename parallel::distributed::Triangulation<1,spacedim> &,
1768  DoFHandler<1,spacedim> &,
1769  const std::vector<::types::global_dof_index> &,
1770  const std::vector<::types::global_dof_index> &)
1771  {
1772  Assert (false, ExcNotImplemented());
1773  }
1774 
1775 
1776 
1777 
1778  template <int dim, int spacedim>
1779  void
1780  set_dofindices_recursively (
1781  const parallel::distributed::Triangulation<dim,spacedim> &tria,
1782  const typename ::internal::p4est::types<dim>::quadrant &p4est_cell,
1783  const typename DoFHandler<dim,spacedim>::level_cell_iterator &dealii_cell,
1784  const typename ::internal::p4est::types<dim>::quadrant &quadrant,
1785  ::types::global_dof_index *dofs)
1786  {
1787  if (internal::p4est::quadrant_is_equal<dim>(p4est_cell, quadrant))
1788  {
1789  Assert(!dealii_cell->has_children(), ExcInternalError());
1790  Assert(dealii_cell->is_ghost(), ExcInternalError());
1791 
1792  // update dof indices of cell
1793  std::vector<::types::global_dof_index>
1794  dof_indices (dealii_cell->get_fe().dofs_per_cell);
1795  dealii_cell->update_cell_dof_indices_cache();
1796  dealii_cell->get_dof_indices(dof_indices);
1797 
1798  bool complete = true;
1799  for (unsigned int i=0; i<dof_indices.size(); ++i)
1800  if (dofs[i] != DoFHandler<dim,spacedim>::invalid_dof_index)
1801  {
1802  Assert((dof_indices[i] ==
1803  (DoFHandler<dim,spacedim>::invalid_dof_index))
1804  ||
1805  (dof_indices[i]==dofs[i]),
1806  ExcInternalError());
1807  dof_indices[i]=dofs[i];
1808  }
1809  else
1810  complete=false;
1811 
1812  if (!complete)
1813  const_cast
1814  <typename DoFHandler<dim,spacedim>::level_cell_iterator &>
1815  (dealii_cell)->set_user_flag();
1816  else
1817  const_cast
1818  <typename DoFHandler<dim,spacedim>::level_cell_iterator &>
1819  (dealii_cell)->clear_user_flag();
1820 
1821  const_cast
1822  <typename DoFHandler<dim,spacedim>::level_cell_iterator &>
1823  (dealii_cell)->set_dof_indices(dof_indices);
1824 
1825  return;
1826  }
1827 
1828  if (! dealii_cell->has_children())
1829  return;
1830 
1831  if (! internal::p4est::quadrant_is_ancestor<dim> (p4est_cell, quadrant))
1832  return;
1833 
1834  typename ::internal::p4est::types<dim>::quadrant
1835  p4est_child[GeometryInfo<dim>::max_children_per_cell];
1836  internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
1837 
1838  for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
1839  set_dofindices_recursively<dim,spacedim> (tria, p4est_child[c],
1840  dealii_cell->child(c),
1841  quadrant, dofs);
1842  }
1843 
1844 
1845 
1846  template <int spacedim>
1847  void
1848  communicate_dof_indices_on_marked_cells
1849  (const DoFHandler<1,spacedim> &,
1850  const std::map<unsigned int, std::set<::types::subdomain_id> > &,
1851  const std::vector<::types::global_dof_index> &,
1852  const std::vector<::types::global_dof_index> &)
1853  {
1854  Assert (false, ExcNotImplemented());
1855  }
1856 
1857 
1858 
1859  template <int dim, int spacedim>
1860  void
1861  communicate_dof_indices_on_marked_cells
1862  (const DoFHandler<dim,spacedim> &dof_handler,
1863  const std::map<unsigned int, std::set<::types::subdomain_id> > &vertices_with_ghost_neighbors,
1864  const std::vector<::types::global_dof_index> &coarse_cell_to_p4est_tree_permutation,
1865  const std::vector<::types::global_dof_index> &p4est_tree_to_coarse_cell_permutation)
1866  {
1867 #ifndef DEAL_II_WITH_P4EST
1868  (void)vertices_with_ghost_neighbors;
1869  Assert (false, ExcNotImplemented());
1870 #else
1871 
1872  const parallel::distributed::Triangulation< dim, spacedim > *tr
1873  = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
1874  (&dof_handler.get_triangulation()));
1875  Assert (tr != 0, ExcInternalError());
1876 
1877  // now collect cells and their
1878  // dof_indices for the
1879  // interested neighbors
1880  typedef
1881  std::map<::types::subdomain_id, typename types<dim>::cellinfo>
1882  cellmap_t;
1883  cellmap_t needs_to_get_cells;
1884 
1885  for (typename DoFHandler<dim,spacedim>::level_cell_iterator
1886  cell = dof_handler.begin(0);
1887  cell != dof_handler.end(0);
1888  ++cell)
1889  {
1890  typename ::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
1891  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
1892 
1893  fill_dofindices_recursively<dim,spacedim>
1894  (*tr,
1895  coarse_cell_to_p4est_tree_permutation[cell->index()],
1896  cell,
1897  p4est_coarse_cell,
1898  vertices_with_ghost_neighbors,
1899  needs_to_get_cells);
1900  }
1901 
1902 
1903  //sending
1904  std::vector<std::vector<char> > sendbuffers (needs_to_get_cells.size());
1905  std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
1906  std::vector<MPI_Request> requests (needs_to_get_cells.size());
1907 
1908  unsigned int idx=0;
1909 
1910  for (typename cellmap_t::iterator it=needs_to_get_cells.begin();
1911  it!=needs_to_get_cells.end();
1912  ++it, ++buffer, ++idx)
1913  {
1914  const unsigned int num_cells = it->second.tree_index.size();
1915  (void)num_cells;
1916 
1917  Assert(num_cells==it->second.quadrants.size(), ExcInternalError());
1918  Assert(num_cells>0, ExcInternalError());
1919 
1920  // pack all the data into
1921  // the buffer for this
1922  // recipient and send
1923  // it. keep data around
1924  // till we can make sure
1925  // that the packet has been
1926  // received
1927  it->second.pack_data (*buffer);
1928  const int ierr = MPI_Isend(&(*buffer)[0], buffer->size(),
1929  MPI_BYTE, it->first,
1930  123, tr->get_communicator(), &requests[idx]);
1931  AssertThrowMPI(ierr);
1932  }
1933 
1934 
1935  // mark all own cells, that miss some
1936  // dof_data and collect the neighbors
1937  // that are going to send stuff to us
1938  std::set<::types::subdomain_id> senders;
1939  {
1940  std::vector<::types::global_dof_index> local_dof_indices;
1941  typename DoFHandler<dim,spacedim>::active_cell_iterator
1942  cell, endc = dof_handler.end();
1943 
1944  for (cell = dof_handler.begin_active(); cell != endc; ++cell)
1945  if (!cell->is_artificial())
1946  {
1947  if (cell->is_ghost())
1948  {
1949  if (cell->user_flag_set())
1950  senders.insert(cell->subdomain_id());
1951  }
1952  else
1953  {
1954  local_dof_indices.resize (cell->get_fe().dofs_per_cell);
1955  cell->get_dof_indices (local_dof_indices);
1956  if (local_dof_indices.end() !=
1957  std::find (local_dof_indices.begin(),
1958  local_dof_indices.end(),
1959  DoFHandler<dim,spacedim>::invalid_dof_index))
1960  cell->set_user_flag();
1961  else
1962  cell->clear_user_flag();
1963  }
1964 
1965  }
1966  }
1967 
1968 
1969  //* 5. receive ghostcelldata
1970  std::vector<char> receive;
1971  for (unsigned int i=0; i<senders.size(); ++i)
1972  {
1973  MPI_Status status;
1974  int len;
1975  int ierr = MPI_Probe(MPI_ANY_SOURCE, 123, tr->get_communicator(), &status);
1976  AssertThrowMPI(ierr);
1977  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
1978  AssertThrowMPI(ierr);
1979  receive.resize(len);
1980 
1981  char *ptr = &receive[0];
1982  ierr = MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG,
1983  tr->get_communicator(), &status);
1984  AssertThrowMPI(ierr);
1985 
1986  typename types<dim>::cellinfo cellinfo;
1987  cellinfo.unpack_data(receive);
1988  unsigned int cells = cellinfo.tree_index.size();
1989  ::types::global_dof_index *dofs = cellinfo.dofs.data();
1990 
1991  // the dofs pointer contains for each cell the number of dofs
1992  // on that cell (dofs[0]) followed by the dof indices itself.
1993  for (unsigned int c=0; c<cells; ++c, dofs+=1+dofs[0])
1994  {
1995  typename DoFHandler<dim,spacedim>::level_cell_iterator
1996  cell (&dof_handler.get_triangulation(),
1997  0,
1998  p4est_tree_to_coarse_cell_permutation[cellinfo.tree_index[c]],
1999  &dof_handler);
2000 
2001  typename ::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
2002  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
2003 
2004  Assert(cell->get_fe().dofs_per_cell==dofs[0], ExcInternalError());
2005 
2006  set_dofindices_recursively<dim,spacedim> (*tr,
2007  p4est_coarse_cell,
2008  cell,
2009  cellinfo.quadrants[c],
2010  (dofs+1));
2011  }
2012  }
2013 
2014  // complete all sends, so that we can
2015  // safely destroy the buffers.
2016  if (requests.size() > 0)
2017  {
2018  const int ierr = MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE);
2019  AssertThrowMPI(ierr);
2020  }
2021 
2022 
2023 #ifdef DEBUG
2024  {
2025  //check all msgs got sent and received
2026  unsigned int sum_send=0;
2027  unsigned int sum_recv=0;
2028  unsigned int sent=needs_to_get_cells.size();
2029  unsigned int recv=senders.size();
2030 
2031  int ierr = MPI_Allreduce(&sent, &sum_send, 1, MPI_UNSIGNED, MPI_SUM, tr->get_communicator());
2032  AssertThrowMPI(ierr);
2033  ierr = MPI_Allreduce(&recv, &sum_recv, 1, MPI_UNSIGNED, MPI_SUM, tr->get_communicator());
2034  AssertThrowMPI(ierr);
2035  Assert(sum_send==sum_recv, ExcInternalError());
2036  }
2037 #endif
2038 
2039  //update dofindices
2040  {
2041  typename DoFHandler<dim,spacedim>::active_cell_iterator
2042  cell, endc = dof_handler.end();
2043 
2044  for (cell = dof_handler.begin_active(); cell != endc; ++cell)
2045  if (!cell->is_artificial())
2046  cell->update_cell_dof_indices_cache();
2047  }
2048 
2049  // important, so that sends between two
2050  // calls to this function are not mixed
2051  // up.
2052  //
2053  // this is necessary because above we
2054  // just see if there are messages and
2055  // then receive them, without
2056  // discriminating where they come from
2057  // and whether they were sent in phase
2058  // 1 or 2. the need for a global
2059  // communication step like this barrier
2060  // could be avoided by receiving
2061  // messages specifically from those
2062  // processors from which we expect
2063  // messages, and by using different
2064  // tags for phase 1 and 2
2065  const int ierr = MPI_Barrier(tr->get_communicator());
2066  AssertThrowMPI(ierr);
2067 #endif
2068  }
2069 
2070 
2071 
2072 
2073 
2074 
2075 
2076  }
2077 
2078 #endif // DEAL_II_WITH_P4EST
2079 
2080 
2081 
2082  template <int dim, int spacedim>
2083  void
2084  ParallelDistributed<dim, spacedim>::
2085  distribute_dofs (DoFHandler<dim,spacedim> &dof_handler,
2086  NumberCache &number_cache_current) const
2087  {
2088  NumberCache number_cache;
2089 
2090 #ifndef DEAL_II_WITH_P4EST
2091  (void)dof_handler;
2092  Assert (false, ExcNotImplemented());
2093 #else
2094 
2095  parallel::distributed::Triangulation< dim, spacedim > *tr
2096  = (dynamic_cast<parallel::distributed::Triangulation<dim,spacedim>*>
2097  (const_cast<::Triangulation< dim, spacedim >*>
2098  (&dof_handler.get_triangulation())));
2099  Assert (tr != 0, ExcInternalError());
2100 
2101  const unsigned int
2102  n_cpus = Utilities::MPI::n_mpi_processes (tr->get_communicator());
2103 
2104  //* 1. distribute on own
2105  //* subdomain
2106  const ::types::global_dof_index n_initial_local_dofs =
2107  Implementation::distribute_dofs (0, tr->locally_owned_subdomain(),
2108  dof_handler);
2109 
2110  //* 2. iterate over ghostcells and
2111  //kill dofs that are not owned
2112  //by us
2113  std::vector<::types::global_dof_index> renumbering(n_initial_local_dofs);
2114  for (unsigned int i=0; i<renumbering.size(); ++i)
2115  renumbering[i] = i;
2116 
2117  {
2118  std::vector<::types::global_dof_index> local_dof_indices;
2119 
2120  typename DoFHandler<dim,spacedim>::active_cell_iterator
2121  cell = dof_handler.begin_active(),
2122  endc = dof_handler.end();
2123 
2124  for (; cell != endc; ++cell)
2125  if (cell->is_ghost() &&
2126  (cell->subdomain_id() < tr->locally_owned_subdomain()))
2127  {
2128  // we found a
2129  // neighboring ghost
2130  // cell whose subdomain
2131  // is "stronger" than
2132  // our own subdomain
2133 
2134  // delete all dofs that
2135  // live there and that
2136  // we have previously
2137  // assigned a number to
2138  // (i.e. the ones on
2139  // the interface)
2140  local_dof_indices.resize (cell->get_fe().dofs_per_cell);
2141  cell->get_dof_indices (local_dof_indices);
2142  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
2143  if (local_dof_indices[i] != DoFHandler<dim,spacedim>::invalid_dof_index)
2144  renumbering[local_dof_indices[i]]
2145  = DoFHandler<dim,spacedim>::invalid_dof_index;
2146  }
2147  }
2148 
2149 
2150  // make indices consecutive
2151  number_cache.n_locally_owned_dofs = 0;
2152  for (std::vector<::types::global_dof_index>::iterator it=renumbering.begin();
2153  it!=renumbering.end(); ++it)
2154  if (*it != DoFHandler<dim,spacedim>::invalid_dof_index)
2155  *it = number_cache.n_locally_owned_dofs++;
2156 
2157  //* 3. communicate local dofcount and
2158  //shift ids to make them unique
2159  number_cache.n_locally_owned_dofs_per_processor.resize(n_cpus);
2160 
2161  const int ierr = MPI_Allgather ( &number_cache.n_locally_owned_dofs,
2162  1, DEAL_II_DOF_INDEX_MPI_TYPE,
2163  &number_cache.n_locally_owned_dofs_per_processor[0],
2164  1, DEAL_II_DOF_INDEX_MPI_TYPE,
2165  tr->get_communicator());
2166  AssertThrowMPI(ierr);
2167 
2168  const ::types::global_dof_index
2169  shift = std::accumulate (number_cache
2170  .n_locally_owned_dofs_per_processor.begin(),
2171  number_cache
2172  .n_locally_owned_dofs_per_processor.begin()
2173  + tr->locally_owned_subdomain(),
2174  static_cast<::types::global_dof_index>(0));
2175  for (std::vector<::types::global_dof_index>::iterator it=renumbering.begin();
2176  it!=renumbering.end(); ++it)
2177  if (*it != DoFHandler<dim,spacedim>::invalid_dof_index)
2178  (*it) += shift;
2179 
2180  // now re-enumerate all dofs to
2181  // this shifted and condensed
2182  // numbering form. we renumber
2183  // some dofs as invalid, so
2184  // choose the nocheck-version.
2185  Implementation::renumber_dofs (renumbering, IndexSet(0),
2186  dof_handler, false);
2187 
2188  // now a little bit of
2189  // housekeeping
2190  number_cache.n_global_dofs
2191  = std::accumulate (number_cache
2192  .n_locally_owned_dofs_per_processor.begin(),
2193  number_cache
2194  .n_locally_owned_dofs_per_processor.end(),
2195  static_cast<::types::global_dof_index>(0));
2196 
2197  number_cache.locally_owned_dofs = IndexSet(number_cache.n_global_dofs);
2198  number_cache.locally_owned_dofs
2199  .add_range(shift,
2200  shift+number_cache.n_locally_owned_dofs);
2201  number_cache.locally_owned_dofs.compress();
2202 
2203  // fill global_dof_indexsets
2204  number_cache.locally_owned_dofs_per_processor.resize(n_cpus);
2205  {
2206  ::types::global_dof_index lshift = 0;
2207  for (unsigned int i=0; i<n_cpus; ++i)
2208  {
2209  number_cache.locally_owned_dofs_per_processor[i]
2210  = IndexSet(number_cache.n_global_dofs);
2211  number_cache.locally_owned_dofs_per_processor[i]
2212  .add_range(lshift,
2213  lshift +
2214  number_cache.n_locally_owned_dofs_per_processor[i]);
2215  lshift += number_cache.n_locally_owned_dofs_per_processor[i];
2216  }
2217  }
2218  Assert(number_cache.locally_owned_dofs_per_processor
2219  [tr->locally_owned_subdomain()].n_elements()
2220  ==
2221  number_cache.n_locally_owned_dofs,
2222  ExcInternalError());
2223  Assert(!number_cache.locally_owned_dofs_per_processor
2224  [tr->locally_owned_subdomain()].n_elements()
2225  ||
2226  number_cache.locally_owned_dofs_per_processor
2227  [tr->locally_owned_subdomain()].nth_index_in_set(0)
2228  == shift,
2229  ExcInternalError());
2230 
2231  //* 4. send dofids of cells that are
2232  //ghostcells on other machines
2233 
2234  std::vector<bool> user_flags;
2235  tr->save_user_flags(user_flags);
2236  tr->clear_user_flags ();
2237 
2238  //mark all own cells for transfer
2239  for (typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof_handler.begin_active();
2240  cell != dof_handler.end(); ++cell)
2241  if (!cell->is_artificial())
2242  cell->set_user_flag();
2243 
2244  // add each ghostcells'
2245  // subdomain to the vertex and
2246  // keep track of interesting
2247  // neighbors
2248  std::map<unsigned int, std::set<::types::subdomain_id> >
2249  vertices_with_ghost_neighbors;
2250 
2251  tr->fill_vertices_with_ghost_neighbors (vertices_with_ghost_neighbors);
2252 
2253 
2254  /* Send and receive cells. After this,
2255  only the local cells are marked,
2256  that received new data. This has to
2257  be communicated in a second
2258  communication step. */
2259  communicate_dof_indices_on_marked_cells (dof_handler,
2260  vertices_with_ghost_neighbors,
2261  tr->coarse_cell_to_p4est_tree_permutation,
2262  tr->p4est_tree_to_coarse_cell_permutation);
2263 
2264  communicate_dof_indices_on_marked_cells (dof_handler,
2265  vertices_with_ghost_neighbors,
2266  tr->coarse_cell_to_p4est_tree_permutation,
2267  tr->p4est_tree_to_coarse_cell_permutation);
2268 
2269  tr->load_user_flags(user_flags);
2270 
2271 #ifdef DEBUG
2272  //check that we are really done
2273  {
2274  std::vector<::types::global_dof_index> local_dof_indices;
2275 
2276  for (typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof_handler.begin_active();
2277  cell != dof_handler.end(); ++cell)
2278  if (!cell->is_artificial())
2279  {
2280  local_dof_indices.resize (cell->get_fe().dofs_per_cell);
2281  cell->get_dof_indices (local_dof_indices);
2282  if (local_dof_indices.end() !=
2283  std::find (local_dof_indices.begin(),
2284  local_dof_indices.end(),
2285  DoFHandler<dim,spacedim>::invalid_dof_index))
2286  {
2287  if (cell->is_ghost())
2288  {
2289  Assert(false, ExcMessage ("Not a ghost cell"));
2290  }
2291  else
2292  {
2293  Assert(false, ExcMessage ("Not one of our own cells"));
2294  }
2295  }
2296  }
2297  }
2298 #endif // DEBUG
2299 #endif // DEAL_II_WITH_P4EST
2300 
2301  number_cache_current = number_cache;
2302  }
2303 
2304 
2305 
2306  template <int dim, int spacedim>
2307  void
2308  ParallelDistributed<dim, spacedim>::
2309  distribute_mg_dofs (DoFHandler<dim,spacedim> &dof_handler,
2310  std::vector<NumberCache> &number_caches) const
2311  {
2312 #ifndef DEAL_II_WITH_P4EST
2313  (void)dof_handler;
2314  (void)number_caches;
2315  Assert (false, ExcNotImplemented());
2316 #else
2317 
2318  parallel::distributed::Triangulation< dim, spacedim > *tr
2319  = (dynamic_cast<parallel::distributed::Triangulation<dim,spacedim>*>
2320  (const_cast<::Triangulation< dim, spacedim >*>
2321  (&dof_handler.get_triangulation())));
2322  Assert (tr != 0, ExcInternalError());
2323 
2324  AssertThrow(
2325  (tr->settings & parallel::distributed::Triangulation< dim, spacedim >::construct_multigrid_hierarchy),
2326  ExcMessage("Multigrid DoFs can only be distributed on a parallel Triangulation if the flag construct_multigrid_hierarchy is set in the constructor."));
2327 
2328 
2329  const unsigned int
2330  n_cpus = Utilities::MPI::n_mpi_processes (tr->get_communicator());
2331 
2332  unsigned int n_levels = Utilities::MPI::max(dof_handler.get_triangulation().n_levels(), tr->get_communicator());
2333 
2334  for (unsigned int level = 0; level < n_levels; ++level)
2335  {
2336  NumberCache &number_cache = number_caches[level];
2337 
2338  //* 1. distribute on own
2339  //* subdomain
2340  const unsigned int n_initial_local_dofs =
2341  Implementation::distribute_dofs_on_level(0, tr->locally_owned_subdomain(), dof_handler, level);
2342 
2343  //* 2. iterate over ghostcells and
2344  //kill dofs that are not owned
2345  //by us
2346  std::vector<::types::global_dof_index> renumbering(n_initial_local_dofs);
2347  for (::types::global_dof_index i=0; i<renumbering.size(); ++i)
2348  renumbering[i] = i;
2349 
2350  if (level<tr->n_levels())
2351  {
2352  std::vector<::types::global_dof_index> local_dof_indices;
2353 
2354  typename DoFHandler<dim,spacedim>::level_cell_iterator
2355  cell = dof_handler.begin(level),
2356  endc = dof_handler.end(level);
2357 
2358  for (; cell != endc; ++cell)
2359  if (cell->level_subdomain_id()!=numbers::artificial_subdomain_id &&
2360  (cell->level_subdomain_id() < tr->locally_owned_subdomain()))
2361  {
2362  // we found a
2363  // neighboring ghost
2364  // cell whose subdomain
2365  // is "stronger" than
2366  // our own subdomain
2367 
2368  // delete all dofs that
2369  // live there and that
2370  // we have previously
2371  // assigned a number to
2372  // (i.e. the ones on
2373  // the interface)
2374  local_dof_indices.resize (cell->get_fe().dofs_per_cell);
2375  cell->get_mg_dof_indices (local_dof_indices);
2376  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
2377  if (local_dof_indices[i] != DoFHandler<dim,spacedim>::invalid_dof_index)
2378  renumbering[local_dof_indices[i]]
2379  = DoFHandler<dim,spacedim>::invalid_dof_index;
2380  }
2381  }
2382 
2383 
2384  // make indices consecutive
2385  number_cache.n_locally_owned_dofs = 0;
2386  for (std::vector<::types::global_dof_index>::iterator it=renumbering.begin();
2387  it!=renumbering.end(); ++it)
2388  if (*it != DoFHandler<dim,spacedim>::invalid_dof_index)
2389  *it = number_cache.n_locally_owned_dofs++;
2390 
2391  //* 3. communicate local dofcount and
2392  //shift ids to make them unique
2393  number_cache.n_locally_owned_dofs_per_processor.resize(n_cpus);
2394 
2395  int ierr = MPI_Allgather ( &number_cache.n_locally_owned_dofs,
2396  1, DEAL_II_DOF_INDEX_MPI_TYPE,
2397  &number_cache.n_locally_owned_dofs_per_processor[0],
2398  1, DEAL_II_DOF_INDEX_MPI_TYPE,
2399  tr->get_communicator());
2400  AssertThrowMPI(ierr);
2401 
2402  const ::types::global_dof_index
2403  shift = std::accumulate (number_cache
2404  .n_locally_owned_dofs_per_processor.begin(),
2405  number_cache
2406  .n_locally_owned_dofs_per_processor.begin()
2407  + tr->locally_owned_subdomain(),
2408  static_cast<::types::global_dof_index>(0));
2409  for (std::vector<::types::global_dof_index>::iterator it=renumbering.begin();
2410  it!=renumbering.end(); ++it)
2411  if (*it != DoFHandler<dim,spacedim>::invalid_dof_index)
2412  (*it) += shift;
2413 
2414  // now re-enumerate all dofs to
2415  // this shifted and condensed
2416  // numbering form. we renumber
2417  // some dofs as invalid, so
2418  // choose the nocheck-version.
2419  Implementation::renumber_mg_dofs (renumbering, IndexSet(0),
2420  dof_handler, level, false);
2421 
2422  // now a little bit of
2423  // housekeeping
2424  number_cache.n_global_dofs
2425  = std::accumulate (number_cache
2426  .n_locally_owned_dofs_per_processor.begin(),
2427  number_cache
2428  .n_locally_owned_dofs_per_processor.end(),
2429  static_cast<::types::global_dof_index>(0));
2430 
2431  number_cache.locally_owned_dofs = IndexSet(number_cache.n_global_dofs);
2432  number_cache.locally_owned_dofs
2433  .add_range(shift,
2434  shift+number_cache.n_locally_owned_dofs);
2435  number_cache.locally_owned_dofs.compress();
2436 
2437  // fill global_dof_indexsets
2438  number_cache.locally_owned_dofs_per_processor.resize(n_cpus);
2439  {
2440  ::types::global_dof_index lshift = 0;
2441  for (unsigned int i=0; i<n_cpus; ++i)
2442  {
2443  number_cache.locally_owned_dofs_per_processor[i]
2444  = IndexSet(number_cache.n_global_dofs);
2445  number_cache.locally_owned_dofs_per_processor[i]
2446  .add_range(lshift,
2447  lshift +
2448  number_cache.n_locally_owned_dofs_per_processor[i]);
2449  lshift += number_cache.n_locally_owned_dofs_per_processor[i];
2450  }
2451  }
2452  Assert(number_cache.locally_owned_dofs_per_processor
2453  [tr->locally_owned_subdomain()].n_elements()
2454  ==
2455  number_cache.n_locally_owned_dofs,
2456  ExcInternalError());
2457  Assert(!number_cache.locally_owned_dofs_per_processor
2458  [tr->locally_owned_subdomain()].n_elements()
2459  ||
2460  number_cache.locally_owned_dofs_per_processor
2461  [tr->locally_owned_subdomain()].nth_index_in_set(0)
2462  == shift,
2463  ExcInternalError());
2464 
2465  }
2466 
2467 
2468  //* communicate ghost DoFs
2469  // We mark all ghost cells by setting the user_flag and then request
2470  // these cells from the corresponding owners. As this information
2471  // can be incomplete,
2472  {
2473  std::vector<bool> user_flags;
2474  tr->save_user_flags(user_flags);
2475  tr->clear_user_flags ();
2476 
2477  // mark all ghost cells for transfer
2478  {
2479  typename DoFHandler<dim,spacedim>::level_cell_iterator
2480  cell, endc = dof_handler.end();
2481  for (cell = dof_handler.begin(); cell != endc; ++cell)
2482  if (cell->level_subdomain_id() != ::numbers::artificial_subdomain_id
2483  && !cell->is_locally_owned_on_level())
2484  cell->set_user_flag();
2485  }
2486 
2487  // Phase 1. Request all marked cells from corresponding owners. If we
2488  // managed to get every DoF, remove the user_flag, otherwise we
2489  // will request them again in the step below.
2490  communicate_mg_ghost_cells(*tr,
2491  dof_handler,
2492  tr->coarse_cell_to_p4est_tree_permutation,
2493  tr->p4est_tree_to_coarse_cell_permutation);
2494 
2495  // This barrier is crucial so that messages between phase 1&2 don't
2496  // mix.
2497  const int ierr = MPI_Barrier(tr->get_communicator());
2498  AssertThrowMPI(ierr);
2499 
2500  // Phase 2, only request the cells that were not completed in Phase
2501  // 1.
2502  communicate_mg_ghost_cells(*tr,
2503  dof_handler,
2504  tr->coarse_cell_to_p4est_tree_permutation,
2505  tr->p4est_tree_to_coarse_cell_permutation);
2506 
2507 #ifdef DEBUG
2508  // make sure we have removed all flags:
2509  {
2510  typename DoFHandler<dim,spacedim>::level_cell_iterator
2511  cell, endc = dof_handler.end();
2512  for (cell = dof_handler.begin(); cell != endc; ++cell)
2513  if (cell->level_subdomain_id() != ::numbers::artificial_subdomain_id
2514  && !cell->is_locally_owned_on_level())
2515  Assert(cell->user_flag_set()==false, ExcInternalError());
2516  }
2517 #endif
2518 
2519  tr->load_user_flags(user_flags);
2520  }
2521 
2522 
2523 
2524 
2525 #ifdef DEBUG
2526  //check that we are really done
2527  {
2528  std::vector<::types::global_dof_index> local_dof_indices;
2529  typename DoFHandler<dim,spacedim>::level_cell_iterator
2530  cell, endc = dof_handler.end();
2531 
2532  for (cell = dof_handler.begin(); cell != endc; ++cell)
2533  if (cell->level_subdomain_id() != ::numbers::artificial_subdomain_id)
2534  {
2535  local_dof_indices.resize (cell->get_fe().dofs_per_cell);
2536  cell->get_mg_dof_indices (local_dof_indices);
2537  if (local_dof_indices.end() !=
2538  std::find (local_dof_indices.begin(),
2539  local_dof_indices.end(),
2540  DoFHandler<dim,spacedim>::invalid_dof_index))
2541  {
2542  Assert(false, ExcMessage ("not all DoFs got distributed!"));
2543  }
2544  }
2545  }
2546 #endif // DEBUG
2547 
2548 #endif // DEAL_II_WITH_P4EST
2549  }
2550 
2551 
2552  template <int dim, int spacedim>
2553  void
2554  ParallelDistributed<dim, spacedim>::
2555  renumber_dofs (const std::vector<::types::global_dof_index> &new_numbers,
2556  ::DoFHandler<dim,spacedim> &dof_handler,
2557  NumberCache &number_cache_current) const
2558  {
2559  (void)new_numbers;
2560  (void)dof_handler;
2561 
2562  Assert (new_numbers.size() == dof_handler.locally_owned_dofs().n_elements(),
2563  ExcInternalError());
2564 
2565  NumberCache number_cache;
2566 
2567 #ifndef DEAL_II_WITH_P4EST
2568  Assert (false, ExcNotImplemented());
2569 #else
2570 
2571 
2572  // calculate new IndexSet. First try to find out if the new indices
2573  // are contiguous blocks. This avoids inserting each index
2574  // individually into the IndexSet, which is slow. If we own no DoFs,
2575  // we still need to go through this function, but we can skip this
2576  // calculation.
2577 
2578  number_cache.locally_owned_dofs = IndexSet (dof_handler.n_dofs());
2579  if (dof_handler.locally_owned_dofs().n_elements()>0)
2580  {
2581  std::vector<::types::global_dof_index> new_numbers_sorted (new_numbers);
2582  std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
2583  std::vector<::types::global_dof_index>::const_iterator it = new_numbers_sorted.begin();
2584  const unsigned int n_blocks = dof_handler.get_fe().n_blocks();
2585  std::vector<std::pair<::types::global_dof_index,unsigned int> > block_indices(n_blocks);
2586  block_indices[0].first = *it++;
2587  block_indices[0].second = 1;
2588  unsigned int current_block = 0, n_filled_blocks = 1;
2589  for ( ; it != new_numbers_sorted.end(); ++it)
2590  {
2591  bool done = false;
2592 
2593  // search from the current block onwards whether the next
2594  // index is shifted by one from the previous one.
2595  for (unsigned int i=0; i<n_filled_blocks; ++i)
2596  if (*it == block_indices[current_block].first
2597  +block_indices[current_block].second)
2598  {
2599  block_indices[current_block].second++;
2600  done = true;
2601  break;
2602  }
2603  else
2604  {
2605  if (current_block == n_filled_blocks-1)
2606  current_block = 0;
2607  else
2608  ++current_block;
2609  }
2610 
2611  // could not find any contiguous range: need to add a new
2612  // block if possible. Abort otherwise, which will add all
2613  // elements individually to the IndexSet.
2614  if (done == false)
2615  {
2616  if (n_filled_blocks < n_blocks)
2617  {
2618  block_indices[n_filled_blocks].first = *it;
2619  block_indices[n_filled_blocks].second = 1;
2620  current_block = n_filled_blocks;
2621  ++n_filled_blocks;
2622  }
2623  else
2624  break;
2625  }
2626  }
2627 
2628  // check whether all indices could be assigned to blocks. If yes,
2629  // we can add the block ranges to the IndexSet, otherwise we need
2630  // to go through the indices once again and add each element
2631  // individually
2632  unsigned int sum = 0;
2633  for (unsigned int i=0; i<n_filled_blocks; ++i)
2634  sum += block_indices[i].second;
2635  if (sum == new_numbers.size())
2636  for (unsigned int i=0; i<n_filled_blocks; ++i)
2637  number_cache.locally_owned_dofs.add_range (block_indices[i].first,
2638  block_indices[i].first+
2639  block_indices[i].second);
2640  else
2641  number_cache.locally_owned_dofs.add_indices(new_numbers_sorted.begin(),
2642  new_numbers_sorted.end());
2643  }
2644 
2645 
2646  number_cache.locally_owned_dofs.compress();
2647  Assert (number_cache.locally_owned_dofs.n_elements() == new_numbers.size(),
2648  ExcInternalError());
2649  // also check with the number of locally owned degrees of freedom that
2650  // the DoFHandler object still stores
2651  Assert (number_cache.locally_owned_dofs.n_elements() ==
2652  dof_handler.n_locally_owned_dofs(),
2653  ExcInternalError());
2654 
2655  // then also set this number in our own copy
2656  number_cache.n_locally_owned_dofs = dof_handler.n_locally_owned_dofs();
2657 
2658  // mark not locally active DoFs as invalid
2659  {
2660  std::vector<::types::global_dof_index> local_dof_indices;
2661 
2662  typename DoFHandler<dim,spacedim>::active_cell_iterator
2663  cell = dof_handler.begin_active(),
2664  endc = dof_handler.end();
2665 
2666  for (; cell != endc; ++cell)
2667  if (!cell->is_artificial())
2668  {
2669  local_dof_indices.resize (cell->get_fe().dofs_per_cell);
2670  cell->get_dof_indices (local_dof_indices);
2671  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
2672  {
2673  if (local_dof_indices[i] == DoFHandler<dim,spacedim>::invalid_dof_index)
2674  continue;
2675 
2676  if (!dof_handler.locally_owned_dofs().is_element(local_dof_indices[i]))
2677  {
2678  //this DoF is not owned by us, so set it to invalid.
2679  local_dof_indices[i]
2680  = DoFHandler<dim,spacedim>::invalid_dof_index;
2681  }
2682  }
2683 
2684  cell->set_dof_indices (local_dof_indices);
2685  }
2686  }
2687 
2688 
2689  // renumber. Skip when there is nothing to do because we own no DoF.
2690  if (dof_handler.locally_owned_dofs().n_elements() > 0)
2691  Implementation::renumber_dofs (new_numbers,
2692  dof_handler.locally_owned_dofs(),
2693  dof_handler,
2694  false);
2695 
2696  // communication
2697  {
2698  parallel::distributed::Triangulation< dim, spacedim > *tr
2699  = (dynamic_cast<parallel::distributed::Triangulation<dim,spacedim>*>
2700  (const_cast<::Triangulation< dim, spacedim >*>
2701  (&dof_handler.get_triangulation())));
2702  Assert (tr != 0, ExcInternalError());
2703 
2704  std::vector<bool> user_flags;
2705  tr->save_user_flags(user_flags);
2706  tr->clear_user_flags ();
2707 
2708  //mark all own cells for transfer
2709  typename DoFHandler<dim,spacedim>::active_cell_iterator
2710  cell, endc = dof_handler.end();
2711  for (cell = dof_handler.begin_active(); cell != endc; ++cell)
2712  if (!cell->is_artificial())
2713  cell->set_user_flag();
2714 
2715  // add each ghostcells' subdomain to the vertex and keep track of
2716  // interesting neighbors
2717  std::map<unsigned int, std::set<::types::subdomain_id> >
2718  vertices_with_ghost_neighbors;
2719 
2720  tr->fill_vertices_with_ghost_neighbors (vertices_with_ghost_neighbors);
2721 
2722  // Send and receive cells. After this, only the local cells are
2723  // marked, that received new data. This has to be communicated in a
2724  // second communication step.
2725  communicate_dof_indices_on_marked_cells (dof_handler,
2726  vertices_with_ghost_neighbors,
2727  tr->coarse_cell_to_p4est_tree_permutation,
2728  tr->p4est_tree_to_coarse_cell_permutation);
2729 
2730  communicate_dof_indices_on_marked_cells (dof_handler,
2731  vertices_with_ghost_neighbors,
2732  tr->coarse_cell_to_p4est_tree_permutation,
2733  tr->p4est_tree_to_coarse_cell_permutation);
2734 
2735 
2736  // * Create global_dof_indexsets by transferring our own owned_dofs
2737  // to every other machine.
2738  const unsigned int n_cpus = Utilities::MPI::n_mpi_processes (tr->get_communicator());
2739 
2740  // Serialize our IndexSet and determine size.
2741  std::ostringstream oss;
2742  number_cache.locally_owned_dofs.block_write(oss);
2743  std::string oss_str=oss.str();
2744  std::vector<char> my_data(oss_str.begin(), oss_str.end());
2745  unsigned int my_size = oss_str.size();
2746 
2747  // determine maximum size of IndexSet
2748  const unsigned int max_size
2749  = Utilities::MPI::max (my_size, tr->get_communicator());
2750 
2751  // as we are reading past the end, we need to increase the size of
2752  // the local buffer. This is filled with zeros.
2753  my_data.resize(max_size);
2754 
2755  std::vector<char> buffer(max_size*n_cpus);
2756  const int ierr = MPI_Allgather(&my_data[0], max_size, MPI_BYTE,
2757  &buffer[0], max_size, MPI_BYTE,
2758  tr->get_communicator());
2759  AssertThrowMPI(ierr);
2760 
2761  number_cache.locally_owned_dofs_per_processor.resize (n_cpus);
2762  number_cache.n_locally_owned_dofs_per_processor.resize (n_cpus);
2763  for (unsigned int i=0; i<n_cpus; ++i)
2764  {
2765  std::stringstream strstr;
2766  strstr.write(&buffer[i*max_size],max_size);
2767  // This does not read the whole buffer, when the size is smaller
2768  // than max_size. Therefore we need to create a new stringstream
2769  // in each iteration (resetting would be fine too).
2770  number_cache.locally_owned_dofs_per_processor[i]
2771  .block_read(strstr);
2772  number_cache.n_locally_owned_dofs_per_processor[i]
2773  = number_cache.locally_owned_dofs_per_processor[i].n_elements();
2774  }
2775 
2776  number_cache.n_global_dofs
2777  = std::accumulate (number_cache
2778  .n_locally_owned_dofs_per_processor.begin(),
2779  number_cache
2780  .n_locally_owned_dofs_per_processor.end(),
2781  static_cast<::types::global_dof_index>(0));
2782 
2783  tr->load_user_flags(user_flags);
2784  }
2785 #endif
2786 
2787  number_cache_current = number_cache;
2788  }
2789  }
2790  }
2791 }
2792 
2793 
2794 
2795 
2796 /*-------------- Explicit Instantiations -------------------------------*/
2797 #include "dof_handler_policy.inst"
2798 
2799 
2800 DEAL_II_NAMESPACE_CLOSE
Triangulation::n_active_cells
unsigned int n_active_cells() const
Definition: tria.cc:11244
DoFHandler::mg_vertex_dofs
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1067
DoFTools::locally_owned_dofs_per_subdomain
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandlerType &dof_handler)
Definition: dof_tools.cc:1166
parallel::distributed::Triangulation::fill_vertices_with_ghost_neighbors
void fill_vertices_with_ghost_neighbors(std::map< unsigned int, std::set<::types::subdomain_id > > &vertices_with_ghost_neighbors)
Definition: tria.cc:3391
numbers::invalid_subdomain_id
const types::subdomain_id invalid_subdomain_id
Definition: types.h:245
internal::DoFHandler::Policy::Sequential::distribute_dofs
virtual void distribute_dofs(::DoFHandler< dim, spacedim > &dof_handler, NumberCache &number_cache) const
Definition: dof_handler_policy.cc:946
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
IndexSet::nth_index_in_set
size_type nth_index_in_set(const unsigned int local_index) const
Definition: index_set.h:1612
parallel::shared::Triangulation::get_true_subdomain_ids_of_cells
const std::vector< types::subdomain_id > & get_true_subdomain_ids_of_cells() const
Definition: shared_tria.cc:102
internal::DoFHandler::Policy::ParallelDistributed::renumber_dofs
virtual void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers, ::DoFHandler< dim, spacedim > &dof_handler, NumberCache &number_cache) const
Definition: dof_handler_policy.cc:2555
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
internal::DoFHandler::Policy::ParallelShared::distribute_dofs
virtual void distribute_dofs(::DoFHandler< dim, spacedim > &dof_handler, NumberCache &number_cache) const
Definition: dof_handler_policy.cc:1044
DoFHandler::selected_fe
SmartPointer< const FiniteElement< dim, spacedim >, DoFHandler< dim, spacedim > > selected_fe
Definition: dof_handler.h:930
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
Triangulation::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:10668
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:369
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe() const
internal::DoFHandler::NumberCache::n_locally_owned_dofs
types::global_dof_index n_locally_owned_dofs
Definition: number_cache.h:63
DoFHandler::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:740
DoFHandler::n_locally_owned_dofs
unsigned int n_locally_owned_dofs() const
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
DoFHandler::cell_iterator
ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:256
types
Definition: types.h:30
types::global_dof_index
unsigned int global_dof_index
Definition: types.h:88
Assert
#define Assert(cond, exc)
Definition: exceptions.h:313
internal::DoFHandler::Policy::Sequential::renumber_dofs
virtual void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers, ::DoFHandler< dim, spacedim > &dof_handler, NumberCache &number_cache) const
Definition: dof_handler_policy.cc:1006
IndexSet::index_within_set
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1653
internal::DoFHandler::Policy::ParallelShared::distribute_mg_dofs
virtual void distribute_mg_dofs(::DoFHandler< dim, spacedim > &dof_handler, std::vector< NumberCache > &number_caches) const
Definition: dof_handler_policy.cc:1091
parallel::Triangulation::get_communicator
virtual MPI_Comm get_communicator() const
Definition: tria_base.cc:146
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
DoFHandler::faces
::internal::DoFHandler::DoFFaces< dim > * faces
Definition: dof_handler.h:1082
types::subdomain_id
unsigned int subdomain_id
Definition: types.h:42
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::DoFHandler::Policy::ParallelDistributed::distribute_dofs
virtual void distribute_dofs(::DoFHandler< dim, spacedim > &dof_handler, NumberCache &number_cache) const
Definition: dof_handler_policy.cc:2085
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:63
parallel::Triangulation::level_ghost_owners
const std::set< types::subdomain_id > & level_ghost_owners() const
Definition: tria_base.cc:246
IndexSet::add_range
void add_range(const size_type begin, const size_type end)
Definition: index_set.cc:88
numbers::artificial_subdomain_id
const types::subdomain_id artificial_subdomain_id
Definition: types.h:261
IndexSet::compress
void compress() const
Definition: index_set.h:1428
internal::DoFHandler::Policy::Sequential::distribute_mg_dofs
virtual void distribute_mg_dofs(::DoFHandler< dim, spacedim > &dof_handler, std::vector< NumberCache > &number_caches) const
Definition: dof_handler_policy.cc:980
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1211
internal::DoFHandler::NumberCache::n_locally_owned_dofs_per_processor
std::vector< types::global_dof_index > n_locally_owned_dofs_per_processor
Definition: number_cache.h:77
internal::DoFHandler::Policy::ParallelDistributed::distribute_mg_dofs
virtual void distribute_mg_dofs(::DoFHandler< dim, spacedim > &dof_handler, std::vector< NumberCache > &number_caches) const
Definition: dof_handler_policy.cc:2309
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:73
Triangulation::save_user_flags
void save_user_flags(std::ostream &out) const
Definition: tria.cc:9899
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
IteratorState::valid
Iterator points to a valid object.
Definition: tria_iterator_base.h:39
internal::DoFHandler::NumberCache::locally_owned_dofs_per_processor
std::vector< IndexSet > locally_owned_dofs_per_processor
Definition: number_cache.h:84
IndexSet::is_element
bool is_element(const size_type index) const
Definition: index_set.h:1489
internal::DoFHandler::NumberCache::n_global_dofs
types::global_dof_index n_global_dofs
Definition: number_cache.h:57
DoFRenumbering::subdomain_wise
void subdomain_wise(DoFHandlerType &dof_handler)
Definition: dof_renumbering.cc:1892
parallel::Triangulation::locally_owned_subdomain
types::subdomain_id locally_owned_subdomain() const
Definition: tria_base.cc:230
DoFHandler::locally_owned_dofs
const IndexSet & locally_owned_dofs() const
IndexSet::n_elements
size_type n_elements() const
Definition: index_set.h:1560
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
internal::DoFHandler::Policy::ParallelShared::renumber_dofs
virtual void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers, ::DoFHandler< dim, spacedim > &dof_handler, NumberCache &number_cache) const
Definition: dof_handler_policy.cc:1104
DoFHandler::vertex_dofs
std::vector< types::global_dof_index > vertex_dofs
Definition: dof_handler.h:1061
parallel::distributed::Triangulation::coarse_cell_to_p4est_tree_permutation
std::vector< types::global_dof_index > coarse_cell_to_p4est_tree_permutation
Definition: tria.h:807
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