Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3104-g8c3bc2e695 2025-04-21 18:40: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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
dof_info.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 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
19
21
24
25#include <deal.II/matrix_free/dof_info.templates.h>
27
28#include <iostream>
29
31
32namespace internal
33{
34 namespace MatrixFreeFunctions
35 {
36 // ensure that the type defined in both dof_info.h and
37 // hanging_nodes_internal.h is consistent
38 static_assert(std::is_same_v<compressed_constraint_kind, std::uint8_t>,
39 "Unexpected type for compressed hanging node indicators!");
40
41
42
44 {
45 clear();
46 }
47
48
49
50 void
52 {
53 row_starts.clear();
54 dof_indices.clear();
56 vector_partitioner.reset();
57 ghost_dofs.clear();
58 dofs_per_cell.clear();
59 dofs_per_face.clear();
63 n_components.clear();
64 start_components.clear();
66 plain_dof_indices.clear();
68 for (unsigned int i = 0; i < 3; ++i)
69 {
70 index_storage_variants[i].clear();
71 dof_indices_contiguous[i].clear();
74 }
75 store_plain_indices = false;
77 max_fe_index = 0;
78 fe_index_conversion.clear();
79 }
80
81
82
83 void
84 DoFInfo::get_dof_indices_on_cell_batch(std::vector<unsigned int> &my_rows,
85 const unsigned int cell,
86 const bool apply_constraints) const
87 {
88 const unsigned int n_fe_components = start_components.back();
89 const unsigned int fe_index =
90 dofs_per_cell.size() == 1 ? 0 : cell_active_fe_index[cell];
91 const unsigned int dofs_this_cell = dofs_per_cell[fe_index];
92
93 const unsigned int n_vectorization = vectorization_length;
94 constexpr auto dof_access_index = dof_access_cell;
96 n_vectorization_lanes_filled[dof_access_index].size());
97 const unsigned int n_vectorization_actual =
98 n_vectorization_lanes_filled[dof_access_index][cell];
99
100 // we might have constraints, so the final number
101 // of indices is not known a priori.
102 // conservatively reserve the maximum without constraints
103 my_rows.reserve(n_vectorization * dofs_this_cell);
104 my_rows.resize(0);
105 unsigned int total_size = 0;
106 for (unsigned int v = 0; v < n_vectorization_actual; ++v)
107 {
108 const unsigned int ib =
109 (cell * n_vectorization + v) * n_fe_components;
110 const unsigned int ie =
111 (cell * n_vectorization + v + 1) * n_fe_components;
112
113 // figure out constraints by comparing constraint_indicator row
114 // shift for this cell within the block as compared to the next
115 // one
116 const bool has_constraints =
117 (hanging_node_constraint_masks.size() != 0 &&
119 hanging_node_constraint_masks[cell * n_vectorization + v] !=
121 hanging_node_constraint_masks_comp[fe_index][0 /*TODO*/]) ||
122 (row_starts[ib].second != row_starts[ib + n_fe_components].second);
123
124 auto do_copy = [&](const unsigned int *begin,
125 const unsigned int *end) {
126 const unsigned int shift = total_size;
127 total_size += (end - begin);
128 my_rows.resize(total_size);
129 std::copy(begin, end, my_rows.begin() + shift);
130 };
131
132 if (!has_constraints || apply_constraints)
133 {
134 const unsigned int *begin =
135 dof_indices.data() + row_starts[ib].first;
136 const unsigned int *end =
137 dof_indices.data() + row_starts[ie].first;
138 do_copy(begin, end);
139 }
140 else
141 {
142 Assert(row_starts_plain_indices[cell * n_vectorization + v] !=
145 const unsigned int *begin =
146 plain_dof_indices.data() +
147 row_starts_plain_indices[cell * n_vectorization + v];
148 const unsigned int *end = begin + dofs_this_cell;
149 do_copy(begin, end);
150 }
151 }
152 }
153
154
155
156 void
157 DoFInfo::assign_ghosts(const std::vector<unsigned int> &boundary_cells,
158 const MPI_Comm communicator_sm,
159 const bool use_vector_data_exchanger_full)
160 {
161 Assert(boundary_cells.size() < row_starts.size(), ExcInternalError());
162
163 // sort ghost dofs and compress out duplicates
164 const unsigned int n_owned = (vector_partitioner->local_range().second -
165 vector_partitioner->local_range().first);
166 const std::size_t n_ghosts = ghost_dofs.size();
167 if constexpr (running_in_debug_mode())
168 {
169 for (const auto dof_index : dof_indices)
170 AssertIndexRange(dof_index, n_owned + n_ghosts);
171 }
172
173 const unsigned int n_components = start_components.back();
174 std::vector<unsigned int> ghost_numbering(n_ghosts);
175 IndexSet ghost_indices(vector_partitioner->size());
176 if (n_ghosts > 0)
177 {
178 unsigned int n_unique_ghosts = 0;
179 // since we need to go back to the local_to_global indices and
180 // replace the temporary numbering of ghosts by the real number in
181 // the index set, we need to store these values
182 std::vector<std::pair<types::global_dof_index, unsigned int>>
183 ghost_origin(n_ghosts);
184 for (std::size_t i = 0; i < n_ghosts; ++i)
185 {
186 ghost_origin[i].first = ghost_dofs[i];
187 ghost_origin[i].second = i;
188 }
189 std::sort(ghost_origin.begin(), ghost_origin.end());
190
191 types::global_dof_index last_contiguous_start = ghost_origin[0].first;
192 ghost_numbering[ghost_origin[0].second] = 0;
193 for (std::size_t i = 1; i < n_ghosts; ++i)
194 {
195 if (ghost_origin[i].first > ghost_origin[i - 1].first + 1)
196 {
197 ghost_indices.add_range(last_contiguous_start,
198 ghost_origin[i - 1].first + 1);
199 last_contiguous_start = ghost_origin[i].first;
200 }
201 if (ghost_origin[i].first > ghost_origin[i - 1].first)
202 ++n_unique_ghosts;
203 ghost_numbering[ghost_origin[i].second] = n_unique_ghosts;
204 }
205 ++n_unique_ghosts;
206 ghost_indices.add_range(last_contiguous_start,
207 ghost_origin.back().first + 1);
208 ghost_indices.compress();
209
210 // make sure that we got the correct local numbering of the ghost
211 // dofs. the ghost index set should store the same number
212 {
213 AssertDimension(n_unique_ghosts, ghost_indices.n_elements());
214 for (std::size_t i = 0; i < n_ghosts; ++i)
215 Assert(ghost_numbering[i] ==
216 ghost_indices.index_within_set(ghost_dofs[i]),
218 }
219
220 // apply correct numbering for ghost indices: We previously just
221 // enumerated them according to their appearance in the
222 // local_to_global structure. Above, we derived a relation between
223 // this enumeration and the actual number
224 const unsigned int n_boundary_cells = boundary_cells.size();
225 for (unsigned int i = 0; i < n_boundary_cells; ++i)
226 {
227 unsigned int *data_ptr =
228 dof_indices.data() +
229 row_starts[boundary_cells[i] * n_components].first;
230 const unsigned int *row_end =
231 dof_indices.data() +
232 row_starts[(boundary_cells[i] + 1) * n_components].first;
233 for (; data_ptr != row_end; ++data_ptr)
234 *data_ptr = ((*data_ptr < n_owned) ?
235 *data_ptr :
236 n_owned + ghost_numbering[*data_ptr - n_owned]);
237
238 // now the same procedure for plain indices
239 if (store_plain_indices == true)
240 {
241 bool has_hanging_nodes = false;
242
243 const unsigned int fe_index =
244 (cell_active_fe_index.empty() ||
245 dofs_per_cell.size() == 1) ?
246 0 :
247 cell_active_fe_index[boundary_cells[i]];
248 AssertIndexRange(fe_index, dofs_per_cell.size());
249
250 if (hanging_node_constraint_masks.size() > 0 &&
252 hanging_node_constraint_masks[boundary_cells[i]] !=
254 for (unsigned int comp = 0; comp < n_components; ++comp)
255 has_hanging_nodes |=
257
258 if (has_hanging_nodes ||
259 row_starts[boundary_cells[i] * n_components].second !=
260 row_starts[(boundary_cells[i] + 1) * n_components]
261 .second)
262 {
263 unsigned int *data_ptr =
264 plain_dof_indices.data() +
265 row_starts_plain_indices[boundary_cells[i]];
266 const unsigned int *row_end =
267 data_ptr + dofs_per_cell[fe_index];
268 for (; data_ptr != row_end; ++data_ptr)
269 *data_ptr =
270 ((*data_ptr < n_owned) ?
271 *data_ptr :
272 n_owned + ghost_numbering[*data_ptr - n_owned]);
273 }
274 }
275 }
276 }
277
278 std::vector<types::global_dof_index> empty;
279 ghost_dofs.swap(empty);
280
281 // set the ghost indices now. need to cast away constness here, but that
282 // is uncritical since we reset the Partitioner in the same initialize
283 // call as this call here.
286 vec_part->set_ghost_indices(ghost_indices);
287
288 if (use_vector_data_exchanger_full == false)
289 vector_exchanger = std::make_shared<
292 else
294 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
295 vector_partitioner, communicator_sm);
296 }
297
298
299
300 void
302 const TaskInfo &task_info,
303 const std::vector<unsigned int> &renumbering,
304 const std::vector<unsigned int> &constraint_pool_row_index,
305 const std::vector<unsigned char> &irregular_cells)
306 {
307 (void)constraint_pool_row_index;
308
309 // first reorder the active FE index.
310 const bool have_hp = dofs_per_cell.size() > 1;
311 if (cell_active_fe_index.size() > 0)
312 {
313 std::vector<unsigned int> new_active_fe_index;
314 new_active_fe_index.reserve(task_info.cell_partition_data.back());
315 unsigned int position_cell = 0;
316 for (unsigned int cell = 0;
317 cell < task_info.cell_partition_data.back();
318 ++cell)
319 {
320 const unsigned int n_comp =
321 (irregular_cells[cell] > 0 ? irregular_cells[cell] :
323
324 // take maximum FE index among the ones present (we might have
325 // lumped some lower indices into higher ones)
326 unsigned int fe_index =
327 cell_active_fe_index[renumbering[position_cell]];
328 for (unsigned int j = 1; j < n_comp; ++j)
329 fe_index = std::max(
330 fe_index,
331 cell_active_fe_index[renumbering[position_cell + j]]);
332
333 new_active_fe_index.push_back(fe_index);
334 position_cell += n_comp;
335 }
336 std::swap(new_active_fe_index, cell_active_fe_index);
337 }
338 if (have_hp)
340 task_info.cell_partition_data.back());
341
342 const unsigned int n_components = start_components.back();
343
344 std::vector<std::pair<unsigned int, unsigned int>> new_row_starts(
346 task_info.cell_partition_data.back() +
347 1);
348 std::vector<unsigned int> new_dof_indices;
349 std::vector<std::pair<unsigned short, unsigned short>>
350 new_constraint_indicator;
351 std::vector<unsigned int> new_plain_indices, new_rowstart_plain;
352 unsigned int position_cell = 0;
353 new_dof_indices.reserve(dof_indices.size());
354 new_constraint_indicator.reserve(constraint_indicator.size());
355
356 std::vector<compressed_constraint_kind> new_hanging_node_constraint_masks;
357 new_hanging_node_constraint_masks.reserve(
359
360 if (store_plain_indices == true)
361 {
362 new_rowstart_plain.resize(vectorization_length *
363 task_info.cell_partition_data.back() +
364 1,
366 new_plain_indices.reserve(plain_dof_indices.size());
367 }
368
369 // copy the indices and the constraint indicators to the new data field,
370 // where we will go through the cells in the renumbered way. in case the
371 // vectorization length does not exactly match up, we fill invalid
372 // numbers to the rowstart data. for contiguous cell indices, we skip
373 // the rowstarts field completely and directly go into the
374 // new_dof_indices field (this layout is used in FEEvaluation).
375 for (unsigned int i = 0; i < task_info.cell_partition_data.back(); ++i)
376 {
377 const unsigned int n_vect =
378 (irregular_cells[i] > 0 ? irregular_cells[i] :
380 const unsigned int dofs_per_cell =
381 have_hp ? this->dofs_per_cell[cell_active_fe_index[i]] :
382 this->dofs_per_cell[0];
383
384 for (unsigned int j = 0; j < n_vect; ++j)
385 {
386 const unsigned int cell_no =
387 renumbering[position_cell + j] * n_components;
388
389 bool has_hanging_nodes = false;
390
391 if (hanging_node_constraint_masks.size() > 0 &&
393 {
394 const auto mask =
395 hanging_node_constraint_masks[renumbering[position_cell +
396 j]];
397 new_hanging_node_constraint_masks.push_back(mask);
398
400 for (unsigned int comp = 0; comp < n_components; ++comp)
401 has_hanging_nodes |= hanging_node_constraint_masks_comp
402 [have_hp ? cell_active_fe_index[i] : 0][comp];
403 }
404
405 for (unsigned int comp = 0; comp < n_components; ++comp)
406 {
407 new_row_starts[(i * vectorization_length + j) * n_components +
408 comp]
409 .first = new_dof_indices.size();
410 new_row_starts[(i * vectorization_length + j) * n_components +
411 comp]
412 .second = new_constraint_indicator.size();
413
414 new_dof_indices.insert(
415 new_dof_indices.end(),
416 dof_indices.data() + row_starts[cell_no + comp].first,
417 dof_indices.data() + row_starts[cell_no + comp + 1].first);
418 for (unsigned int index = row_starts[cell_no + comp].second;
419 index != row_starts[cell_no + comp + 1].second;
420 ++index)
421 new_constraint_indicator.push_back(
423 }
425 ((row_starts[cell_no].second !=
426 row_starts[cell_no + n_components].second) ||
427 has_hanging_nodes))
428 {
429 new_rowstart_plain[i * vectorization_length + j] =
430 new_plain_indices.size();
431 new_plain_indices.insert(
432 new_plain_indices.end(),
433 plain_dof_indices.data() +
435 plain_dof_indices.data() +
438 }
439 }
440 for (unsigned int j = n_vect; j < vectorization_length; ++j)
441 for (unsigned int comp = 0; comp < n_components; ++comp)
442 {
443 new_row_starts[(i * vectorization_length + j) * n_components +
444 comp]
445 .first = new_dof_indices.size();
446 new_row_starts[(i * vectorization_length + j) * n_components +
447 comp]
448 .second = new_constraint_indicator.size();
449 }
450
451 for (unsigned int j = n_vect; j < vectorization_length; ++j)
452 if (hanging_node_constraint_masks.size() > 0)
453 new_hanging_node_constraint_masks.push_back(
455
456 position_cell += n_vect;
457 }
458 AssertDimension(position_cell * n_components + 1, row_starts.size());
459
460 AssertDimension(dof_indices.size(), new_dof_indices.size());
461 new_row_starts[task_info.cell_partition_data.back() *
463 .first = new_dof_indices.size();
464 new_row_starts[task_info.cell_partition_data.back() *
466 .second = new_constraint_indicator.size();
467
469 new_constraint_indicator.size());
470
471 new_row_starts.swap(row_starts);
472 new_dof_indices.swap(dof_indices);
473 new_constraint_indicator.swap(constraint_indicator);
474 new_plain_indices.swap(plain_dof_indices);
475 new_rowstart_plain.swap(row_starts_plain_indices);
476 new_hanging_node_constraint_masks.swap(hanging_node_constraint_masks);
477
478 if constexpr (running_in_debug_mode())
479 {
480 // sanity check 1: all indices should be smaller than the number of
481 // dofs locally owned plus the number of ghosts
482 const unsigned int index_range =
483 (vector_partitioner->local_range().second -
484 vector_partitioner->local_range().first) +
485 vector_partitioner->ghost_indices().n_elements();
486 for (const auto dof_index : dof_indices)
487 AssertIndexRange(dof_index, index_range);
488
489 // sanity check 2: for the constraint indicators, the first index
490 // should be smaller than the number of indices in the row, and the
491 // second index should be smaller than the number of constraints in
492 // the constraint pool.
493 for (unsigned int row = 0; row < task_info.cell_partition_data.back();
494 ++row)
495 {
496 const unsigned int row_length_ind =
498 .first -
502 .second,
503 constraint_indicator.size() + 1);
504 const std::pair<unsigned short, unsigned short>
505 *con_it =
506 constraint_indicator.data() +
508 *end_con =
509 constraint_indicator.data() +
511 .second;
512 for (; con_it != end_con; ++con_it)
513 {
514 AssertIndexRange(con_it->first, row_length_ind + 1);
515 AssertIndexRange(con_it->second,
516 constraint_pool_row_index.size() - 1);
517 }
518 }
519
520 // sanity check 3: check the number of cells once again
521 unsigned int n_active_cells = 0;
522 for (unsigned int c = 0;
523 c < *(task_info.cell_partition_data.end() - 2);
524 ++c)
525 if (irregular_cells[c] > 0)
526 n_active_cells += irregular_cells[c];
527 else
528 n_active_cells += vectorization_length;
529 AssertDimension(n_active_cells, task_info.n_active_cells);
530 }
531
532 compute_cell_index_compression(irregular_cells);
533 }
534
535
536
537 void
539 const std::vector<unsigned char> &irregular_cells)
540 {
541 const bool have_hp = dofs_per_cell.size() > 1;
542 const unsigned int n_components = start_components.back();
543
545 row_starts.size() % vectorization_length == 1,
547 if (vectorization_length > 1)
549 irregular_cells.size());
551 irregular_cells.size(), IndexStorageVariants::full);
553 irregular_cells.size());
554 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
555 if (irregular_cells[i] > 0)
556 n_vectorization_lanes_filled[dof_access_cell][i] = irregular_cells[i];
557 else
560
562 irregular_cells.size() * vectorization_length,
567 irregular_cells.size() * vectorization_length,
569
570 std::vector<unsigned int> index_kinds(
571 static_cast<unsigned int>(
573 1);
574 std::vector<unsigned int> offsets(vectorization_length);
575 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
576 {
577 const unsigned int ndofs =
578 dofs_per_cell[have_hp ? cell_active_fe_index[i] : 0];
579 const unsigned int n_comp =
581
582 // check 1: Check if there are constraints -> no compression possible
583 bool has_constraints = false;
584 for (unsigned int j = 0; j < n_comp; ++j)
585 {
586 const unsigned int cell_no = i * vectorization_length + j;
587 if (row_starts[cell_no * n_components].second !=
588 row_starts[(cell_no + 1) * n_components].second)
589 {
590 has_constraints = true;
591 break;
592 }
593 }
594 if (has_constraints)
597 else
598 {
599 bool indices_are_contiguous = (ndofs > 0);
600 for (unsigned int j = 0; j < n_comp; ++j)
601 {
602 const unsigned int cell_no = i * vectorization_length + j;
603 const unsigned int *dof_indices =
604 this->dof_indices.data() +
605 row_starts[cell_no * n_components].first;
607 ndofs,
608 row_starts[(cell_no + 1) * n_components].first -
609 row_starts[cell_no * n_components].first);
610 for (unsigned int i = 1; i < ndofs; ++i)
611 if (dof_indices[i] != dof_indices[0] + i)
612 {
613 indices_are_contiguous = false;
614 break;
615 }
616 }
617
618 bool indices_are_interleaved_and_contiguous =
619 (ndofs > 1 && n_comp == vectorization_length);
620
621 {
622 const unsigned int *dof_indices =
623 this->dof_indices.data() +
625 for (unsigned int k = 0;
626 k < ndofs && indices_are_interleaved_and_contiguous;
627 ++k)
628 for (unsigned int j = 0; j < n_comp; ++j)
629 if (dof_indices[j * ndofs + k] !=
630 dof_indices[0] + k * n_comp + j)
631 {
632 indices_are_interleaved_and_contiguous = false;
633 break;
634 }
635 }
636
637 if (indices_are_contiguous ||
638 indices_are_interleaved_and_contiguous)
639 {
640 for (unsigned int j = 0; j < n_comp; ++j)
641 {
642 const unsigned int start_index =
645 .first;
646 AssertIndexRange(start_index, dof_indices.size());
648 [i * vectorization_length + j] =
649 this->dof_indices.empty() ?
650 0 :
651 this->dof_indices[start_index];
652 }
653 }
654
655 if (indices_are_interleaved_and_contiguous)
656 {
660 for (unsigned int j = 0; j < n_comp; ++j)
662 j] = n_comp;
663 }
664 else if (indices_are_contiguous)
665 {
668 for (unsigned int j = 0; j < n_comp; ++j)
670 j] = 1;
671 }
672 else if (ndofs > 0)
673 {
674 int indices_are_interleaved_and_mixed = 2;
675 const unsigned int *dof_indices =
676 &this->dof_indices[row_starts[i * vectorization_length *
678 .first];
679 for (unsigned int j = 0; j < n_comp; ++j)
680 offsets[j] =
681 dof_indices[j * ndofs + 1] - dof_indices[j * ndofs];
682 for (unsigned int k = 0;
683 k < ndofs && indices_are_interleaved_and_mixed != 0;
684 ++k)
685 for (unsigned int j = 0; j < n_comp; ++j)
686 // the first if case is to avoid negative offsets
687 // (invalid)
688 if (dof_indices[j * ndofs + 1] < dof_indices[j * ndofs] ||
689 dof_indices[j * ndofs + k] !=
690 dof_indices[j * ndofs] + k * offsets[j])
691 {
692 indices_are_interleaved_and_mixed = 0;
693 break;
694 }
695 if (indices_are_interleaved_and_mixed == 2)
696 {
697 for (unsigned int j = 0; j < n_comp; ++j)
700 offsets[j];
701 for (unsigned int j = 0; j < n_comp; ++j)
703 [i * vectorization_length + j] =
704 dof_indices[j * ndofs];
705 for (unsigned int j = 0; j < n_comp; ++j)
706 if (offsets[j] != vectorization_length)
707 {
708 indices_are_interleaved_and_mixed = 1;
709 break;
710 }
711 if (indices_are_interleaved_and_mixed == 1 ||
712 n_comp != vectorization_length)
716 else
719 }
720 else
721 {
722 const unsigned int *dof_indices =
723 this->dof_indices.data() +
725 .first;
726 if (n_comp == vectorization_length)
729 else
732
733 // do not use interleaved storage if two entries within
734 // vectorized array point to the same index (scatter not
735 // possible due to race condition)
736 for (const unsigned int *indices = dof_indices;
737 indices != dof_indices + ndofs;
738 ++indices)
739 {
740 bool is_sorted = true;
741 unsigned int previous = indices[0];
742 for (unsigned int l = 1; l < n_comp; ++l)
743 {
744 const unsigned int current = indices[l * ndofs];
745 if (current <= previous)
746 is_sorted = false;
747
748 // the simple check failed, must compare all
749 // indices manually - due to short sizes this
750 // O(n^2) algorithm is better than sorting
751 if (!is_sorted)
752 for (unsigned int j = 0; j < l; ++j)
753 if (indices[j * ndofs] == current)
754 {
756 [dof_access_cell][i] =
758 break;
759 }
760 previous = current;
761 }
762 }
763 }
764 }
765 else // ndofs == 0
768 }
769 index_kinds[static_cast<unsigned int>(
771 }
772
773 // Cleanup phase: we want to avoid single cells with different properties
774 // than the bulk of the domain in order to avoid extra checks in the face
775 // identification.
776
777 // Step 1: check whether the interleaved indices were only assigned to
778 // the single cell within a vectorized array.
779 auto fix_single_interleaved_indices =
780 [&](const IndexStorageVariants variant) {
781 if (index_kinds[static_cast<unsigned int>(
783 0 &&
784 index_kinds[static_cast<unsigned int>(variant)] > 0)
785 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
786 {
793 [i * vectorization_length] ==
794 1))
795 {
797 index_kinds[static_cast<unsigned int>(
800 index_kinds[static_cast<unsigned int>(variant)]++;
801 }
802 }
803 };
804
805 fix_single_interleaved_indices(IndexStorageVariants::full);
806 fix_single_interleaved_indices(IndexStorageVariants::contiguous);
807 fix_single_interleaved_indices(IndexStorageVariants::interleaved);
808
809 unsigned int n_interleaved =
810 index_kinds[static_cast<unsigned int>(
812 index_kinds[static_cast<unsigned int>(
814 index_kinds[static_cast<unsigned int>(
816
817 // Step 2: fix single contiguous cell among others with interleaved
818 // storage
819 if (n_interleaved > 0 && index_kinds[static_cast<unsigned int>(
821 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
824 {
827 index_kinds[static_cast<unsigned int>(
829 index_kinds[static_cast<unsigned int>(
831 }
832
833 // Step 3: Interleaved cells are left but also some non-contiguous ones
834 // -> revert all to full storage
835 if (n_interleaved > 0 &&
836 index_kinds[static_cast<unsigned int>(IndexStorageVariants::full)] +
837 index_kinds[static_cast<unsigned int>(
839 0)
840 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
843 {
844 index_kinds[static_cast<unsigned int>(
845 index_storage_variants[2][i])]--;
850 else
853 index_kinds[static_cast<unsigned int>(
855 }
856
857 // Step 4: Copy the interleaved indices into their own data structure
858 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
861 {
864 {
867 continue;
868 }
869 const unsigned int ndofs =
870 dofs_per_cell[have_hp ? cell_active_fe_index[i] : 0];
871 const unsigned int *dof_indices =
872 &this->dof_indices
874 unsigned int *interleaved_dof_indices =
876 [row_starts[i * vectorization_length * n_components].first];
877 AssertDimension(this->dof_indices.size(),
878 this->dof_indices_interleaved.size());
883 this->dof_indices_interleaved.size() + 1);
885 row_starts[i * vectorization_length * n_components].first +
886 ndofs * vectorization_length,
887 this->dof_indices_interleaved.size() + 1);
888 for (unsigned int k = 0; k < ndofs; ++k)
889 {
890 const unsigned int *my_dof_indices = dof_indices + k;
891 const unsigned int *end =
892 interleaved_dof_indices + vectorization_length;
893 for (; interleaved_dof_indices != end;
894 ++interleaved_dof_indices, my_dof_indices += ndofs)
895 *interleaved_dof_indices = *my_dof_indices;
896 }
897 }
898 }
899
900
901
902 void
904 const Table<2, ShapeInfo<double>> &shape_info,
905 const unsigned int n_owned_cells,
906 const unsigned int n_lanes,
907 const std::vector<FaceToCellTopology<1>> &inner_faces,
908 const std::vector<FaceToCellTopology<1>> &ghosted_faces,
909 const bool fill_cell_centric,
910 const MPI_Comm communicator_sm,
911 const bool use_vector_data_exchanger_full)
912 {
914
915 // partitioner 0: no face integrals, simply use the indices present
916 // on the cells
917 std::vector<types::global_dof_index> ghost_indices;
918 {
919 const unsigned int n_components = start_components.back();
920 for (unsigned int cell = 0; cell < n_owned_cells; ++cell)
921 {
922 for (unsigned int i = row_starts[cell * n_components].first;
923 i < row_starts[(cell + 1) * n_components].first;
924 ++i)
925 if (dof_indices[i] >= part.locally_owned_size())
926 ghost_indices.push_back(part.local_to_global(dof_indices[i]));
927
928 const unsigned int fe_index =
929 dofs_per_cell.size() == 1 ? 0 :
930 cell_active_fe_index[cell / n_lanes];
931 const unsigned int dofs_this_cell = dofs_per_cell[fe_index];
932
933 for (unsigned int i = row_starts_plain_indices[cell];
934 i < row_starts_plain_indices[cell] + dofs_this_cell;
935 ++i)
936 if (plain_dof_indices[i] >= part.locally_owned_size())
937 ghost_indices.push_back(
939 }
940 std::sort(ghost_indices.begin(), ghost_indices.end());
941 IndexSet compressed_set(part.size());
942 compressed_set.add_indices(ghost_indices.begin(), ghost_indices.end());
943 compressed_set.subtract_set(part.locally_owned_range());
944 const bool all_ghosts_equal =
945 Utilities::MPI::min<int>(static_cast<int>(
946 compressed_set.n_elements() ==
947 part.ghost_indices().n_elements()),
948 part.get_mpi_communicator()) != 0;
949
950 std::shared_ptr<const Utilities::MPI::Partitioner> temp_0;
951
952 if (all_ghosts_equal)
953 temp_0 = vector_partitioner;
954 else
955 {
956 temp_0 = std::make_shared<Utilities::MPI::Partitioner>(
958 const_cast<Utilities::MPI::Partitioner *>(temp_0.get())
959 ->set_ghost_indices(compressed_set, part.ghost_indices());
960 }
961
962 if (use_vector_data_exchanger_full == false)
963 vector_exchanger_face_variants[0] = std::make_shared<
965 temp_0);
966 else
968 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
969 temp_0, communicator_sm);
970 }
971
972 // construct a numbering of faces
973 std::vector<FaceToCellTopology<1>> all_faces(inner_faces);
974 all_faces.insert(all_faces.end(),
975 ghosted_faces.begin(),
976 ghosted_faces.end());
977 Table<2, unsigned int> cell_and_face_to_faces(
978 (row_starts.size() - 1) / start_components.back(),
979 2 * shape_info(0, 0).n_dimensions);
980 cell_and_face_to_faces.fill(numbers::invalid_unsigned_int);
981 for (unsigned int f = 0; f < all_faces.size(); ++f)
982 {
983 cell_and_face_to_faces(all_faces[f].cells_interior[0],
984 all_faces[f].interior_face_no) = f;
985 Assert(all_faces[f].cells_exterior[0] !=
988 cell_and_face_to_faces(all_faces[f].cells_exterior[0],
989 all_faces[f].exterior_face_no) = f;
990 }
991
992 // lambda function to detect objects on face pairs
993 const auto loop_over_faces =
994 [&](const std::function<
995 void(const unsigned int, const unsigned int, const bool)> &fu) {
996 for (const auto &face : inner_faces)
997 {
998 AssertIndexRange(face.cells_interior[0], n_owned_cells);
999 fu(face.cells_exterior[0], face.exterior_face_no, false /*flag*/);
1000 }
1001 };
1002
1003 const auto loop_over_all_faces =
1004 [&](const std::function<
1005 void(const unsigned int, const unsigned int, const bool)> &fu) {
1006 for (unsigned int c = 0; c < cell_and_face_to_faces.size(0); ++c)
1007 for (unsigned int d = 0; d < cell_and_face_to_faces.size(1); ++d)
1008 {
1009 const unsigned int f = cell_and_face_to_faces(c, d);
1011 continue;
1012
1013 const unsigned int cell_m = all_faces[f].cells_interior[0];
1014 const unsigned int cell_p = all_faces[f].cells_exterior[0];
1015
1016 const bool ext = c == cell_m;
1017
1018 if (ext && cell_p == numbers::invalid_unsigned_int)
1019 continue;
1020
1021 const unsigned int p = ext ? cell_p : cell_m;
1022 const unsigned int face_no = ext ?
1023 all_faces[f].exterior_face_no :
1024 all_faces[f].interior_face_no;
1025
1026 fu(p, face_no, true);
1027 }
1028 };
1029
1030 const auto process_values =
1031 [&](
1032 std::shared_ptr<const Utilities::MPI::Partitioner>
1033 &vector_partitioner_values,
1034 const std::function<void(
1035 const std::function<void(
1036 const unsigned int, const unsigned int, const bool)> &)> &loop) {
1037 bool all_nodal_and_tensorial = shape_info.size(1) == 1;
1038
1039 if (all_nodal_and_tensorial)
1040 for (unsigned int c = 0; c < n_base_elements; ++c)
1041 {
1042 const auto &si =
1043 shape_info(global_base_element_offset + c, 0).data.front();
1044 if (!si.nodal_at_cell_boundaries ||
1045 (si.element_type ==
1047 all_nodal_and_tensorial = false;
1048 }
1049
1050 if (all_nodal_and_tensorial == false)
1051 vector_partitioner_values = vector_partitioner;
1052 else
1053 {
1054 bool has_noncontiguous_cell = false;
1055
1056 loop([&](const unsigned int cell_no,
1057 const unsigned int face_no,
1058 const bool flag) {
1059 const unsigned int index =
1061 if (flag || (index != numbers::invalid_unsigned_int &&
1062 index >= part.locally_owned_size()))
1063 {
1064 const unsigned int stride =
1065 dof_indices_interleave_strides[dof_access_cell][cell_no];
1066 unsigned int i = 0;
1067 for (unsigned int e = 0; e < n_base_elements; ++e)
1068 for (unsigned int c = 0; c < n_components[e]; ++c)
1069 {
1070 const ShapeInfo<double> &shape =
1071 shape_info(global_base_element_offset + e, 0);
1072 for (unsigned int j = 0;
1073 j < shape.dofs_per_component_on_face;
1074 ++j)
1075 ghost_indices.push_back(part.local_to_global(
1076 index + i +
1077 shape.face_to_cell_index_nodal(face_no, j) *
1078 stride));
1079 i += shape.dofs_per_component_on_cell * stride;
1080 }
1081 AssertDimension(i, dofs_per_cell[0] * stride);
1082 }
1084 has_noncontiguous_cell = true;
1085 });
1086 has_noncontiguous_cell =
1087 Utilities::MPI::min<int>(static_cast<int>(
1088 has_noncontiguous_cell),
1089 part.get_mpi_communicator()) != 0;
1090
1091 std::sort(ghost_indices.begin(), ghost_indices.end());
1092 IndexSet compressed_set(part.size());
1093 compressed_set.add_indices(ghost_indices.begin(),
1094 ghost_indices.end());
1095 compressed_set.subtract_set(part.locally_owned_range());
1096 const bool all_ghosts_equal =
1097 Utilities::MPI::min<int>(static_cast<int>(
1098 compressed_set.n_elements() ==
1099 part.ghost_indices().n_elements()),
1100 part.get_mpi_communicator()) != 0;
1101 if (all_ghosts_equal || has_noncontiguous_cell)
1102 vector_partitioner_values = vector_partitioner;
1103 else
1104 {
1105 vector_partitioner_values =
1106 std::make_shared<Utilities::MPI::Partitioner>(
1108 const_cast<Utilities::MPI::Partitioner *>(
1109 vector_partitioner_values.get())
1110 ->set_ghost_indices(compressed_set, part.ghost_indices());
1111 }
1112 }
1113 };
1114
1115
1116 const auto process_gradients =
1117 [&](
1118 const std::shared_ptr<const Utilities::MPI::Partitioner>
1119 &vector_partitoner_values,
1120 std::shared_ptr<const Utilities::MPI::Partitioner>
1121 &vector_partitioner_gradients,
1122 const std::function<void(
1123 const std::function<void(
1124 const unsigned int, const unsigned int, const bool)> &)> &loop) {
1125 bool all_hermite = shape_info.size(1) == 1;
1126
1127 if (all_hermite)
1128 for (unsigned int c = 0; c < n_base_elements; ++c)
1129 if (shape_info(global_base_element_offset + c, 0).element_type !=
1131 all_hermite = false;
1132 if (all_hermite == false ||
1133 vector_partitoner_values.get() == vector_partitioner.get())
1134 vector_partitioner_gradients = vector_partitioner;
1135 else
1136 {
1137 loop([&](const unsigned int cell_no,
1138 const unsigned int face_no,
1139 const bool flag) {
1140 const unsigned int index =
1141 dof_indices_contiguous[dof_access_cell][cell_no];
1142 if (flag || (index != numbers::invalid_unsigned_int &&
1143 index >= part.locally_owned_size()))
1144 {
1145 const unsigned int stride =
1146 dof_indices_interleave_strides[dof_access_cell][cell_no];
1147 unsigned int i = 0;
1148 for (unsigned int e = 0; e < n_base_elements; ++e)
1149 for (unsigned int c = 0; c < n_components[e]; ++c)
1150 {
1151 const ShapeInfo<double> &shape =
1152 shape_info(global_base_element_offset + e, 0);
1153 for (unsigned int j = 0;
1154 j < 2 * shape.dofs_per_component_on_face;
1155 ++j)
1156 ghost_indices.push_back(part.local_to_global(
1157 index + i +
1158 shape.face_to_cell_index_hermite(face_no, j) *
1159 stride));
1160 i += shape.dofs_per_component_on_cell * stride;
1161 }
1162 AssertDimension(i, dofs_per_cell[0] * stride);
1163 }
1164 });
1165 std::sort(ghost_indices.begin(), ghost_indices.end());
1166 IndexSet compressed_set(part.size());
1167 compressed_set.add_indices(ghost_indices.begin(),
1168 ghost_indices.end());
1169 compressed_set.subtract_set(part.locally_owned_range());
1170 const bool all_ghosts_equal =
1171 Utilities::MPI::min<int>(static_cast<int>(
1172 compressed_set.n_elements() ==
1173 part.ghost_indices().n_elements()),
1174 part.get_mpi_communicator()) != 0;
1175 if (all_ghosts_equal)
1176 vector_partitioner_gradients = vector_partitioner;
1177 else
1178 {
1179 vector_partitioner_gradients =
1180 std::make_shared<Utilities::MPI::Partitioner>(
1181 part.locally_owned_range(), part.get_mpi_communicator());
1182 const_cast<Utilities::MPI::Partitioner *>(
1183 vector_partitioner_gradients.get())
1184 ->set_ghost_indices(compressed_set, part.ghost_indices());
1185 }
1186 }
1187 };
1188
1189 std::shared_ptr<const Utilities::MPI::Partitioner> temp_1, temp_2, temp_3,
1190 temp_4;
1191
1192 // partitioner 1: values on faces
1193 process_values(temp_1, loop_over_faces);
1194
1195 // partitioner 2: values and gradients on faces
1196 process_gradients(temp_1, temp_2, loop_over_faces);
1197
1198 if (fill_cell_centric)
1199 {
1200 ghost_indices.clear();
1201 // partitioner 3: values on all faces
1202 process_values(temp_3, loop_over_all_faces);
1203 // partitioner 4: values and gradients on faces
1204 process_gradients(temp_3, temp_4, loop_over_all_faces);
1205 }
1206 else
1207 {
1208 temp_3 = std::make_shared<Utilities::MPI::Partitioner>(
1209 part.locally_owned_range(), part.get_mpi_communicator());
1210 temp_4 = std::make_shared<Utilities::MPI::Partitioner>(
1211 part.locally_owned_range(), part.get_mpi_communicator());
1212 }
1213
1214 if (use_vector_data_exchanger_full == false)
1215 {
1216 vector_exchanger_face_variants[1] = std::make_shared<
1217 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1218 temp_1);
1219 vector_exchanger_face_variants[2] = std::make_shared<
1220 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1221 temp_2);
1222 vector_exchanger_face_variants[3] = std::make_shared<
1223 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1224 temp_3);
1225 vector_exchanger_face_variants[4] = std::make_shared<
1226 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1227 temp_4);
1228 }
1229 else
1230 {
1231 vector_exchanger_face_variants[1] =
1232 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1233 temp_1, communicator_sm);
1234 vector_exchanger_face_variants[2] =
1235 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1236 temp_2, communicator_sm);
1237 vector_exchanger_face_variants[3] =
1238 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1239 temp_3, communicator_sm);
1240 vector_exchanger_face_variants[4] =
1241 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1242 temp_4, communicator_sm);
1243 }
1244 }
1245
1246
1247
1248 void
1249 DoFInfo::compute_shared_memory_contiguous_indices(
1250 std::array<std::vector<std::pair<unsigned int, unsigned int>>, 3>
1251 &cell_indices_contiguous_sm)
1252 {
1253 AssertDimension(dofs_per_cell.size(), 1);
1254
1255 for (unsigned int i = 0; i < 3; ++i)
1256 {
1257 dof_indices_contiguous_sm[i].resize(
1258 cell_indices_contiguous_sm[i].size());
1259
1260 for (unsigned int j = 0; j < cell_indices_contiguous_sm[i].size();
1261 ++j)
1262 if (cell_indices_contiguous_sm[i][j].first !=
1264 dof_indices_contiguous_sm[i][j] = {
1265 cell_indices_contiguous_sm[i][j].first,
1266 cell_indices_contiguous_sm[i][j].second * dofs_per_cell[0]};
1267 else
1268 dof_indices_contiguous_sm[i][j] = {numbers::invalid_unsigned_int,
1270 }
1271 }
1272
1273
1274
1275 namespace internal
1276 {
1277 // We construct the connectivity graph in parallel. we use one lock for
1278 // 256 degrees of freedom to keep the number of locks down to a
1279 // reasonable level and reduce the cost of locking to some extent.
1280 static constexpr unsigned int bucket_size_threading = 256;
1281
1282
1283
1284 void
1285 compute_row_lengths(const unsigned int begin,
1286 const unsigned int end,
1287 const DoFInfo &dof_info,
1288 std::vector<std::mutex> &mutexes,
1289 std::vector<unsigned int> &row_lengths)
1290 {
1291 std::vector<unsigned int> scratch;
1292 const unsigned int n_components = dof_info.start_components.back();
1293 for (unsigned int block = begin; block < end; ++block)
1294 {
1295 scratch.assign(
1296 dof_info.dof_indices.data() +
1297 dof_info.row_starts[block * n_components].first,
1298 dof_info.dof_indices.data() +
1299 dof_info.row_starts[(block + 1) * n_components].first);
1300 std::sort(scratch.begin(), scratch.end());
1301
1302 const std::vector<unsigned int>::const_iterator end_unique =
1303 std::unique(scratch.begin(), scratch.end());
1304 for (std::vector<unsigned int>::const_iterator it = scratch.begin();
1305 it != end_unique;
1306 /* update in loop body */)
1307 {
1308 // In this code, the procedure is that we insert all elements
1309 // that are within the range of one lock at once
1310 const unsigned int next_bucket =
1312
1313 std::lock_guard<std::mutex> lock(
1314 mutexes[*it / bucket_size_threading]);
1315 for (; it != end_unique && *it < next_bucket; ++it)
1316 {
1317 AssertIndexRange(*it, row_lengths.size());
1318 ++row_lengths[*it];
1319 }
1320 }
1321 }
1322 }
1323
1324 void
1325 fill_connectivity_dofs(const unsigned int begin,
1326 const unsigned int end,
1327 const DoFInfo &dof_info,
1328 const std::vector<unsigned int> &row_lengths,
1329 std::vector<std::mutex> &mutexes,
1330 ::SparsityPattern &connectivity_dof)
1331 {
1332 std::vector<unsigned int> scratch;
1333 const unsigned int n_components = dof_info.start_components.back();
1334 for (unsigned int block = begin; block < end; ++block)
1335 {
1336 scratch.assign(
1337 dof_info.dof_indices.data() +
1338 dof_info.row_starts[block * n_components].first,
1339 dof_info.dof_indices.data() +
1340 dof_info.row_starts[(block + 1) * n_components].first);
1341 std::sort(scratch.begin(), scratch.end());
1342
1343 const std::vector<unsigned int>::const_iterator end_unique =
1344 std::unique(scratch.begin(), scratch.end());
1345 for (std::vector<unsigned int>::const_iterator it = scratch.begin();
1346 it != end_unique;
1347 /* update in loop body */)
1348 {
1349 const unsigned int next_bucket =
1351
1352 std::lock_guard<std::mutex> lock(
1353 mutexes[*it / bucket_size_threading]);
1354 for (; it != end_unique && *it < next_bucket; ++it)
1355 if (row_lengths[*it] > 0)
1356 connectivity_dof.add(*it, block);
1357 }
1358 }
1359 }
1360
1361
1362
1363 void
1364 fill_connectivity(const unsigned int begin,
1365 const unsigned int end,
1366 const DoFInfo &dof_info,
1367 const std::vector<unsigned int> &renumbering,
1368 const ::SparsityPattern &connectivity_dof,
1369 DynamicSparsityPattern &connectivity)
1370 {
1371 ordered_vector row_entries;
1372 const unsigned int n_components = dof_info.start_components.back();
1373 for (unsigned int block = begin; block < end; ++block)
1374 {
1375 row_entries.clear();
1376
1377 const unsigned int
1378 *it = dof_info.dof_indices.data() +
1379 dof_info.row_starts[block * n_components].first,
1380 *end_cell = dof_info.dof_indices.data() +
1381 dof_info.row_starts[(block + 1) * n_components].first;
1382 for (; it != end_cell; ++it)
1383 {
1384 SparsityPattern::iterator sp = connectivity_dof.begin(*it);
1385 std::vector<types::global_dof_index>::iterator insert_pos =
1386 row_entries.begin();
1387 for (; sp != connectivity_dof.end(*it); ++sp)
1388 if (sp->column() != block)
1389 row_entries.insert(renumbering[sp->column()], insert_pos);
1390 }
1391 connectivity.add_entries(renumbering[block],
1392 row_entries.begin(),
1393 row_entries.end());
1394 }
1395 }
1396
1397 } // namespace internal
1398
1399 void
1400 DoFInfo::make_connectivity_graph(
1401 const TaskInfo &task_info,
1402 const std::vector<unsigned int> &renumbering,
1403 DynamicSparsityPattern &connectivity) const
1404 {
1405 unsigned int n_rows = (vector_partitioner->local_range().second -
1406 vector_partitioner->local_range().first) +
1407 vector_partitioner->ghost_indices().n_elements();
1408
1409 // Avoid square sparsity patterns that allocate the diagonal entry
1410 if (n_rows == task_info.n_active_cells)
1411 ++n_rows;
1412
1413 // first determine row lengths
1414 std::vector<unsigned int> row_lengths(n_rows);
1415 std::vector<std::mutex> mutexes(n_rows / internal::bucket_size_threading +
1416 1);
1418 0,
1419 task_info.n_active_cells,
1420 [this, &mutexes, &row_lengths](const unsigned int begin,
1421 const unsigned int end) {
1422 internal::compute_row_lengths(
1423 begin, end, *this, mutexes, row_lengths);
1424 },
1425 20);
1426
1427 // disregard dofs that only sit on a single cell because they cannot
1428 // couple
1429 for (unsigned int row = 0; row < n_rows; ++row)
1430 if (row_lengths[row] <= 1)
1431 row_lengths[row] = 0;
1432
1433 // Create a temporary sparsity pattern that holds to each degree of
1434 // freedom on which cells it appears, i.e., store the connectivity
1435 // between cells and dofs
1436 SparsityPattern connectivity_dof(n_rows,
1437 task_info.n_active_cells,
1438 row_lengths);
1440 0,
1441 task_info.n_active_cells,
1442 [this, &row_lengths, &mutexes, &connectivity_dof](
1443 const unsigned int begin, const unsigned int end) {
1444 internal::fill_connectivity_dofs(
1445 begin, end, *this, row_lengths, mutexes, connectivity_dof);
1446 },
1447 20);
1448 connectivity_dof.compress();
1449
1450
1451 // Invert renumbering for use in fill_connectivity.
1452 std::vector<unsigned int> reverse_numbering(task_info.n_active_cells);
1453 reverse_numbering = Utilities::invert_permutation(renumbering);
1454
1455 // From the above connectivity between dofs and cells, we can finally
1456 // create a connectivity list between cells. The connectivity graph
1457 // should apply the renumbering, i.e., the entry for cell j is the entry
1458 // for cell renumbering[j] in the original ordering.
1460 0,
1461 task_info.n_active_cells,
1462 [this, &reverse_numbering, &connectivity_dof, &connectivity](
1463 const unsigned int begin, const unsigned int end) {
1464 internal::fill_connectivity(begin,
1465 end,
1466 *this,
1467 reverse_numbering,
1468 connectivity_dof,
1469 connectivity);
1470 },
1471 20);
1472 }
1473
1474
1475
1476 void
1477 DoFInfo::compute_dof_renumbering(
1478 std::vector<types::global_dof_index> &renumbering)
1479 {
1480 const unsigned int locally_owned_size =
1481 vector_partitioner->locally_owned_size();
1482 renumbering.resize(0);
1484
1485 types::global_dof_index counter = 0;
1486 const unsigned int n_components = start_components.back();
1487 const unsigned int n_cell_batches =
1488 n_vectorization_lanes_filled[dof_access_cell].size();
1489 Assert(n_cell_batches <=
1490 (row_starts.size() - 1) / vectorization_length / n_components,
1492 for (unsigned int cell_no = 0; cell_no < n_cell_batches; ++cell_no)
1493 {
1494 // do not renumber in case we have constraints
1495 if (row_starts[cell_no * n_components * vectorization_length]
1496 .second ==
1497 row_starts[(cell_no + 1) * n_components * vectorization_length]
1498 .second)
1499 {
1500 const unsigned int ndofs =
1501 dofs_per_cell.size() == 1 ?
1502 dofs_per_cell[0] :
1503 (dofs_per_cell[cell_active_fe_index.size() > 0 ?
1504 cell_active_fe_index[cell_no] :
1505 0]);
1506 const unsigned int *dof_ind =
1507 dof_indices.data() +
1508 row_starts[cell_no * n_components * vectorization_length].first;
1509 for (unsigned int i = 0; i < ndofs; ++i)
1510 for (unsigned int j = 0;
1511 j < n_vectorization_lanes_filled[dof_access_cell][cell_no];
1512 ++j)
1513 if (dof_ind[j * ndofs + i] < locally_owned_size)
1514 if (renumbering[dof_ind[j * ndofs + i]] ==
1516 renumbering[dof_ind[j * ndofs + i]] = counter++;
1517 }
1518 }
1519
1521 for (types::global_dof_index &dof_index : renumbering)
1522 if (dof_index == numbers::invalid_dof_index)
1523 dof_index = counter++;
1524
1525 // transform indices to global index space
1526 for (types::global_dof_index &dof_index : renumbering)
1527 dof_index = vector_partitioner->local_to_global(dof_index);
1528
1529 AssertDimension(counter, renumbering.size());
1530 }
1531
1532
1533
1534 std::size_t
1535 DoFInfo::memory_consumption() const
1536 {
1537 std::size_t memory = sizeof(*this);
1538 for (const auto &storage : index_storage_variants)
1539 memory += storage.capacity() * sizeof(storage[0]);
1540 memory +=
1541 (row_starts.capacity() * sizeof(std::pair<unsigned int, unsigned int>));
1542 memory += MemoryConsumption::memory_consumption(dof_indices);
1543 memory += MemoryConsumption::memory_consumption(dof_indices_interleaved);
1544 memory += MemoryConsumption::memory_consumption(dof_indices_contiguous);
1545 memory +=
1546 MemoryConsumption::memory_consumption(dof_indices_contiguous_sm);
1547 memory +=
1548 MemoryConsumption::memory_consumption(dof_indices_interleave_strides);
1549 memory +=
1550 MemoryConsumption::memory_consumption(n_vectorization_lanes_filled);
1552 hanging_node_constraint_masks_comp);
1553 memory +=
1554 MemoryConsumption::memory_consumption(hanging_node_constraint_masks);
1555 memory += MemoryConsumption::memory_consumption(constrained_dofs);
1556 memory += MemoryConsumption::memory_consumption(row_starts_plain_indices);
1557 memory += MemoryConsumption::memory_consumption(plain_dof_indices);
1558 memory += MemoryConsumption::memory_consumption(constraint_indicator);
1559 memory += MemoryConsumption::memory_consumption(*vector_partitioner);
1560 memory += MemoryConsumption::memory_consumption(n_components);
1561 memory += MemoryConsumption::memory_consumption(start_components);
1562 memory += MemoryConsumption::memory_consumption(component_to_base_index);
1563 memory +=
1564 MemoryConsumption::memory_consumption(component_dof_indices_offset);
1565 memory += MemoryConsumption::memory_consumption(dofs_per_cell);
1566 memory += MemoryConsumption::memory_consumption(dofs_per_face);
1567 memory += MemoryConsumption::memory_consumption(cell_active_fe_index);
1568 memory += MemoryConsumption::memory_consumption(fe_index_conversion);
1569 memory +=
1570 MemoryConsumption::memory_consumption(vector_zero_range_list_index);
1571 memory += MemoryConsumption::memory_consumption(vector_zero_range_list);
1572 memory += MemoryConsumption::memory_consumption(cell_loop_pre_list_index);
1573 memory += MemoryConsumption::memory_consumption(cell_loop_pre_list);
1574 memory +=
1575 MemoryConsumption::memory_consumption(cell_loop_post_list_index);
1576 memory += MemoryConsumption::memory_consumption(cell_loop_post_list);
1577 return memory;
1578 }
1579 } // namespace MatrixFreeFunctions
1580} // namespace internal
1581
1582namespace internal
1583{
1584 namespace MatrixFreeFunctions
1585 {
1586 template void
1587 DoFInfo::read_dof_indices<double>(
1588 const std::vector<types::global_dof_index> &,
1589 const std::vector<types::global_dof_index> &,
1590 const bool,
1591 const ::AffineConstraints<double> &,
1592 const unsigned int,
1593 ConstraintValues<double> &,
1594 bool &);
1595
1596 template void
1597 DoFInfo::read_dof_indices<float>(
1598 const std::vector<types::global_dof_index> &,
1599 const std::vector<types::global_dof_index> &,
1600 const bool,
1601 const ::AffineConstraints<float> &,
1602 const unsigned int,
1603 ConstraintValues<double> &,
1604 bool &);
1605
1606 template bool
1607 DoFInfo::process_hanging_node_constraints<1>(
1608 const HangingNodes<1> &,
1609 const std::vector<std::vector<unsigned int>> &,
1610 const unsigned int,
1612 std::vector<types::global_dof_index> &);
1613 template bool
1614 DoFInfo::process_hanging_node_constraints<2>(
1615 const HangingNodes<2> &,
1616 const std::vector<std::vector<unsigned int>> &,
1617 const unsigned int,
1619 std::vector<types::global_dof_index> &);
1620 template bool
1621 DoFInfo::process_hanging_node_constraints<3>(
1622 const HangingNodes<3> &,
1623 const std::vector<std::vector<unsigned int>> &,
1624 const unsigned int,
1626 std::vector<types::global_dof_index> &);
1627
1628 template void
1629 DoFInfo::compute_face_index_compression<1>(
1630 const std::vector<FaceToCellTopology<1>> &,
1631 bool);
1632 template void
1633 DoFInfo::compute_face_index_compression<2>(
1634 const std::vector<FaceToCellTopology<2>> &,
1635 bool);
1636 template void
1637 DoFInfo::compute_face_index_compression<4>(
1638 const std::vector<FaceToCellTopology<4>> &,
1639 bool);
1640 template void
1641 DoFInfo::compute_face_index_compression<8>(
1642 const std::vector<FaceToCellTopology<8>> &,
1643 bool);
1644 template void
1645 DoFInfo::compute_face_index_compression<16>(
1646 const std::vector<FaceToCellTopology<16>> &,
1647 bool);
1648
1649 template void
1650 DoFInfo::compute_vector_zero_access_pattern<1>(
1651 const TaskInfo &,
1652 const std::vector<FaceToCellTopology<1>> &);
1653 template void
1654 DoFInfo::compute_vector_zero_access_pattern<2>(
1655 const TaskInfo &,
1656 const std::vector<FaceToCellTopology<2>> &);
1657 template void
1658 DoFInfo::compute_vector_zero_access_pattern<4>(
1659 const TaskInfo &,
1660 const std::vector<FaceToCellTopology<4>> &);
1661 template void
1662 DoFInfo::compute_vector_zero_access_pattern<8>(
1663 const TaskInfo &,
1664 const std::vector<FaceToCellTopology<8>> &);
1665 template void
1666 DoFInfo::compute_vector_zero_access_pattern<16>(
1667 const TaskInfo &,
1668 const std::vector<FaceToCellTopology<16>> &);
1669
1670 template void
1671 DoFInfo::print_memory_consumption<std::ostream>(std::ostream &,
1672 const TaskInfo &) const;
1673 template void
1674 DoFInfo::print_memory_consumption<ConditionalOStream>(
1676 const TaskInfo &) const;
1677
1678 template void
1679 DoFInfo::print<double>(const std::vector<double> &,
1680 const std::vector<unsigned int> &,
1681 std::ostream &) const;
1682
1683 template void
1684 DoFInfo::print<float>(const std::vector<float> &,
1685 const std::vector<unsigned int> &,
1686 std::ostream &) const;
1687 } // namespace MatrixFreeFunctions
1688} // namespace internal
1689
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
size_type index_within_set(const size_type global_index) const
Definition index_set.h:2013
size_type n_elements() const
Definition index_set.h:1945
void subtract_set(const IndexSet &other)
Definition index_set.cc:478
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1814
void compress() const
Definition index_set.h:1795
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1842
void add(const size_type i, const size_type j)
const IndexSet & locally_owned_range() const
const IndexSet & ghost_indices() const
unsigned int locally_owned_size() const
types::global_dof_index local_to_global(const unsigned int local_index) const
virtual MPI_Comm get_mpi_communicator() const override
void set_ghost_indices(const IndexSet &ghost_indices, const IndexSet &larger_ghost_index_set=IndexSet())
types::global_dof_index size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
Point< 2 > second
Definition grid_out.cc:4633
Point< 2 > first
Definition grid_out.cc:4632
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
std::size_t size
Definition mpi.cc:745
types::global_dof_index locally_owned_size
Definition mpi.cc:833
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition utilities.h:1700
void fill_connectivity_dofs(const unsigned int begin, const unsigned int end, const DoFInfo &dof_info, const std::vector< unsigned int > &row_lengths, std::vector< std::mutex > &mutexes, ::SparsityPattern &connectivity_dof)
Definition dof_info.cc:1325
void fill_connectivity(const unsigned int begin, const unsigned int end, const DoFInfo &dof_info, const std::vector< unsigned int > &renumbering, const ::SparsityPattern &connectivity_dof, DynamicSparsityPattern &connectivity)
Definition dof_info.cc:1364
void compute_row_lengths(const unsigned int begin, const unsigned int end, const DoFInfo &dof_info, std::vector< std::mutex > &mutexes, std::vector< unsigned int > &row_lengths)
Definition dof_info.cc:1285
static constexpr unsigned int bucket_size_threading
Definition dof_info.cc:1280
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
constexpr types::global_dof_index invalid_dof_index
Definition types.h:269
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
void apply_to_subranges(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, const Function &f, const unsigned int grainsize)
Definition parallel.h:540
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
void reorder_cells(const TaskInfo &task_info, const std::vector< unsigned int > &renumbering, const std::vector< unsigned int > &constraint_pool_row_index, const std::vector< unsigned char > &irregular_cells)
Definition dof_info.cc:301
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
Definition dof_info.h:538
std::vector< std::pair< unsigned int, unsigned int > > row_starts
Definition dof_info.h:497
void compute_tight_partitioners(const Table< 2, ShapeInfo< double > > &shape_info, const unsigned int n_owned_cells, const unsigned int n_lanes, const std::vector< FaceToCellTopology< 1 > > &inner_faces, const std::vector< FaceToCellTopology< 1 > > &ghosted_faces, const bool fill_cell_centric, const MPI_Comm communicator_sm, const bool use_vector_data_exchanger_full)
Definition dof_info.cc:903
std::vector< std::vector< bool > > hanging_node_constraint_masks_comp
Definition dof_info.h:520
void get_dof_indices_on_cell_batch(std::vector< unsigned int > &local_indices, const unsigned int cell_batch, const bool with_constraints=true) const
Definition dof_info.cc:84
void compute_cell_index_compression(const std::vector< unsigned char > &irregular_cells)
Definition dof_info.cc:538
std::vector< unsigned int > dofs_per_cell
Definition dof_info.h:695
void assign_ghosts(const std::vector< unsigned int > &boundary_cells, const MPI_Comm communicator_sm, const bool use_vector_data_exchanger_full)
Definition dof_info.cc:157
std::shared_ptr< const internal::MatrixFreeFunctions::VectorDataExchange::Base > vector_exchanger
Definition dof_info.h:600
std::vector< std::vector< unsigned int > > fe_index_conversion
Definition dof_info.h:723
std::vector< unsigned int > dof_indices
Definition dof_info.h:514
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:593
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
Definition dof_info.h:526
std::array< std::vector< unsigned int >, 3 > dof_indices_interleave_strides
Definition dof_info.h:574
std::vector< unsigned int > row_starts_plain_indices
Definition dof_info.h:637
std::vector< unsigned int > dofs_per_face
Definition dof_info.h:700
std::vector< unsigned int > n_components
Definition dof_info.h:665
std::vector< types::global_dof_index > ghost_dofs
Definition dof_info.h:730
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
Definition dof_info.h:553
std::vector< unsigned int > cell_active_fe_index
Definition dof_info.h:710
std::array< std::shared_ptr< const internal::MatrixFreeFunctions::VectorDataExchange::Base >, 5 > vector_exchanger_face_variants
Definition dof_info.h:625
std::vector< unsigned int > plain_dof_indices
Definition dof_info.h:647
std::array< std::vector< unsigned char >, 3 > n_vectorization_lanes_filled
Definition dof_info.h:585
std::vector< unsigned int > start_components
Definition dof_info.h:671
std::vector< unsigned int > dof_indices_interleaved
Definition dof_info.h:543
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
Definition dof_info.h:489
std::vector< unsigned int > cell_partition_data
Definition task_info.h:474