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