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