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