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