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