Reference documentation for deal.II version GIT 01a9543f1b 2023-12-05 20:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
mg_transfer_internal.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2023 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 
18 
19 #include <deal.II/dofs/dof_tools.h>
20 
21 #include <deal.II/fe/fe_tools.h>
22 
24 
26 
27 #include <memory>
28 
30 
31 namespace internal
32 {
33  namespace MGTransfer
34  {
35  // Internal data structure that is used in the MPI communication in
36  // fill_copy_indices(). It represents an entry in the copy_indices* map,
37  // that associates a level dof index with a global dof index.
38  struct DoFPair
39  {
40  unsigned int level;
43 
44  DoFPair(const unsigned int level,
47  : level(level)
50  {}
51 
56  {}
57  };
58 
59  template <int dim, int spacedim>
60  void
62  const DoFHandler<dim, spacedim> &dof_handler,
63  const MGConstrainedDoFs *mg_constrained_dofs,
64  std::vector<std::vector<
65  std::pair<types::global_dof_index, types::global_dof_index>>>
66  &copy_indices,
67  std::vector<std::vector<
68  std::pair<types::global_dof_index, types::global_dof_index>>>
69  &copy_indices_global_mine,
70  std::vector<std::vector<
71  std::pair<types::global_dof_index, types::global_dof_index>>>
72  &copy_indices_level_mine,
73  const bool skip_interface_dofs)
74  {
75  // Now we are filling the variables copy_indices*, which are essentially
76  // maps from global to mgdof for each level stored as a std::vector of
77  // pairs. We need to split this map on each level depending on the
78  // ownership of the global and mgdof, so that we later do not access
79  // non-local elements in copy_to/from_mg.
80  // We keep track in the bitfield dof_touched which global dof has been
81  // processed already (on the current level). This is the same as the
82  // multigrid running in serial.
83 
84  // map cpu_index -> vector of data
85  // that will be copied into copy_indices_level_mine
86  std::vector<DoFPair> send_data_temp;
87 
88  const unsigned int n_levels =
89  dof_handler.get_triangulation().n_global_levels();
90  copy_indices.resize(n_levels);
91  copy_indices_global_mine.resize(n_levels);
92  copy_indices_level_mine.resize(n_levels);
93  const IndexSet &owned_dofs = dof_handler.locally_owned_dofs();
94 
95  const unsigned int dofs_per_cell = dof_handler.get_fe().n_dofs_per_cell();
96  std::vector<types::global_dof_index> global_dof_indices(dofs_per_cell);
97  std::vector<types::global_dof_index> level_dof_indices(dofs_per_cell);
98 
99  for (unsigned int level = 0; level < n_levels; ++level)
100  {
101  std::vector<bool> dof_touched(owned_dofs.n_elements(), false);
102  const IndexSet &owned_level_dofs =
103  dof_handler.locally_owned_mg_dofs(level);
104 
105  // for the most common case where copy_indices are locally owned
106  // both globally and on the level, we want to skip collecting pairs
107  // and later sorting them. instead, we insert these indices into a
108  // plain vector
109  std::vector<types::global_dof_index> unrolled_copy_indices;
110 
111  copy_indices_level_mine[level].clear();
112  copy_indices_global_mine[level].clear();
113 
114  for (const auto &level_cell :
116  {
117  if (dof_handler.get_triangulation().locally_owned_subdomain() !=
119  (level_cell->level_subdomain_id() ==
121  level_cell->subdomain_id() ==
123  continue;
124 
125  unrolled_copy_indices.resize(owned_dofs.n_elements(),
127 
128  // get the dof numbers of this cell for the global and the
129  // level-wise numbering
130  level_cell->get_dof_indices(global_dof_indices);
131  level_cell->get_mg_dof_indices(level_dof_indices);
132 
133  for (unsigned int i = 0; i < dofs_per_cell; ++i)
134  {
135  // we need to ignore if the DoF is on a refinement edge
136  // (hanging node)
137  if (skip_interface_dofs && mg_constrained_dofs != nullptr &&
138  mg_constrained_dofs->at_refinement_edge(
139  level, level_dof_indices[i]))
140  continue;
141 
142  // First check whether we own any of the active dof index
143  // and the level one. This check involves locally owned
144  // indices which often consist only of a single range, so
145  // they are cheap to look up.
146  const types::global_dof_index global_index_in_set =
147  owned_dofs.index_within_set(global_dof_indices[i]);
148  bool global_mine =
149  global_index_in_set != numbers::invalid_dof_index;
150  bool level_mine =
151  owned_level_dofs.is_element(level_dof_indices[i]);
152 
153  if (global_mine && level_mine)
154  {
155  // we own both the active dof index and the level one ->
156  // set them into the vector, indexed by the local index
157  // range of the active dof
158  unrolled_copy_indices[global_index_in_set] =
159  level_dof_indices[i];
160  }
161  else if (global_mine &&
162  dof_touched[global_index_in_set] == false)
163  {
164  copy_indices_global_mine[level].emplace_back(
165  global_dof_indices[i], level_dof_indices[i]);
166 
167  // send this to the owner of the level_dof:
168  send_data_temp.emplace_back(level,
169  global_dof_indices[i],
170  level_dof_indices[i]);
171  dof_touched[global_index_in_set] = true;
172  }
173  else
174  {
175  // somebody will send those to me
176  }
177  }
178  }
179 
180  // we now translate the plain vector for the copy_indices field into
181  // the expected format of a pair of indices
182  if (!unrolled_copy_indices.empty())
183  {
184  copy_indices[level].clear();
185 
186  // reserve the full length in case we did not hit global-mine
187  // indices, so we expect all indices to come into copy_indices
188  if (copy_indices_global_mine[level].empty())
189  copy_indices[level].reserve(unrolled_copy_indices.size());
190 
191  // owned_dofs.nth_index_in_set(i) in this query is
192  // usually cheap to look up as there are few ranges in
193  // the locally owned part
194  for (unsigned int i = 0; i < unrolled_copy_indices.size(); ++i)
195  if (unrolled_copy_indices[i] != numbers::invalid_dof_index)
196  copy_indices[level].emplace_back(
197  owned_dofs.nth_index_in_set(i), unrolled_copy_indices[i]);
198  }
199  }
200 
201  const ::parallel::TriangulationBase<dim, spacedim> *tria =
202  (dynamic_cast<const ::parallel::TriangulationBase<dim, spacedim>
203  *>(&dof_handler.get_triangulation()));
204  AssertThrow(
205  send_data_temp.empty() || tria != nullptr,
206  ExcMessage(
207  "We should only be sending information with a parallel Triangulation!"));
208 
209 #ifdef DEAL_II_WITH_MPI
210  if (tria && Utilities::MPI::sum(send_data_temp.size(),
211  tria->get_communicator()) > 0)
212  {
213  const std::set<types::subdomain_id> &neighbors =
214  tria->level_ghost_owners();
215  std::map<int, std::vector<DoFPair>> send_data;
216 
217  std::sort(send_data_temp.begin(),
218  send_data_temp.end(),
219  [](const DoFPair &lhs, const DoFPair &rhs) {
220  if (lhs.level < rhs.level)
221  return true;
222  if (lhs.level > rhs.level)
223  return false;
224 
225  if (lhs.level_dof_index < rhs.level_dof_index)
226  return true;
227  if (lhs.level_dof_index > rhs.level_dof_index)
228  return false;
229 
230  if (lhs.global_dof_index < rhs.global_dof_index)
231  return true;
232  else
233  return false;
234  });
235  send_data_temp.erase(
236  std::unique(send_data_temp.begin(),
237  send_data_temp.end(),
238  [](const DoFPair &lhs, const DoFPair &rhs) {
239  return (lhs.level == rhs.level) &&
240  (lhs.level_dof_index == rhs.level_dof_index) &&
241  (lhs.global_dof_index == rhs.global_dof_index);
242  }),
243  send_data_temp.end());
244 
245  for (unsigned int level = 0; level < n_levels; ++level)
246  {
247  const IndexSet &owned_level_dofs =
248  dof_handler.locally_owned_mg_dofs(level);
249 
250  std::vector<types::global_dof_index> level_dof_indices;
251  std::vector<types::global_dof_index> global_dof_indices;
252  for (const auto &dofpair : send_data_temp)
253  if (dofpair.level == level)
254  {
255  level_dof_indices.push_back(dofpair.level_dof_index);
256  global_dof_indices.push_back(dofpair.global_dof_index);
257  }
258 
259  IndexSet is_ghost(owned_level_dofs.size());
260  is_ghost.add_indices(level_dof_indices.begin(),
261  level_dof_indices.end());
262 
263  AssertThrow(level_dof_indices.size() == is_ghost.n_elements(),
264  ExcMessage("Size does not match!"));
265 
266  const auto index_owner =
267  Utilities::MPI::compute_index_owner(owned_level_dofs,
268  is_ghost,
270 
271  AssertThrow(level_dof_indices.size() == index_owner.size(),
272  ExcMessage("Size does not match!"));
273 
274  for (unsigned int i = 0; i < index_owner.size(); ++i)
275  send_data[index_owner[i]].emplace_back(level,
276  global_dof_indices[i],
277  level_dof_indices[i]);
278  }
279 
280 
281  // Protect the send/recv logic with a mutex:
282  static Utilities::MPI::CollectiveMutex mutex;
284  mutex, tria->get_communicator());
285 
286  const int mpi_tag =
288 
289  // * send
290  std::vector<MPI_Request> requests;
291  {
292  for (const auto dest : neighbors)
293  {
294  requests.push_back(MPI_Request());
295  std::vector<DoFPair> &data = send_data[dest];
296 
297  const int ierr =
298  MPI_Isend(data.data(),
299  data.size() * sizeof(decltype(*data.data())),
300  MPI_BYTE,
301  dest,
302  mpi_tag,
304  &*requests.rbegin());
305  AssertThrowMPI(ierr);
306  }
307  }
308 
309  // * receive
310  {
311  // We should get one message from each of our neighbors
312  std::vector<DoFPair> receive_buffer;
313  for (unsigned int counter = 0; counter < neighbors.size();
314  ++counter)
315  {
316  MPI_Status status;
317  int ierr = MPI_Probe(MPI_ANY_SOURCE,
318  mpi_tag,
320  &status);
321  AssertThrowMPI(ierr);
322  int len;
323  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
324  AssertThrowMPI(ierr);
325 
326  if (len == 0)
327  {
328  ierr = MPI_Recv(nullptr,
329  0,
330  MPI_BYTE,
331  status.MPI_SOURCE,
332  status.MPI_TAG,
334  &status);
335  AssertThrowMPI(ierr);
336  continue;
337  }
338 
339  int count = len / sizeof(DoFPair);
340  Assert(static_cast<int>(count * sizeof(DoFPair)) == len,
341  ExcInternalError());
342  receive_buffer.resize(count);
343 
344  void *ptr = receive_buffer.data();
345  ierr = MPI_Recv(ptr,
346  len,
347  MPI_BYTE,
348  status.MPI_SOURCE,
349  status.MPI_TAG,
351  &status);
352  AssertThrowMPI(ierr);
353 
354  for (const auto &dof_pair : receive_buffer)
355  {
356  copy_indices_level_mine[dof_pair.level].emplace_back(
357  dof_pair.global_dof_index, dof_pair.level_dof_index);
358  }
359  }
360  }
361 
362  // * wait for all MPI_Isend to complete
363  if (requests.size() > 0)
364  {
365  const int ierr = MPI_Waitall(requests.size(),
366  requests.data(),
367  MPI_STATUSES_IGNORE);
368  AssertThrowMPI(ierr);
369  requests.clear();
370  }
371 # ifdef DEBUG
372  // Make sure in debug mode, that everybody sent/received all packages
373  // on this level. If a deadlock occurs here, the list of expected
374  // senders is not computed correctly.
375  const int ierr = MPI_Barrier(tria->get_communicator());
376  AssertThrowMPI(ierr);
377 # endif
378  }
379 #endif
380 
381  // Sort the indices, except the copy_indices which already are
382  // sorted. This will produce more reliable debug output for regression
383  // tests and won't hurt performance even in release mode because the
384  // non-owned indices are a small subset of all unknowns.
385  std::less<std::pair<types::global_dof_index, types::global_dof_index>>
386  compare;
387  for (auto &level_indices : copy_indices_level_mine)
388  std::sort(level_indices.begin(), level_indices.end(), compare);
389  for (auto &level_indices : copy_indices_global_mine)
390  std::sort(level_indices.begin(), level_indices.end(), compare);
391  }
392 
393 
394 
395  // initialize the vectors needed for the transfer (and merge with the
396  // content in copy_indices_global_mine)
397  void
399  const IndexSet &locally_owned,
400  std::vector<types::global_dof_index> &ghosted_level_dofs,
401  const std::shared_ptr<const Utilities::MPI::Partitioner>
402  &external_partitioner,
403  const MPI_Comm communicator,
404  std::shared_ptr<const Utilities::MPI::Partitioner> &target_partitioner,
405  Table<2, unsigned int> &copy_indices_global_mine)
406  {
407  std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end());
408  IndexSet ghosted_dofs(locally_owned.size());
409  ghosted_dofs.add_indices(ghosted_level_dofs.begin(),
410  ghosted_level_dofs.end());
411  ghosted_dofs.compress();
412 
413  // Add possible ghosts from the previous content in the vector
414  if (target_partitioner.get() != nullptr &&
415  target_partitioner->size() == locally_owned.size())
416  {
417  ghosted_dofs.add_indices(target_partitioner->ghost_indices());
418  }
419 
420  // check if the given partitioner's ghosts represent a superset of the
421  // ghosts we require in this function
422  const int ghosts_locally_contained =
423  (external_partitioner.get() != nullptr &&
424  (external_partitioner->ghost_indices() & ghosted_dofs) ==
425  ghosted_dofs) ?
426  1 :
427  0;
428  if (external_partitioner.get() != nullptr &&
429  Utilities::MPI::min(ghosts_locally_contained, communicator) == 1)
430  {
431  // shift the local number of the copy indices according to the new
432  // partitioner that we are going to use during the access to the
433  // entries
434  if (target_partitioner.get() != nullptr &&
435  target_partitioner->size() == locally_owned.size())
436  for (unsigned int i = 0; i < copy_indices_global_mine.n_cols(); ++i)
437  copy_indices_global_mine(1, i) =
438  external_partitioner->global_to_local(
439  target_partitioner->local_to_global(
440  copy_indices_global_mine(1, i)));
441  target_partitioner = external_partitioner;
442  }
443  else
444  {
445  if (target_partitioner.get() != nullptr &&
446  target_partitioner->size() == locally_owned.size())
447  for (unsigned int i = 0; i < copy_indices_global_mine.n_cols(); ++i)
448  copy_indices_global_mine(1, i) =
449  locally_owned.n_elements() +
450  ghosted_dofs.index_within_set(
451  target_partitioner->local_to_global(
452  copy_indices_global_mine(1, i)));
453  target_partitioner =
454  std::make_shared<Utilities::MPI::Partitioner>(locally_owned,
455  ghosted_dofs,
456  communicator);
457  }
458  }
459 
460 
461 
462  // Transform the ghost indices to local index space for the vector
463  inline void
465  const Utilities::MPI::Partitioner &part,
466  const std::vector<types::global_dof_index> &mine,
467  const std::vector<types::global_dof_index> &remote,
468  std::vector<unsigned int> &localized_indices)
469  {
470  localized_indices.resize(mine.size() + remote.size(),
472  for (unsigned int i = 0; i < mine.size(); ++i)
473  if (mine[i] != numbers::invalid_dof_index)
474  localized_indices[i] = part.global_to_local(mine[i]);
475 
476  for (unsigned int i = 0; i < remote.size(); ++i)
477  if (remote[i] != numbers::invalid_dof_index)
478  localized_indices[i + mine.size()] = part.global_to_local(remote[i]);
479  }
480 
481 
482 
483  // given the collection of child cells in lexicographic ordering as seen
484  // from the parent, compute the first index of the given child
485  template <int dim>
486  unsigned int
487  compute_shift_within_children(const unsigned int child,
488  const unsigned int fe_shift_1d,
489  const unsigned int fe_degree)
490  {
491  // we put the degrees of freedom of all child cells in lexicographic
492  // ordering
493  unsigned int c_tensor_index[dim];
494  unsigned int tmp = child;
495  for (unsigned int d = 0; d < dim; ++d)
496  {
497  c_tensor_index[d] = tmp % 2;
498  tmp /= 2;
499  }
500  const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
501  unsigned int factor = 1;
502  unsigned int shift = fe_shift_1d * c_tensor_index[0];
503  for (unsigned int d = 1; d < dim; ++d)
504  {
505  factor *= n_child_dofs_1d;
506  shift = shift + factor * fe_shift_1d * c_tensor_index[d];
507  }
508  return shift;
509  }
510 
511 
512 
513  // puts the indices on the given child cell in lexicographic ordering with
514  // respect to the collection of all child cells as seen from the parent
515  template <int dim>
516  void
518  const unsigned int child,
519  const unsigned int fe_shift_1d,
520  const unsigned int fe_degree,
521  const std::vector<unsigned int> &lexicographic_numbering,
522  const std::vector<types::global_dof_index> &local_dof_indices,
523  types::global_dof_index *target_indices)
524  {
525  const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
526  const unsigned int shift =
527  compute_shift_within_children<dim>(child, fe_shift_1d, fe_degree);
528  const unsigned int n_components =
529  local_dof_indices.size() / Utilities::fixed_power<dim>(fe_degree + 1);
530  types::global_dof_index *indices = target_indices + shift;
531  const unsigned int n_scalar_cell_dofs =
532  Utilities::fixed_power<dim>(n_child_dofs_1d);
533  for (unsigned int c = 0, m = 0; c < n_components; ++c)
534  for (unsigned int k = 0; k < (dim > 2 ? (fe_degree + 1) : 1); ++k)
535  for (unsigned int j = 0; j < (dim > 1 ? (fe_degree + 1) : 1); ++j)
536  for (unsigned int i = 0; i < (fe_degree + 1); ++i, ++m)
537  {
538  const unsigned int index =
539  c * n_scalar_cell_dofs +
540  k * n_child_dofs_1d * n_child_dofs_1d + j * n_child_dofs_1d +
541  i;
543  indices[index] ==
544  local_dof_indices[lexicographic_numbering[m]],
545  ExcInternalError());
546  indices[index] = local_dof_indices[lexicographic_numbering[m]];
547  }
548  }
549 
550 
551 
552  template <int dim, typename Number>
553  void
555  const FiniteElement<1> &fe,
556  const DoFHandler<dim> &dof_handler)
557  {
558  // currently, we have only FE_Q and FE_DGQ type elements implemented
559  elem_info.n_components = dof_handler.get_fe().element_multiplicity(0);
560  AssertDimension(Utilities::fixed_power<dim>(fe.n_dofs_per_cell()) *
561  elem_info.n_components,
562  dof_handler.get_fe().n_dofs_per_cell());
563  AssertDimension(fe.degree, dof_handler.get_fe().degree);
564  elem_info.fe_degree = fe.degree;
565  elem_info.element_is_continuous = fe.n_dofs_per_vertex() > 0;
567 
568  // step 1.2: get renumbering of 1d basis functions to lexicographic
569  // numbers. The distinction according to fe.n_dofs_per_vertex() is to
570  // support both continuous and discontinuous bases.
571  std::vector<unsigned int> renumbering(fe.n_dofs_per_cell());
572  {
574  renumbering[0] = 0;
575  for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
576  renumbering[i + fe.n_dofs_per_vertex()] =
578  if (fe.n_dofs_per_vertex() > 0)
579  renumbering[fe.n_dofs_per_cell() - fe.n_dofs_per_vertex()] =
580  fe.n_dofs_per_vertex();
581  }
582 
583  // step 1.3: create a dummy 1d quadrature formula to extract the
584  // lexicographic numbering for the elements
585  Assert(fe.n_dofs_per_vertex() == 0 || fe.n_dofs_per_vertex() == 1,
587  const unsigned int shift = fe.n_dofs_per_cell() - fe.n_dofs_per_vertex();
588  const unsigned int n_child_dofs_1d =
589  (fe.n_dofs_per_vertex() > 0 ? (2 * fe.n_dofs_per_cell() - 1) :
590  (2 * fe.n_dofs_per_cell()));
591 
592  elem_info.n_child_cell_dofs =
593  elem_info.n_components * Utilities::fixed_power<dim>(n_child_dofs_1d);
594  const Quadrature<1> dummy_quadrature(
595  std::vector<Point<1>>(1, Point<1>()));
597  shape_info.reinit(dummy_quadrature, dof_handler.get_fe(), 0);
598  elem_info.lexicographic_numbering = shape_info.lexicographic_numbering;
599 
600  // step 1.4: get the 1d prolongation matrix and combine from both children
601  elem_info.prolongation_matrix_1d.resize(fe.n_dofs_per_cell() *
602  n_child_dofs_1d);
603 
604  for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; ++c)
605  for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
606  for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
607  elem_info
608  .prolongation_matrix_1d[i * n_child_dofs_1d + j + c * shift] =
609  fe.get_prolongation_matrix(c)(renumbering[j], renumbering[i]);
610  }
611 
612 
613 
614  // Sets up most of the internal data structures of the MGTransferMatrixFree
615  // class
616  template <int dim, typename Number>
617  void
619  const DoFHandler<dim> &dof_handler,
620  const MGConstrainedDoFs *mg_constrained_dofs,
621  const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
622  &external_partitioners,
623  ElementInfo<Number> &elem_info,
624  std::vector<std::vector<unsigned int>> &level_dof_indices,
625  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
626  &parent_child_connect,
627  std::vector<unsigned int> &n_owned_level_cells,
628  std::vector<std::vector<std::vector<unsigned short>>> &dirichlet_indices,
629  std::vector<std::vector<Number>> &weights_on_refined,
630  std::vector<Table<2, unsigned int>> &copy_indices_global_mine,
631  MGLevelObject<std::shared_ptr<const Utilities::MPI::Partitioner>>
632  &target_partitioners)
633  {
634  level_dof_indices.clear();
635  parent_child_connect.clear();
636  n_owned_level_cells.clear();
637  dirichlet_indices.clear();
638  weights_on_refined.clear();
639 
640 #ifdef DEBUG
641  if (mg_constrained_dofs)
642  {
643  const unsigned int n_levels =
644  dof_handler.get_triangulation().n_global_levels();
645 
646  for (unsigned int l = 0; l < n_levels; ++l)
647  {
648  const auto &constraints =
649  mg_constrained_dofs->get_user_constraint_matrix(l);
650 
651  // no inhomogeneities are supported
652  AssertDimension(constraints.n_inhomogeneities(), 0);
653 
654  for (const auto dof : constraints.get_local_lines())
655  {
656  const auto *entries_ptr =
657  constraints.get_constraint_entries(dof);
658 
659  if (entries_ptr == nullptr)
660  continue;
661 
662  // only homogeneous or identity constraints are supported
663  Assert((entries_ptr->size() == 0) ||
664  ((entries_ptr->size() == 1) &&
665  (std::abs((*entries_ptr)[0].second - 1.) <
668  }
669  }
670  }
671 #endif
672 
673  // we collect all child DoFs of a mother cell together. For faster
674  // tensorized operations, we align the degrees of freedom
675  // lexicographically. We distinguish FE_Q elements and FE_DGQ elements
676 
677  const ::Triangulation<dim> &tria = dof_handler.get_triangulation();
678 
679  // ---------------------------- 1. Extract 1d info about the finite
680  // element step 1.1: create a 1d copy of the finite element from FETools
681  // where we substitute the template argument
682  AssertDimension(dof_handler.get_fe().n_base_elements(), 1);
683  std::string fe_name = dof_handler.get_fe().base_element(0).get_name();
684  {
685  const std::size_t template_starts = fe_name.find_first_of('<');
686  Assert(fe_name[template_starts + 1] ==
687  (dim == 1 ? '1' : (dim == 2 ? '2' : '3')),
688  ExcInternalError());
689  fe_name[template_starts + 1] = '1';
690  }
691  const std::unique_ptr<FiniteElement<1>> fe(
692  FETools::get_fe_by_name<1, 1>(fe_name));
693 
694  setup_element_info(elem_info, *fe, dof_handler);
695 
696 
697  // ---------- 2. Extract and match dof indices between child and parent
698  const unsigned int n_levels = tria.n_global_levels();
699  level_dof_indices.resize(n_levels);
700  parent_child_connect.resize(n_levels - 1);
701  n_owned_level_cells.resize(n_levels - 1);
702  std::vector<std::vector<unsigned int>> coarse_level_indices(n_levels - 1);
703  for (unsigned int level = 0;
704  level < std::min(tria.n_levels(), n_levels - 1);
705  ++level)
706  coarse_level_indices[level].resize(tria.n_raw_cells(level),
708  std::vector<types::global_dof_index> local_dof_indices(
709  dof_handler.get_fe().n_dofs_per_cell());
710  dirichlet_indices.resize(n_levels - 1);
711 
712  AssertDimension(target_partitioners.max_level(), n_levels - 1);
713  Assert(external_partitioners.empty() ||
714  external_partitioners.size() == n_levels,
715  ExcDimensionMismatch(external_partitioners.size(), n_levels));
716 
717  for (unsigned int level = n_levels - 1; level > 0; --level)
718  {
719  unsigned int counter = 0;
720  std::vector<types::global_dof_index> global_level_dof_indices;
721  std::vector<types::global_dof_index> global_level_dof_indices_remote;
722  std::vector<types::global_dof_index> ghosted_level_dofs;
723  std::vector<types::global_dof_index> global_level_dof_indices_l0;
724  std::vector<types::global_dof_index> ghosted_level_dofs_l0;
725 
726  // step 2.1: loop over the cells on the coarse side
727  typename DoFHandler<dim>::cell_iterator cell,
728  endc = dof_handler.end(level - 1);
729  for (cell = dof_handler.begin(level - 1); cell != endc; ++cell)
730  {
731  // need to look into a cell if it has children and it is locally
732  // owned
733  if (!cell->has_children())
734  continue;
735 
736  bool consider_cell =
739  cell->level_subdomain_id() == tria.locally_owned_subdomain());
740 
741  // due to the particular way we store DoF indices (via children),
742  // we also need to add the DoF indices for coarse cells where we
743  // own at least one child
744  const bool cell_is_remote = !consider_cell;
745  for (unsigned int c = 0;
746  c < GeometryInfo<dim>::max_children_per_cell;
747  ++c)
748  if (cell->child(c)->level_subdomain_id() ==
750  {
751  consider_cell = true;
752  break;
753  }
754 
755  if (!consider_cell)
756  continue;
757 
758  // step 2.2: loop through children and append the dof indices to
759  // the appropriate list. We need separate lists for the owned
760  // coarse cell case (which will be part of
761  // restriction/prolongation between level-1 and level) and the
762  // remote case (which needs to store DoF indices for the
763  // operations between level and level+1).
764  AssertDimension(cell->n_children(),
766  std::vector<types::global_dof_index> &next_indices =
767  cell_is_remote ? global_level_dof_indices_remote :
768  global_level_dof_indices;
769  const std::size_t start_index = next_indices.size();
770  next_indices.resize(start_index + elem_info.n_child_cell_dofs,
772  for (unsigned int c = 0;
773  c < GeometryInfo<dim>::max_children_per_cell;
774  ++c)
775  {
776  if (cell_is_remote && cell->child(c)->level_subdomain_id() !=
778  continue;
779  cell->child(c)->get_mg_dof_indices(local_dof_indices);
780 
781  resolve_identity_constraints(mg_constrained_dofs,
782  level,
783  local_dof_indices);
784 
785  const IndexSet &owned_level_dofs =
786  dof_handler.locally_owned_mg_dofs(level);
787  for (const auto local_dof_index : local_dof_indices)
788  if (!owned_level_dofs.is_element(local_dof_index))
789  ghosted_level_dofs.push_back(local_dof_index);
790 
791  add_child_indices<dim>(c,
792  fe->n_dofs_per_cell() -
793  fe->n_dofs_per_vertex(),
794  fe->degree,
795  elem_info.lexicographic_numbering,
796  local_dof_indices,
797  &next_indices[start_index]);
798 
799  // step 2.3 store the connectivity to the parent
800  if (cell->child(c)->has_children() &&
803  cell->child(c)->level_subdomain_id() ==
805  {
806  const unsigned int child_index =
807  coarse_level_indices[level][cell->child(c)->index()];
808  AssertIndexRange(child_index,
809  parent_child_connect[level].size());
810  unsigned int parent_index = counter;
811  // remote cells, i.e., cells where we work on a further
812  // level but are not treated on the current level, need to
813  // be placed at the end of the list; however, we do not
814  // yet know the exact position in the array, so shift
815  // their parent index by the number of cells so we can set
816  // the correct number after the end of this loop
817  if (cell_is_remote)
818  parent_index =
819  start_index / elem_info.n_child_cell_dofs +
820  tria.n_cells(level);
821  parent_child_connect[level][child_index] =
822  std::make_pair(parent_index, c);
823  AssertIndexRange(dof_handler.get_fe().n_dofs_per_cell(),
824  static_cast<unsigned short>(-1));
825 
826  // set Dirichlet boundary conditions (as a list of
827  // constrained DoFs) for the child
828  if (mg_constrained_dofs != nullptr)
829  for (unsigned int i = 0;
830  i < dof_handler.get_fe().n_dofs_per_cell();
831  ++i)
832  if (mg_constrained_dofs->is_boundary_index(
833  level,
834  local_dof_indices
835  [elem_info.lexicographic_numbering[i]]))
836  dirichlet_indices[level][child_index].push_back(i);
837  }
838  }
839  if (!cell_is_remote)
840  {
841  AssertIndexRange(static_cast<unsigned int>(cell->index()),
842  coarse_level_indices[level - 1].size());
843  coarse_level_indices[level - 1][cell->index()] = counter++;
844  }
845 
846  // step 2.4: include indices for the coarsest cells. we still
847  // insert the indices as if they were from a child in order to use
848  // the same code (the coarsest level does not matter much in terms
849  // of memory, so we gain in code simplicity)
850  if (level == 1 && !cell_is_remote)
851  {
852  cell->get_mg_dof_indices(local_dof_indices);
853 
854  resolve_identity_constraints(mg_constrained_dofs,
855  level - 1,
856  local_dof_indices);
857 
858  const IndexSet &owned_level_dofs_l0 =
859  dof_handler.locally_owned_mg_dofs(0);
860  for (const auto local_dof_index : local_dof_indices)
861  if (!owned_level_dofs_l0.is_element(local_dof_index))
862  ghosted_level_dofs_l0.push_back(local_dof_index);
863 
864  const std::size_t start_index =
865  global_level_dof_indices_l0.size();
866  global_level_dof_indices_l0.resize(
867  start_index + elem_info.n_child_cell_dofs,
869  add_child_indices<dim>(
870  0,
871  fe->n_dofs_per_cell() - fe->n_dofs_per_vertex(),
872  fe->degree,
873  elem_info.lexicographic_numbering,
874  local_dof_indices,
875  &global_level_dof_indices_l0[start_index]);
876 
877  dirichlet_indices[0].emplace_back();
878  if (mg_constrained_dofs != nullptr)
879  for (unsigned int i = 0;
880  i < dof_handler.get_fe().n_dofs_per_cell();
881  ++i)
882  if (mg_constrained_dofs->is_boundary_index(
883  0,
884  local_dof_indices[elem_info
886  dirichlet_indices[0].back().push_back(i);
887  }
888  }
889 
890  // step 2.5: store information about the current level and prepare the
891  // Dirichlet indices and parent-child relationship for the next
892  // coarser level
893  AssertDimension(counter * elem_info.n_child_cell_dofs,
894  global_level_dof_indices.size());
895  n_owned_level_cells[level - 1] = counter;
896  dirichlet_indices[level - 1].resize(counter);
897  parent_child_connect[level - 1].resize(
898  counter,
899  std::make_pair(numbers::invalid_unsigned_int,
901 
902  // step 2.6: put the cells with remotely owned parent to the end of
903  // the list (these are needed for the transfer from level to level+1
904  // but not for the transfer from level-1 to level).
905  if (level < n_levels - 1)
906  for (std::vector<std::pair<unsigned int, unsigned int>>::iterator
907  i = parent_child_connect[level].begin();
908  i != parent_child_connect[level].end();
909  ++i)
910  if (i->first >= tria.n_cells(level))
911  {
912  i->first -= tria.n_cells(level);
913  i->first += counter;
914  }
915 
916  // step 2.7: Initialize the partitioner for the ghosted vector
917  //
918  // We use a vector based on the target partitioner handed in also in
919  // the base class for keeping ghosted transfer indices. To avoid
920  // keeping two very similar vectors, we keep one single ghosted
921  // vector that is augmented/filled here.
923  ghosted_level_dofs,
924  external_partitioners.empty() ?
925  nullptr :
926  external_partitioners[level],
928  target_partitioners[level],
929  copy_indices_global_mine[level]);
930 
931  copy_indices_to_mpi_local_numbers(*target_partitioners[level],
932  global_level_dof_indices,
933  global_level_dof_indices_remote,
934  level_dof_indices[level]);
935  // step 2.8: Initialize the ghosted vector for level 0
936  if (level == 1)
937  {
938  for (unsigned int i = 0; i < parent_child_connect[0].size(); ++i)
939  parent_child_connect[0][i] = std::make_pair(i, 0U);
940 
942  ghosted_level_dofs_l0,
943  external_partitioners.empty() ?
944  nullptr :
945  external_partitioners[0],
947  target_partitioners[0],
948  copy_indices_global_mine[0]);
949 
951  *target_partitioners[0],
952  global_level_dof_indices_l0,
953  std::vector<types::global_dof_index>(),
954  level_dof_indices[0]);
955  }
956  }
957 
958  // ---------------------- 3. compute weights to make restriction additive
959 
960  const unsigned int n_child_dofs_1d =
961  fe->degree + 1 + fe->n_dofs_per_cell() - fe->n_dofs_per_vertex();
962 
963  // get the valence of the individual components and compute the weights as
964  // the inverse of the valence
965  weights_on_refined.resize(n_levels - 1);
966  for (unsigned int level = 1; level < n_levels; ++level)
967  {
969  target_partitioners[level]);
970  for (unsigned int c = 0; c < n_owned_level_cells[level - 1]; ++c)
971  for (unsigned int j = 0; j < elem_info.n_child_cell_dofs; ++j)
972  touch_count.local_element(
973  level_dof_indices[level][elem_info.n_child_cell_dofs * c +
974  j]) += Number(1.);
975  touch_count.compress(VectorOperation::add);
976  touch_count.update_ghost_values();
977 
978  std::vector<unsigned int> degree_to_3(n_child_dofs_1d);
979  degree_to_3[0] = 0;
980  for (unsigned int i = 1; i < n_child_dofs_1d - 1; ++i)
981  degree_to_3[i] = 1;
982  degree_to_3.back() = 2;
983 
984  // we only store 3^dim weights because all dofs on a line have the
985  // same valence, and all dofs on a quad have the same valence.
986  weights_on_refined[level - 1].resize(n_owned_level_cells[level - 1] *
987  Utilities::fixed_power<dim>(3));
988  for (unsigned int c = 0; c < n_owned_level_cells[level - 1]; ++c)
989  for (unsigned int k = 0, m = 0; k < (dim > 2 ? n_child_dofs_1d : 1);
990  ++k)
991  for (unsigned int j = 0; j < (dim > 1 ? n_child_dofs_1d : 1); ++j)
992  {
993  unsigned int shift = 9 * degree_to_3[k] + 3 * degree_to_3[j];
994  for (unsigned int i = 0; i < n_child_dofs_1d; ++i, ++m)
995  weights_on_refined[level -
996  1][c * Utilities::fixed_power<dim>(3) +
997  shift + degree_to_3[i]] =
998  Number(1.) /
999  touch_count.local_element(
1000  level_dof_indices[level]
1001  [elem_info.n_child_cell_dofs * c + m]);
1002  }
1003  }
1004  }
1005 
1006 
1007 
1008  void
1010  const MGConstrainedDoFs *mg_constrained_dofs,
1011  const unsigned int level,
1012  std::vector<types::global_dof_index> &dof_indices)
1013  {
1014  if (mg_constrained_dofs != nullptr &&
1015  mg_constrained_dofs->get_level_constraints(level).n_constraints() > 0)
1016  for (auto &ind : dof_indices)
1017  if (mg_constrained_dofs->get_level_constraints(level)
1019  {
1020  Assert(mg_constrained_dofs->get_level_constraints(level)
1022  ->size() == 1,
1023  ExcInternalError());
1024  ind = mg_constrained_dofs->get_level_constraints(level)
1026  ->front()
1027  .first;
1028  }
1029  }
1030 
1031  } // namespace MGTransfer
1032 } // namespace internal
1033 
1034 // Explicit instantiations
1035 
1036 #include "mg_transfer_internal.inst"
1037 
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
bool is_identity_constrained(const size_type line_n) const
size_type n_constraints() const
cell_iterator end() const
const Triangulation< dim, spacedim > & get_triangulation() const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
cell_iterator begin(const unsigned int level=0) const
const IndexSet & locally_owned_dofs() const
unsigned int n_dofs_per_vertex() const
const unsigned int degree
Definition: fe_data.h:453
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
virtual std::string get_name() const =0
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
unsigned int element_multiplicity(const unsigned int index) const
unsigned int n_base_elements() const
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
size_type size() const
Definition: index_set.h:1761
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1986
size_type n_elements() const
Definition: index_set.h:1919
bool is_element(const size_type index) const
Definition: index_set.h:1879
void clear()
Definition: index_set.h:1737
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1967
void compress() const
Definition: index_set.h:1769
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1816
Number local_element(const size_type local_index) const
void compress(VectorOperation::values operation)
bool at_refinement_edge(const unsigned int level, const types::global_dof_index index) const
bool is_boundary_index(const unsigned int level, const types::global_dof_index index) const
const AffineConstraints< double > & get_user_constraint_matrix(const unsigned int level) const
const AffineConstraints< double > & get_level_constraints(const unsigned int level) const
virtual MPI_Comm get_communicator() const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_levels() const
unsigned int n_raw_cells(const unsigned int level) const
virtual unsigned int n_global_levels() const
unsigned int n_cells() const
unsigned int global_to_local(const types::global_dof_index global_index) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
unsigned int level
Definition: grid_out.cc:4617
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1947
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:469
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:191
static const char U
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
@ mg_transfer_fill_copy_indices
mg_transfer_internal.cc: fill_copy_indices()
Definition: mpi_tags.h:72
T sum(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
std::vector< unsigned int > compute_index_owner(const IndexSet &owned_indices, const IndexSet &indices_to_look_up, const MPI_Comm comm)
Definition: mpi.cc:1089
void setup_element_info(ElementInfo< Number > &elem_info, const FiniteElement< 1 > &fe, const DoFHandler< dim > &dof_handler)
void resolve_identity_constraints(const MGConstrainedDoFs *mg_constrained_dofs, const unsigned int level, std::vector< types::global_dof_index > &dof_indices)
unsigned int compute_shift_within_children(const unsigned int child, const unsigned int fe_shift_1d, const unsigned int fe_degree)
void add_child_indices(const unsigned int child, const unsigned int fe_shift_1d, const unsigned int fe_degree, const std::vector< unsigned int > &lexicographic_numbering, const std::vector< types::global_dof_index > &local_dof_indices, types::global_dof_index *target_indices)
void setup_transfer(const DoFHandler< dim > &dof_handler, const MGConstrainedDoFs *mg_constrained_dofs, const std::vector< std::shared_ptr< const Utilities::MPI::Partitioner >> &external_partitioners, ElementInfo< Number > &elem_info, std::vector< std::vector< unsigned int >> &level_dof_indices, std::vector< std::vector< std::pair< unsigned int, unsigned int >>> &parent_child_connect, std::vector< unsigned int > &n_owned_level_cells, std::vector< std::vector< std::vector< unsigned short >>> &dirichlet_indices, std::vector< std::vector< Number >> &weights_on_refined, std::vector< Table< 2, unsigned int >> &copy_indices_global_mine, MGLevelObject< std::shared_ptr< const Utilities::MPI::Partitioner >> &vector_partitioners)
void reinit_level_partitioner(const IndexSet &locally_owned, std::vector< types::global_dof_index > &ghosted_level_dofs, const std::shared_ptr< const Utilities::MPI::Partitioner > &external_partitioner, const MPI_Comm communicator, std::shared_ptr< const Utilities::MPI::Partitioner > &target_partitioner, Table< 2, unsigned int > &copy_indices_global_mine)
void copy_indices_to_mpi_local_numbers(const Utilities::MPI::Partitioner &part, const std::vector< types::global_dof_index > &mine, const std::vector< types::global_dof_index > &remote, std::vector< unsigned int > &localized_indices)
void fill_copy_indices(const DoFHandler< dim, spacedim > &dof_handler, const MGConstrainedDoFs *mg_constrained_dofs, std::vector< std::vector< std::pair< types::global_dof_index, types::global_dof_index >>> &copy_indices, std::vector< std::vector< std::pair< types::global_dof_index, types::global_dof_index >>> &copy_indices_global_mine, std::vector< std::vector< std::pair< types::global_dof_index, types::global_dof_index >>> &copy_indices_level_mine, const bool skip_interface_dofs=true)
const types::subdomain_id artificial_subdomain_id
Definition: types.h:363
const types::subdomain_id invalid_subdomain_id
Definition: types.h:342
static const unsigned int invalid_unsigned_int
Definition: types.h:221
const types::global_dof_index invalid_dof_index
Definition: types.h:253
unsigned int global_dof_index
Definition: types.h:82
types::global_dof_index global_dof_index
types::global_dof_index level_dof_index
DoFPair(const unsigned int level, const types::global_dof_index global_dof_index, const types::global_dof_index level_dof_index)
std::vector< Number > prolongation_matrix_1d
std::vector< unsigned int > lexicographic_numbering
void reinit(const Quadrature< dim_q > &quad, const FiniteElement< dim, spacedim > &fe_dim, const unsigned int base_element=0)
std::vector< unsigned int > lexicographic_numbering
Definition: shape_info.h:480
const ::Triangulation< dim, spacedim > & tria