Reference documentation for deal.II version GIT 7bd1c95dbb 2022-06-29 20:50:01+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\}}\)
tools.h
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 #ifndef dealii_matrix_free_tools_h
17 #define dealii_matrix_free_tools_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/grid/tria.h>
22 
26 
27 
29 
34 namespace MatrixFreeTools
35 {
41  template <int dim, typename AdditionalData>
42  void
44  AdditionalData & additional_data);
45 
54  template <int dim,
55  int fe_degree,
56  int n_q_points_1d,
57  int n_components,
58  typename Number,
59  typename VectorizedArrayType,
60  typename VectorType>
61  void
64  VectorType & diagonal_global,
65  const std::function<void(FEEvaluation<dim,
66  fe_degree,
67  n_q_points_1d,
68  n_components,
69  Number,
70  VectorizedArrayType> &)> &local_vmult,
71  const unsigned int dof_no = 0,
72  const unsigned int quad_no = 0,
73  const unsigned int first_selected_component = 0);
74 
78  template <typename CLASS,
79  int dim,
80  int fe_degree,
81  int n_q_points_1d,
82  int n_components,
83  typename Number,
84  typename VectorizedArrayType,
85  typename VectorType>
86  void
89  VectorType & diagonal_global,
90  void (CLASS::*cell_operation)(FEEvaluation<dim,
91  fe_degree,
92  n_q_points_1d,
93  n_components,
94  Number,
95  VectorizedArrayType> &) const,
96  const CLASS * owning_class,
97  const unsigned int dof_no = 0,
98  const unsigned int quad_no = 0,
99  const unsigned int first_selected_component = 0);
100 
101 
110  template <int dim,
111  int fe_degree,
112  int n_q_points_1d,
113  int n_components,
114  typename Number,
115  typename VectorizedArrayType,
116  typename MatrixType>
117  void
120  const AffineConstraints<Number> & constraints,
121  MatrixType & matrix,
122  const std::function<void(FEEvaluation<dim,
123  fe_degree,
124  n_q_points_1d,
125  n_components,
126  Number,
127  VectorizedArrayType> &)> &local_vmult,
128  const unsigned int dof_no = 0,
129  const unsigned int quad_no = 0,
130  const unsigned int first_selected_component = 0);
131 
135  template <typename CLASS,
136  int dim,
137  int fe_degree,
138  int n_q_points_1d,
139  int n_components,
140  typename Number,
141  typename VectorizedArrayType,
142  typename MatrixType>
143  void
146  const AffineConstraints<Number> & constraints,
147  MatrixType & matrix,
148  void (CLASS::*cell_operation)(FEEvaluation<dim,
149  fe_degree,
150  n_q_points_1d,
151  n_components,
152  Number,
153  VectorizedArrayType> &) const,
154  const CLASS * owning_class,
155  const unsigned int dof_no = 0,
156  const unsigned int quad_no = 0,
157  const unsigned int first_selected_component = 0);
158 
159 
160  // implementations
161 
162 #ifndef DOXYGEN
163 
164  template <int dim, typename AdditionalData>
165  void
167  AdditionalData & additional_data)
168  {
169  // ... determine if we are on an active or a multigrid level
170  const unsigned int level = additional_data.mg_level;
171  const bool is_mg = (level != numbers::invalid_unsigned_int);
172 
173  // ... create empty list for the category of each cell
174  if (is_mg)
175  additional_data.cell_vectorization_category.assign(
176  std::distance(tria.begin(level), tria.end(level)), 0);
177  else
178  additional_data.cell_vectorization_category.assign(tria.n_active_cells(),
179  0);
180 
181  // ... set up scaling factor
182  std::vector<unsigned int> factors(GeometryInfo<dim>::faces_per_cell);
183 
184  auto bids = tria.get_boundary_ids();
185  std::sort(bids.begin(), bids.end());
186 
187  {
188  unsigned int n_bids = bids.size() + 1;
189  int offset = 1;
190  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
191  i++, offset = offset * n_bids)
192  factors[i] = offset;
193  }
194 
195  const auto to_category = [&](const auto &cell) {
196  unsigned int c_num = 0;
197  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
198  {
199  auto &face = *cell->face(i);
200  if (face.at_boundary() && !cell->has_periodic_neighbor(i))
201  c_num +=
202  factors[i] * (1 + std::distance(bids.begin(),
203  std::find(bids.begin(),
204  bids.end(),
205  face.boundary_id())));
206  }
207  return c_num;
208  };
209 
210  if (!is_mg)
211  {
212  for (auto cell = tria.begin_active(); cell != tria.end(); ++cell)
213  {
214  if (cell->is_locally_owned())
215  additional_data
216  .cell_vectorization_category[cell->active_cell_index()] =
217  to_category(cell);
218  }
219  }
220  else
221  {
222  for (auto cell = tria.begin(level); cell != tria.end(level); ++cell)
223  {
224  if (cell->is_locally_owned_on_level())
225  additional_data.cell_vectorization_category[cell->index()] =
226  to_category(cell);
227  }
228  }
229 
230  // ... finalize set up of matrix_free
231  additional_data.hold_all_faces_to_owned_cells = true;
232  additional_data.cell_vectorization_categories_strict = true;
233  additional_data.mapping_update_flags_faces_by_cells =
234  additional_data.mapping_update_flags_inner_faces |
235  additional_data.mapping_update_flags_boundary_faces;
236  }
237 
238  namespace internal
239  {
240  template <typename Number>
241  struct LocalCSR
242  {
243  LocalCSR()
244  : row{0}
245  {}
246 
247  std::vector<unsigned int> row_lid_to_gid;
248  std::vector<unsigned int> row;
249  std::vector<unsigned int> col;
250  std::vector<Number> val;
251  };
252 
253  template <int dim,
254  int fe_degree,
255  int n_q_points_1d,
256  int n_components,
257  typename Number,
258  typename VectorizedArrayType>
259  class ComputeDiagonalHelper
260  {
261  public:
262  static const unsigned int n_lanes = VectorizedArrayType::size();
263 
264  ComputeDiagonalHelper(FEEvaluation<dim,
265  fe_degree,
266  n_q_points_1d,
267  n_components,
268  Number,
269  VectorizedArrayType> &phi)
270  : phi(phi)
271  {}
272 
273  void
274  reinit(const unsigned int cell)
275  {
276  this->phi.reinit(cell);
277  // STEP 1: get relevant information from FEEvaluation
278  const auto & dof_info = phi.get_dof_info();
279  const unsigned int n_fe_components = dof_info.start_components.back();
280  const unsigned int dofs_per_component = phi.dofs_per_component;
281  const auto & matrix_free = phi.get_matrix_free();
282 
283  // if we have a block vector with components with the same DoFHandler,
284  // each component is described with same set of constraints and
285  // we consider the shift in components only during access of the global
286  // vector
287  const unsigned int first_selected_component =
288  n_fe_components == 1 ? 0 : phi.get_first_selected_component();
289 
290  const unsigned int n_lanes_filled =
291  matrix_free.n_active_entries_per_cell_batch(cell);
292 
293  std::array<const unsigned int *, n_lanes> dof_indices{};
294  {
295  for (unsigned int v = 0; v < n_lanes_filled; ++v)
296  dof_indices[v] =
297  dof_info.dof_indices.data() +
298  dof_info
299  .row_starts[(cell * n_lanes + v) * n_fe_components +
300  first_selected_component]
301  .first;
302  }
303 
304  // STEP 2: setup CSR storage of transposed locally-relevant
305  // constraint matrix
306  c_pools = std::array<internal::LocalCSR<Number>, n_lanes>();
307 
308  for (unsigned int v = 0; v < n_lanes_filled; ++v)
309  {
310  unsigned int index_indicators, next_index_indicators;
311 
312  index_indicators =
313  dof_info
314  .row_starts[(cell * n_lanes + v) * n_fe_components +
315  first_selected_component]
316  .second;
317  next_index_indicators =
318  dof_info
319  .row_starts[(cell * n_lanes + v) * n_fe_components +
320  first_selected_component + 1]
321  .second;
322 
323  // STEP 2a: setup locally-relevant constraint matrix in a
324  // coordinate list (COO)
325  std::vector<std::tuple<unsigned int, unsigned int, Number>>
326  locally_relevant_constrains; // (constrained local index,
327  // global index of dof which
328  // constrains, weight)
329 
330  if (n_components == 1 || n_fe_components == 1)
331  {
332  unsigned int ind_local = 0;
333  for (; index_indicators != next_index_indicators;
334  ++index_indicators, ++ind_local)
335  {
336  const std::pair<unsigned short, unsigned short> indicator =
337  dof_info.constraint_indicator[index_indicators];
338 
339  for (unsigned int j = 0; j < indicator.first;
340  ++j, ++ind_local)
341  locally_relevant_constrains.emplace_back(
342  ind_local, dof_indices[v][j], 1.0);
343 
344  dof_indices[v] += indicator.first;
345 
346  const Number *data_val =
347  matrix_free.constraint_pool_begin(indicator.second);
348  const Number *end_pool =
349  matrix_free.constraint_pool_end(indicator.second);
350 
351  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
352  locally_relevant_constrains.emplace_back(ind_local,
353  *dof_indices[v],
354  *data_val);
355  }
356 
357  AssertIndexRange(ind_local, dofs_per_component + 1);
358 
359  for (; ind_local < dofs_per_component;
360  ++dof_indices[v], ++ind_local)
361  locally_relevant_constrains.emplace_back(ind_local,
362  *dof_indices[v],
363  1.0);
364  }
365  else
366  {
367  // case with vector-valued finite elements where all
368  // components are included in one single vector. Assumption:
369  // first come all entries to the first component, then all
370  // entries to the second one, and so on. This is ensured by
371  // the way MatrixFree reads out the indices.
372  for (unsigned int comp = 0; comp < n_components; ++comp)
373  {
374  unsigned int ind_local = 0;
375 
376  // check whether there is any constraint on the current
377  // cell
378  for (; index_indicators != next_index_indicators;
379  ++index_indicators, ++ind_local)
380  {
381  const std::pair<unsigned short, unsigned short>
382  indicator =
383  dof_info.constraint_indicator[index_indicators];
384 
385  // run through values up to next constraint
386  for (unsigned int j = 0; j < indicator.first;
387  ++j, ++ind_local)
388  locally_relevant_constrains.emplace_back(
389  comp * dofs_per_component + ind_local,
390  dof_indices[v][j],
391  1.0);
392  dof_indices[v] += indicator.first;
393 
394  const Number *data_val =
395  matrix_free.constraint_pool_begin(indicator.second);
396  const Number *end_pool =
397  matrix_free.constraint_pool_end(indicator.second);
398 
399  for (; data_val != end_pool;
400  ++data_val, ++dof_indices[v])
401  locally_relevant_constrains.emplace_back(
402  comp * dofs_per_component + ind_local,
403  *dof_indices[v],
404  *data_val);
405  }
406 
407  AssertIndexRange(ind_local, dofs_per_component + 1);
408 
409  // get the dof values past the last constraint
410  for (; ind_local < dofs_per_component;
411  ++dof_indices[v], ++ind_local)
412  locally_relevant_constrains.emplace_back(
413  comp * dofs_per_component + ind_local,
414  *dof_indices[v],
415  1.0);
416 
417  if (comp + 1 < n_components)
418  {
419  next_index_indicators =
420  dof_info
421  .row_starts[(cell * n_lanes + v) * n_fe_components +
422  first_selected_component + comp + 2]
423  .second;
424  }
425  }
426  }
427  // STEP 2b: sort and make unique
428 
429  // sort vector
430  std::sort(locally_relevant_constrains.begin(),
431  locally_relevant_constrains.end(),
432  [](const auto &a, const auto &b) {
433  if (std::get<0>(a) < std::get<0>(b))
434  return true;
435  return (std::get<0>(a) == std::get<0>(b)) &&
436  (std::get<1>(a) < std::get<1>(b));
437  });
438 
439  // make sure that all entries are unique
440  locally_relevant_constrains.erase(
441  unique(locally_relevant_constrains.begin(),
442  locally_relevant_constrains.end(),
443  [](const auto &a, const auto &b) {
444  return (std::get<1>(a) == std::get<1>(b)) &&
445  (std::get<0>(a) == std::get<0>(b));
446  }),
447  locally_relevant_constrains.end());
448 
449  // STEP 2c: apply hanging-node constraints
450  if (dof_info.hanging_node_constraint_masks.size() > 0 &&
451  dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
452  dof_info
453  .hanging_node_constraint_masks_comp[phi.get_active_fe_index()]
454  [first_selected_component])
455  {
456  const auto mask =
457  dof_info.hanging_node_constraint_masks[cell * n_lanes + v];
458 
459  // cell has hanging nodes
460  if (mask != ::internal::MatrixFreeFunctions::
462  {
463  // check if hanging node internpolation matrix has been set
464  // up
465  if (locally_relevant_constrains_hn_map.find(mask) ==
466  locally_relevant_constrains_hn_map.end())
467  {
468  // 1) collect hanging-node constraints for cell assuming
469  // scalar finite element
471  dofs_per_component);
472 
475  VectorizedArrayType::size()>
476  constraint_mask;
477  constraint_mask.fill(
480 
481  constraint_mask[0] = mask;
482 
483  std::vector<
484  std::tuple<unsigned int, unsigned int, Number>>
485  locally_relevant_constrains_hn;
486 
487  for (unsigned int i = 0; i < dofs_per_component; ++i)
488  {
489  for (unsigned int j = 0; j < dofs_per_component;
490  ++j)
491  values_dofs[j][0] = static_cast<Number>(i == j);
492 
494  dim,
495  Number,
496  VectorizedArrayType>::apply(1,
497  phi.get_shape_info()
498  .data.front()
499  .fe_degree,
500  phi.get_shape_info(),
501  false,
502  constraint_mask,
503  values_dofs.data());
504 
505  for (unsigned int j = 0; j < dofs_per_component;
506  ++j)
507  if (1e-10 < std::abs(values_dofs[j][0]) &&
508  (j != i ||
509  1e-10 < std::abs(values_dofs[j][0] - 1.0)))
510  locally_relevant_constrains_hn.emplace_back(
511  j, i, values_dofs[j][0]);
512  }
513 
514 
515  std::sort(locally_relevant_constrains_hn.begin(),
516  locally_relevant_constrains_hn.end(),
517  [](const auto &a, const auto &b) {
518  if (std::get<0>(a) < std::get<0>(b))
519  return true;
520  return (std::get<0>(a) == std::get<0>(b)) &&
521  (std::get<1>(a) < std::get<1>(b));
522  });
523 
524  // 1b) extend for multiple components
525  const unsigned int n_hn_constraints =
526  locally_relevant_constrains_hn.size();
527  locally_relevant_constrains_hn.resize(n_hn_constraints *
528  n_components);
529 
530  for (unsigned int c = 0; c < n_components; ++c)
531  for (unsigned int i = 0; i < n_hn_constraints; ++i)
532  locally_relevant_constrains_hn[c *
533  n_hn_constraints +
534  i] =
535  std::tuple<unsigned int, unsigned int, Number>{
536  std::get<0>(locally_relevant_constrains_hn[i]) +
537  c * dofs_per_component,
538  std::get<1>(locally_relevant_constrains_hn[i]) +
539  c * dofs_per_component,
540  std::get<2>(locally_relevant_constrains_hn[i])};
541 
542 
543  locally_relevant_constrains_hn_map[mask] =
544  locally_relevant_constrains_hn;
545  }
546 
547  const auto &locally_relevant_constrains_hn =
548  locally_relevant_constrains_hn_map[mask];
549 
550  // 2) perform vmult with other constraints
551  std::vector<std::tuple<unsigned int, unsigned int, Number>>
552  locally_relevant_constrains_temp;
553 
554  for (unsigned int i = 0;
555  i < dofs_per_component * n_components;
556  ++i)
557  {
558  const auto lower_bound_fu = [](const auto &a,
559  const auto &b) {
560  return std::get<0>(a) < b;
561  };
562 
563  const auto upper_bound_fu = [](const auto &a,
564  const auto &b) {
565  return a < std::get<0>(b);
566  };
567 
568  const auto i_begin = std::lower_bound(
569  locally_relevant_constrains_hn.begin(),
570  locally_relevant_constrains_hn.end(),
571  i,
572  lower_bound_fu);
573  const auto i_end = std::upper_bound(
574  locally_relevant_constrains_hn.begin(),
575  locally_relevant_constrains_hn.end(),
576  i,
577  upper_bound_fu);
578 
579  if (i_begin == i_end)
580  {
581  // dof is not constrained by hanging-node constraint
582  // (identity matrix): simply copy constraints
583  const auto j_begin = std::lower_bound(
584  locally_relevant_constrains.begin(),
585  locally_relevant_constrains.end(),
586  i,
587  lower_bound_fu);
588  const auto j_end = std::upper_bound(
589  locally_relevant_constrains.begin(),
590  locally_relevant_constrains.end(),
591  i,
592  upper_bound_fu);
593 
594  for (auto v = j_begin; v != j_end; ++v)
595  locally_relevant_constrains_temp.emplace_back(*v);
596  }
597  else
598  {
599  // dof is constrained: build transitive closure
600  for (auto v0 = i_begin; v0 != i_end; ++v0)
601  {
602  const auto j_begin = std::lower_bound(
603  locally_relevant_constrains.begin(),
604  locally_relevant_constrains.end(),
605  std::get<1>(*v0),
606  lower_bound_fu);
607  const auto j_end = std::upper_bound(
608  locally_relevant_constrains.begin(),
609  locally_relevant_constrains.end(),
610  std::get<1>(*v0),
611  upper_bound_fu);
612 
613  for (auto v1 = j_begin; v1 != j_end; ++v1)
614  locally_relevant_constrains_temp.emplace_back(
615  std::get<0>(*v0),
616  std::get<1>(*v1),
617  std::get<2>(*v0) * std::get<2>(*v1));
618  }
619  }
620  }
621 
622  locally_relevant_constrains =
623  locally_relevant_constrains_temp;
624  }
625  }
626 
627  // STEP 2d: transpose COO
628  std::sort(locally_relevant_constrains.begin(),
629  locally_relevant_constrains.end(),
630  [](const auto &a, const auto &b) {
631  if (std::get<1>(a) < std::get<1>(b))
632  return true;
633  return (std::get<1>(a) == std::get<1>(b)) &&
634  (std::get<0>(a) < std::get<0>(b));
635  });
636 
637  // STEP 2e: translate COO to CRS
638  auto &c_pool = c_pools[v];
639  {
640  if (locally_relevant_constrains.size() > 0)
641  c_pool.row_lid_to_gid.emplace_back(
642  std::get<1>(locally_relevant_constrains.front()));
643  for (const auto &j : locally_relevant_constrains)
644  {
645  if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
646  {
647  c_pool.row_lid_to_gid.push_back(std::get<1>(j));
648  c_pool.row.push_back(c_pool.val.size());
649  }
650 
651  c_pool.col.emplace_back(std::get<0>(j));
652  c_pool.val.emplace_back(std::get<2>(j));
653  }
654 
655  if (c_pool.val.size() > 0)
656  c_pool.row.push_back(c_pool.val.size());
657  }
658  }
659  // STEP 3: compute element matrix A_e, apply
660  // locally-relevant constraints C_e^T * A_e * C_e, and get the
661  // the diagonal entry
662  // (C_e^T * A_e * C_e)(i,i)
663  // or
664  // C_e^T(i,:) * A_e * C_e(:,i).
665  //
666  // Since, we compute the element matrix column-by-column and as a
667  // result never actually have the full element matrix, we actually
668  // perform following steps:
669  // 1) loop over all columns of the element matrix
670  // a) compute column i
671  // b) compute for each j (rows of C_e^T):
672  // (C_e^T(j,:) * A_e(:,i)) * C_e(i,j)
673  // or
674  // (C_e^T(j,:) * A_e(:,i)) * C_e^T(j,i)
675  // This gives a contribution the j-th entry of the
676  // locally-relevant diagonal and comprises the multiplication
677  // by the locally-relevant constraint matrix from the left and
678  // the right. There is no contribution to the j-th vector
679  // entry if the j-th row of C_e^T is empty or C_e^T(j,i) is
680  // zero.
681 
682  // set size locally-relevant diagonal
683  for (unsigned int v = 0; v < n_lanes_filled; ++v)
684  diagonals_local_constrained[v].assign(
685  c_pools[v].row_lid_to_gid.size() *
686  (n_fe_components == 1 ? n_components : 1),
687  Number(0.0));
688  }
689 
690  void
691  prepare_basis_vector(const unsigned int i)
692  {
693  this->i = i;
694 
695  // compute i-th column of element stiffness matrix:
696  // this could be simply performed as done at the moment with
697  // matrix-free operator evaluation applied to a ith-basis vector
698  for (unsigned int j = 0; j < phi.dofs_per_cell; ++j)
699  phi.begin_dof_values()[j] = static_cast<Number>(i == j);
700  }
701 
702  void
703  submit()
704  {
705  // if we have a block vector with components with the same DoFHandler,
706  // we need to figure out which component and which DoF within the
707  // component are we currently considering
708  const unsigned int n_fe_components =
709  phi.get_dof_info().start_components.back();
710  const unsigned int comp =
711  n_fe_components == 1 ? i / phi.dofs_per_component : 0;
712  const unsigned int i_comp =
713  n_fe_components == 1 ? (i % phi.dofs_per_component) : i;
714 
715  // apply local constraint matrix from left and from right:
716  // loop over all rows of transposed constrained matrix
717  for (unsigned int v = 0;
718  v < phi.get_matrix_free().n_active_entries_per_cell_batch(
719  phi.get_current_cell_index());
720  ++v)
721  {
722  const auto &c_pool = c_pools[v];
723 
724  for (unsigned int j = 0; j < c_pool.row.size() - 1; ++j)
725  {
726  // check if the result will be zero, so that we can skip
727  // the following computations -> binary search
728  const auto scale_iterator =
729  std::lower_bound(c_pool.col.begin() + c_pool.row[j],
730  c_pool.col.begin() + c_pool.row[j + 1],
731  i_comp);
732 
733  // explanation: j-th row of C_e^T is empty (see above)
734  if (scale_iterator == c_pool.col.begin() + c_pool.row[j + 1])
735  continue;
736 
737  // explanation: C_e^T(j,i) is zero (see above)
738  if (*scale_iterator != i_comp)
739  continue;
740 
741  // apply constraint matrix from the left
742  Number temp = 0.0;
743  for (unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
744  temp += c_pool.val[k] *
745  phi.begin_dof_values()[comp * phi.dofs_per_component +
746  c_pool.col[k]][v];
747 
748  // apply constraint matrix from the right
749  diagonals_local_constrained
750  [v][j + comp * c_pools[v].row_lid_to_gid.size()] +=
751  temp *
752  c_pool.val[std::distance(c_pool.col.begin(), scale_iterator)];
753  }
754  }
755  }
756 
757  template <typename VectorType>
758  inline void
759  distribute_local_to_global(
760  std::array<VectorType *, n_components> &diagonal_global)
761  {
762  // STEP 4: assembly results: add into global vector
763  const unsigned int n_fe_components =
764  phi.get_dof_info().start_components.back();
765 
766  for (unsigned int v = 0;
767  v < phi.get_matrix_free().n_active_entries_per_cell_batch(
768  phi.get_current_cell_index());
769  ++v)
770  // if we have a block vector with components with the same
771  // DoFHandler, we need to loop over all components manually and
772  // need to apply the correct shift
773  for (unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
774  for (unsigned int comp = 0;
775  comp < (n_fe_components == 1 ?
776  static_cast<unsigned int>(n_components) :
777  1);
778  ++comp)
780  *diagonal_global[n_fe_components == 1 ? comp : 0],
781  c_pools[v].row_lid_to_gid[j],
782  diagonals_local_constrained
783  [v][j + comp * c_pools[v].row_lid_to_gid.size()]);
784  }
785 
786  private:
787  FEEvaluation<dim,
788  fe_degree,
789  n_q_points_1d,
790  n_components,
791  Number,
792  VectorizedArrayType> &phi;
793 
794  unsigned int i;
795 
796  std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
797 
798  // local storage: buffer so that we access the global vector once
799  // note: may be larger then dofs_per_cell in the presence of
800  // constraints!
801  std::array<std::vector<Number>, n_lanes> diagonals_local_constrained;
802 
803  std::map<
805  std::vector<std::tuple<unsigned int, unsigned int, Number>>>
806  locally_relevant_constrains_hn_map;
807  };
808 
809  } // namespace internal
810 
811  template <int dim,
812  int fe_degree,
813  int n_q_points_1d,
814  int n_components,
815  typename Number,
816  typename VectorizedArrayType,
817  typename VectorType>
818  void
821  VectorType & diagonal_global,
822  const std::function<void(FEEvaluation<dim,
823  fe_degree,
824  n_q_points_1d,
825  n_components,
826  Number,
827  VectorizedArrayType> &)> &local_vmult,
828  const unsigned int dof_no,
829  const unsigned int quad_no,
830  const unsigned int first_selected_component)
831  {
832  int dummy = 0;
833 
834  std::array<typename ::internal::BlockVectorSelector<
835  VectorType,
836  IsBlockVector<VectorType>::value>::BaseVectorType *,
837  n_components>
838  diagonal_global_components;
839 
840  for (unsigned int d = 0; d < n_components; ++d)
841  diagonal_global_components[d] = ::internal::
842  BlockVectorSelector<VectorType, IsBlockVector<VectorType>::value>::
843  get_vector_component(diagonal_global, d + first_selected_component);
844 
845  const auto &dof_info = matrix_free.get_dof_info(dof_no);
846 
847  if (dof_info.start_components.back() == 1)
848  for (unsigned int comp = 0; comp < n_components; ++comp)
849  {
850  Assert(diagonal_global_components[comp] != nullptr,
851  ExcMessage("The finite element underlying this FEEvaluation "
852  "object is scalar, but you requested " +
853  std::to_string(n_components) +
854  " components via the template argument in "
855  "FEEvaluation. In that case, you must pass an "
856  "std::vector<VectorType> or a BlockVector to " +
857  "read_dof_values and distribute_local_to_global."));
859  *diagonal_global_components[comp], matrix_free, dof_info);
860  }
861  else
862  {
864  *diagonal_global_components[0], matrix_free, dof_info);
865  }
866 
867  matrix_free.template cell_loop<VectorType, int>(
868  [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
869  VectorType &,
870  const int &,
871  const std::pair<unsigned int, unsigned int> &range) mutable {
872  FEEvaluation<dim,
873  fe_degree,
874  n_q_points_1d,
875  n_components,
876  Number,
877  VectorizedArrayType>
878  phi(matrix_free, range, dof_no, quad_no, first_selected_component);
879 
880  internal::ComputeDiagonalHelper<dim,
881  fe_degree,
882  n_q_points_1d,
883  n_components,
884  Number,
885  VectorizedArrayType>
886  helper(phi);
887 
888  for (unsigned int cell = range.first; cell < range.second; ++cell)
889  {
890  helper.reinit(cell);
891 
892  for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
893  {
894  helper.prepare_basis_vector(i);
895  local_vmult(phi);
896  helper.submit();
897  }
898 
899  helper.distribute_local_to_global(diagonal_global_components);
900  }
901  },
902  diagonal_global,
903  dummy,
904  false);
905  }
906 
907  template <typename CLASS,
908  int dim,
909  int fe_degree,
910  int n_q_points_1d,
911  int n_components,
912  typename Number,
913  typename VectorizedArrayType,
914  typename VectorType>
915  void
918  VectorType & diagonal_global,
919  void (CLASS::*cell_operation)(FEEvaluation<dim,
920  fe_degree,
921  n_q_points_1d,
922  n_components,
923  Number,
924  VectorizedArrayType> &) const,
925  const CLASS * owning_class,
926  const unsigned int dof_no,
927  const unsigned int quad_no,
928  const unsigned int first_selected_component)
929  {
930  compute_diagonal<dim,
931  fe_degree,
932  n_q_points_1d,
933  n_components,
934  Number,
935  VectorizedArrayType,
936  VectorType>(
937  matrix_free,
938  diagonal_global,
939  [&](auto &feeval) { (owning_class->*cell_operation)(feeval); },
940  dof_no,
941  quad_no,
942  first_selected_component);
943  }
944 
945  namespace internal
946  {
951  template <typename MatrixType,
952  typename Number,
953  typename std::enable_if<std::is_same<
954  typename std::remove_const<typename std::remove_reference<
955  typename MatrixType::value_type>::type>::type,
956  typename std::remove_const<typename std::remove_reference<
957  Number>::type>::type>::value>::type * = nullptr>
959  create_new_affine_constraints_if_needed(
960  const MatrixType &,
961  const AffineConstraints<Number> &constraints,
963  {
964  return constraints;
965  }
966 
972  template <typename MatrixType,
973  typename Number,
974  typename std::enable_if<!std::is_same<
975  typename std::remove_const<typename std::remove_reference<
976  typename MatrixType::value_type>::type>::type,
977  typename std::remove_const<typename std::remove_reference<
978  Number>::type>::type>::value>::type * = nullptr>
980  create_new_affine_constraints_if_needed(
981  const MatrixType &,
982  const AffineConstraints<Number> &constraints,
984  &new_constraints)
985  {
986  new_constraints =
987  std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
988  new_constraints->copy_from(constraints);
989 
990  return *new_constraints;
991  }
992  } // namespace internal
993 
994  template <int dim,
995  int fe_degree,
996  int n_q_points_1d,
997  int n_components,
998  typename Number,
999  typename VectorizedArrayType,
1000  typename MatrixType>
1001  void
1004  const AffineConstraints<Number> & constraints_in,
1005  MatrixType & matrix,
1006  const std::function<void(FEEvaluation<dim,
1007  fe_degree,
1008  n_q_points_1d,
1009  n_components,
1010  Number,
1011  VectorizedArrayType> &)> &local_vmult,
1012  const unsigned int dof_no,
1013  const unsigned int quad_no,
1014  const unsigned int first_selected_component)
1015  {
1016  std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
1017  constraints_for_matrix;
1019  internal::create_new_affine_constraints_if_needed(matrix,
1020  constraints_in,
1021  constraints_for_matrix);
1022 
1023  matrix_free.template cell_loop<MatrixType, MatrixType>(
1024  [&](const auto &, auto &dst, const auto &, const auto range) {
1025  FEEvaluation<dim,
1026  fe_degree,
1027  n_q_points_1d,
1028  n_components,
1029  Number,
1030  VectorizedArrayType>
1031  integrator(
1032  matrix_free, range, dof_no, quad_no, first_selected_component);
1033 
1034  const unsigned int dofs_per_cell = integrator.dofs_per_cell;
1035 
1036  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
1037  std::vector<types::global_dof_index> dof_indices_mf(dofs_per_cell);
1038 
1039  std::array<FullMatrix<typename MatrixType::value_type>,
1040  VectorizedArrayType::size()>
1041  matrices;
1042 
1043  std::fill_n(matrices.begin(),
1044  VectorizedArrayType::size(),
1046  dofs_per_cell));
1047 
1048  const auto lexicographic_numbering =
1049  matrix_free
1050  .get_shape_info(dof_no,
1051  quad_no,
1052  first_selected_component,
1053  integrator.get_active_fe_index(),
1054  integrator.get_active_quadrature_index())
1055  .lexicographic_numbering;
1056 
1057  for (auto cell = range.first; cell < range.second; ++cell)
1058  {
1059  integrator.reinit(cell);
1060 
1061  const unsigned int n_filled_lanes =
1062  matrix_free.n_active_entries_per_cell_batch(cell);
1063 
1064  for (unsigned int v = 0; v < n_filled_lanes; ++v)
1065  matrices[v] = 0.0;
1066 
1067  for (unsigned int j = 0; j < dofs_per_cell; ++j)
1068  {
1069  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1070  integrator.begin_dof_values()[i] =
1071  static_cast<Number>(i == j);
1072 
1073  local_vmult(integrator);
1074 
1075  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1076  for (unsigned int v = 0; v < n_filled_lanes; ++v)
1077  matrices[v](i, j) = integrator.begin_dof_values()[i][v];
1078  }
1079 
1080  for (unsigned int v = 0; v < n_filled_lanes; ++v)
1081  {
1082  const auto cell_v =
1083  matrix_free.get_cell_iterator(cell, v, dof_no);
1084 
1085  if (matrix_free.get_mg_level() != numbers::invalid_unsigned_int)
1086  cell_v->get_mg_dof_indices(dof_indices);
1087  else
1088  cell_v->get_dof_indices(dof_indices);
1089 
1090  for (unsigned int j = 0; j < dof_indices.size(); ++j)
1091  dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]];
1092 
1093  constraints.distribute_local_to_global(matrices[v],
1094  dof_indices_mf,
1095  dst);
1096  }
1097  }
1098  },
1099  matrix,
1100  matrix);
1101 
1102  matrix.compress(VectorOperation::add);
1103  }
1104 
1105  template <typename CLASS,
1106  int dim,
1107  int fe_degree,
1108  int n_q_points_1d,
1109  int n_components,
1110  typename Number,
1111  typename VectorizedArrayType,
1112  typename MatrixType>
1113  void
1116  const AffineConstraints<Number> & constraints,
1117  MatrixType & matrix,
1118  void (CLASS::*cell_operation)(FEEvaluation<dim,
1119  fe_degree,
1120  n_q_points_1d,
1121  n_components,
1122  Number,
1123  VectorizedArrayType> &) const,
1124  const CLASS * owning_class,
1125  const unsigned int dof_no,
1126  const unsigned int quad_no,
1127  const unsigned int first_selected_component)
1128  {
1129  compute_matrix<dim,
1130  fe_degree,
1131  n_q_points_1d,
1132  n_components,
1133  Number,
1134  VectorizedArrayType,
1135  MatrixType>(
1136  matrix_free,
1137  constraints,
1138  matrix,
1139  [&](auto &feeval) { (owning_class->*cell_operation)(feeval); },
1140  dof_no,
1141  quad_no,
1142  first_selected_component);
1143  }
1144 
1145 #endif // DOXYGEN
1146 
1147 } // namespace MatrixFreeTools
1148 
1150 
1151 
1152 #endif
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void copy_from(const AffineConstraints< other_number > &other)
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
unsigned int get_mg_level() const
const Number * constraint_pool_begin(const unsigned int pool_index) const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
const Number * constraint_pool_end(const unsigned int pool_index) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
virtual std::vector< types::boundary_id > get_boundary_ids() const
cell_iterator begin(const unsigned int level=0) const
unsigned int n_active_cells() const
cell_iterator end() const
active_cell_iterator begin_active(const unsigned int level=0) 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
unsigned int level
Definition: grid_out.cc:4607
const unsigned int v0
Definition: grid_tools.cc:1000
const unsigned int v1
Definition: grid_tools.cc:1000
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::string to_string(const T &t)
Definition: patterns.h:2403
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
@ matrix
Contents is actually a matrix.
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, VectorType &diagonal_global, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &local_vmult, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
void categorize_by_boundary_ids(const Triangulation< dim > &tria, AdditionalData &additional_data)
void compute_matrix(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const AffineConstraints< Number > &constraints, MatrixType &matrix, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &local_vmult, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1153
std::uint8_t compressed_constraint_kind
Definition: dof_info.h:83
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
void vector_access_add(VectorType &vec, const unsigned int entry, const typename VectorType::value_type &val)
void check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
static const unsigned int invalid_unsigned_int
Definition: types.h:201
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:36
const ::Triangulation< dim, spacedim > & tria