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