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