Reference documentation for deal.II version GIT 75f1417c0d 2023-02-03 16:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
169  template <int dim,
170  typename Number,
171  typename VectorizedArrayType = VectorizedArray<Number>>
173  {
174  public:
180  {
184  AdditionalData(const unsigned int dof_index = 0)
186  {}
187 
191  unsigned int dof_index;
192  };
193 
203  void
205  const AdditionalData &additional_data = AdditionalData())
206  {
207  this->matrix_free = &matrix_free;
208 
209  std::vector<unsigned int> valid_fe_indices;
210 
211  const auto &fe_collection =
212  matrix_free.get_dof_handler(additional_data.dof_index)
213  .get_fe_collection();
214 
215  for (unsigned int i = 0; i < fe_collection.size(); ++i)
216  if (fe_collection[i].n_dofs_per_cell() > 0)
217  valid_fe_indices.push_back(i);
218 
219  // TODO: relax this so that arbitrary number of valid
220  // FEs can be accepted
221  AssertDimension(valid_fe_indices.size(), 1);
222 
223  fe_index_valid = *valid_fe_indices.begin();
224  }
225 
231  template <typename VectorTypeOut, typename VectorTypeIn>
232  void
233  cell_loop(const std::function<void(
235  VectorTypeOut &,
236  const VectorTypeIn &,
237  const std::pair<unsigned int, unsigned int>)> &cell_operation,
238  VectorTypeOut & dst,
239  const VectorTypeIn & src,
240  const bool zero_dst_vector = false)
241  {
242  const auto ebd_cell_operation = [&](const auto &matrix_free,
243  auto & dst,
244  const auto &src,
245  const auto range) {
246  const auto category = matrix_free.get_cell_range_category(range);
247 
248  if (category != fe_index_valid)
249  return;
250 
251  cell_operation(matrix_free, dst, src, range);
252  };
253 
254  matrix_free->template cell_loop<VectorTypeOut, VectorTypeIn>(
255  ebd_cell_operation, dst, src, zero_dst_vector);
256  }
257 
265  template <typename VectorTypeOut, typename VectorTypeIn>
266  void
267  loop(const std::function<
269  VectorTypeOut &,
270  const VectorTypeIn &,
271  const std::pair<unsigned int, unsigned int>)> &cell_operation,
272  const std::function<
274  VectorTypeOut &,
275  const VectorTypeIn &,
276  const std::pair<unsigned int, unsigned int>)> &face_operation,
277  const std::function<
279  VectorTypeOut &,
280  const VectorTypeIn &,
281  const std::pair<unsigned int, unsigned int>,
282  const bool)> &boundary_operation,
283  VectorTypeOut & dst,
284  const VectorTypeIn & src,
285  const bool zero_dst_vector = false)
286  {
287  const auto ebd_cell_operation = [&](const auto &matrix_free,
288  auto & dst,
289  const auto &src,
290  const auto range) {
291  const auto category = matrix_free.get_cell_range_category(range);
292 
293  if (category != fe_index_valid)
294  return;
295 
296  cell_operation(matrix_free, dst, src, range);
297  };
298 
299  const auto ebd_internal_or_boundary_face_operation =
300  [&](const auto &matrix_free,
301  auto & dst,
302  const auto &src,
303  const auto range) {
304  const auto category = matrix_free.get_face_range_category(range);
305 
306  const unsigned int type =
307  static_cast<unsigned int>(category.first == fe_index_valid) +
308  static_cast<unsigned int>(category.second == fe_index_valid);
309 
310  if (type == 0)
311  return; // deactivated face -> nothing to do
312 
313  if (type == 1) // boundary face
314  boundary_operation(
315  matrix_free, dst, src, range, category.first == fe_index_valid);
316  else if (type == 2) // internal face
317  face_operation(matrix_free, dst, src, range);
318  };
319 
320  matrix_free->template loop<VectorTypeOut, VectorTypeIn>(
321  ebd_cell_operation,
322  ebd_internal_or_boundary_face_operation,
323  ebd_internal_or_boundary_face_operation,
324  dst,
325  src,
326  zero_dst_vector);
327  }
328 
329  private:
335 
339  unsigned int fe_index_valid;
340  };
341 
342  // implementations
343 
344 #ifndef DOXYGEN
345 
346  template <int dim, typename AdditionalData>
347  void
349  AdditionalData & additional_data)
350  {
351  // ... determine if we are on an active or a multigrid level
352  const unsigned int level = additional_data.mg_level;
353  const bool is_mg = (level != numbers::invalid_unsigned_int);
354 
355  // ... create empty list for the category of each cell
356  if (is_mg)
357  additional_data.cell_vectorization_category.assign(
358  std::distance(tria.begin(level), tria.end(level)), 0);
359  else
360  additional_data.cell_vectorization_category.assign(tria.n_active_cells(),
361  0);
362 
363  // ... set up scaling factor
364  std::vector<unsigned int> factors(GeometryInfo<dim>::faces_per_cell);
365 
366  auto bids = tria.get_boundary_ids();
367  std::sort(bids.begin(), bids.end());
368 
369  {
370  unsigned int n_bids = bids.size() + 1;
371  int offset = 1;
372  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
373  i++, offset = offset * n_bids)
374  factors[i] = offset;
375  }
376 
377  const auto to_category = [&](const auto &cell) {
378  unsigned int c_num = 0;
379  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
380  {
381  auto &face = *cell->face(i);
382  if (face.at_boundary() && !cell->has_periodic_neighbor(i))
383  c_num +=
384  factors[i] * (1 + std::distance(bids.begin(),
385  std::find(bids.begin(),
386  bids.end(),
387  face.boundary_id())));
388  }
389  return c_num;
390  };
391 
392  if (!is_mg)
393  {
394  for (auto cell = tria.begin_active(); cell != tria.end(); ++cell)
395  {
396  if (cell->is_locally_owned())
397  additional_data
398  .cell_vectorization_category[cell->active_cell_index()] =
399  to_category(cell);
400  }
401  }
402  else
403  {
404  for (auto cell = tria.begin(level); cell != tria.end(level); ++cell)
405  {
406  if (cell->is_locally_owned_on_level())
407  additional_data.cell_vectorization_category[cell->index()] =
408  to_category(cell);
409  }
410  }
411 
412  // ... finalize set up of matrix_free
413  additional_data.hold_all_faces_to_owned_cells = true;
414  additional_data.cell_vectorization_categories_strict = true;
415  additional_data.mapping_update_flags_faces_by_cells =
416  additional_data.mapping_update_flags_inner_faces |
417  additional_data.mapping_update_flags_boundary_faces;
418  }
419 
420  namespace internal
421  {
422  template <typename Number>
423  struct LocalCSR
424  {
425  LocalCSR()
426  : row{0}
427  {}
428 
429  std::vector<unsigned int> row_lid_to_gid;
430  std::vector<unsigned int> row;
431  std::vector<unsigned int> col;
432  std::vector<Number> val;
433  };
434 
435  template <int dim,
436  int fe_degree,
437  int n_q_points_1d,
438  int n_components,
439  typename Number,
440  typename VectorizedArrayType>
441  class ComputeDiagonalHelper
442  {
443  public:
444  static const unsigned int n_lanes = VectorizedArrayType::size();
445 
446  ComputeDiagonalHelper(FEEvaluation<dim,
447  fe_degree,
448  n_q_points_1d,
449  n_components,
450  Number,
451  VectorizedArrayType> &phi)
452  : phi(phi)
453  {}
454 
455  void
456  reinit(const unsigned int cell)
457  {
458  this->phi.reinit(cell);
459  // STEP 1: get relevant information from FEEvaluation
460  const auto & dof_info = phi.get_dof_info();
461  const unsigned int n_fe_components = dof_info.start_components.back();
462  const unsigned int dofs_per_component = phi.dofs_per_component;
463  const auto & matrix_free = phi.get_matrix_free();
464 
465  // if we have a block vector with components with the same DoFHandler,
466  // each component is described with same set of constraints and
467  // we consider the shift in components only during access of the global
468  // vector
469  const unsigned int first_selected_component =
470  n_fe_components == 1 ? 0 : phi.get_first_selected_component();
471 
472  const unsigned int n_lanes_filled =
473  matrix_free.n_active_entries_per_cell_batch(cell);
474 
475  std::array<const unsigned int *, n_lanes> dof_indices{};
476  {
477  for (unsigned int v = 0; v < n_lanes_filled; ++v)
478  dof_indices[v] =
479  dof_info.dof_indices.data() +
480  dof_info
481  .row_starts[(cell * n_lanes + v) * n_fe_components +
482  first_selected_component]
483  .first;
484  }
485 
486  // STEP 2: setup CSR storage of transposed locally-relevant
487  // constraint matrix
488  c_pools = std::array<internal::LocalCSR<Number>, n_lanes>();
489 
490  for (unsigned int v = 0; v < n_lanes_filled; ++v)
491  {
492  unsigned int index_indicators, next_index_indicators;
493 
494  index_indicators =
495  dof_info
496  .row_starts[(cell * n_lanes + v) * n_fe_components +
497  first_selected_component]
498  .second;
499  next_index_indicators =
500  dof_info
501  .row_starts[(cell * n_lanes + v) * n_fe_components +
502  first_selected_component + 1]
503  .second;
504 
505  // STEP 2a: setup locally-relevant constraint matrix in a
506  // coordinate list (COO)
507  std::vector<std::tuple<unsigned int, unsigned int, Number>>
508  locally_relevant_constrains; // (constrained local index,
509  // global index of dof which
510  // constrains, weight)
511 
512  if (n_components == 1 || n_fe_components == 1)
513  {
514  unsigned int ind_local = 0;
515  for (; index_indicators != next_index_indicators;
516  ++index_indicators, ++ind_local)
517  {
518  const std::pair<unsigned short, unsigned short> indicator =
519  dof_info.constraint_indicator[index_indicators];
520 
521  for (unsigned int j = 0; j < indicator.first;
522  ++j, ++ind_local)
523  locally_relevant_constrains.emplace_back(
524  ind_local, dof_indices[v][j], 1.0);
525 
526  dof_indices[v] += indicator.first;
527 
528  const Number *data_val =
529  matrix_free.constraint_pool_begin(indicator.second);
530  const Number *end_pool =
531  matrix_free.constraint_pool_end(indicator.second);
532 
533  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
534  locally_relevant_constrains.emplace_back(ind_local,
535  *dof_indices[v],
536  *data_val);
537  }
538 
539  AssertIndexRange(ind_local, dofs_per_component + 1);
540 
541  for (; ind_local < dofs_per_component;
542  ++dof_indices[v], ++ind_local)
543  locally_relevant_constrains.emplace_back(ind_local,
544  *dof_indices[v],
545  1.0);
546  }
547  else
548  {
549  // case with vector-valued finite elements where all
550  // components are included in one single vector. Assumption:
551  // first come all entries to the first component, then all
552  // entries to the second one, and so on. This is ensured by
553  // the way MatrixFree reads out the indices.
554  for (unsigned int comp = 0; comp < n_components; ++comp)
555  {
556  unsigned int ind_local = 0;
557 
558  // check whether there is any constraint on the current
559  // cell
560  for (; index_indicators != next_index_indicators;
561  ++index_indicators, ++ind_local)
562  {
563  const std::pair<unsigned short, unsigned short>
564  indicator =
565  dof_info.constraint_indicator[index_indicators];
566 
567  // run through values up to next constraint
568  for (unsigned int j = 0; j < indicator.first;
569  ++j, ++ind_local)
570  locally_relevant_constrains.emplace_back(
571  comp * dofs_per_component + ind_local,
572  dof_indices[v][j],
573  1.0);
574  dof_indices[v] += indicator.first;
575 
576  const Number *data_val =
577  matrix_free.constraint_pool_begin(indicator.second);
578  const Number *end_pool =
579  matrix_free.constraint_pool_end(indicator.second);
580 
581  for (; data_val != end_pool;
582  ++data_val, ++dof_indices[v])
583  locally_relevant_constrains.emplace_back(
584  comp * dofs_per_component + ind_local,
585  *dof_indices[v],
586  *data_val);
587  }
588 
589  AssertIndexRange(ind_local, dofs_per_component + 1);
590 
591  // get the dof values past the last constraint
592  for (; ind_local < dofs_per_component;
593  ++dof_indices[v], ++ind_local)
594  locally_relevant_constrains.emplace_back(
595  comp * dofs_per_component + ind_local,
596  *dof_indices[v],
597  1.0);
598 
599  if (comp + 1 < n_components)
600  {
601  next_index_indicators =
602  dof_info
603  .row_starts[(cell * n_lanes + v) * n_fe_components +
604  first_selected_component + comp + 2]
605  .second;
606  }
607  }
608  }
609  // STEP 2b: sort and make unique
610 
611  // sort vector
612  std::sort(locally_relevant_constrains.begin(),
613  locally_relevant_constrains.end(),
614  [](const auto &a, const auto &b) {
615  if (std::get<0>(a) < std::get<0>(b))
616  return true;
617  return (std::get<0>(a) == std::get<0>(b)) &&
618  (std::get<1>(a) < std::get<1>(b));
619  });
620 
621  // make sure that all entries are unique
622  locally_relevant_constrains.erase(
623  unique(locally_relevant_constrains.begin(),
624  locally_relevant_constrains.end(),
625  [](const auto &a, const auto &b) {
626  return (std::get<1>(a) == std::get<1>(b)) &&
627  (std::get<0>(a) == std::get<0>(b));
628  }),
629  locally_relevant_constrains.end());
630 
631  // STEP 2c: apply hanging-node constraints
632  if (dof_info.hanging_node_constraint_masks.size() > 0 &&
633  dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
634  dof_info
635  .hanging_node_constraint_masks_comp[phi.get_active_fe_index()]
636  [first_selected_component])
637  {
638  const auto mask =
639  dof_info.hanging_node_constraint_masks[cell * n_lanes + v];
640 
641  // cell has hanging nodes
642  if (mask != ::internal::MatrixFreeFunctions::
644  {
645  // check if hanging node internpolation matrix has been set
646  // up
647  if (locally_relevant_constrains_hn_map.find(mask) ==
648  locally_relevant_constrains_hn_map.end())
649  {
650  // 1) collect hanging-node constraints for cell assuming
651  // scalar finite element
653  dofs_per_component);
654 
657  VectorizedArrayType::size()>
658  constraint_mask;
659  constraint_mask.fill(
662 
663  constraint_mask[0] = mask;
664 
665  std::vector<
666  std::tuple<unsigned int, unsigned int, Number>>
667  locally_relevant_constrains_hn;
668 
669  for (unsigned int i = 0; i < dofs_per_component; ++i)
670  {
671  for (unsigned int j = 0; j < dofs_per_component;
672  ++j)
673  values_dofs[j][0] = static_cast<Number>(i == j);
674 
676  dim,
677  Number,
678  VectorizedArrayType>::apply(1,
679  phi.get_shape_info()
680  .data.front()
681  .fe_degree,
682  phi.get_shape_info(),
683  false,
684  constraint_mask,
685  values_dofs.data());
686 
687  for (unsigned int j = 0; j < dofs_per_component;
688  ++j)
689  if (1e-10 < std::abs(values_dofs[j][0]) &&
690  (j != i ||
691  1e-10 < std::abs(values_dofs[j][0] - 1.0)))
692  locally_relevant_constrains_hn.emplace_back(
693  j, i, values_dofs[j][0]);
694  }
695 
696 
697  std::sort(locally_relevant_constrains_hn.begin(),
698  locally_relevant_constrains_hn.end(),
699  [](const auto &a, const auto &b) {
700  if (std::get<0>(a) < std::get<0>(b))
701  return true;
702  return (std::get<0>(a) == std::get<0>(b)) &&
703  (std::get<1>(a) < std::get<1>(b));
704  });
705 
706  // 1b) extend for multiple components
707  const unsigned int n_hn_constraints =
708  locally_relevant_constrains_hn.size();
709  locally_relevant_constrains_hn.resize(n_hn_constraints *
710  n_components);
711 
712  for (unsigned int c = 0; c < n_components; ++c)
713  for (unsigned int i = 0; i < n_hn_constraints; ++i)
714  locally_relevant_constrains_hn[c *
715  n_hn_constraints +
716  i] =
717  std::tuple<unsigned int, unsigned int, Number>{
718  std::get<0>(locally_relevant_constrains_hn[i]) +
719  c * dofs_per_component,
720  std::get<1>(locally_relevant_constrains_hn[i]) +
721  c * dofs_per_component,
722  std::get<2>(locally_relevant_constrains_hn[i])};
723 
724 
725  locally_relevant_constrains_hn_map[mask] =
726  locally_relevant_constrains_hn;
727  }
728 
729  const auto &locally_relevant_constrains_hn =
730  locally_relevant_constrains_hn_map[mask];
731 
732  // 2) perform vmult with other constraints
733  std::vector<std::tuple<unsigned int, unsigned int, Number>>
734  locally_relevant_constrains_temp;
735 
736  for (unsigned int i = 0;
737  i < dofs_per_component * n_components;
738  ++i)
739  {
740  const auto lower_bound_fu = [](const auto &a,
741  const auto &b) {
742  return std::get<0>(a) < b;
743  };
744 
745  const auto upper_bound_fu = [](const auto &a,
746  const auto &b) {
747  return a < std::get<0>(b);
748  };
749 
750  const auto i_begin = std::lower_bound(
751  locally_relevant_constrains_hn.begin(),
752  locally_relevant_constrains_hn.end(),
753  i,
754  lower_bound_fu);
755  const auto i_end = std::upper_bound(
756  locally_relevant_constrains_hn.begin(),
757  locally_relevant_constrains_hn.end(),
758  i,
759  upper_bound_fu);
760 
761  if (i_begin == i_end)
762  {
763  // dof is not constrained by hanging-node constraint
764  // (identity matrix): simply copy constraints
765  const auto j_begin = std::lower_bound(
766  locally_relevant_constrains.begin(),
767  locally_relevant_constrains.end(),
768  i,
769  lower_bound_fu);
770  const auto j_end = std::upper_bound(
771  locally_relevant_constrains.begin(),
772  locally_relevant_constrains.end(),
773  i,
774  upper_bound_fu);
775 
776  for (auto v = j_begin; v != j_end; ++v)
777  locally_relevant_constrains_temp.emplace_back(*v);
778  }
779  else
780  {
781  // dof is constrained: build transitive closure
782  for (auto v0 = i_begin; v0 != i_end; ++v0)
783  {
784  const auto j_begin = std::lower_bound(
785  locally_relevant_constrains.begin(),
786  locally_relevant_constrains.end(),
787  std::get<1>(*v0),
788  lower_bound_fu);
789  const auto j_end = std::upper_bound(
790  locally_relevant_constrains.begin(),
791  locally_relevant_constrains.end(),
792  std::get<1>(*v0),
793  upper_bound_fu);
794 
795  for (auto v1 = j_begin; v1 != j_end; ++v1)
796  locally_relevant_constrains_temp.emplace_back(
797  std::get<0>(*v0),
798  std::get<1>(*v1),
799  std::get<2>(*v0) * std::get<2>(*v1));
800  }
801  }
802  }
803 
804  locally_relevant_constrains =
805  locally_relevant_constrains_temp;
806  }
807  }
808 
809  // STEP 2d: transpose COO
810  std::sort(locally_relevant_constrains.begin(),
811  locally_relevant_constrains.end(),
812  [](const auto &a, const auto &b) {
813  if (std::get<1>(a) < std::get<1>(b))
814  return true;
815  return (std::get<1>(a) == std::get<1>(b)) &&
816  (std::get<0>(a) < std::get<0>(b));
817  });
818 
819  // STEP 2e: translate COO to CRS
820  auto &c_pool = c_pools[v];
821  {
822  if (locally_relevant_constrains.size() > 0)
823  c_pool.row_lid_to_gid.emplace_back(
824  std::get<1>(locally_relevant_constrains.front()));
825  for (const auto &j : locally_relevant_constrains)
826  {
827  if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
828  {
829  c_pool.row_lid_to_gid.push_back(std::get<1>(j));
830  c_pool.row.push_back(c_pool.val.size());
831  }
832 
833  c_pool.col.emplace_back(std::get<0>(j));
834  c_pool.val.emplace_back(std::get<2>(j));
835  }
836 
837  if (c_pool.val.size() > 0)
838  c_pool.row.push_back(c_pool.val.size());
839  }
840  }
841  // STEP 3: compute element matrix A_e, apply
842  // locally-relevant constraints C_e^T * A_e * C_e, and get the
843  // the diagonal entry
844  // (C_e^T * A_e * C_e)(i,i)
845  // or
846  // C_e^T(i,:) * A_e * C_e(:,i).
847  //
848  // Since, we compute the element matrix column-by-column and as a
849  // result never actually have the full element matrix, we actually
850  // perform following steps:
851  // 1) loop over all columns of the element matrix
852  // a) compute column i
853  // b) compute for each j (rows of C_e^T):
854  // (C_e^T(j,:) * A_e(:,i)) * C_e(i,j)
855  // or
856  // (C_e^T(j,:) * A_e(:,i)) * C_e^T(j,i)
857  // This gives a contribution the j-th entry of the
858  // locally-relevant diagonal and comprises the multiplication
859  // by the locally-relevant constraint matrix from the left and
860  // the right. There is no contribution to the j-th vector
861  // entry if the j-th row of C_e^T is empty or C_e^T(j,i) is
862  // zero.
863 
864  // set size locally-relevant diagonal
865  for (unsigned int v = 0; v < n_lanes_filled; ++v)
866  diagonals_local_constrained[v].assign(
867  c_pools[v].row_lid_to_gid.size() *
868  (n_fe_components == 1 ? n_components : 1),
869  Number(0.0));
870  }
871 
872  void
873  prepare_basis_vector(const unsigned int i)
874  {
875  this->i = i;
876 
877  // compute i-th column of element stiffness matrix:
878  // this could be simply performed as done at the moment with
879  // matrix-free operator evaluation applied to a ith-basis vector
880  for (unsigned int j = 0; j < phi.dofs_per_cell; ++j)
881  phi.begin_dof_values()[j] = static_cast<Number>(i == j);
882  }
883 
884  void
885  submit()
886  {
887  // if we have a block vector with components with the same DoFHandler,
888  // we need to figure out which component and which DoF within the
889  // component are we currently considering
890  const unsigned int n_fe_components =
891  phi.get_dof_info().start_components.back();
892  const unsigned int comp =
893  n_fe_components == 1 ? i / phi.dofs_per_component : 0;
894  const unsigned int i_comp =
895  n_fe_components == 1 ? (i % phi.dofs_per_component) : i;
896 
897  // apply local constraint matrix from left and from right:
898  // loop over all rows of transposed constrained matrix
899  for (unsigned int v = 0;
900  v < phi.get_matrix_free().n_active_entries_per_cell_batch(
901  phi.get_current_cell_index());
902  ++v)
903  {
904  const auto &c_pool = c_pools[v];
905 
906  for (unsigned int j = 0; j < c_pool.row.size() - 1; ++j)
907  {
908  // check if the result will be zero, so that we can skip
909  // the following computations -> binary search
910  const auto scale_iterator =
911  std::lower_bound(c_pool.col.begin() + c_pool.row[j],
912  c_pool.col.begin() + c_pool.row[j + 1],
913  i_comp);
914 
915  // explanation: j-th row of C_e^T is empty (see above)
916  if (scale_iterator == c_pool.col.begin() + c_pool.row[j + 1])
917  continue;
918 
919  // explanation: C_e^T(j,i) is zero (see above)
920  if (*scale_iterator != i_comp)
921  continue;
922 
923  // apply constraint matrix from the left
924  Number temp = 0.0;
925  for (unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
926  temp += c_pool.val[k] *
927  phi.begin_dof_values()[comp * phi.dofs_per_component +
928  c_pool.col[k]][v];
929 
930  // apply constraint matrix from the right
931  diagonals_local_constrained
932  [v][j + comp * c_pools[v].row_lid_to_gid.size()] +=
933  temp *
934  c_pool.val[std::distance(c_pool.col.begin(), scale_iterator)];
935  }
936  }
937  }
938 
939  template <typename VectorType>
940  inline void
941  distribute_local_to_global(
942  std::array<VectorType *, n_components> &diagonal_global)
943  {
944  // STEP 4: assembly results: add into global vector
945  const unsigned int n_fe_components =
946  phi.get_dof_info().start_components.back();
947 
948  for (unsigned int v = 0;
949  v < phi.get_matrix_free().n_active_entries_per_cell_batch(
950  phi.get_current_cell_index());
951  ++v)
952  // if we have a block vector with components with the same
953  // DoFHandler, we need to loop over all components manually and
954  // need to apply the correct shift
955  for (unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
956  for (unsigned int comp = 0;
957  comp < (n_fe_components == 1 ?
958  static_cast<unsigned int>(n_components) :
959  1);
960  ++comp)
962  *diagonal_global[n_fe_components == 1 ? comp : 0],
963  c_pools[v].row_lid_to_gid[j],
964  diagonals_local_constrained
965  [v][j + comp * c_pools[v].row_lid_to_gid.size()]);
966  }
967 
968  private:
969  FEEvaluation<dim,
970  fe_degree,
971  n_q_points_1d,
972  n_components,
973  Number,
974  VectorizedArrayType> &phi;
975 
976  unsigned int i;
977 
978  std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
979 
980  // local storage: buffer so that we access the global vector once
981  // note: may be larger then dofs_per_cell in the presence of
982  // constraints!
983  std::array<std::vector<Number>, n_lanes> diagonals_local_constrained;
984 
985  std::map<
987  std::vector<std::tuple<unsigned int, unsigned int, Number>>>
988  locally_relevant_constrains_hn_map;
989  };
990 
991  } // namespace internal
992 
993  template <int dim,
994  int fe_degree,
995  int n_q_points_1d,
996  int n_components,
997  typename Number,
998  typename VectorizedArrayType,
999  typename VectorType>
1000  void
1003  VectorType & diagonal_global,
1004  const std::function<void(FEEvaluation<dim,
1005  fe_degree,
1006  n_q_points_1d,
1007  n_components,
1008  Number,
1009  VectorizedArrayType> &)> &local_vmult,
1010  const unsigned int dof_no,
1011  const unsigned int quad_no,
1012  const unsigned int first_selected_component)
1013  {
1014  int dummy = 0;
1015 
1016  std::array<typename ::internal::BlockVectorSelector<
1017  VectorType,
1018  IsBlockVector<VectorType>::value>::BaseVectorType *,
1019  n_components>
1020  diagonal_global_components;
1021 
1022  for (unsigned int d = 0; d < n_components; ++d)
1023  diagonal_global_components[d] = ::internal::
1024  BlockVectorSelector<VectorType, IsBlockVector<VectorType>::value>::
1025  get_vector_component(diagonal_global, d + first_selected_component);
1026 
1027  const auto &dof_info = matrix_free.get_dof_info(dof_no);
1028 
1029  if (dof_info.start_components.back() == 1)
1030  for (unsigned int comp = 0; comp < n_components; ++comp)
1031  {
1032  Assert(diagonal_global_components[comp] != nullptr,
1033  ExcMessage("The finite element underlying this FEEvaluation "
1034  "object is scalar, but you requested " +
1035  std::to_string(n_components) +
1036  " components via the template argument in "
1037  "FEEvaluation. In that case, you must pass an "
1038  "std::vector<VectorType> or a BlockVector to " +
1039  "read_dof_values and distribute_local_to_global."));
1041  *diagonal_global_components[comp], matrix_free, dof_info);
1042  }
1043  else
1044  {
1046  *diagonal_global_components[0], matrix_free, dof_info);
1047  }
1048 
1049  matrix_free.template cell_loop<VectorType, int>(
1050  [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1051  VectorType &,
1052  const int &,
1053  const std::pair<unsigned int, unsigned int> &range) mutable {
1054  FEEvaluation<dim,
1055  fe_degree,
1056  n_q_points_1d,
1057  n_components,
1058  Number,
1059  VectorizedArrayType>
1060  phi(matrix_free, range, dof_no, quad_no, first_selected_component);
1061 
1062  internal::ComputeDiagonalHelper<dim,
1063  fe_degree,
1064  n_q_points_1d,
1065  n_components,
1066  Number,
1067  VectorizedArrayType>
1068  helper(phi);
1069 
1070  for (unsigned int cell = range.first; cell < range.second; ++cell)
1071  {
1072  helper.reinit(cell);
1073 
1074  for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
1075  {
1076  helper.prepare_basis_vector(i);
1077  local_vmult(phi);
1078  helper.submit();
1079  }
1080 
1081  helper.distribute_local_to_global(diagonal_global_components);
1082  }
1083  },
1084  diagonal_global,
1085  dummy,
1086  false);
1087  }
1088 
1089  template <typename CLASS,
1090  int dim,
1091  int fe_degree,
1092  int n_q_points_1d,
1093  int n_components,
1094  typename Number,
1095  typename VectorizedArrayType,
1096  typename VectorType>
1097  void
1100  VectorType & diagonal_global,
1101  void (CLASS::*cell_operation)(FEEvaluation<dim,
1102  fe_degree,
1103  n_q_points_1d,
1104  n_components,
1105  Number,
1106  VectorizedArrayType> &) const,
1107  const CLASS * owning_class,
1108  const unsigned int dof_no,
1109  const unsigned int quad_no,
1110  const unsigned int first_selected_component)
1111  {
1112  compute_diagonal<dim,
1113  fe_degree,
1114  n_q_points_1d,
1115  n_components,
1116  Number,
1117  VectorizedArrayType,
1118  VectorType>(
1119  matrix_free,
1120  diagonal_global,
1121  [&](auto &feeval) { (owning_class->*cell_operation)(feeval); },
1122  dof_no,
1123  quad_no,
1124  first_selected_component);
1125  }
1126 
1127  namespace internal
1128  {
1133  template <typename MatrixType,
1134  typename Number,
1135  std::enable_if_t<std::is_same<
1136  typename std::remove_const<typename std::remove_reference<
1137  typename MatrixType::value_type>::type>::type,
1138  typename std::remove_const<typename std::remove_reference<
1139  Number>::type>::type>::value> * = nullptr>
1141  create_new_affine_constraints_if_needed(
1142  const MatrixType &,
1143  const AffineConstraints<Number> &constraints,
1145  {
1146  return constraints;
1147  }
1148 
1154  template <typename MatrixType,
1155  typename Number,
1156  std::enable_if_t<!std::is_same<
1157  typename std::remove_const<typename std::remove_reference<
1158  typename MatrixType::value_type>::type>::type,
1159  typename std::remove_const<typename std::remove_reference<
1160  Number>::type>::type>::value> * = nullptr>
1162  create_new_affine_constraints_if_needed(
1163  const MatrixType &,
1164  const AffineConstraints<Number> &constraints,
1166  &new_constraints)
1167  {
1168  new_constraints =
1169  std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
1170  new_constraints->copy_from(constraints);
1171 
1172  return *new_constraints;
1173  }
1174  } // namespace internal
1175 
1176  template <int dim,
1177  int fe_degree,
1178  int n_q_points_1d,
1179  int n_components,
1180  typename Number,
1181  typename VectorizedArrayType,
1182  typename MatrixType>
1183  void
1186  const AffineConstraints<Number> & constraints_in,
1187  MatrixType & matrix,
1188  const std::function<void(FEEvaluation<dim,
1189  fe_degree,
1190  n_q_points_1d,
1191  n_components,
1192  Number,
1193  VectorizedArrayType> &)> &local_vmult,
1194  const unsigned int dof_no,
1195  const unsigned int quad_no,
1196  const unsigned int first_selected_component)
1197  {
1198  std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
1199  constraints_for_matrix;
1201  internal::create_new_affine_constraints_if_needed(matrix,
1202  constraints_in,
1203  constraints_for_matrix);
1204 
1205  matrix_free.template cell_loop<MatrixType, MatrixType>(
1206  [&](const auto &, auto &dst, const auto &, const auto range) {
1207  FEEvaluation<dim,
1208  fe_degree,
1209  n_q_points_1d,
1210  n_components,
1211  Number,
1212  VectorizedArrayType>
1213  integrator(
1214  matrix_free, range, dof_no, quad_no, first_selected_component);
1215 
1216  const unsigned int dofs_per_cell = integrator.dofs_per_cell;
1217 
1218  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
1219  std::vector<types::global_dof_index> dof_indices_mf(dofs_per_cell);
1220 
1221  std::array<FullMatrix<typename MatrixType::value_type>,
1222  VectorizedArrayType::size()>
1223  matrices;
1224 
1225  std::fill_n(matrices.begin(),
1226  VectorizedArrayType::size(),
1228  dofs_per_cell));
1229 
1230  const auto lexicographic_numbering =
1231  matrix_free
1232  .get_shape_info(dof_no,
1233  quad_no,
1234  first_selected_component,
1235  integrator.get_active_fe_index(),
1236  integrator.get_active_quadrature_index())
1237  .lexicographic_numbering;
1238 
1239  for (auto cell = range.first; cell < range.second; ++cell)
1240  {
1241  integrator.reinit(cell);
1242 
1243  const unsigned int n_filled_lanes =
1244  matrix_free.n_active_entries_per_cell_batch(cell);
1245 
1246  for (unsigned int v = 0; v < n_filled_lanes; ++v)
1247  matrices[v] = 0.0;
1248 
1249  for (unsigned int j = 0; j < dofs_per_cell; ++j)
1250  {
1251  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1252  integrator.begin_dof_values()[i] =
1253  static_cast<Number>(i == j);
1254 
1255  local_vmult(integrator);
1256 
1257  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1258  for (unsigned int v = 0; v < n_filled_lanes; ++v)
1259  matrices[v](i, j) = integrator.begin_dof_values()[i][v];
1260  }
1261 
1262  for (unsigned int v = 0; v < n_filled_lanes; ++v)
1263  {
1264  const auto cell_v =
1265  matrix_free.get_cell_iterator(cell, v, dof_no);
1266 
1267  if (matrix_free.get_mg_level() != numbers::invalid_unsigned_int)
1268  cell_v->get_mg_dof_indices(dof_indices);
1269  else
1270  cell_v->get_dof_indices(dof_indices);
1271 
1272  for (unsigned int j = 0; j < dof_indices.size(); ++j)
1273  dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]];
1274 
1275  constraints.distribute_local_to_global(matrices[v],
1276  dof_indices_mf,
1277  dst);
1278  }
1279  }
1280  },
1281  matrix,
1282  matrix);
1283 
1284  matrix.compress(VectorOperation::add);
1285  }
1286 
1287  template <typename CLASS,
1288  int dim,
1289  int fe_degree,
1290  int n_q_points_1d,
1291  int n_components,
1292  typename Number,
1293  typename VectorizedArrayType,
1294  typename MatrixType>
1295  void
1298  const AffineConstraints<Number> & constraints,
1299  MatrixType & matrix,
1300  void (CLASS::*cell_operation)(FEEvaluation<dim,
1301  fe_degree,
1302  n_q_points_1d,
1303  n_components,
1304  Number,
1305  VectorizedArrayType> &) const,
1306  const CLASS * owning_class,
1307  const unsigned int dof_no,
1308  const unsigned int quad_no,
1309  const unsigned int first_selected_component)
1310  {
1311  compute_matrix<dim,
1312  fe_degree,
1313  n_q_points_1d,
1314  n_components,
1315  Number,
1316  VectorizedArrayType,
1317  MatrixType>(
1318  matrix_free,
1319  constraints,
1320  matrix,
1321  [&](auto &feeval) { (owning_class->*cell_operation)(feeval); },
1322  dof_no,
1323  quad_no,
1324  first_selected_component);
1325  }
1326 
1327 #endif // DOXYGEN
1328 
1329 } // namespace MatrixFreeTools
1330 
1332 
1333 
1334 #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)
SmartPointer< const MatrixFree< dim, Number, VectorizedArrayType > > matrix_free
Definition: tools.h:334
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, VectorTypeOut &, const VectorTypeIn &, const std::pair< unsigned int, unsigned int >)> &cell_operation, VectorTypeOut &dst, const VectorTypeIn &src, const bool zero_dst_vector=false)
Definition: tools.h:233
void reinit(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const AdditionalData &additional_data=AdditionalData())
Definition: tools.h:204
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, VectorTypeOut &, const VectorTypeIn &, const std::pair< unsigned int, unsigned int >)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, VectorTypeOut &, const VectorTypeIn &, const std::pair< unsigned int, unsigned int >)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, VectorTypeOut &, const VectorTypeIn &, const std::pair< unsigned int, unsigned int >, const bool)> &boundary_operation, VectorTypeOut &dst, const VectorTypeIn &src, const bool zero_dst_vector=false)
Definition: tools.h:267
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:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
Point< 2 > second
Definition: grid_out.cc:4606
Point< 2 > first
Definition: grid_out.cc:4605
unsigned int level
Definition: grid_out.cc:4608
const unsigned int v0
Definition: grid_tools.cc:1048
const unsigned int v1
Definition: grid_tools.cc:1048
#define Assert(cond, exc)
Definition: exceptions.h:1583
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1756
#define AssertIndexRange(index, range)
Definition: exceptions.h:1821
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)
std::string to_string(const T &t)
Definition: patterns.h:2394
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:999
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:206
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