Reference documentation for deal.II version GIT c9976103bc 2022-12-09 17:30:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria_base.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/logstream.h>
19 #include <deal.II/base/mpi.templates.h>
21 #include <deal.II/base/utilities.h>
22 
25 
28 #include <deal.II/grid/tria.h>
31 
33 
34 #include <algorithm>
35 #include <cstdint>
36 #include <fstream>
37 #include <iostream>
38 #include <limits>
39 #include <numeric>
40 
41 
43 
44 namespace parallel
45 {
46  template <int dim, int spacedim>
48  const MPI_Comm &mpi_communicator,
49  const typename ::Triangulation<dim, spacedim>::MeshSmoothing
50  smooth_grid,
51  const bool check_for_distorted_cells)
52  : ::Triangulation<dim, spacedim>(smooth_grid,
53  check_for_distorted_cells)
54  , mpi_communicator(mpi_communicator)
55  , my_subdomain(Utilities::MPI::this_mpi_process(this->mpi_communicator))
56  , n_subdomains(Utilities::MPI::n_mpi_processes(this->mpi_communicator))
57  {
58 #ifndef DEAL_II_WITH_MPI
59  Assert(false, ExcNeedsMPI());
60 #endif
61  }
62 
63 
64 
65  template <int dim, int spacedim>
66  void
68  const ::Triangulation<dim, spacedim> &other_tria)
69  {
70 #ifndef DEAL_II_WITH_MPI
71  (void)other_tria;
72  Assert(false, ExcNeedsMPI());
73 #else
75 
76  if (const ::parallel::TriangulationBase<dim, spacedim> *other_tria_x =
77  dynamic_cast<const ::parallel::TriangulationBase<dim, spacedim>
78  *>(&other_tria))
79  {
80  // release unused vector memory because we will have very different
81  // vectors now
84  }
85 #endif
86  }
87 
88 
89 
90  template <int dim, int spacedim>
91  std::size_t
93  {
94  std::size_t mem =
96  MemoryConsumption::memory_consumption(this->mpi_communicator) +
99  number_cache.n_global_active_cells) +
100  MemoryConsumption::memory_consumption(number_cache.n_global_levels);
101  return mem;
102  }
103 
104  template <int dim, int spacedim>
106  {
107  // release unused vector memory because the vector layout is going to look
108  // very different now
111  }
112 
113  template <int dim, int spacedim>
115  : n_locally_owned_active_cells(0)
116  , number_of_global_coarse_cells(0)
117  , n_global_levels(0)
118  {}
119 
120  template <int dim, int spacedim>
121  unsigned int
123  {
124  return number_cache.n_locally_owned_active_cells;
125  }
126 
127  template <int dim, int spacedim>
128  unsigned int
130  {
131  return number_cache.n_global_levels;
132  }
133 
134  template <int dim, int spacedim>
137  {
138  return number_cache.n_global_active_cells;
139  }
140 
141  template <int dim, int spacedim>
142  MPI_Comm
144  {
145  return mpi_communicator;
146  }
147 
148 #ifdef DEAL_II_WITH_MPI
149  template <int dim, int spacedim>
150  void
152  {
153  number_cache.ghost_owners.clear();
154  number_cache.level_ghost_owners.clear();
155  number_cache.n_locally_owned_active_cells = 0;
156 
157  if (this->n_levels() == 0)
158  {
159  // Skip communication done below if we do not have any cells
160  // (meaning the Triangulation is empty on all processors). This will
161  // happen when called from the destructor of Triangulation, which
162  // can get called during exception handling causing a hang in this
163  // function.
164  number_cache.n_global_active_cells = 0;
165  number_cache.n_global_levels = 0;
166  return;
167  }
168 
169 
170  {
171  // find ghost owners
172  for (const auto &cell : this->active_cell_iterators())
173  if (cell->is_ghost())
174  number_cache.ghost_owners.insert(cell->subdomain_id());
175 
176  Assert(number_cache.ghost_owners.size() <
177  Utilities::MPI::n_mpi_processes(this->mpi_communicator),
178  ExcInternalError());
179  }
180 
181  if (this->n_levels() > 0)
182  number_cache.n_locally_owned_active_cells = std::count_if(
183  this->begin_active(),
185  this->end()),
186  [](const auto &i) { return i.is_locally_owned(); });
187  else
188  number_cache.n_locally_owned_active_cells = 0;
189 
190  // Potentially cast to a 64 bit type before accumulating to avoid
191  // overflow:
192  number_cache.n_global_active_cells =
194  number_cache.n_locally_owned_active_cells),
195  this->mpi_communicator);
196 
197  number_cache.n_global_levels =
198  Utilities::MPI::max(this->n_levels(), this->mpi_communicator);
199 
200  // Store MPI ranks of level ghost owners of this processor on all
201  // levels.
202  if (this->is_multilevel_hierarchy_constructed() == true)
203  {
204  number_cache.level_ghost_owners.clear();
205 
206  // if there is nothing to do, then do nothing
207  if (this->n_levels() == 0)
208  return;
209 
210  // find level ghost owners
211  for (const auto &cell : this->cell_iterators())
212  if (cell->level_subdomain_id() != numbers::artificial_subdomain_id &&
213  cell->level_subdomain_id() != this->locally_owned_subdomain())
214  this->number_cache.level_ghost_owners.insert(
215  cell->level_subdomain_id());
216 
217 # ifdef DEBUG
218  // Check that level_ghost_owners is symmetric by sending a message
219  // to everyone
220  {
221  int ierr = MPI_Barrier(this->mpi_communicator);
222  AssertThrowMPI(ierr);
223 
224  const int mpi_tag = Utilities::MPI::internal::Tags::
226 
227  // important: preallocate to avoid (re)allocation:
228  std::vector<MPI_Request> requests(
229  this->number_cache.level_ghost_owners.size());
230  unsigned int dummy = 0;
231  unsigned int req_counter = 0;
232 
233  for (const auto &it : this->number_cache.level_ghost_owners)
234  {
235  ierr = MPI_Isend(&dummy,
236  1,
237  MPI_UNSIGNED,
238  it,
239  mpi_tag,
240  this->mpi_communicator,
241  &requests[req_counter]);
242  AssertThrowMPI(ierr);
243  ++req_counter;
244  }
245 
246  for (const auto &it : this->number_cache.level_ghost_owners)
247  {
248  unsigned int dummy;
249  ierr = MPI_Recv(&dummy,
250  1,
251  MPI_UNSIGNED,
252  it,
253  mpi_tag,
254  this->mpi_communicator,
255  MPI_STATUS_IGNORE);
256  AssertThrowMPI(ierr);
257  }
258 
259  if (requests.size() > 0)
260  {
261  ierr = MPI_Waitall(requests.size(),
262  requests.data(),
263  MPI_STATUSES_IGNORE);
264  AssertThrowMPI(ierr);
265  }
266 
267  ierr = MPI_Barrier(this->mpi_communicator);
268  AssertThrowMPI(ierr);
269  }
270 # endif
271 
272  Assert(this->number_cache.level_ghost_owners.size() <
273  Utilities::MPI::n_mpi_processes(this->mpi_communicator),
274  ExcInternalError());
275  }
276 
277  this->number_cache.number_of_global_coarse_cells = this->n_cells(0);
278 
279  // reset global cell ids
280  this->reset_global_cell_indices();
281  }
282 
283 #else
284 
285  template <int dim, int spacedim>
286  void
288  {
289  Assert(false, ExcNeedsMPI());
290  }
291 
292 #endif
293 
294  template <int dim, int spacedim>
295  void
297  {
298  // run algorithm for locally-owned cells
300 
301  // translate ReferenceCell to unsigned int (needed by
302  // Utilities::MPI::compute_set_union)
303  std::vector<unsigned int> reference_cells_ui;
304 
305  for (const auto &i : this->reference_cells)
306  reference_cells_ui.push_back(static_cast<unsigned int>(i));
307 
308  // create union
309  reference_cells_ui =
310  Utilities::MPI::compute_set_union(reference_cells_ui,
311  this->mpi_communicator);
312 
313  // transform back and store result
314  this->reference_cells.clear();
315  for (const auto &i : reference_cells_ui)
316  this->reference_cells.emplace_back(
318  }
319 
320 
321 
322  template <int dim, int spacedim>
325  {
326  return my_subdomain;
327  }
328 
329 
330 
331  template <int dim, int spacedim>
332  const std::set<types::subdomain_id> &
334  {
335  return number_cache.ghost_owners;
336  }
337 
338 
339 
340  template <int dim, int spacedim>
341  const std::set<types::subdomain_id> &
343  {
344  return number_cache.level_ghost_owners;
345  }
346 
347 
348 
349  template <int dim, int spacedim>
350  std::vector<types::boundary_id>
352  {
355  this->mpi_communicator);
356  }
357 
358 
359 
360  template <int dim, int spacedim>
361  std::vector<types::manifold_id>
363  {
366  this->mpi_communicator);
367  }
368 
369 
370 
371  template <int dim, int spacedim>
372  void
374  {
375 #ifndef DEAL_II_WITH_MPI
376  Assert(false, ExcNeedsMPI());
377 #else
378  if (const auto pst =
380  this))
381  if (pst->with_artificial_cells() == false)
382  {
383  // Specialization for parallel::shared::Triangulation without
384  // artificial cells. The code below only works if a halo of a single
385  // ghost cells is needed.
386 
387  std::vector<unsigned int> cell_counter(n_subdomains + 1);
388 
389  // count number of cells of each process
390  for (const auto &cell : this->active_cell_iterators())
391  cell_counter[cell->subdomain_id() + 1]++;
392 
393  // take prefix sum to obtain offset of each process
394  for (unsigned int i = 0; i < n_subdomains; ++i)
395  cell_counter[i + 1] += cell_counter[i];
396 
397  AssertDimension(cell_counter.back(), this->n_active_cells());
398 
399  // create partitioners
400  IndexSet is_local(this->n_active_cells());
401  is_local.add_range(cell_counter[my_subdomain],
402  cell_counter[my_subdomain + 1]);
403  number_cache.active_cell_index_partitioner =
404  std::make_shared<const Utilities::MPI::Partitioner>(
405  is_local,
407  this->mpi_communicator);
408 
409  // set global active cell indices and increment process-local counters
410  for (const auto &cell : this->active_cell_iterators())
411  cell->set_global_active_cell_index(
412  cell_counter[cell->subdomain_id()]++);
413 
414  Assert(this->is_multilevel_hierarchy_constructed() == false,
416 
417  return;
418  }
419 
420  // 1) determine number of active locally-owned cells
421  const types::global_cell_index n_locally_owned_cells =
422  this->n_locally_owned_active_cells();
423 
424  // 2) determine the offset of each process
426 
427  const int ierr = MPI_Exscan(
428  &n_locally_owned_cells,
429  &cell_index,
430  1,
431  Utilities::MPI::mpi_type_id_for_type<decltype(n_locally_owned_cells)>,
432  MPI_SUM,
433  this->mpi_communicator);
434  AssertThrowMPI(ierr);
435 
436  // 3) give global indices to locally-owned cells and mark all other cells as
437  // invalid
438  std::pair<types::global_cell_index, types::global_cell_index> my_range;
439  my_range.first = cell_index;
440 
441  for (const auto &cell : this->active_cell_iterators())
442  if (cell->is_locally_owned())
443  cell->set_global_active_cell_index(cell_index++);
444  else
445  cell->set_global_active_cell_index(numbers::invalid_dof_index);
446 
447  my_range.second = cell_index;
448 
449  // 4) determine the global indices of ghost cells
450  std::vector<types::global_dof_index> is_ghost_vector;
451  GridTools::exchange_cell_data_to_ghosts<types::global_cell_index>(
452  *this,
453  [](const auto &cell) { return cell->global_active_cell_index(); },
454  [&is_ghost_vector](const auto &cell, const auto &id) {
455  cell->set_global_active_cell_index(id);
456  is_ghost_vector.push_back(id);
457  });
458 
459  // 5) set up new partitioner
460  IndexSet is_local(this->n_global_active_cells());
461  is_local.add_range(my_range.first, my_range.second);
462 
463  std::sort(is_ghost_vector.begin(), is_ghost_vector.end());
464  IndexSet is_ghost(this->n_global_active_cells());
465  is_ghost.add_indices(is_ghost_vector.begin(), is_ghost_vector.end());
466 
467  number_cache.active_cell_index_partitioner =
468  std::make_shared<const Utilities::MPI::Partitioner>(
469  is_local, is_ghost, this->mpi_communicator);
470 
471  // 6) proceed with multigrid levels if requested
472  if (this->is_multilevel_hierarchy_constructed() == true)
473  {
474  // 1) determine number of locally-owned cells on levels
475  std::vector<types::global_cell_index> n_cells_level(
476  this->n_global_levels(), 0);
477 
478  for (auto cell : this->cell_iterators())
479  if (cell->level_subdomain_id() == this->locally_owned_subdomain())
480  n_cells_level[cell->level()]++;
481 
482  // 2) determine the offset of each process
483  std::vector<types::global_cell_index> cell_index(
484  this->n_global_levels(), 0);
485 
486  int ierr = MPI_Exscan(
487  n_cells_level.data(),
488  cell_index.data(),
489  this->n_global_levels(),
490  Utilities::MPI::mpi_type_id_for_type<decltype(*n_cells_level.data())>,
491  MPI_SUM,
492  this->mpi_communicator);
493  AssertThrowMPI(ierr);
494 
495  // 3) determine global number of "active" cells on each level
496  Utilities::MPI::sum(n_cells_level,
497  this->mpi_communicator,
498  n_cells_level);
499 
500  // 4) give global indices to locally-owned cells on level and mark
501  // all other cells as invalid
502  std::vector<
503  std::pair<types::global_cell_index, types::global_cell_index>>
504  my_ranges(this->n_global_levels());
505  for (unsigned int l = 0; l < this->n_global_levels(); ++l)
506  my_ranges[l].first = cell_index[l];
507 
508  for (auto cell : this->cell_iterators())
509  if (cell->level_subdomain_id() == this->locally_owned_subdomain())
510  cell->set_global_level_cell_index(cell_index[cell->level()]++);
511  else
512  cell->set_global_level_cell_index(numbers::invalid_dof_index);
513 
514  for (unsigned int l = 0; l < this->n_global_levels(); ++l)
515  my_ranges[l].second = cell_index[l];
516 
517  // 5) update the numbers of ghost level cells
518  std::vector<std::vector<types::global_dof_index>> is_ghost_vectors(
519  this->n_global_levels());
523  *this,
524  [](const auto &cell) { return cell->global_level_cell_index(); },
525  [&is_ghost_vectors](const auto &cell, const auto &id) {
526  cell->set_global_level_cell_index(id);
527  is_ghost_vectors[cell->level()].push_back(id);
528  });
529 
530  number_cache.level_cell_index_partitioners.resize(
531  this->n_global_levels());
532 
533  // 6) set up cell partitioners for each level
534  for (unsigned int l = 0; l < this->n_global_levels(); ++l)
535  {
536  IndexSet is_local(n_cells_level[l]);
537  is_local.add_range(my_ranges[l].first, my_ranges[l].second);
538 
539  IndexSet is_ghost(n_cells_level[l]);
540  std::sort(is_ghost_vectors[l].begin(), is_ghost_vectors[l].end());
541  is_ghost.add_indices(is_ghost_vectors[l].begin(),
542  is_ghost_vectors[l].end());
543 
544  number_cache.level_cell_index_partitioners[l] =
545  std::make_shared<const Utilities::MPI::Partitioner>(
546  is_local, is_ghost, this->mpi_communicator);
547  }
548  }
549 
550 #endif
551  }
552 
553 
554 
555  template <int dim, int spacedim>
556  void
558  const std::vector<bool> &vertex_locally_moved)
559  {
560  AssertDimension(vertex_locally_moved.size(), this->n_vertices());
561 #ifdef DEBUG
562  {
563  const std::vector<bool> locally_owned_vertices =
565  for (unsigned int i = 0; i < locally_owned_vertices.size(); ++i)
566  Assert((vertex_locally_moved[i] == false) ||
567  (locally_owned_vertices[i] == true),
568  ExcMessage("The vertex_locally_moved argument must not "
569  "contain vertices that are not locally owned"));
570  }
571 #endif
572 
573  Point<spacedim> invalid_point;
574  for (unsigned int d = 0; d < spacedim; ++d)
575  invalid_point[d] = std::numeric_limits<double>::quiet_NaN();
576 
577  const auto pack = [&](const auto &cell) {
578  std::vector<Point<spacedim>> vertices(cell->n_vertices());
579 
580  for (const auto v : cell->vertex_indices())
581  if (vertex_locally_moved[cell->vertex_index(v)])
582  vertices[v] = cell->vertex(v);
583  else
584  vertices[v] = invalid_point;
585 
586  return vertices;
587  };
588 
589  const auto unpack = [&](const auto &cell, const auto &vertices) {
590  for (const auto v : cell->vertex_indices())
591  if (std::isnan(vertices[v][0]) == false)
592  cell->vertex(v) = vertices[v];
593  };
594 
595  if (this->is_multilevel_hierarchy_constructed())
597  std::vector<Point<spacedim>>>(*this, pack, unpack);
598  else
599  GridTools::exchange_cell_data_to_ghosts<std::vector<Point<spacedim>>>(
600  *this, pack, unpack);
601  }
602 
603 
604 
605  template <int dim, int spacedim>
606  const std::weak_ptr<const Utilities::MPI::Partitioner>
608  {
609  return number_cache.active_cell_index_partitioner;
610  }
611 
612 
613 
614  template <int dim, int spacedim>
615  const std::weak_ptr<const Utilities::MPI::Partitioner>
617  const unsigned int level) const
618  {
619  Assert(this->is_multilevel_hierarchy_constructed(), ExcNotImplemented());
620  AssertIndexRange(level, this->n_global_levels());
621 
622  return number_cache.level_cell_index_partitioners[level];
623  }
624 
625 
626 
627  template <int dim, int spacedim>
630  {
631  return number_cache.number_of_global_coarse_cells;
632  }
633 
634 
635 
636  template <int dim, int spacedim>
638  const MPI_Comm &mpi_communicator,
639  const typename ::Triangulation<dim, spacedim>::MeshSmoothing
640  smooth_grid,
641  const bool check_for_distorted_cells)
642  : ::parallel::TriangulationBase<dim, spacedim>(
643  mpi_communicator,
644  smooth_grid,
645  check_for_distorted_cells)
646  , cell_attached_data({0, 0, {}, {}})
647  , data_transfer(mpi_communicator)
648  {}
649 
650 
651 
652  template <int dim, int spacedim>
653  void
655  {
656  cell_attached_data = {0, 0, {}, {}};
657  data_transfer.clear();
658 
660  }
661 
662 
663 
664  template <int dim, int spacedim>
665  void
667  const unsigned int global_first_cell,
668  const unsigned int global_num_cells,
669  const std::string &filename) const
670  {
671  // cast away constness
672  auto tria = const_cast<
674 
675  if (this->cell_attached_data.n_attached_data_sets > 0)
676  {
677  // pack attached data first
678  tria->data_transfer.pack_data(
679  tria->local_cell_relations,
680  tria->cell_attached_data.pack_callbacks_fixed,
681  tria->cell_attached_data.pack_callbacks_variable);
682 
683  // then store buffers in file
684  tria->data_transfer.save(global_first_cell, global_num_cells, filename);
685 
686  // and release the memory afterwards
687  tria->data_transfer.clear();
688  }
689 
690  // clear all of the callback data, as explained in the documentation of
691  // register_data_attach()
692  {
693  tria->cell_attached_data.n_attached_data_sets = 0;
694  tria->cell_attached_data.pack_callbacks_fixed.clear();
695  tria->cell_attached_data.pack_callbacks_variable.clear();
696  }
697  }
698 
699 
700 
701  template <int dim, int spacedim>
702  bool
704  {
705  if (this->n_global_levels() <= 1)
706  return false; // can not have hanging nodes without refined cells
707 
708  // if there are any active cells with level less than n_global_levels()-1,
709  // then there is obviously also one with level n_global_levels()-1, and
710  // consequently there must be a hanging node somewhere.
711  //
712  // The problem is that we cannot just ask for the first active cell, but
713  // instead need to filter over locally owned cells.
714  const bool have_coarser_cell =
715  std::any_of(this->begin_active(this->n_global_levels() - 2),
716  this->end_active(this->n_global_levels() - 2),
717  [](const CellAccessor<dim, spacedim> &cell) {
718  return cell.is_locally_owned();
719  });
720 
721  // return true if at least one process has a coarser cell
722  return Utilities::MPI::max(have_coarser_cell ? 1 : 0,
723  this->mpi_communicator) != 0;
724  }
725 
726 
727 
728  template <int dim, int spacedim>
729  void
731  const unsigned int global_first_cell,
732  const unsigned int global_num_cells,
733  const unsigned int local_num_cells,
734  const std::string &filename,
735  const unsigned int n_attached_deserialize_fixed,
736  const unsigned int n_attached_deserialize_variable)
737  {
738  // load saved data, if any was stored
739  if (this->cell_attached_data.n_attached_deserialize > 0)
740  {
741  this->data_transfer.load(global_first_cell,
742  global_num_cells,
743  local_num_cells,
744  filename,
745  n_attached_deserialize_fixed,
746  n_attached_deserialize_variable);
747 
748  this->data_transfer.unpack_cell_status(this->local_cell_relations);
749 
750  // the CellStatus of all stored cells should always be CELL_PERSIST.
751  for (const auto &cell_rel : this->local_cell_relations)
752  {
753  (void)cell_rel;
754  Assert(
755  (cell_rel.second == // cell_status
757  spacedim>::CELL_PERSIST),
758  ExcInternalError());
759  }
760  }
761  }
762 
763 
764 
765  template <int dim, int spacedim>
766  unsigned int
768  const std::function<std::vector<char>(const cell_iterator &,
769  const CellStatus)> &pack_callback,
770  const bool returns_variable_size_data)
771  {
772  unsigned int handle = numbers::invalid_unsigned_int;
773 
774  // Add new callback function to the corresponding register.
775  // Encode handles according to returns_variable_size_data.
776  if (returns_variable_size_data)
777  {
778  handle = 2 * cell_attached_data.pack_callbacks_variable.size();
779  cell_attached_data.pack_callbacks_variable.push_back(pack_callback);
780  }
781  else
782  {
783  handle = 2 * cell_attached_data.pack_callbacks_fixed.size() + 1;
784  cell_attached_data.pack_callbacks_fixed.push_back(pack_callback);
785  }
786 
787  // Increase overall counter.
788  ++cell_attached_data.n_attached_data_sets;
789 
790  return handle;
791  }
792 
793 
794 
795  template <int dim, int spacedim>
796  void
798  const unsigned int handle,
799  const std::function<
800  void(const cell_iterator &,
801  const CellStatus,
802  const boost::iterator_range<std::vector<char>::const_iterator> &)>
803  &unpack_callback)
804  {
805  // perform unpacking
806  data_transfer.unpack_data(local_cell_relations, handle, unpack_callback);
807 
808  // decrease counters
809  --cell_attached_data.n_attached_data_sets;
810  if (cell_attached_data.n_attached_deserialize > 0)
811  --cell_attached_data.n_attached_deserialize;
812 
813  // important: only remove data if we are not in the deserialization
814  // process. There, each SolutionTransfer registers and unpacks before
815  // the next one does this, so n_attached_data_sets is only 1 here. This
816  // would destroy the saved data before the second SolutionTransfer can
817  // get it. This created a bug that is documented in
818  // tests/mpi/p4est_save_03 with more than one SolutionTransfer.
819 
820  if (cell_attached_data.n_attached_data_sets == 0 &&
821  cell_attached_data.n_attached_deserialize == 0)
822  {
823  // everybody got their data, time for cleanup!
824  cell_attached_data.pack_callbacks_fixed.clear();
825  cell_attached_data.pack_callbacks_variable.clear();
826  data_transfer.clear();
827 
828  // reset all cell_status entries after coarsening/refinement
829  for (auto &cell_rel : local_cell_relations)
830  cell_rel.second =
832  }
833  }
834 
835 
836 
837  /* ------------------ class DataTransfer<dim,spacedim> ----------------- */
838 
839 
840  template <int dim, int spacedim>
842  const MPI_Comm &mpi_communicator)
843  : variable_size_data_stored(false)
844  , mpi_communicator(mpi_communicator)
845  {}
846 
847 
848 
849  template <int dim, int spacedim>
850  void
852  const std::vector<cell_relation_t> &cell_relations,
853  const std::vector<typename CellAttachedData::pack_callback_t>
854  &pack_callbacks_fixed,
855  const std::vector<typename CellAttachedData::pack_callback_t>
856  &pack_callbacks_variable)
857  {
858  Assert(src_data_fixed.size() == 0,
859  ExcMessage("Previously packed data has not been released yet!"));
861 
862  const unsigned int n_callbacks_fixed = pack_callbacks_fixed.size();
863  const unsigned int n_callbacks_variable = pack_callbacks_variable.size();
864 
865  // Store information that we packed variable size data in
866  // a member variable for later.
867  variable_size_data_stored = (n_callbacks_variable > 0);
868 
869  // If variable transfer is scheduled, we will store the data size that
870  // each variable size callback function writes in this auxiliary
871  // container. The information will be stored by each cell in this vector
872  // temporarily.
873  std::vector<unsigned int> cell_sizes_variable_cumulative(
874  n_callbacks_variable);
875 
876  // Prepare the buffer structure, in which each callback function will
877  // store its data for each active cell.
878  // The outmost shell in this container construct corresponds to the
879  // data packed per cell. The next layer resembles the data that
880  // each callback function packs on the corresponding cell. These
881  // buffers are chains of chars stored in an std::vector<char>.
882  // A visualisation of the data structure:
883  /* clang-format off */
884  // | cell_1 | | cell_2 | ...
885  // || callback_1 || callback_2 |...| || callback_1 || callback_2 |...| ...
886  // |||char|char|...|||char|char|...|...| |||char|char|...|||char|char|...|...| ...
887  /* clang-format on */
888  std::vector<std::vector<std::vector<char>>> packed_fixed_size_data(
889  cell_relations.size());
890  std::vector<std::vector<std::vector<char>>> packed_variable_size_data(
891  variable_size_data_stored ? cell_relations.size() : 0);
892 
893  //
894  // --------- Pack data for fixed and variable size transfer ---------
895  //
896  // Iterate over all cells, call all callback functions on each cell,
897  // and store their data in the corresponding buffer scope.
898  {
899  auto cell_rel_it = cell_relations.cbegin();
900  auto data_cell_fixed_it = packed_fixed_size_data.begin();
901  auto data_cell_variable_it = packed_variable_size_data.begin();
902  for (; cell_rel_it != cell_relations.cend(); ++cell_rel_it)
903  {
904  const auto &dealii_cell = cell_rel_it->first;
905  const auto &cell_status = cell_rel_it->second;
906 
907  // Assertions about the tree structure.
908  switch (cell_status)
909  {
914  // double check the condition that we will only ever attach
915  // data to active cells when we get here
916  Assert(dealii_cell->is_active(), ExcInternalError());
917  break;
918 
921  // double check the condition that we will only ever attach
922  // data to cells with children when we get here. however, we
923  // can only tolerate one level of coarsening at a time, so
924  // check that the children are all active
925  Assert(dealii_cell->is_active() == false, ExcInternalError());
926  for (unsigned int c = 0;
927  c < GeometryInfo<dim>::max_children_per_cell;
928  ++c)
929  Assert(dealii_cell->child(c)->is_active(),
930  ExcInternalError());
931  break;
932 
935  // do nothing on invalid cells
936  break;
937 
938  default:
939  Assert(false, ExcInternalError());
940  break;
941  }
942 
943  // Reserve memory corresponding to the number of callback
944  // functions that will be called.
945  // If variable size transfer is scheduled, we need to leave
946  // room for an array that holds information about how many
947  // bytes each of the variable size callback functions will
948  // write.
949  // On cells flagged with CELL_INVALID, only its CellStatus
950  // will be stored.
951  const unsigned int n_fixed_size_data_sets_on_cell =
952  1 +
953  ((cell_status ==
955  spacedim>::CELL_INVALID) ?
956  0 :
957  ((variable_size_data_stored ? 1 : 0) + n_callbacks_fixed));
958  data_cell_fixed_it->resize(n_fixed_size_data_sets_on_cell);
959 
960  // We continue with packing all data on this specific cell.
961  auto data_fixed_it = data_cell_fixed_it->begin();
962 
963  // First, we pack the CellStatus information.
964  // to get consistent data sizes on each cell for the fixed size
965  // transfer, we won't allow compression
966  *data_fixed_it =
967  Utilities::pack(cell_status, /*allow_compression=*/false);
968  ++data_fixed_it;
969 
970  // Proceed with all registered callback functions.
971  // Skip cells with the CELL_INVALID flag.
972  if (cell_status !=
974  spacedim>::CELL_INVALID)
975  {
976  // Pack fixed size data.
977  for (auto callback_it = pack_callbacks_fixed.cbegin();
978  callback_it != pack_callbacks_fixed.cend();
979  ++callback_it, ++data_fixed_it)
980  {
981  *data_fixed_it = (*callback_it)(dealii_cell, cell_status);
982  }
983 
984  // Pack variable size data.
985  // If we store variable size data, we need to transfer
986  // the sizes of each corresponding callback function
987  // via fixed size transfer as well.
989  {
990  const unsigned int n_variable_size_data_sets_on_cell =
991  ((cell_status ==
993  CELL_INVALID) ?
994  0 :
995  n_callbacks_variable);
996  data_cell_variable_it->resize(
997  n_variable_size_data_sets_on_cell);
998 
999  auto callback_it = pack_callbacks_variable.cbegin();
1000  auto data_variable_it = data_cell_variable_it->begin();
1001  auto sizes_variable_it =
1002  cell_sizes_variable_cumulative.begin();
1003  for (; callback_it != pack_callbacks_variable.cend();
1004  ++callback_it, ++data_variable_it, ++sizes_variable_it)
1005  {
1006  *data_variable_it =
1007  (*callback_it)(dealii_cell, cell_status);
1008 
1009  // Store data sizes for each callback function first.
1010  // Make it cumulative below.
1011  *sizes_variable_it = data_variable_it->size();
1012  }
1013 
1014  // Turn size vector into its cumulative representation.
1015  std::partial_sum(cell_sizes_variable_cumulative.begin(),
1016  cell_sizes_variable_cumulative.end(),
1017  cell_sizes_variable_cumulative.begin());
1018 
1019  // Serialize cumulative variable size vector
1020  // value-by-value. This way we can circumvent the overhead
1021  // of storing the container object as a whole, since we
1022  // know its size by the number of registered callback
1023  // functions.
1024  data_fixed_it->resize(n_callbacks_variable *
1025  sizeof(unsigned int));
1026  for (unsigned int i = 0; i < n_callbacks_variable; ++i)
1027  std::memcpy(&(data_fixed_it->at(i * sizeof(unsigned int))),
1028  &(cell_sizes_variable_cumulative.at(i)),
1029  sizeof(unsigned int));
1030 
1031  ++data_fixed_it;
1032  }
1033 
1034  // Double check that we packed everything we wanted
1035  // in the fixed size buffers.
1036  Assert(data_fixed_it == data_cell_fixed_it->end(),
1037  ExcInternalError());
1038  }
1039 
1040  ++data_cell_fixed_it;
1041 
1042  // Increment the variable size data iterator
1043  // only if we actually pack this kind of data
1044  // to avoid getting out of bounds.
1046  ++data_cell_variable_it;
1047  } // loop over cell_relations
1048  }
1049 
1050  //
1051  // ----------- Gather data sizes for fixed size transfer ------------
1052  //
1053  // Generate a vector which stores the sizes of each callback function,
1054  // including the packed CellStatus transfer.
1055  // Find the very first cell that we wrote to with all callback
1056  // functions (i.e. a cell that was not flagged with CELL_INVALID)
1057  // and store the sizes of each buffer.
1058  //
1059  // To deal with the case that at least one of the processors does not
1060  // own any cell at all, we will exchange the information about the data
1061  // sizes among them later. The code in between is still well-defined,
1062  // since the following loops will be skipped.
1063  std::vector<unsigned int> local_sizes_fixed(
1064  1 + n_callbacks_fixed + (variable_size_data_stored ? 1 : 0));
1065  for (const auto &data_cell : packed_fixed_size_data)
1066  {
1067  if (data_cell.size() == local_sizes_fixed.size())
1068  {
1069  auto sizes_fixed_it = local_sizes_fixed.begin();
1070  auto data_fixed_it = data_cell.cbegin();
1071  for (; data_fixed_it != data_cell.cend();
1072  ++data_fixed_it, ++sizes_fixed_it)
1073  {
1074  *sizes_fixed_it = data_fixed_it->size();
1075  }
1076 
1077  break;
1078  }
1079  }
1080 
1081  // Check if all cells have valid sizes.
1082  for (auto data_cell_fixed_it = packed_fixed_size_data.cbegin();
1083  data_cell_fixed_it != packed_fixed_size_data.cend();
1084  ++data_cell_fixed_it)
1085  {
1086  Assert((data_cell_fixed_it->size() == 1) ||
1087  (data_cell_fixed_it->size() == local_sizes_fixed.size()),
1088  ExcInternalError());
1089  }
1090 
1091  // Share information about the packed data sizes
1092  // of all callback functions across all processors, in case one
1093  // of them does not own any cells at all.
1094  std::vector<unsigned int> global_sizes_fixed(local_sizes_fixed.size());
1095  Utilities::MPI::max(local_sizes_fixed,
1096  this->mpi_communicator,
1097  global_sizes_fixed);
1098 
1099  // Construct cumulative sizes, since this is the only information
1100  // we need from now on.
1101  sizes_fixed_cumulative.resize(global_sizes_fixed.size());
1102  std::partial_sum(global_sizes_fixed.begin(),
1103  global_sizes_fixed.end(),
1104  sizes_fixed_cumulative.begin());
1105 
1106  //
1107  // ---------- Gather data sizes for variable size transfer ----------
1108  //
1110  {
1111  src_sizes_variable.reserve(packed_variable_size_data.size());
1112  for (const auto &data_cell : packed_variable_size_data)
1113  {
1114  int variable_data_size_on_cell = 0;
1115 
1116  for (const auto &data : data_cell)
1117  variable_data_size_on_cell += data.size();
1118 
1119  src_sizes_variable.push_back(variable_data_size_on_cell);
1120  }
1121  }
1122 
1123  //
1124  // ------------------------ Build buffers ---------------------------
1125  //
1126  const unsigned int expected_size_fixed =
1127  cell_relations.size() * sizes_fixed_cumulative.back();
1128  const unsigned int expected_size_variable =
1129  std::accumulate(src_sizes_variable.begin(),
1130  src_sizes_variable.end(),
1132 
1133  // Move every piece of packed fixed size data into the consecutive
1134  // buffer.
1135  src_data_fixed.reserve(expected_size_fixed);
1136  for (const auto &data_cell_fixed : packed_fixed_size_data)
1137  {
1138  // Move every fraction of packed data into the buffer
1139  // reserved for this particular cell.
1140  for (const auto &data_fixed : data_cell_fixed)
1141  std::move(data_fixed.begin(),
1142  data_fixed.end(),
1143  std::back_inserter(src_data_fixed));
1144 
1145  // If we only packed the CellStatus information
1146  // (i.e. encountered a cell flagged CELL_INVALID),
1147  // fill the remaining space with invalid entries.
1148  // We can skip this if there is nothing else to pack.
1149  if ((data_cell_fixed.size() == 1) &&
1150  (sizes_fixed_cumulative.size() > 1))
1151  {
1152  const std::size_t bytes_skipped =
1154 
1155  src_data_fixed.insert(src_data_fixed.end(),
1156  bytes_skipped,
1157  static_cast<char>(-1)); // invalid_char
1158  }
1159  }
1160 
1161  // Move every piece of packed variable size data into the consecutive
1162  // buffer.
1164  {
1165  src_data_variable.reserve(expected_size_variable);
1166  for (const auto &data_cell : packed_variable_size_data)
1167  {
1168  // Move every fraction of packed data into the buffer
1169  // reserved for this particular cell.
1170  for (const auto &data : data_cell)
1171  std::move(data.begin(),
1172  data.end(),
1173  std::back_inserter(src_data_variable));
1174  }
1175  }
1176 
1177  // Double check that we packed everything correctly.
1178  Assert(src_data_fixed.size() == expected_size_fixed, ExcInternalError());
1179  Assert(src_data_variable.size() == expected_size_variable,
1180  ExcInternalError());
1181  }
1182 
1183 
1184 
1185  template <int dim, int spacedim>
1186  void
1188  std::vector<cell_relation_t> &cell_relations) const
1189  {
1190  Assert(sizes_fixed_cumulative.size() > 0,
1191  ExcMessage("No data has been packed!"));
1192  if (cell_relations.size() > 0)
1193  {
1194  Assert(dest_data_fixed.size() > 0,
1195  ExcMessage("No data has been received!"));
1196  }
1197 
1198  // Size of CellStatus object that will be unpacked on each cell.
1199  const unsigned int size = sizes_fixed_cumulative.front();
1200 
1201  // Iterate over all cells and overwrite the CellStatus
1202  // information from the transferred data.
1203  // Proceed buffer iterator position to next cell after
1204  // each iteration.
1205  auto cell_rel_it = cell_relations.begin();
1206  auto dest_fixed_it = dest_data_fixed.cbegin();
1207  for (; cell_rel_it != cell_relations.end();
1208  ++cell_rel_it, dest_fixed_it += sizes_fixed_cumulative.back())
1209  {
1210  cell_rel_it->second = // cell_status
1212  dim,
1213  spacedim>::CellStatus>(dest_fixed_it,
1214  dest_fixed_it + size,
1215  /*allow_compression=*/false);
1216  }
1217  }
1218 
1219 
1220 
1221  template <int dim, int spacedim>
1222  void
1224  const std::vector<cell_relation_t> &cell_relations,
1225  const unsigned int handle,
1226  const std::function<
1227  void(const typename ::Triangulation<dim, spacedim>::cell_iterator &,
1228  const typename ::Triangulation<dim, spacedim>::CellStatus &,
1229  const boost::iterator_range<std::vector<char>::const_iterator> &)>
1230  &unpack_callback) const
1231  {
1232  // We decode the handle returned by register_data_attach() back into
1233  // a format we can use. All even handles belong to those callback
1234  // functions which write/read variable size data, all odd handles
1235  // interact with fixed size buffers.
1236  const bool callback_variable_transfer = (handle % 2 == 0);
1237  const unsigned int callback_index = handle / 2;
1238 
1239  // Cells will always receive fixed size data (i.e., CellStatus
1240  // information), but not necessarily variable size data (e.g., with a
1241  // ParticleHandler a cell might not contain any particle at all).
1242  // Thus it is sufficient to check if fixed size data has been received.
1243  Assert(sizes_fixed_cumulative.size() > 0,
1244  ExcMessage("No data has been packed!"));
1245  if (cell_relations.size() > 0)
1246  {
1247  Assert(dest_data_fixed.size() > 0,
1248  ExcMessage("No data has been received!"));
1249  }
1250 
1251  std::vector<char>::const_iterator dest_data_it;
1252  std::vector<char>::const_iterator dest_sizes_cell_it;
1253 
1254  // Depending on whether our callback function unpacks fixed or
1255  // variable size data, we have to pursue different approaches
1256  // to localize the correct fraction of the buffer from which
1257  // we are allowed to read.
1258  unsigned int offset = numbers::invalid_unsigned_int;
1259  unsigned int size = numbers::invalid_unsigned_int;
1260  unsigned int data_increment = numbers::invalid_unsigned_int;
1261 
1262  if (callback_variable_transfer)
1263  {
1264  // For the variable size data, we need to extract the
1265  // data size from the fixed size buffer on each cell.
1266  //
1267  // We packed this information last, so the last packed
1268  // object in the fixed size buffer corresponds to the
1269  // variable data sizes.
1270  //
1271  // The last entry of sizes_fixed_cumulative corresponds
1272  // to the size of all fixed size data packed on the cell.
1273  // To get the offset for the last packed object, we need
1274  // to get the next-to-last entry.
1275  const unsigned int offset_variable_data_sizes =
1276  sizes_fixed_cumulative[sizes_fixed_cumulative.size() - 2];
1277 
1278  // This iterator points to the data size that the
1279  // callback_function packed for each specific cell.
1280  // Adjust buffer iterator to the offset of the callback
1281  // function so that we only have to advance its position
1282  // to the next cell after each iteration.
1283  dest_sizes_cell_it = dest_data_fixed.cbegin() +
1284  offset_variable_data_sizes +
1285  callback_index * sizeof(unsigned int);
1286 
1287  // Let the data iterator point to the correct buffer.
1288  dest_data_it = dest_data_variable.cbegin();
1289  }
1290  else
1291  {
1292  // For the fixed size data, we can get the information about
1293  // the buffer location on each cell directly from the
1294  // sizes_fixed_cumulative vector.
1295  offset = sizes_fixed_cumulative[callback_index];
1296  size = sizes_fixed_cumulative[callback_index + 1] - offset;
1297  data_increment = sizes_fixed_cumulative.back();
1298 
1299  // Let the data iterator point to the correct buffer.
1300  // Adjust buffer iterator to the offset of the callback
1301  // function so that we only have to advance its position
1302  // to the next cell after each iteration.
1303  if (cell_relations.begin() != cell_relations.end())
1304  dest_data_it = dest_data_fixed.cbegin() + offset;
1305  }
1306 
1307  // Iterate over all cells and unpack the transferred data.
1308  auto cell_rel_it = cell_relations.begin();
1309  auto dest_sizes_it = dest_sizes_variable.cbegin();
1310  for (; cell_rel_it != cell_relations.end(); ++cell_rel_it)
1311  {
1312  const auto &dealii_cell = cell_rel_it->first;
1313  const auto &cell_status = cell_rel_it->second;
1314 
1315  if (callback_variable_transfer)
1316  {
1317  // Update the increment according to the whole data size
1318  // of the current cell.
1319  data_increment = *dest_sizes_it;
1320 
1321  if (cell_status !=
1323  spacedim>::CELL_INVALID)
1324  {
1325  // Extract the corresponding values for offset and size from
1326  // the cumulative sizes array stored in the fixed size
1327  // buffer.
1328  if (callback_index == 0)
1329  offset = 0;
1330  else
1331  std::memcpy(&offset,
1332  &(*(dest_sizes_cell_it - sizeof(unsigned int))),
1333  sizeof(unsigned int));
1334 
1335  std::memcpy(&size,
1336  &(*dest_sizes_cell_it),
1337  sizeof(unsigned int));
1338 
1339  size -= offset;
1340 
1341  // Move the data iterator to the corresponding position
1342  // of the callback function and adjust the increment
1343  // accordingly.
1344  dest_data_it += offset;
1345  data_increment -= offset;
1346  }
1347 
1348  // Advance data size iterators to the next cell, avoid iterating
1349  // past the end of dest_sizes_cell_it
1350  if (cell_rel_it != cell_relations.end() - 1)
1351  dest_sizes_cell_it += sizes_fixed_cumulative.back();
1352  ++dest_sizes_it;
1353  }
1354 
1355  switch (cell_status)
1356  {
1358  spacedim>::CELL_PERSIST:
1360  spacedim>::CELL_COARSEN:
1361  unpack_callback(dealii_cell,
1362  cell_status,
1363  boost::make_iterator_range(dest_data_it,
1364  dest_data_it + size));
1365  break;
1366 
1368  spacedim>::CELL_REFINE:
1369  unpack_callback(dealii_cell->parent(),
1370  cell_status,
1371  boost::make_iterator_range(dest_data_it,
1372  dest_data_it + size));
1373  break;
1374 
1376  spacedim>::CELL_INVALID:
1377  // Skip this cell.
1378  break;
1379 
1380  default:
1381  Assert(false, ExcInternalError());
1382  break;
1383  }
1384 
1385  if (cell_rel_it != cell_relations.end() - 1)
1386  dest_data_it += data_increment;
1387  }
1388  }
1389 
1390 
1391 
1392  template <int dim, int spacedim>
1393  void
1395  const unsigned int global_first_cell,
1396  const unsigned int global_num_cells,
1397  const std::string &filename) const
1398  {
1399 #ifdef DEAL_II_WITH_MPI
1400  // Large fractions of this function have been copied from
1401  // DataOutInterface::write_vtu_in_parallel.
1402  // TODO: Write general MPIIO interface.
1403 
1404  Assert(sizes_fixed_cumulative.size() > 0,
1405  ExcMessage("No data has been packed!"));
1406 
1408 
1409  const unsigned int bytes_per_cell = sizes_fixed_cumulative.back();
1410 
1411  //
1412  // ---------- Fixed size data ----------
1413  //
1414  {
1415  const std::string fname_fixed = std::string(filename) + "_fixed.data";
1416 
1417  MPI_Info info;
1418  int ierr = MPI_Info_create(&info);
1419  AssertThrowMPI(ierr);
1420 
1421  MPI_File fh;
1422  ierr = MPI_File_open(mpi_communicator,
1423  fname_fixed.c_str(),
1424  MPI_MODE_CREATE | MPI_MODE_WRONLY,
1425  info,
1426  &fh);
1427  AssertThrowMPI(ierr);
1428 
1429  ierr = MPI_File_set_size(fh, 0); // delete the file contents
1430  AssertThrowMPI(ierr);
1431  // this barrier is necessary, because otherwise others might already
1432  // write while one core is still setting the size to zero.
1433  ierr = MPI_Barrier(mpi_communicator);
1434  AssertThrowMPI(ierr);
1435  ierr = MPI_Info_free(&info);
1436  AssertThrowMPI(ierr);
1437  // ------------------
1438 
1439  // Write cumulative sizes to file.
1440  // Since each processor owns the same information about the data
1441  // sizes, it is sufficient to let only the first processor perform
1442  // this task.
1443  if (myrank == 0)
1444  {
1446  fh,
1447  0,
1448  sizes_fixed_cumulative.data(),
1449  sizes_fixed_cumulative.size(),
1450  MPI_UNSIGNED,
1451  MPI_STATUS_IGNORE);
1452  AssertThrowMPI(ierr);
1453  }
1454 
1455  // Write packed data to file simultaneously.
1456  const MPI_Offset size_header =
1457  sizes_fixed_cumulative.size() * sizeof(unsigned int);
1458 
1459  // Make sure we do the following computation in 64bit integers to be
1460  // able to handle 4GB+ files:
1461  const MPI_Offset my_global_file_position =
1462  size_header +
1463  static_cast<MPI_Offset>(global_first_cell) * bytes_per_cell;
1464 
1465  ierr =
1467  my_global_file_position,
1468  src_data_fixed.data(),
1469  src_data_fixed.size(),
1470  MPI_BYTE,
1471  MPI_STATUS_IGNORE);
1472  AssertThrowMPI(ierr);
1473 
1474  ierr = MPI_File_close(&fh);
1475  AssertThrowMPI(ierr);
1476  }
1477 
1478 
1479 
1480  //
1481  // ---------- Variable size data ----------
1482  //
1483  if (variable_size_data_stored)
1484  {
1485  const std::string fname_variable =
1486  std::string(filename) + "_variable.data";
1487 
1488  MPI_Info info;
1489  int ierr = MPI_Info_create(&info);
1490  AssertThrowMPI(ierr);
1491 
1492  MPI_File fh;
1493  ierr = MPI_File_open(mpi_communicator,
1494  fname_variable.c_str(),
1495  MPI_MODE_CREATE | MPI_MODE_WRONLY,
1496  info,
1497  &fh);
1498  AssertThrowMPI(ierr);
1499 
1500  ierr = MPI_File_set_size(fh, 0); // delete the file contents
1501  AssertThrowMPI(ierr);
1502  // this barrier is necessary, because otherwise others might already
1503  // write while one core is still setting the size to zero.
1504  ierr = MPI_Barrier(mpi_communicator);
1505  AssertThrowMPI(ierr);
1506  ierr = MPI_Info_free(&info);
1507  AssertThrowMPI(ierr);
1508 
1509  // Write sizes of each cell into file simultaneously.
1510  {
1511  const MPI_Offset my_global_file_position =
1512  static_cast<MPI_Offset>(global_first_cell) * sizeof(unsigned int);
1513 
1514  // It is very unlikely that a single process has more than
1515  // 2 billion cells, but we might as well check.
1516  AssertThrow(src_sizes_variable.size() <
1517  static_cast<std::size_t>(
1519  ExcNotImplemented());
1520 
1522  fh,
1523  my_global_file_position,
1524  src_sizes_variable.data(),
1525  src_sizes_variable.size(),
1526  MPI_INT,
1527  MPI_STATUS_IGNORE);
1528  AssertThrowMPI(ierr);
1529  }
1530 
1531  // Gather size of data in bytes we want to store from this
1532  // processor and compute the prefix sum. We do this in 64 bit
1533  // to avoid overflow for files larger than 4GB:
1534  const std::uint64_t size_on_proc = src_data_variable.size();
1535  std::uint64_t prefix_sum = 0;
1536  ierr = MPI_Exscan(&size_on_proc,
1537  &prefix_sum,
1538  1,
1539  MPI_UINT64_T,
1540  MPI_SUM,
1542  AssertThrowMPI(ierr);
1543 
1544  const MPI_Offset my_global_file_position =
1545  static_cast<MPI_Offset>(global_num_cells) * sizeof(unsigned int) +
1546  prefix_sum;
1547 
1548  // Write data consecutively into file.
1549  ierr =
1551  my_global_file_position,
1552  src_data_variable.data(),
1553  src_data_variable.size(),
1554  MPI_BYTE,
1555  MPI_STATUS_IGNORE);
1556  AssertThrowMPI(ierr);
1557 
1558 
1559  ierr = MPI_File_close(&fh);
1560  AssertThrowMPI(ierr);
1561  }
1562 #else
1563  (void)global_first_cell;
1564  (void)global_num_cells;
1565  (void)filename;
1566 
1567  AssertThrow(false, ExcNeedsMPI());
1568 #endif
1569  }
1570 
1571 
1572 
1573  template <int dim, int spacedim>
1574  void
1576  const unsigned int global_first_cell,
1577  const unsigned int global_num_cells,
1578  const unsigned int local_num_cells,
1579  const std::string &filename,
1580  const unsigned int n_attached_deserialize_fixed,
1581  const unsigned int n_attached_deserialize_variable)
1582  {
1583 #ifdef DEAL_II_WITH_MPI
1584  // Large fractions of this function have been copied from
1585  // DataOutInterface::write_vtu_in_parallel.
1586  // TODO: Write general MPIIO interface.
1587 
1588  Assert(dest_data_fixed.size() == 0,
1589  ExcMessage("Previously loaded data has not been released yet!"));
1590 
1591  variable_size_data_stored = (n_attached_deserialize_variable > 0);
1592 
1593  //
1594  // ---------- Fixed size data ----------
1595  //
1596  {
1597  const std::string fname_fixed = std::string(filename) + "_fixed.data";
1598 
1599  MPI_Info info;
1600  int ierr = MPI_Info_create(&info);
1601  AssertThrowMPI(ierr);
1602 
1603  MPI_File fh;
1604  ierr = MPI_File_open(
1605  mpi_communicator, fname_fixed.c_str(), MPI_MODE_RDONLY, info, &fh);
1606  AssertThrowMPI(ierr);
1607 
1608  ierr = MPI_Info_free(&info);
1609  AssertThrowMPI(ierr);
1610 
1611  // Read cumulative sizes from file.
1612  // Since all processors need the same information about the data
1613  // sizes, let each of them retrieve it by reading from the same
1614  // location in the file.
1615  sizes_fixed_cumulative.resize(1 + n_attached_deserialize_fixed +
1616  (variable_size_data_stored ? 1 : 0));
1618  fh,
1619  0,
1620  sizes_fixed_cumulative.data(),
1621  sizes_fixed_cumulative.size(),
1622  MPI_UNSIGNED,
1623  MPI_STATUS_IGNORE);
1624  AssertThrowMPI(ierr);
1625 
1626  // Allocate sufficient memory.
1627  const unsigned int bytes_per_cell = sizes_fixed_cumulative.back();
1628  dest_data_fixed.resize(static_cast<size_t>(local_num_cells) *
1629  bytes_per_cell);
1630 
1631  // Read packed data from file simultaneously.
1632  const MPI_Offset size_header =
1633  sizes_fixed_cumulative.size() * sizeof(unsigned int);
1634 
1635  // Make sure we do the following computation in 64bit integers to be
1636  // able to handle 4GB+ files:
1637  const MPI_Offset my_global_file_position =
1638  size_header +
1639  static_cast<MPI_Offset>(global_first_cell) * bytes_per_cell;
1640 
1642  my_global_file_position,
1643  dest_data_fixed.data(),
1644  dest_data_fixed.size(),
1645  MPI_BYTE,
1646  MPI_STATUS_IGNORE);
1647  AssertThrowMPI(ierr);
1648 
1649 
1650  ierr = MPI_File_close(&fh);
1651  AssertThrowMPI(ierr);
1652  }
1653 
1654  //
1655  // ---------- Variable size data ----------
1656  //
1657  if (variable_size_data_stored)
1658  {
1659  const std::string fname_variable =
1660  std::string(filename) + "_variable.data";
1661 
1662  MPI_Info info;
1663  int ierr = MPI_Info_create(&info);
1664  AssertThrowMPI(ierr);
1665 
1666  MPI_File fh;
1667  ierr = MPI_File_open(
1668  mpi_communicator, fname_variable.c_str(), MPI_MODE_RDONLY, info, &fh);
1669  AssertThrowMPI(ierr);
1670 
1671  ierr = MPI_Info_free(&info);
1672  AssertThrowMPI(ierr);
1673 
1674  // Read sizes of all locally owned cells.
1675  dest_sizes_variable.resize(local_num_cells);
1676 
1677  const MPI_Offset my_global_file_position_sizes =
1678  static_cast<MPI_Offset>(global_first_cell) * sizeof(unsigned int);
1679 
1681  fh,
1682  my_global_file_position_sizes,
1683  dest_sizes_variable.data(),
1684  dest_sizes_variable.size(),
1685  MPI_INT,
1686  MPI_STATUS_IGNORE);
1687  AssertThrowMPI(ierr);
1688 
1689 
1690  // Compute my data size in bytes and compute prefix sum. We do this
1691  // in 64 bit to avoid overflow for files larger than 4 GB:
1692  const std::uint64_t size_on_proc =
1693  std::accumulate(dest_sizes_variable.begin(),
1694  dest_sizes_variable.end(),
1695  0ULL);
1696 
1697  std::uint64_t prefix_sum = 0;
1698  ierr = MPI_Exscan(&size_on_proc,
1699  &prefix_sum,
1700  1,
1701  MPI_UINT64_T,
1702  MPI_SUM,
1704  AssertThrowMPI(ierr);
1705 
1706  const MPI_Offset my_global_file_position =
1707  static_cast<MPI_Offset>(global_num_cells) * sizeof(unsigned int) +
1708  prefix_sum;
1709 
1710  dest_data_variable.resize(size_on_proc);
1711 
1712  ierr =
1714  my_global_file_position,
1715  dest_data_variable.data(),
1716  dest_data_variable.size(),
1717  MPI_BYTE,
1718  MPI_STATUS_IGNORE);
1719  AssertThrowMPI(ierr);
1720 
1721  ierr = MPI_File_close(&fh);
1722  AssertThrowMPI(ierr);
1723  }
1724 #else
1725  (void)global_first_cell;
1726  (void)global_num_cells;
1727  (void)local_num_cells;
1728  (void)filename;
1729  (void)n_attached_deserialize_fixed;
1730  (void)n_attached_deserialize_variable;
1731 
1732  AssertThrow(false, ExcNeedsMPI());
1733 #endif
1734  }
1735 
1736 
1737 
1738  template <int dim, int spacedim>
1739  void
1741  {
1742  variable_size_data_stored = false;
1743 
1744  // free information about data sizes
1745  sizes_fixed_cumulative.clear();
1746  sizes_fixed_cumulative.shrink_to_fit();
1747 
1748  // free fixed size transfer data
1749  src_data_fixed.clear();
1750  src_data_fixed.shrink_to_fit();
1751 
1752  dest_data_fixed.clear();
1753  dest_data_fixed.shrink_to_fit();
1754 
1755  // free variable size transfer data
1756  src_sizes_variable.clear();
1757  src_sizes_variable.shrink_to_fit();
1758 
1759  src_data_variable.clear();
1760  src_data_variable.shrink_to_fit();
1761 
1762  dest_sizes_variable.clear();
1763  dest_sizes_variable.shrink_to_fit();
1764 
1765  dest_data_variable.clear();
1766  dest_data_variable.shrink_to_fit();
1767  }
1768 
1769 } // end namespace parallel
1770 
1771 
1772 
1773 /*-------------- Explicit Instantiations -------------------------------*/
1774 #include "tria_base.inst"
1775 
bool is_locally_owned() const
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1054
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1662
virtual void clear()
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
void save(Archive &ar, const unsigned int version) const
virtual std::size_t memory_consumption() const
virtual void update_reference_cells()
void pack_data(const std::vector< cell_relation_t > &cell_relations, const std::vector< typename CellAttachedData::pack_callback_t > &pack_callbacks_fixed, const std::vector< typename CellAttachedData::pack_callback_t > &pack_callbacks_variable)
Definition: tria_base.cc:851
void save(const unsigned int global_first_cell, const unsigned int global_num_cells, const std::string &filename) const
Definition: tria_base.cc:1394
void unpack_cell_status(std::vector< cell_relation_t > &cell_relations) const
Definition: tria_base.cc:1187
void unpack_data(const std::vector< cell_relation_t > &cell_relations, const unsigned int handle, const std::function< void(const typename ::Triangulation< dim, spacedim >::cell_iterator &, const typename ::Triangulation< dim, spacedim >::CellStatus &, const boost::iterator_range< std::vector< char >::const_iterator > &)> &unpack_callback) const
Definition: tria_base.cc:1223
DataTransfer(const MPI_Comm &mpi_communicator)
Definition: tria_base.cc:841
std::vector< unsigned int > sizes_fixed_cumulative
Definition: tria_base.h:897
void load(const unsigned int global_first_cell, const unsigned int global_num_cells, const unsigned int local_num_cells, const std::string &filename, const unsigned int n_attached_deserialize_fixed, const unsigned int n_attached_deserialize_variable)
Definition: tria_base.cc:1575
virtual void clear() override
Definition: tria_base.cc:654
void save_attached_data(const unsigned int global_first_cell, const unsigned int global_num_cells, const std::string &filename) const
Definition: tria_base.cc:666
virtual bool has_hanging_nodes() const override
Definition: tria_base.cc:703
void notify_ready_to_unpack(const unsigned int handle, const std::function< void(const cell_iterator &, const CellStatus, const boost::iterator_range< std::vector< char >::const_iterator > &)> &unpack_callback)
Definition: tria_base.cc:797
typename ::Triangulation< dim, spacedim >::CellStatus CellStatus
Definition: tria_base.h:459
unsigned int register_data_attach(const std::function< std::vector< char >(const cell_iterator &, const CellStatus)> &pack_callback, const bool returns_variable_size_data)
Definition: tria_base.cc:767
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition: tria_base.h:456
void load_attached_data(const unsigned int global_first_cell, const unsigned int global_num_cells, const unsigned int local_num_cells, const std::string &filename, const unsigned int n_attached_deserialize_fixed, const unsigned int n_attached_deserialize_variable)
Definition: tria_base.cc:730
const std::set< types::subdomain_id > & level_ghost_owners() const
Definition: tria_base.cc:342
virtual types::global_cell_index n_global_active_cells() const override
Definition: tria_base.cc:136
const std::set< types::subdomain_id > & ghost_owners() const
Definition: tria_base.cc:333
TriangulationBase(const MPI_Comm &mpi_communicator, const typename ::Triangulation< dim, spacedim >::MeshSmoothing smooth_grid=(::Triangulation< dim, spacedim >::none), const bool check_for_distorted_cells=false)
Definition: tria_base.cc:47
virtual unsigned int n_global_levels() const override
Definition: tria_base.cc:129
unsigned int n_locally_owned_active_cells() const
Definition: tria_base.cc:122
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
Definition: tria_base.cc:67
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
Point< 3 > vertices[4]
Point< 2 > second
Definition: grid_out.cc:4606
Point< 2 > first
Definition: grid_out.cc:4605
unsigned int level
Definition: grid_out.cc:4608
unsigned int cell_index
Definition: grid_tools.cc:1177
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsMPI()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1818
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
IteratorRange< BaseIterator > make_iterator_range(const BaseIterator &begin, const typename identity< BaseIterator >::type &end)
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4458
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
types::global_dof_index size_type
Definition: cuda_kernels.h:45
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
int File_write_at_c(MPI_File fh, MPI_Offset offset, const void *buf, MPI_Count count, MPI_Datatype datatype, MPI_Status *status)
int File_read_at_c(MPI_File fh, MPI_Offset offset, void *buf, MPI_Count count, MPI_Datatype datatype, MPI_Status *status)
@ triangulation_base_fill_level_ghost_owners
TriangulationBase<dim, spacedim>::fill_level_ghost_owners()
Definition: mpi_tags.h:102
std::vector< T > compute_set_union(const std::vector< T > &vec, const MPI_Comm &comm)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:157
T sum(const T &t, const MPI_Comm &mpi_communicator)
const MPI_Datatype mpi_type_id_for_type
Definition: mpi.h:1728
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:146
T max(const T &t, const MPI_Comm &mpi_communicator)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1334
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1499
constexpr ::ReferenceCell make_reference_cell_from_int(const std::uint8_t kind)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13757
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13750
const types::subdomain_id artificial_subdomain_id
Definition: types.h:308
static const unsigned int invalid_unsigned_int
Definition: types.h:206
const types::global_dof_index invalid_dof_index
Definition: types.h:226
global_cell_index coarse_cell_id
Definition: types.h:119
unsigned int subdomain_id
Definition: types.h:43
unsigned int global_cell_index
Definition: types.h:110
const ::Triangulation< dim, spacedim > & tria