Reference documentation for deal.II version 9.3.3
\(\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\}}\)
dof_info.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 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
17#include <deal.II/base/conditional_ostream.h>
18#include <deal.II/base/utilities.h>
19#include <deal.II/base/vectorization.h>
20
21#include <deal.II/matrix_free/dof_info.templates.h>
22
23#include <iostream>
24
25DEAL_II_NAMESPACE_OPEN
26
27namespace internal
28{
29 namespace MatrixFreeFunctions
30 {
31 DoFInfo::DoFInfo()
32 {
33 clear();
34 }
35
36
37
38 void
39 DoFInfo::clear()
40 {
41 row_starts.clear();
42 dof_indices.clear();
43 constraint_indicator.clear();
44 vector_partitioner.reset();
45 ghost_dofs.clear();
46 dofs_per_cell.clear();
47 dofs_per_face.clear();
48 vectorization_length = 1;
49 dimension = 2;
50 global_base_element_offset = 0;
51 n_base_elements = 0;
52 n_components.clear();
53 start_components.clear();
54 row_starts_plain_indices.clear();
55 plain_dof_indices.clear();
56 dof_indices_interleaved.clear();
57 for (unsigned int i = 0; i < 3; ++i)
58 {
59 index_storage_variants[i].clear();
60 dof_indices_contiguous[i].clear();
61 dof_indices_interleave_strides[i].clear();
62 n_vectorization_lanes_filled[i].clear();
63 }
64 store_plain_indices = false;
65 cell_active_fe_index.clear();
66 max_fe_index = 0;
67 fe_index_conversion.clear();
68 }
69
70
71
72 void
73 DoFInfo::get_dof_indices_on_cell_batch(std::vector<unsigned int> &my_rows,
74 const unsigned int cell,
75 const bool apply_constraints) const
76 {
77 const unsigned int n_fe_components = start_components.back();
78 const unsigned int fe_index =
79 dofs_per_cell.size() == 1 ? 0 : cell_active_fe_index[cell];
80 const unsigned int dofs_this_cell = dofs_per_cell[fe_index];
81
82 const unsigned int n_vectorization = vectorization_length;
83 constexpr auto dof_access_index = dof_access_cell;
84 AssertIndexRange(cell,
85 n_vectorization_lanes_filled[dof_access_index].size());
86 const unsigned int n_vectorization_actual =
87 n_vectorization_lanes_filled[dof_access_index][cell];
88
89 // we might have constraints, so the final number
90 // of indices is not known a priori.
91 // conservatively reserve the maximum without constraints
92 my_rows.reserve(n_vectorization * dofs_this_cell);
93 my_rows.resize(0);
94 unsigned int total_size = 0;
95 for (unsigned int v = 0; v < n_vectorization_actual; ++v)
96 {
97 const unsigned int ib =
98 (cell * n_vectorization + v) * n_fe_components;
99 const unsigned int ie =
100 (cell * n_vectorization + v + 1) * n_fe_components;
101
102 // figure out constraints by comparing constraint_indicator row
103 // shift for this cell within the block as compared to the next
104 // one
105 const bool has_constraints =
106 row_starts[ib].second != row_starts[ib + n_fe_components].second;
107
108 auto do_copy = [&](const unsigned int *begin,
109 const unsigned int *end) {
110 const unsigned int shift = total_size;
111 total_size += (end - begin);
112 my_rows.resize(total_size);
113 std::copy(begin, end, my_rows.begin() + shift);
114 };
115
116 if (!has_constraints || apply_constraints)
117 {
118 const unsigned int *begin =
119 dof_indices.data() + row_starts[ib].first;
120 const unsigned int *end =
121 dof_indices.data() + row_starts[ie].first;
122 do_copy(begin, end);
123 }
124 else
125 {
126 Assert(row_starts_plain_indices[cell * n_vectorization + v] !=
127 numbers::invalid_unsigned_int,
128 ExcNotInitialized());
129 const unsigned int *begin =
130 plain_dof_indices.data() +
131 row_starts_plain_indices[cell * n_vectorization + v];
132 const unsigned int *end = begin + dofs_this_cell;
133 do_copy(begin, end);
134 }
135 }
136 }
137
138
139
140 void
141 DoFInfo::assign_ghosts(const std::vector<unsigned int> &boundary_cells,
142 const MPI_Comm & communicator_sm,
143 const bool use_vector_data_exchanger_full)
144 {
145 Assert(boundary_cells.size() < row_starts.size(), ExcInternalError());
146
147 // sort ghost dofs and compress out duplicates
148 const unsigned int n_owned = (vector_partitioner->local_range().second -
149 vector_partitioner->local_range().first);
150 const std::size_t n_ghosts = ghost_dofs.size();
151#ifdef DEBUG
152 for (const auto dof_index : dof_indices)
153 AssertIndexRange(dof_index, n_owned + n_ghosts);
154#endif
155
156 const unsigned int n_components = start_components.back();
157 std::vector<unsigned int> ghost_numbering(n_ghosts);
158 IndexSet ghost_indices(vector_partitioner->size());
159 if (n_ghosts > 0)
160 {
161 unsigned int n_unique_ghosts = 0;
162 // since we need to go back to the local_to_global indices and
163 // replace the temporary numbering of ghosts by the real number in
164 // the index set, we need to store these values
165 std::vector<std::pair<types::global_dof_index, unsigned int>>
166 ghost_origin(n_ghosts);
167 for (std::size_t i = 0; i < n_ghosts; ++i)
168 {
169 ghost_origin[i].first = ghost_dofs[i];
170 ghost_origin[i].second = i;
171 }
172 std::sort(ghost_origin.begin(), ghost_origin.end());
173
174 types::global_dof_index last_contiguous_start = ghost_origin[0].first;
175 ghost_numbering[ghost_origin[0].second] = 0;
176 for (std::size_t i = 1; i < n_ghosts; i++)
177 {
178 if (ghost_origin[i].first > ghost_origin[i - 1].first + 1)
179 {
180 ghost_indices.add_range(last_contiguous_start,
181 ghost_origin[i - 1].first + 1);
182 last_contiguous_start = ghost_origin[i].first;
183 }
184 if (ghost_origin[i].first > ghost_origin[i - 1].first)
185 ++n_unique_ghosts;
186 ghost_numbering[ghost_origin[i].second] = n_unique_ghosts;
187 }
188 ++n_unique_ghosts;
189 ghost_indices.add_range(last_contiguous_start,
190 ghost_origin.back().first + 1);
191 ghost_indices.compress();
192
193 // make sure that we got the correct local numbering of the ghost
194 // dofs. the ghost index set should store the same number
195 {
196 AssertDimension(n_unique_ghosts, ghost_indices.n_elements());
197 for (std::size_t i = 0; i < n_ghosts; ++i)
198 Assert(ghost_numbering[i] ==
199 ghost_indices.index_within_set(ghost_dofs[i]),
200 ExcInternalError());
201 }
202
203 // apply correct numbering for ghost indices: We previously just
204 // enumerated them according to their appearance in the
205 // local_to_global structure. Above, we derived a relation between
206 // this enumeration and the actual number
207 const unsigned int n_boundary_cells = boundary_cells.size();
208 for (unsigned int i = 0; i < n_boundary_cells; ++i)
209 {
210 unsigned int *data_ptr =
211 dof_indices.data() +
212 row_starts[boundary_cells[i] * n_components].first;
213 const unsigned int *row_end =
214 dof_indices.data() +
215 row_starts[(boundary_cells[i] + 1) * n_components].first;
216 for (; data_ptr != row_end; ++data_ptr)
217 *data_ptr = ((*data_ptr < n_owned) ?
218 *data_ptr :
219 n_owned + ghost_numbering[*data_ptr - n_owned]);
220
221 // now the same procedure for plain indices
222 if (store_plain_indices == true)
223 {
224 if (row_starts[boundary_cells[i] * n_components].second !=
225 row_starts[(boundary_cells[i] + 1) * n_components].second)
226 {
227 unsigned int *data_ptr =
228 plain_dof_indices.data() +
229 row_starts_plain_indices[boundary_cells[i]];
230 const unsigned int fe_index =
231 (cell_active_fe_index.size() == 0 ||
232 dofs_per_cell.size() == 1) ?
233 0 :
234 cell_active_fe_index[boundary_cells[i]];
235 AssertIndexRange(fe_index, dofs_per_cell.size());
236 const unsigned int *row_end =
237 data_ptr + dofs_per_cell[fe_index];
238 for (; data_ptr != row_end; ++data_ptr)
239 *data_ptr =
240 ((*data_ptr < n_owned) ?
241 *data_ptr :
242 n_owned + ghost_numbering[*data_ptr - n_owned]);
243 }
244 }
245 }
246 }
247
248 std::vector<types::global_dof_index> empty;
249 ghost_dofs.swap(empty);
250
251 // set the ghost indices now. need to cast away constness here, but that
252 // is uncritical since we reset the Partitioner in the same initialize
253 // call as this call here.
254 Utilities::MPI::Partitioner *vec_part =
255 const_cast<Utilities::MPI::Partitioner *>(vector_partitioner.get());
256 vec_part->set_ghost_indices(ghost_indices);
257
258 if (use_vector_data_exchanger_full == false)
259 vector_exchanger = std::make_shared<
260 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
261 vector_partitioner);
262 else
263 vector_exchanger =
264 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
265 vector_partitioner, communicator_sm);
266 }
267
268
269
270 void
271 DoFInfo::reorder_cells(
272 const TaskInfo & task_info,
273 const std::vector<unsigned int> & renumbering,
274 const std::vector<unsigned int> & constraint_pool_row_index,
275 const std::vector<unsigned char> &irregular_cells)
276 {
277 (void)constraint_pool_row_index;
278
279 // first reorder the active FE index.
280 const bool have_hp = dofs_per_cell.size() > 1;
281 if (cell_active_fe_index.size() > 0)
282 {
283 std::vector<unsigned int> new_active_fe_index;
284 new_active_fe_index.reserve(task_info.cell_partition_data.back());
285 unsigned int position_cell = 0;
286 for (unsigned int cell = 0;
287 cell < task_info.cell_partition_data.back();
288 ++cell)
289 {
290 const unsigned int n_comp =
291 (irregular_cells[cell] > 0 ? irregular_cells[cell] :
292 vectorization_length);
293
294 // take maximum FE index among the ones present (we might have
295 // lumped some lower indices into higher ones)
296 unsigned int fe_index =
297 cell_active_fe_index[renumbering[position_cell]];
298 for (unsigned int j = 1; j < n_comp; ++j)
299 fe_index = std::max(
300 fe_index,
301 cell_active_fe_index[renumbering[position_cell + j]]);
302
303 new_active_fe_index.push_back(fe_index);
304 position_cell += n_comp;
305 }
306 std::swap(new_active_fe_index, cell_active_fe_index);
307 }
308 if (have_hp)
309 AssertDimension(cell_active_fe_index.size(),
310 task_info.cell_partition_data.back());
311
312 const unsigned int n_components = start_components.back();
313
314 std::vector<std::pair<unsigned int, unsigned int>> new_row_starts(
315 vectorization_length * n_components *
316 task_info.cell_partition_data.back() +
317 1);
318 std::vector<unsigned int> new_dof_indices;
319 std::vector<std::pair<unsigned short, unsigned short>>
320 new_constraint_indicator;
321 std::vector<unsigned int> new_plain_indices, new_rowstart_plain;
322 unsigned int position_cell = 0;
323 new_dof_indices.reserve(dof_indices.size());
324 new_constraint_indicator.reserve(constraint_indicator.size());
325 if (store_plain_indices == true)
326 {
327 new_rowstart_plain.resize(vectorization_length *
328 task_info.cell_partition_data.back() +
329 1,
330 numbers::invalid_unsigned_int);
331 new_plain_indices.reserve(plain_dof_indices.size());
332 }
333
334 // copy the indices and the constraint indicators to the new data field,
335 // where we will go through the cells in the renumbered way. in case the
336 // vectorization length does not exactly match up, we fill invalid
337 // numbers to the rowstart data. for contiguous cell indices, we skip
338 // the rowstarts field completely and directly go into the
339 // new_dof_indices field (this layout is used in FEEvaluation).
340 for (unsigned int i = 0; i < task_info.cell_partition_data.back(); ++i)
341 {
342 const unsigned int n_vect =
343 (irregular_cells[i] > 0 ? irregular_cells[i] :
344 vectorization_length);
345 const unsigned int dofs_per_cell =
346 have_hp ? this->dofs_per_cell[cell_active_fe_index[i]] :
347 this->dofs_per_cell[0];
348
349 for (unsigned int j = 0; j < n_vect; ++j)
350 {
351 const unsigned int cell_no =
352 renumbering[position_cell + j] * n_components;
353 for (unsigned int comp = 0; comp < n_components; ++comp)
354 {
355 new_row_starts[(i * vectorization_length + j) * n_components +
356 comp]
357 .first = new_dof_indices.size();
358 new_row_starts[(i * vectorization_length + j) * n_components +
359 comp]
360 .second = new_constraint_indicator.size();
361
362 new_dof_indices.insert(
363 new_dof_indices.end(),
364 dof_indices.data() + row_starts[cell_no + comp].first,
365 dof_indices.data() + row_starts[cell_no + comp + 1].first);
366 for (unsigned int index = row_starts[cell_no + comp].second;
367 index != row_starts[cell_no + comp + 1].second;
368 ++index)
369 new_constraint_indicator.push_back(
370 constraint_indicator[index]);
371 }
372 if (store_plain_indices &&
373 row_starts[cell_no].second !=
374 row_starts[cell_no + n_components].second)
375 {
376 new_rowstart_plain[i * vectorization_length + j] =
377 new_plain_indices.size();
378 new_plain_indices.insert(
379 new_plain_indices.end(),
380 plain_dof_indices.data() +
381 row_starts_plain_indices[cell_no / n_components],
382 plain_dof_indices.data() +
383 row_starts_plain_indices[cell_no / n_components] +
384 dofs_per_cell);
385 }
386 }
387 for (unsigned int j = n_vect; j < vectorization_length; ++j)
388 for (unsigned int comp = 0; comp < n_components; ++comp)
389 {
390 new_row_starts[(i * vectorization_length + j) * n_components +
391 comp]
392 .first = new_dof_indices.size();
393 new_row_starts[(i * vectorization_length + j) * n_components +
394 comp]
395 .second = new_constraint_indicator.size();
396 }
397 position_cell += n_vect;
398 }
399 AssertDimension(position_cell * n_components + 1, row_starts.size());
400
401 AssertDimension(dof_indices.size(), new_dof_indices.size());
402 new_row_starts[task_info.cell_partition_data.back() *
403 vectorization_length * n_components]
404 .first = new_dof_indices.size();
405 new_row_starts[task_info.cell_partition_data.back() *
406 vectorization_length * n_components]
407 .second = new_constraint_indicator.size();
408
409 AssertDimension(constraint_indicator.size(),
410 new_constraint_indicator.size());
411
412 new_row_starts.swap(row_starts);
413 new_dof_indices.swap(dof_indices);
414 new_constraint_indicator.swap(constraint_indicator);
415 new_plain_indices.swap(plain_dof_indices);
416 new_rowstart_plain.swap(row_starts_plain_indices);
417
418#ifdef DEBUG
419 // sanity check 1: all indices should be smaller than the number of dofs
420 // locally owned plus the number of ghosts
421 const unsigned int index_range =
422 (vector_partitioner->local_range().second -
423 vector_partitioner->local_range().first) +
424 vector_partitioner->ghost_indices().n_elements();
425 for (const auto dof_index : dof_indices)
426 AssertIndexRange(dof_index, index_range);
427
428 // sanity check 2: for the constraint indicators, the first index should
429 // be smaller than the number of indices in the row, and the second
430 // index should be smaller than the number of constraints in the
431 // constraint pool.
432 for (unsigned int row = 0; row < task_info.cell_partition_data.back();
433 ++row)
434 {
435 const unsigned int row_length_ind =
436 row_starts[(row * vectorization_length + 1) * n_components].first -
437 row_starts[row * vectorization_length * n_components].first;
438 AssertIndexRange(
439 row_starts[(row * vectorization_length + 1) * n_components].second,
440 constraint_indicator.size() + 1);
441 const std::pair<unsigned short, unsigned short> *
442 con_it =
443 constraint_indicator.data() +
444 row_starts[row * vectorization_length * n_components].second,
445 *end_con =
446 constraint_indicator.data() +
447 row_starts[(row * vectorization_length + 1) * n_components].second;
448 for (; con_it != end_con; ++con_it)
449 {
450 AssertIndexRange(con_it->first, row_length_ind + 1);
451 AssertIndexRange(con_it->second,
452 constraint_pool_row_index.size() - 1);
453 }
454 }
455
456 // sanity check 3: check the number of cells once again
457 unsigned int n_active_cells = 0;
458 for (unsigned int c = 0; c < *(task_info.cell_partition_data.end() - 2);
459 ++c)
460 if (irregular_cells[c] > 0)
461 n_active_cells += irregular_cells[c];
462 else
463 n_active_cells += vectorization_length;
464 AssertDimension(n_active_cells, task_info.n_active_cells);
465#endif
466
467 compute_cell_index_compression(irregular_cells);
468 }
469
470
471
472 void
473 DoFInfo::compute_cell_index_compression(
474 const std::vector<unsigned char> &irregular_cells)
475 {
476 const bool have_hp = dofs_per_cell.size() > 1;
477 const unsigned int n_components = start_components.back();
478
479 Assert(vectorization_length == 1 ||
480 row_starts.size() % vectorization_length == 1,
481 ExcInternalError());
482 if (vectorization_length > 1)
483 AssertDimension(row_starts.size() / vectorization_length / n_components,
484 irregular_cells.size());
485 index_storage_variants[dof_access_cell].resize(
486 irregular_cells.size(), IndexStorageVariants::full);
487 n_vectorization_lanes_filled[dof_access_cell].resize(
488 irregular_cells.size());
489 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
490 if (irregular_cells[i] > 0)
491 n_vectorization_lanes_filled[dof_access_cell][i] = irregular_cells[i];
492 else
493 n_vectorization_lanes_filled[dof_access_cell][i] =
494 vectorization_length;
495
496 dof_indices_contiguous[dof_access_cell].resize(
497 irregular_cells.size() * vectorization_length,
498 numbers::invalid_unsigned_int);
499 dof_indices_interleaved.resize(dof_indices.size(),
500 numbers::invalid_unsigned_int);
501 dof_indices_interleave_strides[dof_access_cell].resize(
502 irregular_cells.size() * vectorization_length,
503 numbers::invalid_unsigned_int);
504
505 std::vector<unsigned int> index_kinds(
506 static_cast<unsigned int>(
507 IndexStorageVariants::interleaved_contiguous_mixed_strides) +
508 1);
509 std::vector<unsigned int> offsets(vectorization_length);
510 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
511 {
512 const unsigned int ndofs =
513 dofs_per_cell[have_hp ? cell_active_fe_index[i] : 0];
514 const unsigned int n_comp =
515 n_vectorization_lanes_filled[dof_access_cell][i];
516
517 // check 1: Check if there are constraints -> no compression possible
518 bool has_constraints = false;
519 for (unsigned int j = 0; j < n_comp; ++j)
520 {
521 const unsigned int cell_no = i * vectorization_length + j;
522 if (row_starts[cell_no * n_components].second !=
523 row_starts[(cell_no + 1) * n_components].second)
524 {
525 has_constraints = true;
526 break;
527 }
528 }
529 if (has_constraints)
530 index_storage_variants[dof_access_cell][i] =
531 IndexStorageVariants::full;
532 else
533 {
534 bool indices_are_contiguous = true;
535 for (unsigned int j = 0; j < n_comp; ++j)
536 {
537 const unsigned int cell_no = i * vectorization_length + j;
538 const unsigned int *dof_indices =
539 this->dof_indices.data() +
540 row_starts[cell_no * n_components].first;
541 AssertDimension(
542 ndofs,
543 row_starts[(cell_no + 1) * n_components].first -
544 row_starts[cell_no * n_components].first);
545 for (unsigned int i = 1; i < ndofs; ++i)
546 if (dof_indices[i] != dof_indices[0] + i)
547 {
548 indices_are_contiguous = false;
549 break;
550 }
551 }
552
553 bool indices_are_interleaved_and_contiguous =
554 (ndofs > 1 && n_comp == vectorization_length);
555
556 {
557 const unsigned int *dof_indices =
558 this->dof_indices.data() +
559 row_starts[i * vectorization_length * n_components].first;
560 for (unsigned int k = 0; k < ndofs; ++k)
561 for (unsigned int j = 0; j < n_comp; ++j)
562 if (dof_indices[j * ndofs + k] !=
563 dof_indices[0] + k * n_comp + j)
564 {
565 indices_are_interleaved_and_contiguous = false;
566 break;
567 }
568 }
569
570 if (indices_are_contiguous ||
571 indices_are_interleaved_and_contiguous)
572 {
573 for (unsigned int j = 0; j < n_comp; ++j)
574 dof_indices_contiguous
575 [dof_access_cell][i * vectorization_length + j] =
576 this->dof_indices[row_starts[(i * vectorization_length +
577 j) *
578 n_components]
579 .first];
580 }
581
582 if (indices_are_interleaved_and_contiguous)
583 {
584 Assert(n_comp == vectorization_length, ExcInternalError());
585 index_storage_variants[dof_access_cell][i] =
586 IndexStorageVariants::interleaved_contiguous;
587 for (unsigned int j = 0; j < n_comp; ++j)
588 dof_indices_interleave_strides[2][i * vectorization_length +
589 j] = n_comp;
590 }
591 else if (indices_are_contiguous)
592 {
593 index_storage_variants[dof_access_cell][i] =
594 IndexStorageVariants::contiguous;
595 for (unsigned int j = 0; j < n_comp; ++j)
596 dof_indices_interleave_strides[2][i * vectorization_length +
597 j] = 1;
598 }
599 else
600 {
601 int indices_are_interleaved_and_mixed = 2;
602 const unsigned int *dof_indices =
603 &this->dof_indices[row_starts[i * vectorization_length *
604 n_components]
605 .first];
606 for (unsigned int j = 0; j < n_comp; ++j)
607 offsets[j] =
608 dof_indices[j * ndofs + 1] - dof_indices[j * ndofs];
609 for (unsigned int k = 0; k < ndofs; ++k)
610 for (unsigned int j = 0; j < n_comp; ++j)
611 // the first if case is to avoid negative offsets
612 // (invalid)
613 if (dof_indices[j * ndofs + 1] < dof_indices[j * ndofs] ||
614 dof_indices[j * ndofs + k] !=
615 dof_indices[j * ndofs] + k * offsets[j])
616 {
617 indices_are_interleaved_and_mixed = 0;
618 break;
619 }
620 if (indices_are_interleaved_and_mixed == 2)
621 {
622 for (unsigned int j = 0; j < n_comp; ++j)
623 dof_indices_interleave_strides
624 [dof_access_cell][i * vectorization_length + j] =
625 offsets[j];
626 for (unsigned int j = 0; j < n_comp; ++j)
627 dof_indices_contiguous[dof_access_cell]
628 [i * vectorization_length + j] =
629 dof_indices[j * ndofs];
630 for (unsigned int j = 0; j < n_comp; ++j)
631 if (offsets[j] != vectorization_length)
632 {
633 indices_are_interleaved_and_mixed = 1;
634 break;
635 }
636 if (indices_are_interleaved_and_mixed == 1 ||
637 n_comp != vectorization_length)
638 index_storage_variants[dof_access_cell][i] =
639 IndexStorageVariants::
640 interleaved_contiguous_mixed_strides;
641 else
642 index_storage_variants[dof_access_cell][i] =
643 IndexStorageVariants::interleaved_contiguous_strided;
644 }
645 else
646 {
647 const unsigned int *dof_indices =
648 this->dof_indices.data() +
649 row_starts[i * vectorization_length * n_components]
650 .first;
651 if (n_comp == vectorization_length)
652 index_storage_variants[dof_access_cell][i] =
653 IndexStorageVariants::interleaved;
654 else
655 index_storage_variants[dof_access_cell][i] =
656 IndexStorageVariants::full;
657
658 // do not use interleaved storage if two vectorized
659 // components point to the same field (scatter not
660 // possible)
661 for (unsigned int k = 0; k < ndofs; ++k)
662 for (unsigned int l = 0; l < n_comp; ++l)
663 for (unsigned int j = l + 1; j < n_comp; ++j)
664 if (dof_indices[j * ndofs + k] ==
665 dof_indices[l * ndofs + k])
666 {
667 index_storage_variants[dof_access_cell][i] =
668 IndexStorageVariants::full;
669 break;
670 }
671 }
672 }
673 }
674 index_kinds[static_cast<unsigned int>(
675 index_storage_variants[dof_access_cell][i])]++;
676 }
677
678 // Cleanup phase: we want to avoid single cells with different properties
679 // than the bulk of the domain in order to avoid extra checks in the face
680 // identification.
681
682 // Step 1: check whether the interleaved indices were only assigned to
683 // the single cell within a vectorized array.
684 auto fix_single_interleaved_indices =
685 [&](const IndexStorageVariants variant) {
686 if (index_kinds[static_cast<unsigned int>(
687 IndexStorageVariants::interleaved_contiguous_mixed_strides)] >
688 0 &&
689 index_kinds[static_cast<unsigned int>(variant)] > 0)
690 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
691 {
692 if (index_storage_variants[dof_access_cell][i] ==
693 IndexStorageVariants::
694 interleaved_contiguous_mixed_strides &&
695 n_vectorization_lanes_filled[dof_access_cell][i] == 1 &&
696 (variant != IndexStorageVariants::contiguous ||
697 dof_indices_interleave_strides[dof_access_cell]
698 [i * vectorization_length] ==
699 1))
700 {
701 index_storage_variants[dof_access_cell][i] = variant;
702 index_kinds[static_cast<unsigned int>(
703 IndexStorageVariants::
704 interleaved_contiguous_mixed_strides)]--;
705 index_kinds[static_cast<unsigned int>(variant)]++;
706 }
707 }
708 };
709
710 fix_single_interleaved_indices(IndexStorageVariants::full);
711 fix_single_interleaved_indices(IndexStorageVariants::contiguous);
712 fix_single_interleaved_indices(IndexStorageVariants::interleaved);
713
714 unsigned int n_interleaved =
715 index_kinds[static_cast<unsigned int>(
716 IndexStorageVariants::interleaved_contiguous)] +
717 index_kinds[static_cast<unsigned int>(
718 IndexStorageVariants::interleaved_contiguous_strided)] +
719 index_kinds[static_cast<unsigned int>(
720 IndexStorageVariants::interleaved_contiguous_mixed_strides)];
721
722 // Step 2: fix single contiguous cell among others with interleaved
723 // storage
724 if (n_interleaved > 0 && index_kinds[static_cast<unsigned int>(
725 IndexStorageVariants::contiguous)] > 0)
726 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
727 if (index_storage_variants[dof_access_cell][i] ==
728 IndexStorageVariants::contiguous)
729 {
730 index_storage_variants[dof_access_cell][i] =
731 IndexStorageVariants::interleaved_contiguous_mixed_strides;
732 index_kinds[static_cast<unsigned int>(
733 IndexStorageVariants::contiguous)]--;
734 index_kinds[static_cast<unsigned int>(
735 IndexStorageVariants::interleaved_contiguous_mixed_strides)]++;
736 }
737
738 // Step 3: Interleaved cells are left but also some non-contiguous ones
739 // -> revert all to full storage
740 if (n_interleaved > 0 &&
741 index_kinds[static_cast<unsigned int>(IndexStorageVariants::full)] +
742 index_kinds[static_cast<unsigned int>(
743 IndexStorageVariants::interleaved)] >
744 0)
745 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
746 if (index_storage_variants[dof_access_cell][i] >
747 IndexStorageVariants::contiguous)
748 {
749 index_kinds[static_cast<unsigned int>(
750 index_storage_variants[2][i])]--;
751 if (n_vectorization_lanes_filled[dof_access_cell][i] ==
752 vectorization_length)
753 index_storage_variants[dof_access_cell][i] =
754 IndexStorageVariants::interleaved;
755 else
756 index_storage_variants[dof_access_cell][i] =
757 IndexStorageVariants::full;
758 index_kinds[static_cast<unsigned int>(
759 index_storage_variants[dof_access_cell][i])]++;
760 }
761
762 // Step 4: Copy the interleaved indices into their own data structure
763 for (unsigned int i = 0; i < irregular_cells.size(); ++i)
764 if (index_storage_variants[dof_access_cell][i] ==
765 IndexStorageVariants::interleaved)
766 {
767 if (n_vectorization_lanes_filled[dof_access_cell][i] <
768 vectorization_length)
769 {
770 index_storage_variants[dof_access_cell][i] =
771 IndexStorageVariants::full;
772 continue;
773 }
774 const unsigned int ndofs =
775 dofs_per_cell[have_hp ? cell_active_fe_index[i] : 0];
776 const unsigned int *dof_indices =
777 &this->dof_indices
778 [row_starts[i * vectorization_length * n_components].first];
779 unsigned int *interleaved_dof_indices =
780 &this->dof_indices_interleaved
781 [row_starts[i * vectorization_length * n_components].first];
782 AssertDimension(this->dof_indices.size(),
783 this->dof_indices_interleaved.size());
784 AssertDimension(n_vectorization_lanes_filled[dof_access_cell][i],
785 vectorization_length);
786 AssertIndexRange(
787 row_starts[i * vectorization_length * n_components].first,
788 this->dof_indices_interleaved.size());
789 AssertIndexRange(
790 row_starts[i * vectorization_length * n_components].first +
791 ndofs * vectorization_length,
792 this->dof_indices_interleaved.size() + 1);
793 for (unsigned int k = 0; k < ndofs; ++k)
794 for (unsigned int j = 0; j < vectorization_length; ++j)
795 interleaved_dof_indices[k * vectorization_length + j] =
796 dof_indices[j * ndofs + k];
797 }
798 }
799
800
801
802 void
803 DoFInfo::compute_tight_partitioners(
804 const Table<2, ShapeInfo<double>> & shape_info,
805 const unsigned int n_owned_cells,
806 const unsigned int n_lanes,
807 const std::vector<FaceToCellTopology<1>> &inner_faces,
808 const std::vector<FaceToCellTopology<1>> &ghosted_faces,
809 const bool fill_cell_centric,
810 const MPI_Comm & communicator_sm,
811 const bool use_vector_data_exchanger_full)
812 {
813 const Utilities::MPI::Partitioner &part = *vector_partitioner;
814
815 // partitioner 0: no face integrals, simply use the indices present
816 // on the cells
817 std::vector<types::global_dof_index> ghost_indices;
818 {
819 const unsigned int n_components = start_components.back();
820 for (unsigned int cell = 0; cell < n_owned_cells; ++cell)
821 {
822 for (unsigned int i = row_starts[cell * n_components].first;
823 i < row_starts[(cell + 1) * n_components].first;
824 ++i)
825 if (dof_indices[i] >= part.local_size())
826 ghost_indices.push_back(part.local_to_global(dof_indices[i]));
827
828 const unsigned int fe_index =
829 dofs_per_cell.size() == 1 ? 0 :
830 cell_active_fe_index[cell / n_lanes];
831 const unsigned int dofs_this_cell = dofs_per_cell[fe_index];
832
833 for (unsigned int i = row_starts_plain_indices[cell];
834 i < row_starts_plain_indices[cell] + dofs_this_cell;
835 ++i)
836 if (plain_dof_indices[i] >= part.local_size())
837 ghost_indices.push_back(
838 part.local_to_global(plain_dof_indices[i]));
839 }
840 std::sort(ghost_indices.begin(), ghost_indices.end());
841 ghost_indices.erase(std::unique(ghost_indices.begin(),
842 ghost_indices.end()),
843 ghost_indices.end());
844 IndexSet compressed_set(part.size());
845 compressed_set.add_indices(ghost_indices.begin(), ghost_indices.end());
846 compressed_set.subtract_set(part.locally_owned_range());
847 const bool all_ghosts_equal =
848 Utilities::MPI::min<int>(static_cast<int>(
849 compressed_set.n_elements() ==
850 part.ghost_indices().n_elements()),
851 part.get_mpi_communicator()) != 0;
852
853 std::shared_ptr<const Utilities::MPI::Partitioner> temp_0;
854
855 if (all_ghosts_equal)
856 temp_0 = vector_partitioner;
857 else
858 {
859 temp_0 = std::make_shared<Utilities::MPI::Partitioner>(
860 part.locally_owned_range(), part.get_mpi_communicator());
861 const_cast<Utilities::MPI::Partitioner *>(temp_0.get())
862 ->set_ghost_indices(compressed_set, part.ghost_indices());
863 }
864
865 if (use_vector_data_exchanger_full == false)
866 vector_exchanger_face_variants[0] = std::make_shared<
867 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
868 temp_0);
869 else
870 vector_exchanger_face_variants[0] =
871 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
872 temp_0, communicator_sm);
873 }
874
875 // construct a numbering of faces
876 std::vector<FaceToCellTopology<1>> all_faces(inner_faces);
877 all_faces.insert(all_faces.end(),
878 ghosted_faces.begin(),
879 ghosted_faces.end());
880 Table<2, unsigned int> cell_and_face_to_faces(
881 (row_starts.size() - 1) / start_components.back(),
882 2 * shape_info(0, 0).n_dimensions);
883 cell_and_face_to_faces.fill(numbers::invalid_unsigned_int);
884 for (unsigned int f = 0; f < all_faces.size(); ++f)
885 {
886 cell_and_face_to_faces(all_faces[f].cells_interior[0],
887 all_faces[f].interior_face_no) = f;
888 Assert(all_faces[f].cells_exterior[0] !=
889 numbers::invalid_unsigned_int,
890 ExcInternalError());
891 cell_and_face_to_faces(all_faces[f].cells_exterior[0],
892 all_faces[f].exterior_face_no) = f;
893 }
894
895 // lambda function to detect objects on face pairs
896 const auto loop_over_faces =
897 [&](const std::function<
898 void(const unsigned int, const unsigned int, const bool)> &fu) {
899 for (const auto &face : inner_faces)
900 {
901 AssertIndexRange(face.cells_interior[0], n_owned_cells);
902 fu(face.cells_exterior[0], face.exterior_face_no, false /*flag*/);
903 }
904 };
905
906 const auto loop_over_all_faces =
907 [&](const std::function<
908 void(const unsigned int, const unsigned int, const bool)> &fu) {
909 for (unsigned int c = 0; c < cell_and_face_to_faces.size(0); ++c)
910 for (unsigned int d = 0; d < cell_and_face_to_faces.size(1); ++d)
911 {
912 const unsigned int f = cell_and_face_to_faces(c, d);
913 if (f == numbers::invalid_unsigned_int)
914 continue;
915
916 const unsigned int cell_m = all_faces[f].cells_interior[0];
917 const unsigned int cell_p = all_faces[f].cells_exterior[0];
918
919 const bool ext = c == cell_m;
920
921 if (ext && cell_p == numbers::invalid_unsigned_int)
922 continue;
923
924 const unsigned int p = ext ? cell_p : cell_m;
925 const unsigned int face_no = ext ?
926 all_faces[f].exterior_face_no :
927 all_faces[f].interior_face_no;
928
929 fu(p, face_no, true);
930 }
931 };
932
933 const auto process_values =
934 [&](
935 std::shared_ptr<const Utilities::MPI::Partitioner>
936 &vector_partitioner_values,
937 const std::function<void(
938 const std::function<void(
939 const unsigned int, const unsigned int, const bool)> &)> &loop) {
940 bool all_nodal_and_tensorial = true;
941 for (unsigned int c = 0; c < n_base_elements; ++c)
942 {
943 const auto &si =
944 shape_info(global_base_element_offset + c, 0).data.front();
945 if (!si.nodal_at_cell_boundaries ||
946 (si.element_type ==
947 MatrixFreeFunctions::ElementType::tensor_none))
948 all_nodal_and_tensorial = false;
949 }
950 if (all_nodal_and_tensorial == false)
951 vector_partitioner_values = vector_partitioner;
952 else
953 {
954 bool has_noncontiguous_cell = false;
955
956 loop([&](const unsigned int cell_no,
957 const unsigned int face_no,
958 const bool flag) {
959 const unsigned int index =
960 dof_indices_contiguous[dof_access_cell][cell_no];
961 if (flag || (index != numbers::invalid_unsigned_int &&
962 index >= part.local_size()))
963 {
964 const unsigned int stride =
965 dof_indices_interleave_strides[dof_access_cell][cell_no];
966 unsigned int i = 0;
967 for (unsigned int e = 0; e < n_base_elements; ++e)
968 for (unsigned int c = 0; c < n_components[e]; ++c)
969 {
970 const ShapeInfo<double> &shape =
971 shape_info(global_base_element_offset + e, 0);
972 for (unsigned int j = 0;
973 j < shape.dofs_per_component_on_face;
974 ++j)
975 ghost_indices.push_back(part.local_to_global(
976 index + i +
977 shape.face_to_cell_index_nodal(face_no, j) *
978 stride));
979 i += shape.dofs_per_component_on_cell * stride;
980 }
981 AssertDimension(i, dofs_per_cell[0] * stride);
982 }
983 else if (index == numbers::invalid_unsigned_int)
984 has_noncontiguous_cell = true;
985 });
986 has_noncontiguous_cell =
987 Utilities::MPI::min<int>(static_cast<int>(
988 has_noncontiguous_cell),
989 part.get_mpi_communicator()) != 0;
990
991 std::sort(ghost_indices.begin(), ghost_indices.end());
992 ghost_indices.erase(std::unique(ghost_indices.begin(),
993 ghost_indices.end()),
994 ghost_indices.end());
995 IndexSet compressed_set(part.size());
996 compressed_set.add_indices(ghost_indices.begin(),
997 ghost_indices.end());
998 compressed_set.subtract_set(part.locally_owned_range());
999 const bool all_ghosts_equal =
1000 Utilities::MPI::min<int>(static_cast<int>(
1001 compressed_set.n_elements() ==
1002 part.ghost_indices().n_elements()),
1003 part.get_mpi_communicator()) != 0;
1004 if (all_ghosts_equal || has_noncontiguous_cell)
1005 vector_partitioner_values = vector_partitioner;
1006 else
1007 {
1008 vector_partitioner_values =
1009 std::make_shared<Utilities::MPI::Partitioner>(
1010 part.locally_owned_range(), part.get_mpi_communicator());
1011 const_cast<Utilities::MPI::Partitioner *>(
1012 vector_partitioner_values.get())
1013 ->set_ghost_indices(compressed_set, part.ghost_indices());
1014 }
1015 }
1016 };
1017
1018
1019 const auto process_gradients =
1020 [&](
1021 const std::shared_ptr<const Utilities::MPI::Partitioner>
1022 &vector_partitoner_values,
1023 std::shared_ptr<const Utilities::MPI::Partitioner>
1024 &vector_partitioner_gradients,
1025 const std::function<void(
1026 const std::function<void(
1027 const unsigned int, const unsigned int, const bool)> &)> &loop) {
1028 bool all_hermite = true;
1029 for (unsigned int c = 0; c < n_base_elements; ++c)
1030 if (shape_info(global_base_element_offset + c, 0).element_type !=
1031 MatrixFreeFunctions::tensor_symmetric_hermite)
1032 all_hermite = false;
1033 if (all_hermite == false ||
1034 vector_partitoner_values.get() == vector_partitioner.get())
1035 vector_partitioner_gradients = vector_partitioner;
1036 else
1037 {
1038 loop([&](const unsigned int cell_no,
1039 const unsigned int face_no,
1040 const bool flag) {
1041 const unsigned int index =
1042 dof_indices_contiguous[dof_access_cell][cell_no];
1043 if (flag || (index != numbers::invalid_unsigned_int &&
1044 index >= part.local_size()))
1045 {
1046 const unsigned int stride =
1047 dof_indices_interleave_strides[dof_access_cell][cell_no];
1048 unsigned int i = 0;
1049 for (unsigned int e = 0; e < n_base_elements; ++e)
1050 for (unsigned int c = 0; c < n_components[e]; ++c)
1051 {
1052 const ShapeInfo<double> &shape =
1053 shape_info(global_base_element_offset + e, 0);
1054 for (unsigned int j = 0;
1055 j < 2 * shape.dofs_per_component_on_face;
1056 ++j)
1057 ghost_indices.push_back(part.local_to_global(
1058 index + i +
1059 shape.face_to_cell_index_hermite(face_no, j) *
1060 stride));
1061 i += shape.dofs_per_component_on_cell * stride;
1062 }
1063 AssertDimension(i, dofs_per_cell[0] * stride);
1064 }
1065 });
1066 std::sort(ghost_indices.begin(), ghost_indices.end());
1067 ghost_indices.erase(std::unique(ghost_indices.begin(),
1068 ghost_indices.end()),
1069 ghost_indices.end());
1070 IndexSet compressed_set(part.size());
1071 compressed_set.add_indices(ghost_indices.begin(),
1072 ghost_indices.end());
1073 compressed_set.subtract_set(part.locally_owned_range());
1074 const bool all_ghosts_equal =
1075 Utilities::MPI::min<int>(static_cast<int>(
1076 compressed_set.n_elements() ==
1077 part.ghost_indices().n_elements()),
1078 part.get_mpi_communicator()) != 0;
1079 if (all_ghosts_equal)
1080 vector_partitioner_gradients = vector_partitioner;
1081 else
1082 {
1083 vector_partitioner_gradients =
1084 std::make_shared<Utilities::MPI::Partitioner>(
1085 part.locally_owned_range(), part.get_mpi_communicator());
1086 const_cast<Utilities::MPI::Partitioner *>(
1087 vector_partitioner_gradients.get())
1088 ->set_ghost_indices(compressed_set, part.ghost_indices());
1089 }
1090 }
1091 };
1092
1093 std::shared_ptr<const Utilities::MPI::Partitioner> temp_1, temp_2, temp_3,
1094 temp_4;
1095
1096 // partitioner 1: values on faces
1097 process_values(temp_1, loop_over_faces);
1098
1099 // partitioner 2: values and gradients on faces
1100 process_gradients(temp_1, temp_2, loop_over_faces);
1101
1102 if (fill_cell_centric)
1103 {
1104 ghost_indices.clear();
1105 // partitioner 3: values on all faces
1106 process_values(temp_3, loop_over_all_faces);
1107 // partitioner 4: values and gradients on faces
1108 process_gradients(temp_3, temp_4, loop_over_all_faces);
1109 }
1110 else
1111 {
1112 temp_3 = std::make_shared<Utilities::MPI::Partitioner>(
1113 part.locally_owned_range(), part.get_mpi_communicator());
1114 temp_4 = std::make_shared<Utilities::MPI::Partitioner>(
1115 part.locally_owned_range(), part.get_mpi_communicator());
1116 }
1117
1118 if (use_vector_data_exchanger_full == false)
1119 {
1120 vector_exchanger_face_variants[1] = std::make_shared<
1121 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1122 temp_1);
1123 vector_exchanger_face_variants[2] = std::make_shared<
1124 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1125 temp_2);
1126 vector_exchanger_face_variants[3] = std::make_shared<
1127 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1128 temp_3);
1129 vector_exchanger_face_variants[4] = std::make_shared<
1130 MatrixFreeFunctions::VectorDataExchange::PartitionerWrapper>(
1131 temp_4);
1132 }
1133 else
1134 {
1135 vector_exchanger_face_variants[1] =
1136 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1137 temp_1, communicator_sm);
1138 vector_exchanger_face_variants[2] =
1139 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1140 temp_2, communicator_sm);
1141 vector_exchanger_face_variants[3] =
1142 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1143 temp_3, communicator_sm);
1144 vector_exchanger_face_variants[4] =
1145 std::make_shared<MatrixFreeFunctions::VectorDataExchange::Full>(
1146 temp_4, communicator_sm);
1147 }
1148 }
1149
1150
1151
1152 void
1153 DoFInfo::compute_shared_memory_contiguous_indices(
1154 std::array<std::vector<std::pair<unsigned int, unsigned int>>, 3>
1155 &cell_indices_contiguous_sm)
1156 {
1157 AssertDimension(dofs_per_cell.size(), 1);
1158
1159 for (unsigned int i = 0; i < 3; ++i)
1160 {
1161 dof_indices_contiguous_sm[i].resize(
1162 cell_indices_contiguous_sm[i].size());
1163
1164 for (unsigned int j = 0; j < cell_indices_contiguous_sm[i].size();
1165 ++j)
1166 if (cell_indices_contiguous_sm[i][j].first !=
1167 numbers::invalid_unsigned_int)
1168 dof_indices_contiguous_sm[i][j] = {
1169 cell_indices_contiguous_sm[i][j].first,
1170 cell_indices_contiguous_sm[i][j].second * dofs_per_cell[0]};
1171 else
1172 dof_indices_contiguous_sm[i][j] = {numbers::invalid_unsigned_int,
1173 numbers::invalid_unsigned_int};
1174 }
1175 }
1176
1177
1178
1179 namespace internal
1180 {
1181 // We construct the connectivity graph in parallel. we use one lock for
1182 // 256 degrees of freedom to keep the number of locks down to a
1183 // reasonable level and reduce the cost of locking to some extent.
1184 static constexpr unsigned int bucket_size_threading = 256;
1185
1186
1187
1188 void
1189 compute_row_lengths(const unsigned int begin,
1190 const unsigned int end,
1191 const DoFInfo & dof_info,
1192 std::vector<std::mutex> & mutexes,
1193 std::vector<unsigned int> &row_lengths)
1194 {
1195 std::vector<unsigned int> scratch;
1196 const unsigned int n_components = dof_info.start_components.back();
1197 for (unsigned int block = begin; block < end; ++block)
1198 {
1199 scratch.clear();
1200 scratch.insert(
1201 scratch.end(),
1202 dof_info.dof_indices.data() +
1203 dof_info.row_starts[block * n_components].first,
1204 dof_info.dof_indices.data() +
1205 dof_info.row_starts[(block + 1) * n_components].first);
1206 std::sort(scratch.begin(), scratch.end());
1207 std::vector<unsigned int>::const_iterator end_unique =
1208 std::unique(scratch.begin(), scratch.end());
1209 std::vector<unsigned int>::const_iterator it = scratch.begin();
1210 while (it != end_unique)
1211 {
1212 // In this code, the procedure is that we insert all elements
1213 // that are within the range of one lock at once
1214 const unsigned int next_bucket =
1215 (*it / bucket_size_threading + 1) * bucket_size_threading;
1216 std::lock_guard<std::mutex> lock(
1217 mutexes[*it / bucket_size_threading]);
1218 for (; it != end_unique && *it < next_bucket; ++it)
1219 {
1220 AssertIndexRange(*it, row_lengths.size());
1221 row_lengths[*it]++;
1222 }
1223 }
1224 }
1225 }
1226
1227 void
1228 fill_connectivity_dofs(const unsigned int begin,
1229 const unsigned int end,
1230 const DoFInfo & dof_info,
1231 const std::vector<unsigned int> &row_lengths,
1232 std::vector<std::mutex> & mutexes,
1233 ::SparsityPattern & connectivity_dof)
1234 {
1235 std::vector<unsigned int> scratch;
1236 const unsigned int n_components = dof_info.start_components.back();
1237 for (unsigned int block = begin; block < end; ++block)
1238 {
1239 scratch.clear();
1240 scratch.insert(
1241 scratch.end(),
1242 dof_info.dof_indices.data() +
1243 dof_info.row_starts[block * n_components].first,
1244 dof_info.dof_indices.data() +
1245 dof_info.row_starts[(block + 1) * n_components].first);
1246 std::sort(scratch.begin(), scratch.end());
1247 std::vector<unsigned int>::const_iterator end_unique =
1248 std::unique(scratch.begin(), scratch.end());
1249 std::vector<unsigned int>::const_iterator it = scratch.begin();
1250 while (it != end_unique)
1251 {
1252 const unsigned int next_bucket =
1253 (*it / bucket_size_threading + 1) * bucket_size_threading;
1254 std::lock_guard<std::mutex> lock(
1255 mutexes[*it / bucket_size_threading]);
1256 for (; it != end_unique && *it < next_bucket; ++it)
1257 if (row_lengths[*it] > 0)
1258 connectivity_dof.add(*it, block);
1259 }
1260 }
1261 }
1262
1263
1264
1265 void
1266 fill_connectivity(const unsigned int begin,
1267 const unsigned int end,
1268 const DoFInfo & dof_info,
1269 const std::vector<unsigned int> &renumbering,
1270 const ::SparsityPattern & connectivity_dof,
1271 DynamicSparsityPattern & connectivity)
1272 {
1273 ordered_vector row_entries;
1274 const unsigned int n_components = dof_info.start_components.back();
1275 for (unsigned int block = begin; block < end; ++block)
1276 {
1277 row_entries.clear();
1278
1279 const unsigned int
1280 *it = dof_info.dof_indices.data() +
1281 dof_info.row_starts[block * n_components].first,
1282 *end_cell = dof_info.dof_indices.data() +
1283 dof_info.row_starts[(block + 1) * n_components].first;
1284 for (; it != end_cell; ++it)
1285 {
1286 SparsityPattern::iterator sp = connectivity_dof.begin(*it);
1287 std::vector<types::global_dof_index>::iterator insert_pos =
1288 row_entries.begin();
1289 for (; sp != connectivity_dof.end(*it); ++sp)
1290 if (sp->column() != block)
1291 row_entries.insert(renumbering[sp->column()], insert_pos);
1292 }
1293 connectivity.add_entries(renumbering[block],
1294 row_entries.begin(),
1295 row_entries.end());
1296 }
1297 }
1298
1299 } // namespace internal
1300
1301 void
1302 DoFInfo::make_connectivity_graph(
1303 const TaskInfo & task_info,
1304 const std::vector<unsigned int> &renumbering,
1305 DynamicSparsityPattern & connectivity) const
1306 {
1307 unsigned int n_rows = (vector_partitioner->local_range().second -
1308 vector_partitioner->local_range().first) +
1309 vector_partitioner->ghost_indices().n_elements();
1310
1311 // Avoid square sparsity patterns that allocate the diagonal entry
1312 if (n_rows == task_info.n_active_cells)
1313 ++n_rows;
1314
1315 // first determine row lengths
1316 std::vector<unsigned int> row_lengths(n_rows);
1317 std::vector<std::mutex> mutexes(n_rows / internal::bucket_size_threading +
1318 1);
1319 ::parallel::apply_to_subranges(
1320 0,
1321 task_info.n_active_cells,
1322 [this, &mutexes, &row_lengths](const unsigned int begin,
1323 const unsigned int end) {
1324 internal::compute_row_lengths(
1325 begin, end, *this, mutexes, row_lengths);
1326 },
1327 20);
1328
1329 // disregard dofs that only sit on a single cell because they cannot
1330 // couple
1331 for (unsigned int row = 0; row < n_rows; ++row)
1332 if (row_lengths[row] <= 1)
1333 row_lengths[row] = 0;
1334
1335 // Create a temporary sparsity pattern that holds to each degree of
1336 // freedom on which cells it appears, i.e., store the connectivity
1337 // between cells and dofs
1338 SparsityPattern connectivity_dof(n_rows,
1339 task_info.n_active_cells,
1340 row_lengths);
1341 ::parallel::apply_to_subranges(
1342 0,
1343 task_info.n_active_cells,
1344 [this, &row_lengths, &mutexes, &connectivity_dof](
1345 const unsigned int begin, const unsigned int end) {
1346 internal::fill_connectivity_dofs(
1347 begin, end, *this, row_lengths, mutexes, connectivity_dof);
1348 },
1349 20);
1350 connectivity_dof.compress();
1351
1352
1353 // Invert renumbering for use in fill_connectivity.
1354 std::vector<unsigned int> reverse_numbering(task_info.n_active_cells);
1355 reverse_numbering = Utilities::invert_permutation(renumbering);
1356
1357 // From the above connectivity between dofs and cells, we can finally
1358 // create a connectivity list between cells. The connectivity graph
1359 // should apply the renumbering, i.e., the entry for cell j is the entry
1360 // for cell renumbering[j] in the original ordering.
1361 ::parallel::apply_to_subranges(
1362 0,
1363 task_info.n_active_cells,
1364 [this, &reverse_numbering, &connectivity_dof, &connectivity](
1365 const unsigned int begin, const unsigned int end) {
1366 internal::fill_connectivity(begin,
1367 end,
1368 *this,
1369 reverse_numbering,
1370 connectivity_dof,
1371 connectivity);
1372 },
1373 20);
1374 }
1375
1376
1377
1378 void
1379 DoFInfo::compute_dof_renumbering(
1380 std::vector<types::global_dof_index> &renumbering)
1381 {
1382 const unsigned int locally_owned_size =
1383 vector_partitioner->locally_owned_size();
1384 renumbering.resize(0);
1385 renumbering.resize(locally_owned_size, numbers::invalid_dof_index);
1386
1387 types::global_dof_index counter = 0;
1388 const unsigned int n_components = start_components.back();
1389 const unsigned int n_cell_batches =
1390 n_vectorization_lanes_filled[dof_access_cell].size();
1391 Assert(n_cell_batches <=
1392 (row_starts.size() - 1) / vectorization_length / n_components,
1393 ExcInternalError());
1394 for (unsigned int cell_no = 0; cell_no < n_cell_batches; ++cell_no)
1395 {
1396 // do not renumber in case we have constraints
1397 if (row_starts[cell_no * n_components * vectorization_length]
1398 .second ==
1399 row_starts[(cell_no + 1) * n_components * vectorization_length]
1400 .second)
1401 {
1402 const unsigned int ndofs =
1403 dofs_per_cell.size() == 1 ?
1404 dofs_per_cell[0] :
1405 (dofs_per_cell[cell_active_fe_index.size() > 0 ?
1406 cell_active_fe_index[cell_no] :
1407 0]);
1408 const unsigned int *dof_ind =
1409 dof_indices.data() +
1410 row_starts[cell_no * n_components * vectorization_length].first;
1411 for (unsigned int i = 0; i < ndofs; ++i)
1412 for (unsigned int j = 0;
1413 j < n_vectorization_lanes_filled[dof_access_cell][cell_no];
1414 ++j)
1415 if (dof_ind[j * ndofs + i] < locally_owned_size)
1416 if (renumbering[dof_ind[j * ndofs + i]] ==
1417 numbers::invalid_dof_index)
1418 renumbering[dof_ind[j * ndofs + i]] = counter++;
1419 }
1420 }
1421
1422 AssertIndexRange(counter, locally_owned_size + 1);
1423 for (types::global_dof_index &dof_index : renumbering)
1424 if (dof_index == numbers::invalid_dof_index)
1425 dof_index = counter++;
1426
1427 // transform indices to global index space
1428 for (types::global_dof_index &dof_index : renumbering)
1429 dof_index = vector_partitioner->local_to_global(dof_index);
1430
1431 AssertDimension(counter, renumbering.size());
1432 }
1433
1434
1435
1436 std::size_t
1437 DoFInfo::memory_consumption() const
1438 {
1439 std::size_t memory = sizeof(*this);
1440 memory +=
1441 (row_starts.capacity() * sizeof(std::pair<unsigned int, unsigned int>));
1442 memory += MemoryConsumption::memory_consumption(dof_indices);
1443 memory += MemoryConsumption::memory_consumption(row_starts_plain_indices);
1444 memory += MemoryConsumption::memory_consumption(plain_dof_indices);
1445 memory += MemoryConsumption::memory_consumption(constraint_indicator);
1446 memory += MemoryConsumption::memory_consumption(*vector_partitioner);
1447 return memory;
1448 }
1449 } // namespace MatrixFreeFunctions
1450} // namespace internal
1451
1452namespace internal
1453{
1454 namespace MatrixFreeFunctions
1455 {
1456 template struct ConstraintValues<double>;
1457 template struct ConstraintValues<float>;
1458
1459 template void
1460 DoFInfo::read_dof_indices<double>(
1461 const std::vector<types::global_dof_index> &,
1462 const std::vector<unsigned int> &,
1463 const ::AffineConstraints<double> &,
1464 const unsigned int,
1465 ConstraintValues<double> &,
1466 bool &);
1467
1468 template void
1469 DoFInfo::read_dof_indices<float>(
1470 const std::vector<types::global_dof_index> &,
1471 const std::vector<unsigned int> &,
1472 const ::AffineConstraints<float> &,
1473 const unsigned int,
1474 ConstraintValues<double> &,
1475 bool &);
1476
1477 template void
1478 DoFInfo::compute_face_index_compression<1>(
1479 const std::vector<FaceToCellTopology<1>> &);
1480 template void
1481 DoFInfo::compute_face_index_compression<2>(
1482 const std::vector<FaceToCellTopology<2>> &);
1483 template void
1484 DoFInfo::compute_face_index_compression<4>(
1485 const std::vector<FaceToCellTopology<4>> &);
1486 template void
1487 DoFInfo::compute_face_index_compression<8>(
1488 const std::vector<FaceToCellTopology<8>> &);
1489 template void
1490 DoFInfo::compute_face_index_compression<16>(
1491 const std::vector<FaceToCellTopology<16>> &);
1492
1493 template void
1494 DoFInfo::compute_vector_zero_access_pattern<1>(
1495 const TaskInfo &,
1496 const std::vector<FaceToCellTopology<1>> &);
1497 template void
1498 DoFInfo::compute_vector_zero_access_pattern<2>(
1499 const TaskInfo &,
1500 const std::vector<FaceToCellTopology<2>> &);
1501 template void
1502 DoFInfo::compute_vector_zero_access_pattern<4>(
1503 const TaskInfo &,
1504 const std::vector<FaceToCellTopology<4>> &);
1505 template void
1506 DoFInfo::compute_vector_zero_access_pattern<8>(
1507 const TaskInfo &,
1508 const std::vector<FaceToCellTopology<8>> &);
1509 template void
1510 DoFInfo::compute_vector_zero_access_pattern<16>(
1511 const TaskInfo &,
1512 const std::vector<FaceToCellTopology<16>> &);
1513
1514 template void
1515 DoFInfo::print_memory_consumption<std::ostream>(std::ostream &,
1516 const TaskInfo &) const;
1517 template void
1518 DoFInfo::print_memory_consumption<ConditionalOStream>(
1519 ConditionalOStream &,
1520 const TaskInfo &) const;
1521
1522 template void
1523 DoFInfo::print<double>(const std::vector<double> &,
1524 const std::vector<unsigned int> &,
1525 std::ostream &) const;
1526
1527 template void
1528 DoFInfo::print<float>(const std::vector<float> &,
1529 const std::vector<unsigned int> &,
1530 std::ostream &) const;
1531 } // namespace MatrixFreeFunctions
1532} // namespace internal
1533
1534DEAL_II_NAMESPACE_CLOSE
IndexSet::index_within_set
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
IndexSet::n_elements
size_type n_elements() const
Definition: index_set.h:1832
IndexSet::subtract_set
void subtract_set(const IndexSet &other)
Definition: index_set.cc:258
IndexSet::add_range
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1673
IndexSet::compress
void compress() const
Definition: index_set.h:1642
IndexSet::add_indices
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1703
MPI_Comm
Utilities::MPI::Partitioner::locally_owned_range
const IndexSet & locally_owned_range() const
Utilities::MPI::Partitioner::local_size
unsigned int local_size() const
Utilities::MPI::Partitioner::ghost_indices
const IndexSet & ghost_indices() const
Utilities::MPI::Partitioner::local_to_global
types::global_dof_index local_to_global(const unsigned int local_index) const
Utilities::MPI::Partitioner::set_ghost_indices
void set_ghost_indices(const IndexSet &ghost_indices, const IndexSet &larger_ghost_index_set=IndexSet())
Utilities::MPI::Partitioner::get_mpi_communicator
virtual const MPI_Comm & get_mpi_communicator() const override
Utilities::MPI::Partitioner::size
types::global_dof_index size() const
conditional_ostream.h
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
second
Point< 2 > second
Definition: grid_out.cc:4588
first
Point< 2 > first
Definition: grid_out.cc:4587
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
SparsityPattern::compress
void compress()
Definition: sparsity_pattern.cc:338
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:439
DynamicSparsityPattern::add_entries
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
Definition: dynamic_sparsity_pattern.h:1046
SparsityPatternBase::add
void add(const size_type i, const size_type j)
Definition: sparsity_pattern.cc:704
internal::MatrixFreeFunctions::tensor_none
@ tensor_none
Definition: shape_info.h:99
internal::MatrixFreeFunctions::tensor_symmetric_hermite
@ tensor_symmetric_hermite
Definition: shape_info.h:69
GridTools::shift
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2022
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:266
MemorySpace::swap
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
Definition: memory_space_data.h:92
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:67
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:53
Utilities::invert_permutation
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1482
internal::MatrixFreeFunctions::internal::fill_connectivity_dofs
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:1228
internal::MatrixFreeFunctions::internal::fill_connectivity
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:1266
internal::MatrixFreeFunctions::internal::compute_row_lengths
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:1189
internal::MatrixFreeFunctions::internal::bucket_size_threading
static constexpr unsigned int bucket_size_threading
Definition: dof_info.cc:1184
internal::TriangulationImplementation::n_active_cells
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12594
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:68
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:196
numbers::invalid_dof_index
const types::global_dof_index invalid_dof_index
Definition: types.h:211
parallel::internal::apply_to_subranges
void apply_to_subranges(const tbb::blocked_range< RangeType > &range, const Function &f)
Definition: parallel.h:367
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5484
internal::MatrixFreeFunctions::DoFInfo::reorder_cells
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:271
internal::MatrixFreeFunctions::DoFInfo::constraint_indicator
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
Definition: dof_info.h:510
internal::MatrixFreeFunctions::DoFInfo::row_starts
std::vector< std::pair< unsigned int, unsigned int > > row_starts
Definition: dof_info.h:481
internal::MatrixFreeFunctions::DoFInfo::compute_cell_index_compression
void compute_cell_index_compression(const std::vector< unsigned char > &irregular_cells)
Definition: dof_info.cc:473
internal::MatrixFreeFunctions::DoFInfo::dofs_per_cell
std::vector< unsigned int > dofs_per_cell
Definition: dof_info.h:667
internal::MatrixFreeFunctions::DoFInfo::vector_exchanger
std::shared_ptr< const internal::MatrixFreeFunctions::VectorDataExchange::Base > vector_exchanger
Definition: dof_info.h:572
internal::MatrixFreeFunctions::DoFInfo::fe_index_conversion
std::vector< std::vector< unsigned int > > fe_index_conversion
Definition: dof_info.h:695
internal::MatrixFreeFunctions::DoFInfo::compute_tight_partitioners
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:803
internal::MatrixFreeFunctions::DoFInfo::dof_indices
std::vector< unsigned int > dof_indices
Definition: dof_info.h:498
internal::MatrixFreeFunctions::DoFInfo::vector_partitioner
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:565
internal::MatrixFreeFunctions::DoFInfo::dof_indices_interleave_strides
std::array< std::vector< unsigned int >, 3 > dof_indices_interleave_strides
Definition: dof_info.h:546
internal::MatrixFreeFunctions::DoFInfo::row_starts_plain_indices
std::vector< unsigned int > row_starts_plain_indices
Definition: dof_info.h:609
internal::MatrixFreeFunctions::DoFInfo::dofs_per_face
std::vector< unsigned int > dofs_per_face
Definition: dof_info.h:672
internal::MatrixFreeFunctions::DoFInfo::n_components
std::vector< unsigned int > n_components
Definition: dof_info.h:637
internal::MatrixFreeFunctions::DoFInfo::assign_ghosts
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:141
internal::MatrixFreeFunctions::DoFInfo::ghost_dofs
std::vector< types::global_dof_index > ghost_dofs
Definition: dof_info.h:702
internal::MatrixFreeFunctions::DoFInfo::dof_indices_contiguous
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
Definition: dof_info.h:525
internal::MatrixFreeFunctions::DoFInfo::get_dof_indices_on_cell_batch
void get_dof_indices_on_cell_batch(std::vector< unsigned int > &locall_indices, const unsigned int cell, const bool with_constraints=true) const
Definition: dof_info.cc:73
internal::MatrixFreeFunctions::DoFInfo::cell_active_fe_index
std::vector< unsigned int > cell_active_fe_index
Definition: dof_info.h:682
internal::MatrixFreeFunctions::DoFInfo::vector_exchanger_face_variants
std::array< std::shared_ptr< const internal::MatrixFreeFunctions::VectorDataExchange::Base >, 5 > vector_exchanger_face_variants
Definition: dof_info.h:597
internal::MatrixFreeFunctions::DoFInfo::plain_dof_indices
std::vector< unsigned int > plain_dof_indices
Definition: dof_info.h:619
internal::MatrixFreeFunctions::DoFInfo::n_vectorization_lanes_filled
std::array< std::vector< unsigned char >, 3 > n_vectorization_lanes_filled
Definition: dof_info.h:557
internal::MatrixFreeFunctions::DoFInfo::start_components
std::vector< unsigned int > start_components
Definition: dof_info.h:643
internal::MatrixFreeFunctions::DoFInfo::dof_indices_interleaved
std::vector< unsigned int > dof_indices_interleaved
Definition: dof_info.h:515
internal::MatrixFreeFunctions::DoFInfo::index_storage_variants
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
Definition: dof_info.h:473
internal::MatrixFreeFunctions::TaskInfo::cell_partition_data
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:474
Generated by doxygen 1.9.4