Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.1.1
\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
mg_transfer_internal.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2019 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/distributed/tria.h>
18 
19 #include <deal.II/dofs/dof_tools.h>
20 
21 #include <deal.II/fe/fe_tools.h>
22 
23 #include <deal.II/matrix_free/shape_info.h>
24 
25 #include <deal.II/multigrid/mg_transfer_internal.h>
26 
27 DEAL_II_NAMESPACE_OPEN
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;
40  types::global_dof_index level_dof_index;
41 
42  DoFPair(const unsigned int level,
43  const types::global_dof_index global_dof_index,
44  const types::global_dof_index level_dof_index)
45  : level(level)
47  , level_dof_index(level_dof_index)
48  {}
49 
50  DoFPair()
51  : level(numbers::invalid_unsigned_int)
53  , level_dof_index(numbers::invalid_dof_index)
54  {}
55  };
56 
57 
58 
59  template <int dim, int spacedim>
60  void
61  fill_copy_indices(
62  const ::DoFHandler<dim, spacedim> &mg_dof,
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  mg_dof.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  IndexSet globally_relevant;
94  DoFTools::extract_locally_relevant_dofs(mg_dof, globally_relevant);
95 
96  const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
97  std::vector<types::global_dof_index> global_dof_indices(dofs_per_cell);
98  std::vector<types::global_dof_index> level_dof_indices(dofs_per_cell);
99 
100  for (unsigned int level = 0; level < n_levels; ++level)
101  {
102  std::vector<bool> dof_touched(globally_relevant.n_elements(), false);
103  copy_indices[level].clear();
104  copy_indices_level_mine[level].clear();
105  copy_indices_global_mine[level].clear();
106 
107  typename ::DoFHandler<dim, spacedim>::active_cell_iterator
108  level_cell = mg_dof.begin_active(level);
109  const typename ::DoFHandler<dim, spacedim>::active_cell_iterator
110  level_end = mg_dof.end_active(level);
111 
112  for (; level_cell != level_end; ++level_cell)
113  {
114  if (mg_dof.get_triangulation().locally_owned_subdomain() !=
116  (level_cell->level_subdomain_id() ==
118  level_cell->subdomain_id() ==
120  continue;
121 
122  // get the dof numbers of this cell for the global and the
123  // level-wise numbering
124  level_cell->get_dof_indices(global_dof_indices);
125  level_cell->get_mg_dof_indices(level_dof_indices);
126 
127  for (unsigned int i = 0; i < dofs_per_cell; ++i)
128  {
129  // we need to ignore if the DoF is on a refinement edge
130  // (hanging node)
131  if (skip_interface_dofs && mg_constrained_dofs != nullptr &&
132  mg_constrained_dofs->at_refinement_edge(
133  level, level_dof_indices[i]))
134  continue;
135 
136  types::global_dof_index global_idx =
137  globally_relevant.index_within_set(global_dof_indices[i]);
138  // skip if we did this global dof already (on this or a
139  // coarser level)
140  if (dof_touched[global_idx])
141  continue;
142  bool global_mine = mg_dof.locally_owned_dofs().is_element(
143  global_dof_indices[i]);
144  bool level_mine =
145  mg_dof.locally_owned_mg_dofs(level).is_element(
146  level_dof_indices[i]);
147 
148 
149  if (global_mine && level_mine)
150  {
151  copy_indices[level].emplace_back(global_dof_indices[i],
152  level_dof_indices[i]);
153  }
154  else if (global_mine)
155  {
156  copy_indices_global_mine[level].emplace_back(
157  global_dof_indices[i], level_dof_indices[i]);
158 
159  // send this to the owner of the level_dof:
160  send_data_temp.emplace_back(level,
161  global_dof_indices[i],
162  level_dof_indices[i]);
163  }
164  else
165  {
166  // somebody will send those to me
167  }
168 
169  dof_touched[global_idx] = true;
170  }
171  }
172  }
173 
174  const ::parallel::Triangulation<dim, spacedim> *tria =
175  (dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
176  &mg_dof.get_triangulation()));
177  AssertThrow(
178  send_data_temp.size() == 0 || tria != nullptr,
179  ExcMessage(
180  "We should only be sending information with a parallel Triangulation!"));
181 
182 #ifdef DEAL_II_WITH_MPI
183  if (tria)
184  {
185  // TODO: Searching the owner for every single DoF becomes quite
186  // inefficient. Please fix this, Timo.
187  // The list of neighbors is symmetric (our neighbors have us as a
188  // neighbor), so we can use it to send and to know how many messages
189  // we will get.
190  const std::set<types::subdomain_id> &neighbors =
191  tria->level_ghost_owners();
192  std::map<int, std::vector<DoFPair>> send_data;
193 
194  // * find owners of the level dofs and insert into send_data
195  // accordingly
196  for (const auto &dofpair : send_data_temp)
197  {
198  std::set<types::subdomain_id>::iterator it;
199  for (it = neighbors.begin(); it != neighbors.end(); ++it)
200  {
201  if (mg_dof
202  .locally_owned_mg_dofs_per_processor(dofpair.level)[*it]
203  .is_element(dofpair.level_dof_index))
204  {
205  send_data[*it].push_back(dofpair);
206  break;
207  }
208  }
209  // Is this level dof not owned by any of our neighbors? That would
210  // certainly be a bug!
211  Assert(it != neighbors.end(),
212  ExcMessage("could not find DoF owner."));
213  }
214 
215  // * send
216  std::vector<MPI_Request> requests;
217  {
218  for (const auto dest : neighbors)
219  {
220  requests.push_back(MPI_Request());
221  std::vector<DoFPair> &data = send_data[dest];
222  // If there is nothing to send, we still need to send a message,
223  // because the receiving end will be waitng. In that case we
224  // just send an empty message.
225  if (data.size())
226  {
227  const int ierr = MPI_Isend(data.data(),
228  data.size() * sizeof(data[0]),
229  MPI_BYTE,
230  dest,
231  71,
232  tria->get_communicator(),
233  &*requests.rbegin());
234  AssertThrowMPI(ierr);
235  }
236  else
237  {
238  const int ierr = MPI_Isend(nullptr,
239  0,
240  MPI_BYTE,
241  dest,
242  71,
243  tria->get_communicator(),
244  &*requests.rbegin());
245  AssertThrowMPI(ierr);
246  }
247  }
248  }
249 
250  // * receive
251  {
252  // We should get one message from each of our neighbors
253  std::vector<DoFPair> receive_buffer;
254  for (unsigned int counter = 0; counter < neighbors.size();
255  ++counter)
256  {
257  MPI_Status status;
258  int len;
259  int ierr = MPI_Probe(MPI_ANY_SOURCE,
260  71,
261  tria->get_communicator(),
262  &status);
263  AssertThrowMPI(ierr);
264  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
265  AssertThrowMPI(ierr);
266 
267  if (len == 0)
268  {
269  ierr = MPI_Recv(nullptr,
270  0,
271  MPI_BYTE,
272  status.MPI_SOURCE,
273  status.MPI_TAG,
274  tria->get_communicator(),
275  &status);
276  AssertThrowMPI(ierr);
277  continue;
278  }
279 
280  int count = len / sizeof(DoFPair);
281  Assert(static_cast<int>(count * sizeof(DoFPair)) == len,
282  ExcInternalError());
283  receive_buffer.resize(count);
284 
285  void *ptr = receive_buffer.data();
286  ierr = MPI_Recv(ptr,
287  len,
288  MPI_BYTE,
289  status.MPI_SOURCE,
290  status.MPI_TAG,
291  tria->get_communicator(),
292  &status);
293  AssertThrowMPI(ierr);
294 
295  for (const auto &dof_pair : receive_buffer)
296  {
297  copy_indices_level_mine[dof_pair.level].emplace_back(
298  dof_pair.global_dof_index, dof_pair.level_dof_index);
299  }
300  }
301  }
302 
303  // * wait for all MPI_Isend to complete
304  if (requests.size() > 0)
305  {
306  const int ierr = MPI_Waitall(requests.size(),
307  requests.data(),
308  MPI_STATUSES_IGNORE);
309  AssertThrowMPI(ierr);
310  requests.clear();
311  }
312 # ifdef DEBUG
313  // Make sure in debug mode, that everybody sent/received all packages
314  // on this level. If a deadlock occurs here, the list of expected
315  // senders is not computed correctly.
316  const int ierr = MPI_Barrier(tria->get_communicator());
317  AssertThrowMPI(ierr);
318 # endif
319  }
320 #endif
321 
322  // Sort the indices. This will produce more reliable debug output for
323  // regression tests and likely won't hurt performance even in release
324  // mode.
325  std::less<std::pair<types::global_dof_index, types::global_dof_index>>
326  compare;
327  for (auto &level_indices : copy_indices)
328  std::sort(level_indices.begin(), level_indices.end(), compare);
329  for (auto &level_indices : copy_indices_level_mine)
330  std::sort(level_indices.begin(), level_indices.end(), compare);
331  for (auto &level_indices : copy_indices_global_mine)
332  std::sort(level_indices.begin(), level_indices.end(), compare);
333  }
334 
335 
336 
337  // initialize the vectors needed for the transfer (and merge with the
338  // content in copy_indices_global_mine)
339  template <typename Number>
340  void
341  reinit_ghosted_vector(
342  const IndexSet & locally_owned,
343  std::vector<types::global_dof_index> & ghosted_level_dofs,
344  const MPI_Comm & communicator,
345  LinearAlgebra::distributed::Vector<Number> &ghosted_level_vector,
346  std::vector<std::pair<unsigned int, unsigned int>>
347  &copy_indices_global_mine)
348  {
349  std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end());
350  IndexSet ghosted_dofs(locally_owned.size());
351  ghosted_dofs.add_indices(ghosted_level_dofs.begin(),
352  std::unique(ghosted_level_dofs.begin(),
353  ghosted_level_dofs.end()));
354  ghosted_dofs.compress();
355 
356  // Add possible ghosts from the previous content in the vector
357  if (ghosted_level_vector.size() == locally_owned.size())
358  {
359  // shift the local number of the copy indices according to the new
360  // partitioner that we are going to use for the vector
361  const auto &part = ghosted_level_vector.get_partitioner();
362  ghosted_dofs.add_indices(part->ghost_indices());
363  for (auto &indices : copy_indices_global_mine)
364  indices.second = locally_owned.n_elements() +
365  ghosted_dofs.index_within_set(
366  part->local_to_global(indices.second));
367  }
368  ghosted_level_vector.reinit(locally_owned, ghosted_dofs, communicator);
369  }
370 
371  // Transform the ghost indices to local index space for the vector
372  inline void
373  copy_indices_to_mpi_local_numbers(
374  const Utilities::MPI::Partitioner & part,
375  const std::vector<types::global_dof_index> &mine,
376  const std::vector<types::global_dof_index> &remote,
377  std::vector<unsigned int> & localized_indices)
378  {
379  localized_indices.resize(mine.size() + remote.size(),
381  for (unsigned int i = 0; i < mine.size(); ++i)
382  if (mine[i] != numbers::invalid_dof_index)
383  localized_indices[i] = part.global_to_local(mine[i]);
384 
385  for (unsigned int i = 0; i < remote.size(); ++i)
386  if (remote[i] != numbers::invalid_dof_index)
387  localized_indices[i + mine.size()] = part.global_to_local(remote[i]);
388  }
389 
390  // given the collection of child cells in lexicographic ordering as seen
391  // from the parent, compute the first index of the given child
392  template <int dim>
393  unsigned int
394  compute_shift_within_children(const unsigned int child,
395  const unsigned int fe_shift_1d,
396  const unsigned int fe_degree)
397  {
398  // we put the degrees of freedom of all child cells in lexicographic
399  // ordering
400  unsigned int c_tensor_index[dim];
401  unsigned int tmp = child;
402  for (unsigned int d = 0; d < dim; ++d)
403  {
404  c_tensor_index[d] = tmp % 2;
405  tmp /= 2;
406  }
407  const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
408  unsigned int factor = 1;
409  unsigned int shift = fe_shift_1d * c_tensor_index[0];
410  for (unsigned int d = 1; d < dim; ++d)
411  {
412  factor *= n_child_dofs_1d;
413  shift = shift + factor * fe_shift_1d * c_tensor_index[d];
414  }
415  return shift;
416  }
417 
418  // puts the indices on the given child cell in lexicographic ordering with
419  // respect to the collection of all child cells as seen from the parent
420  template <int dim>
421  void
422  add_child_indices(
423  const unsigned int child,
424  const unsigned int fe_shift_1d,
425  const unsigned int fe_degree,
426  const std::vector<unsigned int> & lexicographic_numbering,
427  const std::vector<types::global_dof_index> &local_dof_indices,
428  types::global_dof_index * target_indices)
429  {
430  const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
431  const unsigned int shift =
432  compute_shift_within_children<dim>(child, fe_shift_1d, fe_degree);
433  const unsigned int n_components =
434  local_dof_indices.size() / Utilities::fixed_power<dim>(fe_degree + 1);
435  types::global_dof_index *indices = target_indices + shift;
436  const unsigned int n_scalar_cell_dofs =
437  Utilities::fixed_power<dim>(n_child_dofs_1d);
438  for (unsigned int c = 0, m = 0; c < n_components; ++c)
439  for (unsigned int k = 0; k < (dim > 2 ? (fe_degree + 1) : 1); ++k)
440  for (unsigned int j = 0; j < (dim > 1 ? (fe_degree + 1) : 1); ++j)
441  for (unsigned int i = 0; i < (fe_degree + 1); ++i, ++m)
442  {
443  const unsigned int index =
444  c * n_scalar_cell_dofs +
445  k * n_child_dofs_1d * n_child_dofs_1d + j * n_child_dofs_1d +
446  i;
447  Assert(indices[index] == numbers::invalid_dof_index ||
448  indices[index] ==
449  local_dof_indices[lexicographic_numbering[m]],
450  ExcInternalError());
451  indices[index] = local_dof_indices[lexicographic_numbering[m]];
452  }
453  }
454 
455  template <int dim, typename Number>
456  void
457  setup_element_info(ElementInfo<Number> & elem_info,
458  const FiniteElement<1> & fe,
459  const ::DoFHandler<dim> &mg_dof)
460  {
461  // currently, we have only FE_Q and FE_DGQ type elements implemented
462  elem_info.n_components = mg_dof.get_fe().element_multiplicity(0);
463  AssertDimension(Utilities::fixed_power<dim>(fe.dofs_per_cell) *
464  elem_info.n_components,
465  mg_dof.get_fe().dofs_per_cell);
466  AssertDimension(fe.degree, mg_dof.get_fe().degree);
467  elem_info.fe_degree = fe.degree;
468  elem_info.element_is_continuous = fe.dofs_per_vertex > 0;
470 
471  // step 1.2: get renumbering of 1D basis functions to lexicographic
472  // numbers. The distinction according to fe.dofs_per_vertex is to support
473  // both continuous and discontinuous bases.
474  std::vector<unsigned int> renumbering(fe.dofs_per_cell);
475  {
477  renumbering[0] = 0;
478  for (unsigned int i = 0; i < fe.dofs_per_line; ++i)
479  renumbering[i + fe.dofs_per_vertex] =
481  if (fe.dofs_per_vertex > 0)
482  renumbering[fe.dofs_per_cell - fe.dofs_per_vertex] =
483  fe.dofs_per_vertex;
484  }
485 
486  // step 1.3: create a dummy 1D quadrature formula to extract the
487  // lexicographic numbering for the elements
488  std::vector<Point<1>> basic_support_points = fe.get_unit_support_points();
489  Assert(fe.dofs_per_vertex == 0 || fe.dofs_per_vertex == 1,
491  const unsigned int shift = fe.dofs_per_cell - fe.dofs_per_vertex;
492  const unsigned int n_child_dofs_1d =
493  (fe.dofs_per_vertex > 0 ? (2 * fe.dofs_per_cell - 1) :
494  (2 * fe.dofs_per_cell));
495 
496  elem_info.n_child_cell_dofs =
497  elem_info.n_components * Utilities::fixed_power<dim>(n_child_dofs_1d);
498  const Quadrature<1> dummy_quadrature(
499  std::vector<Point<1>>(1, Point<1>()));
501  shape_info.reinit(dummy_quadrature, mg_dof.get_fe(), 0);
502  elem_info.lexicographic_numbering = shape_info.lexicographic_numbering;
503 
504  // step 1.4: get the 1d prolongation matrix and combine from both children
505  elem_info.prolongation_matrix_1d.resize(fe.dofs_per_cell *
506  n_child_dofs_1d);
507 
508  for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; ++c)
509  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
510  for (unsigned int j = 0; j < fe.dofs_per_cell; ++j)
511  elem_info
512  .prolongation_matrix_1d[i * n_child_dofs_1d + j + c * shift] =
513  fe.get_prolongation_matrix(c)(renumbering[j], renumbering[i]);
514  }
515 
516  namespace
517  {
523  void
524  replace(const MGConstrainedDoFs * mg_constrained_dofs,
525  const unsigned int level,
526  std::vector<types::global_dof_index> &dof_indices)
527  {
528  if (mg_constrained_dofs != nullptr &&
529  mg_constrained_dofs->get_level_constraints(level).n_constraints() >
530  0)
531  for (auto &ind : dof_indices)
532  if (mg_constrained_dofs->get_level_constraints(level)
534  {
535  Assert(mg_constrained_dofs->get_level_constraints(level)
537  ->size() == 1,
538  ExcInternalError());
539  ind = mg_constrained_dofs->get_level_constraints(level)
541  ->front()
542  .first;
543  }
544  }
545  } // namespace
546 
547  // Sets up most of the internal data structures of the MGTransferMatrixFree
548  // class
549  template <int dim, typename Number>
550  void
551  setup_transfer(
552  const ::DoFHandler<dim> & mg_dof,
553  const MGConstrainedDoFs * mg_constrained_dofs,
554  ElementInfo<Number> & elem_info,
555  std::vector<std::vector<unsigned int>> &level_dof_indices,
556  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
557  & parent_child_connect,
558  std::vector<unsigned int> &n_owned_level_cells,
559  std::vector<std::vector<std::vector<unsigned short>>> &dirichlet_indices,
560  std::vector<std::vector<Number>> & weights_on_refined,
561  std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
562  &copy_indices_global_mine,
564  &ghosted_level_vector)
565  {
566  level_dof_indices.clear();
567  parent_child_connect.clear();
568  n_owned_level_cells.clear();
569  dirichlet_indices.clear();
570  weights_on_refined.clear();
571 
572  // we collect all child DoFs of a mother cell together. For faster
573  // tensorized operations, we align the degrees of freedom
574  // lexicographically. We distinguish FE_Q elements and FE_DGQ elements
575 
576  const ::Triangulation<dim> &tria = mg_dof.get_triangulation();
577 
578  // ---------------------------- 1. Extract 1D info about the finite
579  // element step 1.1: create a 1D copy of the finite element from FETools
580  // where we substitute the template argument
581  AssertDimension(mg_dof.get_fe().n_base_elements(), 1);
582  std::string fe_name = mg_dof.get_fe().base_element(0).get_name();
583  {
584  const std::size_t template_starts = fe_name.find_first_of('<');
585  Assert(fe_name[template_starts + 1] ==
586  (dim == 1 ? '1' : (dim == 2 ? '2' : '3')),
587  ExcInternalError());
588  fe_name[template_starts + 1] = '1';
589  }
590  const std::unique_ptr<FiniteElement<1>> fe(
591  FETools::get_fe_by_name<1, 1>(fe_name));
592 
593  setup_element_info(elem_info, *fe, mg_dof);
594 
595 
596  // -------------- 2. Extract and match dof indices between child and
597  // parent
598  const unsigned int n_levels = tria.n_global_levels();
599  level_dof_indices.resize(n_levels);
600  parent_child_connect.resize(n_levels - 1);
601  n_owned_level_cells.resize(n_levels - 1);
602  std::vector<std::vector<unsigned int>> coarse_level_indices(n_levels - 1);
603  for (unsigned int level = 0;
604  level < std::min(tria.n_levels(), n_levels - 1);
605  ++level)
606  coarse_level_indices[level].resize(tria.n_raw_cells(level),
608  std::vector<types::global_dof_index> local_dof_indices(
609  mg_dof.get_fe().dofs_per_cell);
610  dirichlet_indices.resize(n_levels - 1);
611 
612  // We use the vectors stored ghosted_level_vector in the base class for
613  // keeping ghosted transfer indices. To avoid keeping two very similar
614  // vectors, we merge them here.
615  if (ghosted_level_vector.max_level() != n_levels - 1)
616  ghosted_level_vector.resize(0, n_levels - 1);
617 
618  for (unsigned int level = n_levels - 1; level > 0; --level)
619  {
620  unsigned int counter = 0;
621  std::vector<types::global_dof_index> global_level_dof_indices;
622  std::vector<types::global_dof_index> global_level_dof_indices_remote;
623  std::vector<types::global_dof_index> ghosted_level_dofs;
624  std::vector<types::global_dof_index> global_level_dof_indices_l0;
625  std::vector<types::global_dof_index> ghosted_level_dofs_l0;
626 
627  // step 2.1: loop over the cells on the coarse side
628  typename ::DoFHandler<dim>::cell_iterator cell,
629  endc = mg_dof.end(level - 1);
630  for (cell = mg_dof.begin(level - 1); cell != endc; ++cell)
631  {
632  // need to look into a cell if it has children and it is locally
633  // owned
634  if (!cell->has_children())
635  continue;
636 
637  bool consider_cell = false;
638  if (tria.locally_owned_subdomain() ==
640  cell->level_subdomain_id() == tria.locally_owned_subdomain())
641  consider_cell = true;
642 
643  // due to the particular way we store DoF indices (via children),
644  // we also need to add the DoF indices for coarse cells where we
645  // own at least one child
646  bool cell_is_remote = !consider_cell;
647  for (unsigned int c = 0;
648  c < GeometryInfo<dim>::max_children_per_cell;
649  ++c)
650  if (cell->child(c)->level_subdomain_id() ==
651  tria.locally_owned_subdomain())
652  {
653  consider_cell = true;
654  break;
655  }
656 
657  if (!consider_cell)
658  continue;
659 
660  // step 2.2: loop through children and append the dof indices to
661  // the appropriate list. We need separate lists for the owned
662  // coarse cell case (which will be part of
663  // restriction/prolongation between level-1 and level) and the
664  // remote case (which needs to store DoF indices for the
665  // operations between level and level+1).
666  AssertDimension(cell->n_children(),
668  std::vector<types::global_dof_index> &next_indices =
669  cell_is_remote ? global_level_dof_indices_remote :
670  global_level_dof_indices;
671  const std::size_t start_index = next_indices.size();
672  next_indices.resize(start_index + elem_info.n_child_cell_dofs,
674  for (unsigned int c = 0;
675  c < GeometryInfo<dim>::max_children_per_cell;
676  ++c)
677  {
678  if (cell_is_remote && cell->child(c)->level_subdomain_id() !=
679  tria.locally_owned_subdomain())
680  continue;
681  cell->child(c)->get_mg_dof_indices(local_dof_indices);
682 
683  replace(mg_constrained_dofs, level, local_dof_indices);
684 
685  const IndexSet &owned_level_dofs =
686  mg_dof.locally_owned_mg_dofs(level);
687  for (const auto local_dof_index : local_dof_indices)
688  if (!owned_level_dofs.is_element(local_dof_index))
689  ghosted_level_dofs.push_back(local_dof_index);
690 
691  add_child_indices<dim>(c,
692  fe->dofs_per_cell -
693  fe->dofs_per_vertex,
694  fe->degree,
695  elem_info.lexicographic_numbering,
696  local_dof_indices,
697  &next_indices[start_index]);
698 
699  // step 2.3 store the connectivity to the parent
700  if (cell->child(c)->has_children() &&
701  (tria.locally_owned_subdomain() ==
703  cell->child(c)->level_subdomain_id() ==
704  tria.locally_owned_subdomain()))
705  {
706  const unsigned int child_index =
707  coarse_level_indices[level][cell->child(c)->index()];
708  AssertIndexRange(child_index,
709  parent_child_connect[level].size());
710  unsigned int parent_index = counter;
711  // remote cells, i.e., cells where we work on a further
712  // level but are not treated on the current level, need to
713  // be placed at the end of the list; however, we do not
714  // yet know the exact position in the array, so shift
715  // their parent index by the number of cells so we can set
716  // the correct number after the end of this loop
717  if (cell_is_remote)
718  parent_index =
719  start_index / elem_info.n_child_cell_dofs +
720  tria.n_cells(level);
721  parent_child_connect[level][child_index] =
722  std::make_pair(parent_index, c);
723  AssertIndexRange(mg_dof.get_fe().dofs_per_cell,
724  static_cast<unsigned short>(-1));
725 
726  // set Dirichlet boundary conditions (as a list of
727  // constrained DoFs) for the child
728  if (mg_constrained_dofs != nullptr)
729  for (unsigned int i = 0;
730  i < mg_dof.get_fe().dofs_per_cell;
731  ++i)
732  if (mg_constrained_dofs->is_boundary_index(
733  level,
734  local_dof_indices
735  [elem_info.lexicographic_numbering[i]]))
736  dirichlet_indices[level][child_index].push_back(i);
737  }
738  }
739  if (!cell_is_remote)
740  {
741  AssertIndexRange(static_cast<unsigned int>(cell->index()),
742  coarse_level_indices[level - 1].size());
743  coarse_level_indices[level - 1][cell->index()] = counter++;
744  }
745 
746  // step 2.4: include indices for the coarsest cells. we still
747  // insert the indices as if they were from a child in order to use
748  // the same code (the coarsest level does not matter much in terms
749  // of memory, so we gain in code simplicity)
750  if (level == 1 && !cell_is_remote)
751  {
752  cell->get_mg_dof_indices(local_dof_indices);
753 
754  replace(mg_constrained_dofs, level - 1, local_dof_indices);
755 
756  const IndexSet &owned_level_dofs_l0 =
757  mg_dof.locally_owned_mg_dofs(0);
758  for (const auto local_dof_index : local_dof_indices)
759  if (!owned_level_dofs_l0.is_element(local_dof_index))
760  ghosted_level_dofs_l0.push_back(local_dof_index);
761 
762  const std::size_t start_index =
763  global_level_dof_indices_l0.size();
764  global_level_dof_indices_l0.resize(
765  start_index + elem_info.n_child_cell_dofs,
767  add_child_indices<dim>(
768  0,
769  fe->dofs_per_cell - fe->dofs_per_vertex,
770  fe->degree,
771  elem_info.lexicographic_numbering,
772  local_dof_indices,
773  &global_level_dof_indices_l0[start_index]);
774 
775  dirichlet_indices[0].emplace_back();
776  if (mg_constrained_dofs != nullptr)
777  for (unsigned int i = 0; i < mg_dof.get_fe().dofs_per_cell;
778  ++i)
779  if (mg_constrained_dofs->is_boundary_index(
780  0,
781  local_dof_indices[elem_info
782  .lexicographic_numbering[i]]))
783  dirichlet_indices[0].back().push_back(i);
784  }
785  }
786 
787  // step 2.5: store information about the current level and prepare the
788  // Dirichlet indices and parent-child relationship for the next
789  // coarser level
790  AssertDimension(counter * elem_info.n_child_cell_dofs,
791  global_level_dof_indices.size());
792  n_owned_level_cells[level - 1] = counter;
793  dirichlet_indices[level - 1].resize(counter);
794  parent_child_connect[level - 1].resize(
795  counter,
796  std::make_pair(numbers::invalid_unsigned_int,
798 
799  // step 2.6: put the cells with remotely owned parent to the end of
800  // the list (these are needed for the transfer from level to level+1
801  // but not for the transfer from level-1 to level).
802  if (level < n_levels - 1)
803  for (std::vector<std::pair<unsigned int, unsigned int>>::iterator
804  i = parent_child_connect[level].begin();
805  i != parent_child_connect[level].end();
806  ++i)
807  if (i->first >= tria.n_cells(level))
808  {
809  i->first -= tria.n_cells(level);
810  i->first += counter;
811  }
812 
813  // step 2.7: Initialize the ghosted vector
814  const parallel::Triangulation<dim, dim> *ptria =
815  (dynamic_cast<const parallel::Triangulation<dim, dim> *>(&tria));
816  const MPI_Comm communicator =
817  ptria != nullptr ? ptria->get_communicator() : MPI_COMM_SELF;
818 
819  reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(level),
820  ghosted_level_dofs,
821  communicator,
822  ghosted_level_vector[level],
823  copy_indices_global_mine[level]);
824 
825  copy_indices_to_mpi_local_numbers(
826  *ghosted_level_vector[level].get_partitioner(),
827  global_level_dof_indices,
828  global_level_dof_indices_remote,
829  level_dof_indices[level]);
830  // step 2.8: Initialize the ghosted vector for level 0
831  if (level == 1)
832  {
833  for (unsigned int i = 0; i < parent_child_connect[0].size(); ++i)
834  parent_child_connect[0][i] = std::make_pair(i, 0U);
835 
836  reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(0),
837  ghosted_level_dofs_l0,
838  communicator,
839  ghosted_level_vector[0],
840  copy_indices_global_mine[0]);
841 
842  copy_indices_to_mpi_local_numbers(
843  *ghosted_level_vector[0].get_partitioner(),
844  global_level_dof_indices_l0,
845  std::vector<types::global_dof_index>(),
846  level_dof_indices[0]);
847  }
848  }
849 
850  // ---------------------- 3. compute weights to make restriction additive
851 
852  const unsigned int n_child_dofs_1d =
853  fe->degree + 1 + fe->dofs_per_cell - fe->dofs_per_vertex;
854 
855  // get the valence of the individual components and compute the weights as
856  // the inverse of the valence
857  weights_on_refined.resize(n_levels - 1);
858  for (unsigned int level = 1; level < n_levels; ++level)
859  {
860  ghosted_level_vector[level] = 0;
861  for (unsigned int c = 0; c < n_owned_level_cells[level - 1]; ++c)
862  for (unsigned int j = 0; j < elem_info.n_child_cell_dofs; ++j)
863  ghosted_level_vector[level].local_element(
864  level_dof_indices[level][elem_info.n_child_cell_dofs * c +
865  j]) += Number(1.);
866  ghosted_level_vector[level].compress(VectorOperation::add);
867  ghosted_level_vector[level].update_ghost_values();
868 
869  std::vector<unsigned int> degree_to_3(n_child_dofs_1d);
870  degree_to_3[0] = 0;
871  for (unsigned int i = 1; i < n_child_dofs_1d - 1; ++i)
872  degree_to_3[i] = 1;
873  degree_to_3.back() = 2;
874 
875  // we only store 3^dim weights because all dofs on a line have the
876  // same valence, and all dofs on a quad have the same valence.
877  weights_on_refined[level - 1].resize(n_owned_level_cells[level - 1] *
878  Utilities::fixed_power<dim>(3));
879  for (unsigned int c = 0; c < n_owned_level_cells[level - 1]; ++c)
880  for (unsigned int k = 0, m = 0; k < (dim > 2 ? n_child_dofs_1d : 1);
881  ++k)
882  for (unsigned int j = 0; j < (dim > 1 ? n_child_dofs_1d : 1); ++j)
883  {
884  unsigned int shift = 9 * degree_to_3[k] + 3 * degree_to_3[j];
885  for (unsigned int i = 0; i < n_child_dofs_1d; ++i, ++m)
886  weights_on_refined[level -
887  1][c * Utilities::fixed_power<dim>(3) +
888  shift + degree_to_3[i]] =
889  Number(1.) /
890  ghosted_level_vector[level].local_element(
891  level_dof_indices[level]
892  [elem_info.n_child_cell_dofs * c + m]);
893  }
894  }
895  }
896 
897  } // namespace MGTransfer
898 } // namespace internal
899 
900 // Explicit instantiations
901 
902 #include "mg_transfer_internal.inst"
903 
904 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
Definition: types.h:173
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
const types::subdomain_id invalid_subdomain_id
Definition: types.h:258
void reinit(const Quadrature< 1 > &quad, const FiniteElement< dim > &fe_dim, const unsigned int base_element=0)
bool is_boundary_index(const unsigned int level, const types::global_dof_index index) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1641
const unsigned int degree
Definition: fe_base.h:298
const AffineConstraints< double > & get_level_constraints(const unsigned int level) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const
Definition: fe.cc:335
unsigned int global_to_local(const types::global_dof_index global_index) const
const unsigned int dofs_per_line
Definition: fe_base.h:231
Definition: point.h:110
size_type size() const
Definition: index_set.h:1600
size_type n_constraints() const
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1821
virtual MPI_Comm get_communicator() const
Definition: tria_base.cc:155
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1137
bool is_identity_constrained(const size_type line_n) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
const unsigned int dofs_per_cell
Definition: fe_base.h:282
const std::vector< Point< dim > > & get_unit_support_points() const
Definition: fe.cc:1009
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
unsigned int global_dof_index
Definition: types.h:89
const types::subdomain_id artificial_subdomain_id
Definition: types.h:275
virtual size_type size() const override
bool at_refinement_edge(const unsigned int level, const types::global_dof_index index) const
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1695
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
static ::ExceptionBase & ExcNotImplemented()
const unsigned int dofs_per_vertex
Definition: fe_base.h:225
bool is_element(const size_type index) const
Definition: index_set.h:1665
std::vector< unsigned int > lexicographic_numbering
Definition: shape_info.h:242
const types::global_dof_index invalid_dof_index
Definition: types.h:188
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
size_type n_elements() const
Definition: index_set.h:1732
static ::ExceptionBase & ExcInternalError()