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