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