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