deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20: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\}}\)
Loading...
Searching...
No Matches
tools.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_matrix_free_tools_h
16#define dealii_matrix_free_tools_h
17
18#include <deal.II/base/config.h>
19
20#include <deal.II/grid/tria.h>
21
27
28
30
36{
37 namespace internal
38 {
39 template <int dim, typename Number, bool is_face_>
41 {
42 public:
44
45 std::vector<unsigned int> dof_numbers;
46 std::vector<unsigned int> quad_numbers;
47 std::vector<unsigned int> n_components;
48 std::vector<unsigned int> first_selected_components;
49 std::vector<unsigned int> batch_type;
50 static const bool is_face = is_face_;
51
52 std::function<std::vector<std::unique_ptr<FEEvalType>>(
53 const std::pair<unsigned int, unsigned int> &)>
55 std::function<void(std::vector<std::unique_ptr<FEEvalType>> &,
56 const unsigned int)>
58 std::function<void(std::vector<std::unique_ptr<FEEvalType>> &)>
60 };
61 } // namespace internal
62
68 template <int dim, typename AdditionalData>
69 void
71 AdditionalData &additional_data);
72
73
74
86 template <int dim,
87 int fe_degree,
88 int n_q_points_1d,
89 int n_components,
90 typename Number,
91 typename VectorizedArrayType,
92 typename VectorType>
93 void
96 VectorType &diagonal_global,
97 const std::function<void(FEEvaluation<dim,
98 fe_degree,
99 n_q_points_1d,
100 n_components,
101 Number,
102 VectorizedArrayType> &)>
104 const unsigned int dof_no = 0,
105 const unsigned int quad_no = 0,
106 const unsigned int first_selected_component = 0,
107 const unsigned int first_vector_component = 0);
108
112 template <int dim,
113 int fe_degree,
114 int n_q_points_1d,
115 int n_components,
116 typename Number,
117 typename MemorySpace,
118 typename QuadOperation>
119 void
121 const Portable::MatrixFree<dim, Number> &matrix_free,
124 EvaluationFlags::EvaluationFlags evaluation_flags,
126 const unsigned int dof_no = 0,
127 const unsigned int quad_no = 0,
128 const unsigned int first_selected_component = 0,
129 const unsigned int first_vector_component = 0);
130
134 template <typename CLASS,
135 int dim,
136 int fe_degree,
137 int n_q_points_1d,
138 int n_components,
139 typename Number,
140 typename VectorizedArrayType,
141 typename VectorType>
142 void
145 VectorType &diagonal_global,
146 void (CLASS::*cell_operation)(FEEvaluation<dim,
147 fe_degree,
148 n_q_points_1d,
149 n_components,
150 Number,
151 VectorizedArrayType> &) const,
153 const unsigned int dof_no = 0,
154 const unsigned int quad_no = 0,
155 const unsigned int first_selected_component = 0,
156 const unsigned int first_vector_component = 0);
157
158
159
170 template <int dim,
171 int fe_degree,
172 int n_q_points_1d,
173 int n_components,
174 typename Number,
175 typename VectorizedArrayType,
176 typename VectorType>
177 void
180 VectorType &diagonal_global,
181 const std::function<void(FEEvaluation<dim,
182 fe_degree,
183 n_q_points_1d,
184 n_components,
185 Number,
186 VectorizedArrayType> &)>
188 const std::function<void(FEFaceEvaluation<dim,
189 fe_degree,
190 n_q_points_1d,
191 n_components,
192 Number,
193 VectorizedArrayType> &,
195 fe_degree,
196 n_q_points_1d,
197 n_components,
198 Number,
199 VectorizedArrayType> &)>
201 const std::function<void(FEFaceEvaluation<dim,
202 fe_degree,
203 n_q_points_1d,
204 n_components,
205 Number,
206 VectorizedArrayType> &)>
208 const unsigned int dof_no = 0,
209 const unsigned int quad_no = 0,
210 const unsigned int first_selected_component = 0,
211 const unsigned int first_vector_component = 0);
212
213
214
218 template <typename CLASS,
219 int dim,
220 int fe_degree,
221 int n_q_points_1d,
222 int n_components,
223 typename Number,
224 typename VectorizedArrayType,
225 typename VectorType>
226 void
229 VectorType &diagonal_global,
230 void (CLASS::*cell_operation)(FEEvaluation<dim,
231 fe_degree,
232 n_q_points_1d,
233 n_components,
234 Number,
235 VectorizedArrayType> &) const,
237 fe_degree,
238 n_q_points_1d,
239 n_components,
240 Number,
241 VectorizedArrayType> &,
243 fe_degree,
244 n_q_points_1d,
245 n_components,
246 Number,
247 VectorizedArrayType> &)
248 const,
250 fe_degree,
251 n_q_points_1d,
252 n_components,
253 Number,
254 VectorizedArrayType> &)
255 const,
257 const unsigned int dof_no = 0,
258 const unsigned int quad_no = 0,
259 const unsigned int first_selected_component = 0,
260 const unsigned int first_vector_component = 0);
261
262
263
272 template <int dim,
273 int fe_degree,
274 int n_q_points_1d,
275 int n_components,
276 typename Number,
277 typename VectorizedArrayType,
278 typename MatrixType>
279 void
282 const AffineConstraints<Number> &constraints,
283 MatrixType &matrix,
284 const std::function<void(FEEvaluation<dim,
285 fe_degree,
286 n_q_points_1d,
287 n_components,
288 Number,
289 VectorizedArrayType> &)>
291 const unsigned int dof_no = 0,
292 const unsigned int quad_no = 0,
293 const unsigned int first_selected_component = 0);
294
295
296
300 template <typename CLASS,
301 int dim,
302 int fe_degree,
303 int n_q_points_1d,
304 int n_components,
305 typename Number,
306 typename VectorizedArrayType,
307 typename MatrixType>
308 void
311 const AffineConstraints<Number> &constraints,
312 MatrixType &matrix,
313 void (CLASS::*cell_operation)(FEEvaluation<dim,
314 fe_degree,
315 n_q_points_1d,
316 n_components,
317 Number,
318 VectorizedArrayType> &) const,
320 const unsigned int dof_no = 0,
321 const unsigned int quad_no = 0,
322 const unsigned int first_selected_component = 0);
323
324
325 namespace internal
326 {
331 template <int dim,
332 typename Number,
333 typename VectorizedArrayType,
334 typename VectorType,
335 typename VectorType2>
336 void
340 &data_cell,
342 &data_face,
345 VectorType &diagonal_global,
346 std::vector<VectorType2 *> &diagonal_global_components);
347
352 template <int dim,
353 typename Number,
354 typename VectorizedArrayType,
355 typename MatrixType>
356 void
359 const AffineConstraints<Number> &constraints,
366 MatrixType &matrix);
367 } // namespace internal
368
369
370
380 template <int dim,
381 int fe_degree,
382 int n_q_points_1d,
383 int n_components,
384 typename Number,
385 typename VectorizedArrayType,
386 typename MatrixType>
387 void
390 const AffineConstraints<Number> &constraints,
391 MatrixType &matrix,
392 const std::function<void(FEEvaluation<dim,
393 fe_degree,
394 n_q_points_1d,
395 n_components,
396 Number,
397 VectorizedArrayType> &)>
399 const std::function<void(FEFaceEvaluation<dim,
400 fe_degree,
401 n_q_points_1d,
402 n_components,
403 Number,
404 VectorizedArrayType> &,
406 fe_degree,
407 n_q_points_1d,
408 n_components,
409 Number,
410 VectorizedArrayType> &)>
412 const std::function<void(FEFaceEvaluation<dim,
413 fe_degree,
414 n_q_points_1d,
415 n_components,
416 Number,
417 VectorizedArrayType> &)>
419 const unsigned int dof_no = 0,
420 const unsigned int quad_no = 0,
421 const unsigned int first_selected_component = 0);
422
423
424
428 template <typename CLASS,
429 int dim,
430 int fe_degree,
431 int n_q_points_1d,
432 int n_components,
433 typename Number,
434 typename VectorizedArrayType,
435 typename MatrixType>
436 void
439 const AffineConstraints<Number> &constraints,
440 MatrixType &matrix,
441 void (CLASS::*cell_operation)(FEEvaluation<dim,
442 fe_degree,
443 n_q_points_1d,
444 n_components,
445 Number,
446 VectorizedArrayType> &) const,
448 fe_degree,
449 n_q_points_1d,
450 n_components,
451 Number,
452 VectorizedArrayType> &,
454 fe_degree,
455 n_q_points_1d,
456 n_components,
457 Number,
458 VectorizedArrayType> &)
459 const,
461 fe_degree,
462 n_q_points_1d,
463 n_components,
464 Number,
465 VectorizedArrayType> &)
466 const,
468 const unsigned int dof_no = 0,
469 const unsigned int quad_no = 0,
470 const unsigned int first_selected_component = 0);
471
472
473
482 template <int dim,
483 typename Number,
484 typename VectorizedArrayType = VectorizedArray<Number>>
486 {
487 public:
493 {
497 AdditionalData(const unsigned int dof_index = 0)
499 {}
500
504 unsigned int dof_index;
505 };
506
516 void
518 const AdditionalData &additional_data = AdditionalData())
519 {
520 this->matrix_free = &matrix_free;
521
522 std::vector<unsigned int> valid_fe_indices;
523
524 const auto &fe_collection =
525 matrix_free.get_dof_handler(additional_data.dof_index)
526 .get_fe_collection();
527
528 for (unsigned int i = 0; i < fe_collection.size(); ++i)
529 if (fe_collection[i].n_dofs_per_cell() > 0)
530 valid_fe_indices.push_back(i);
531
532 // TODO: relax this so that arbitrary number of valid
533 // FEs can be accepted
535
537 }
538
544 template <typename VectorTypeOut, typename VectorTypeIn>
545 void
546 cell_loop(const std::function<void(
549 const VectorTypeIn &,
550 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
551 VectorTypeOut &dst,
552 const VectorTypeIn &src,
553 const bool zero_dst_vector = false) const
554 {
555 const auto ebd_cell_operation = [&](const auto &matrix_free,
556 auto &dst,
557 const auto &src,
558 const auto &range) {
559 const auto category = matrix_free.get_cell_range_category(range);
560
562 return;
563
564 cell_operation(matrix_free, dst, src, range);
565 };
566
569 }
570
578 template <typename VectorTypeOut, typename VectorTypeIn>
579 void
580 loop(const std::function<
583 const VectorTypeIn &,
584 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
585 const std::function<
588 const VectorTypeIn &,
589 const std::pair<unsigned int, unsigned int> &)> &face_operation,
590 const std::function<
593 const VectorTypeIn &,
594 const std::pair<unsigned int, unsigned int> &,
595 const bool)> &boundary_operation,
596 VectorTypeOut &dst,
597 const VectorTypeIn &src,
598 const bool zero_dst_vector = false) const
599 {
600 const auto ebd_cell_operation = [&](const auto &matrix_free,
601 auto &dst,
602 const auto &src,
603 const auto &range) {
604 const auto category = matrix_free.get_cell_range_category(range);
605
607 return;
608
609 cell_operation(matrix_free, dst, src, range);
610 };
611
613 [&](const auto &matrix_free,
614 auto &dst,
615 const auto &src,
616 const auto &range) {
617 const auto category = matrix_free.get_face_range_category(range);
618
619 const unsigned int type =
620 static_cast<unsigned int>(category.first == fe_index_valid) +
621 static_cast<unsigned int>(category.second == fe_index_valid);
622
623 if (type == 0)
624 return; // deactivated face -> nothing to do
625
626 if (type == 1) // boundary face
628 matrix_free, dst, src, range, category.first == fe_index_valid);
629 else if (type == 2) // internal face
630 face_operation(matrix_free, dst, src, range);
631 };
632
637 dst,
638 src,
640 }
641
642 private:
648
652 unsigned int fe_index_valid;
653 };
654
655 // implementations
656
657#ifndef DOXYGEN
658
659 template <int dim, typename AdditionalData>
660 void
662 AdditionalData &additional_data)
663 {
664 // ... determine if we are on an active or a multigrid level
665 const unsigned int level = additional_data.mg_level;
666 const bool is_mg = (level != numbers::invalid_unsigned_int);
667
668 // ... create empty list for the category of each cell
669 if (is_mg)
670 additional_data.cell_vectorization_category.assign(
671 std::distance(tria.begin(level), tria.end(level)), 0);
672 else
673 additional_data.cell_vectorization_category.assign(tria.n_active_cells(),
674 0);
675
676 // ... set up scaling factor
677 std::vector<unsigned int> factors(GeometryInfo<dim>::faces_per_cell);
678
679 auto bids = tria.get_boundary_ids();
680 std::sort(bids.begin(), bids.end());
681
682 {
683 unsigned int n_bids = bids.size() + 1;
684 int offset = 1;
685 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
686 i++, offset = offset * n_bids)
687 factors[i] = offset;
688 }
689
690 const auto to_category = [&](const auto &cell) {
691 unsigned int c_num = 0;
692 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
693 {
694 const auto face = cell->face(i);
695 if (face->at_boundary() && !cell->has_periodic_neighbor(i))
696 c_num +=
697 factors[i] * (1 + std::distance(bids.begin(),
698 std::find(bids.begin(),
699 bids.end(),
700 face->boundary_id())));
701 }
702 return c_num;
703 };
704
705 if (!is_mg)
706 {
707 for (auto cell = tria.begin_active(); cell != tria.end(); ++cell)
708 {
709 if (cell->is_locally_owned())
710 additional_data
711 .cell_vectorization_category[cell->active_cell_index()] =
712 to_category(cell);
713 }
714 }
715 else
716 {
717 for (auto cell = tria.begin(level); cell != tria.end(level); ++cell)
718 {
719 if (cell->is_locally_owned_on_level())
720 additional_data.cell_vectorization_category[cell->index()] =
721 to_category(cell);
722 }
723 }
724
725 // ... finalize set up of matrix_free
726 additional_data.hold_all_faces_to_owned_cells = true;
727 additional_data.cell_vectorization_categories_strict = true;
728 additional_data.mapping_update_flags_faces_by_cells =
729 additional_data.mapping_update_flags_inner_faces |
730 additional_data.mapping_update_flags_boundary_faces;
731 }
732
733 namespace internal
734 {
735 template <typename Number>
736 struct LocalCSR
737 {
738 std::vector<unsigned int> row_lid_to_gid;
739 std::vector<unsigned int> row;
740 std::vector<unsigned int> col;
741 std::vector<Number> val;
742
743 std::vector<unsigned int> inverse_lookup_rows;
744 std::vector<std::pair<unsigned int, unsigned int>> inverse_lookup_origins;
745 };
746
747 template <int dim, typename VectorizedArrayType, bool is_face>
749 {
750 public:
751 using FEEvaluationType =
753
754 using Number = typename VectorizedArrayType::value_type;
755 static const unsigned int n_lanes = VectorizedArrayType::size();
756
758 : phi(nullptr)
759 , matrix_free(nullptr)
760 , dofs_per_component(0)
761 , n_components(0)
762 {}
763
765 : phi(nullptr)
766 , matrix_free(nullptr)
767 , dofs_per_component(0)
768 , n_components(0)
769 {}
770
771 void
772 initialize(
775 const unsigned int n_components)
776 {
777 // if we are in hp mode and the number of unknowns changed, we must
778 // clear the map of entries
779 if (dofs_per_component !=
780 phi.get_shape_info().dofs_per_component_on_cell)
781 {
783 dofs_per_component =
784 phi.get_shape_info().dofs_per_component_on_cell;
785 }
786 this->n_components = n_components;
787 this->dofs_per_cell = n_components * dofs_per_component;
788 this->phi = &phi;
789 this->matrix_free = &matrix_free;
790 }
791
792 void
793 reinit(const unsigned int cell)
794 {
795 // STEP 1: get relevant information from FEEvaluation
796 const auto &matrix_free = *this->matrix_free;
797 const auto &dof_info = phi->get_dof_info();
798 const unsigned int n_fe_components = dof_info.start_components.back();
799
800 // if we have a block vector with components with the same DoFHandler,
801 // each component is described with same set of constraints and
802 // we consider the shift in components only during access of the global
803 // vector
804 const unsigned int first_selected_component =
805 n_fe_components == 1 ? 0 : phi->get_first_selected_component();
806
807 this->n_lanes_filled =
808 is_face ? matrix_free.n_active_entries_per_face_batch(cell) :
809 matrix_free.n_active_entries_per_cell_batch(cell);
810
811 // STEP 2: setup CSR storage of transposed locally-relevant
812 // constraint matrix
813
814 // (constrained local index, global index of dof
815 // constraints, weight)
816 std::vector<std::tuple<unsigned int, unsigned int, Number>>
818 locally_relevant_constraints.reserve(dofs_per_cell);
819 std::vector<unsigned int> constraint_position;
820 std::vector<unsigned char> is_constrained_hn;
821
822 const std::array<unsigned int, n_lanes> &cells =
823 this->phi->get_cell_ids();
824
825 for (unsigned int v = 0; v < n_lanes_filled; ++v)
826 {
829
830 const unsigned int *dof_indices;
832
833 const unsigned int start =
834 cells[v] * n_fe_components + first_selected_component;
835 dof_indices =
836 dof_info.dof_indices.data() + dof_info.row_starts[start].first;
837 index_indicators = dof_info.row_starts[start].second;
838 next_index_indicators = dof_info.row_starts[start + 1].second;
839
840 // STEP 2a: setup locally-relevant constraint matrix in a
841 // coordinate list (COO)
843
844 if (n_components == 1 || n_fe_components == 1)
845 {
846 unsigned int ind_local = 0;
849 {
850 const std::pair<unsigned short, unsigned short> indicator =
851 dof_info.constraint_indicator[index_indicators];
852
853 for (unsigned int j = 0; j < indicator.first;
854 ++j, ++ind_local)
856 dof_indices[j],
857 1.0);
858
859 dof_indices += indicator.first;
860
861 const Number *data_val =
862 matrix_free.constraint_pool_begin(indicator.second);
863 const Number *end_pool =
864 matrix_free.constraint_pool_end(indicator.second);
865
866 for (; data_val != end_pool; ++data_val, ++dof_indices)
868 *dof_indices,
869 *data_val);
870 }
871
872 AssertIndexRange(ind_local, dofs_per_component + 1);
873
874 for (; ind_local < dofs_per_component;
875 ++dof_indices, ++ind_local)
877 *dof_indices,
878 1.0);
879 }
880 else
881 {
882 // case with vector-valued finite elements where all
883 // components are included in one single vector. Assumption:
884 // first come all entries to the first component, then all
885 // entries to the second one, and so on. This is ensured by
886 // the way MatrixFree reads out the indices.
887 for (unsigned int comp = 0; comp < n_components; ++comp)
888 {
889 unsigned int ind_local = 0;
890
891 // check whether there is any constraint on the current
892 // cell
895 {
896 const std::pair<unsigned short, unsigned short>
897 indicator =
898 dof_info.constraint_indicator[index_indicators];
899
900 // run through values up to next constraint
901 for (unsigned int j = 0; j < indicator.first;
902 ++j, ++ind_local)
903 locally_relevant_constraints.emplace_back(
904 comp * dofs_per_component + ind_local,
905 dof_indices[j],
906 1.0);
907 dof_indices += indicator.first;
908
909 const Number *data_val =
910 matrix_free.constraint_pool_begin(indicator.second);
911 const Number *end_pool =
912 matrix_free.constraint_pool_end(indicator.second);
913
914 for (; data_val != end_pool; ++data_val, ++dof_indices)
915 locally_relevant_constraints.emplace_back(
916 comp * dofs_per_component + ind_local,
917 *dof_indices,
918 *data_val);
919 }
920
921 AssertIndexRange(ind_local, dofs_per_component + 1);
922
923 // get the dof values past the last constraint
924 for (; ind_local < dofs_per_component;
925 ++dof_indices, ++ind_local)
926 locally_relevant_constraints.emplace_back(
927 comp * dofs_per_component + ind_local,
928 *dof_indices,
929 1.0);
930
931 if (comp + 1 < n_components)
933 dof_info.row_starts[start + comp + 2].second;
934 }
935 }
936
937 // we only need partial sortedness for the algorithm below in that
938 // all entries for a particular row must be adjacent. this is
939 // ensured by the way we fill the field, but check it again
940 for (unsigned int i = 1; i < locally_relevant_constraints.size();
941 ++i)
942 Assert(std::get<0>(locally_relevant_constraints[i]) >=
943 std::get<0>(locally_relevant_constraints[i - 1]),
945
946 // STEP 2c: apply hanging-node constraints
947 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
948 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
949 dof_info.hanging_node_constraint_masks_comp
950 [phi->get_active_fe_index()][first_selected_component])
951 {
952 const auto mask =
953 dof_info.hanging_node_constraint_masks[cells[v]];
954
955 // cell has hanging nodes
956 if (mask != ::internal::MatrixFreeFunctions::
958 {
959 // check if hanging node internpolation matrix has been set
960 // up
964
967
969 if (locally_relevant_constraints_tmp.capacity() <
974
975 // combine with other constraints: to avoid binary
976 // searches, we first build a list of where constraints
977 // are pointing to, and then merge the two lists
978 constraint_position.assign(dofs_per_cell,
980 for (auto &a : locally_relevant_constraints)
981 if (constraint_position[std::get<0>(a)] ==
983 constraint_position[std::get<0>(a)] =
985 &a);
986 is_constrained_hn.assign(dofs_per_cell, false);
988 is_constrained_hn[std::get<0>(hn)] = 1;
989
990 // not constrained from hanging nodes
991 for (const auto &a : locally_relevant_constraints)
992 if (is_constrained_hn[std::get<0>(a)] == 0)
994
995 // dof is constrained by hanging nodes: build transitive
996 // closure
997 for (const auto &hn : locally_relevant_constraints_hn)
998 if (constraint_position[std::get<1>(hn)] !=
1000 {
1003 auto other = locally_relevant_constraints.begin() +
1004 constraint_position[std::get<1>(hn)];
1005 AssertDimension(std::get<0>(*other), std::get<1>(hn));
1006
1007 for (; other != locally_relevant_constraints.end() &&
1008 std::get<0>(*other) == std::get<1>(hn);
1009 ++other)
1011 std::get<0>(hn),
1012 std::get<1>(*other),
1013 std::get<2>(hn) * std::get<2>(*other));
1014 }
1015
1018 }
1019 }
1020
1021 // STEP 2d: transpose COO
1022 std::sort(locally_relevant_constraints.begin(),
1024 [](const auto &a, const auto &b) {
1025 if (std::get<1>(a) < std::get<1>(b))
1026 return true;
1027 return (std::get<1>(a) == std::get<1>(b)) &&
1028 (std::get<0>(a) < std::get<0>(b));
1029 });
1030
1031 // STEP 2e: translate COO to CRS
1032 auto &c_pool = c_pools[v];
1033 {
1034 c_pool.row_lid_to_gid.clear();
1035 c_pool.row.clear();
1036 c_pool.row.push_back(0);
1037 c_pool.col.clear();
1038 c_pool.val.clear();
1039
1040 if (locally_relevant_constraints.size() > 0)
1041 c_pool.row_lid_to_gid.emplace_back(
1042 std::get<1>(locally_relevant_constraints.front()));
1043 for (const auto &j : locally_relevant_constraints)
1044 {
1045 if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
1046 {
1047 c_pool.row_lid_to_gid.push_back(std::get<1>(j));
1048 c_pool.row.push_back(c_pool.val.size());
1049 }
1050
1051 c_pool.col.emplace_back(std::get<0>(j));
1052 c_pool.val.emplace_back(std::get<2>(j));
1053 }
1054
1055 if (c_pool.val.size() > 0)
1056 c_pool.row.push_back(c_pool.val.size());
1057
1058 c_pool.inverse_lookup_rows.clear();
1059 c_pool.inverse_lookup_rows.resize(1 + dofs_per_cell);
1060 for (const unsigned int i : c_pool.col)
1061 c_pool.inverse_lookup_rows[1 + i]++;
1062 // transform to offsets
1063 std::partial_sum(c_pool.inverse_lookup_rows.begin(),
1064 c_pool.inverse_lookup_rows.end(),
1065 c_pool.inverse_lookup_rows.begin());
1066 AssertDimension(c_pool.inverse_lookup_rows.back(),
1067 c_pool.col.size());
1068
1069 c_pool.inverse_lookup_origins.resize(c_pool.col.size());
1070 std::vector<unsigned int> inverse_lookup_count(dofs_per_cell);
1071 for (unsigned int row = 0; row < c_pool.row.size() - 1; ++row)
1072 for (unsigned int col = c_pool.row[row];
1073 col < c_pool.row[row + 1];
1074 ++col)
1075 {
1076 const unsigned int index = c_pool.col[col];
1077 c_pool.inverse_lookup_origins
1078 [c_pool.inverse_lookup_rows[index] +
1079 inverse_lookup_count[index]] = std::make_pair(row, col);
1081 }
1082 }
1083 }
1084
1085 // STEP 3: compute element matrix A_e, apply
1086 // locally-relevant constraints C_e^T * A_e * C_e, and get the
1087 // the diagonal entry
1088 // (C_e^T * A_e * C_e)(i,i)
1089 // or
1090 // C_e^T(i,:) * A_e * C_e(:,i).
1091 //
1092 // Since, we compute the element matrix column-by-column and as a
1093 // result never actually have the full element matrix, we actually
1094 // perform following steps:
1095 // 1) loop over all columns of the element matrix
1096 // a) compute column i
1097 // b) compute for each j (rows of C_e^T):
1098 // (C_e^T(j,:) * A_e(:,i)) * C_e(i,j)
1099 // or
1100 // (C_e^T(j,:) * A_e(:,i)) * C_e^T(j,i)
1101 // This gives a contribution the j-th entry of the
1102 // locally-relevant diagonal and comprises the multiplication
1103 // by the locally-relevant constraint matrix from the left and
1104 // the right. There is no contribution to the j-th vector
1105 // entry if the j-th row of C_e^T is empty or C_e^T(j,i) is
1106 // zero.
1107
1108 // set size locally-relevant diagonal
1109 for (unsigned int v = 0; v < n_lanes_filled; ++v)
1111 c_pools[v].row_lid_to_gid.size() *
1112 (n_fe_components == 1 ? n_components : 1),
1113 Number(0.0));
1114
1115 // check if fast path can be taken via FEEvaluation
1116 bool use_fast_path = true;
1117
1118 for (unsigned int v = 0; v < n_lanes_filled; ++v)
1119 {
1120 auto &c_pool = c_pools[v];
1121
1122 for (unsigned int i = 0; i < c_pool.row.size() - 1; ++i)
1123 {
1124 if ((c_pool.row[i + 1] - c_pool.row[i]) > 1)
1125 {
1126 use_fast_path = false;
1127 break;
1128 }
1129 else if (((c_pool.row[i + 1] - c_pool.row[i]) == 1) &&
1130 (c_pool.val[c_pool.row[i]] != 1.0))
1131 {
1132 use_fast_path = false;
1133 break;
1134 }
1135 }
1136
1137 if (use_fast_path == false)
1138 break;
1139 }
1140
1142 }
1143
1144 void
1146 const ::internal::MatrixFreeFunctions::compressed_constraint_kind
1147 mask)
1148 {
1150
1151 // assume that we constrain one face, i.e., (fe_degree + 1)^(dim-1)
1152 // unknowns - we might have more or less entries, but this is a good
1153 // first guess
1154 const unsigned int degree =
1155 phi->get_shape_info().data.front().fe_degree;
1156 constraints_hn.reserve(Utilities::pow(degree + 1, dim - 1));
1157
1158 // 1) collect hanging-node constraints for cell assuming
1159 // scalar finite element
1160 values_dofs.resize(dofs_per_component);
1161 std::array<
1163 VectorizedArrayType::size()>
1164 constraint_mask;
1165 constraint_mask.fill(::internal::MatrixFreeFunctions::
1167 constraint_mask[0] = mask;
1168
1169 for (unsigned int i = 0; i < dofs_per_component; ++i)
1170 {
1171 for (unsigned int j = 0; j < dofs_per_component; ++j)
1172 values_dofs[j] = VectorizedArrayType();
1173 values_dofs[i] = Number(1);
1174
1176 dim,
1177 Number,
1178 VectorizedArrayType>::apply(1,
1179 degree,
1180 phi->get_shape_info(),
1181 false,
1182 constraint_mask,
1183 values_dofs.data());
1184
1185 const Number tolerance =
1186 std::max(Number(1e-12),
1187 std::numeric_limits<Number>::epsilon() * 16);
1188 for (unsigned int j = 0; j < dofs_per_component; ++j)
1189 if (std::abs(values_dofs[j][0]) > tolerance &&
1190 (j != i ||
1191 std::abs(values_dofs[j][0] - Number(1)) > tolerance))
1192 constraints_hn.emplace_back(j, i, values_dofs[j][0]);
1193 }
1194
1195 // 2) extend for multiple components
1196 const unsigned int n_hn_constraints = constraints_hn.size();
1197 constraints_hn.resize(n_hn_constraints * n_components);
1198
1199 for (unsigned int c = 1; c < n_components; ++c)
1200 for (unsigned int i = 0; i < n_hn_constraints; ++i)
1201 constraints_hn[c * n_hn_constraints + i] = std::make_tuple(
1202 std::get<0>(constraints_hn[i]) + c * dofs_per_component,
1203 std::get<1>(constraints_hn[i]) + c * dofs_per_component,
1204 std::get<2>(constraints_hn[i]));
1205 }
1206
1207 void
1208 prepare_basis_vector(const unsigned int i)
1209 {
1210 this->i = i;
1211
1212 // compute i-th column of element stiffness matrix:
1213 // this could be simply performed as done at the moment with
1214 // matrix-free operator evaluation applied to a ith-basis vector
1215 VectorizedArrayType *dof_values = phi->begin_dof_values();
1216 for (unsigned int j = 0; j < dofs_per_cell; ++j)
1217 dof_values[j] = VectorizedArrayType();
1218 dof_values[i] = Number(1);
1219 }
1220
1221 void
1223 {
1224 VectorizedArrayType *dof_values = phi->begin_dof_values();
1225 for (unsigned int j = 0; j < dofs_per_cell; ++j)
1226 dof_values[j] = VectorizedArrayType();
1227 }
1228
1229 void
1230 submit()
1231 {
1232 // if we have a block vector with components with the same DoFHandler,
1233 // we need to figure out which component and which DoF within the
1234 // component are we currently considering
1235 const unsigned int n_fe_components =
1236 phi->get_dof_info().start_components.back();
1237 const unsigned int comp =
1238 n_fe_components == 1 ? i / dofs_per_component : 0;
1239 const unsigned int i_comp =
1240 n_fe_components == 1 ? (i % dofs_per_component) : i;
1241
1242 // apply local constraint matrix from left and from right:
1243 // loop over all rows of transposed constrained matrix
1244 for (unsigned int v = 0; v < n_lanes_filled; ++v)
1245 {
1246 const auto &c_pool = c_pools[v];
1247
1248 for (unsigned int jj = c_pool.inverse_lookup_rows[i_comp];
1249 jj < c_pool.inverse_lookup_rows[i_comp + 1];
1250 ++jj)
1251 {
1252 const unsigned int j = c_pool.inverse_lookup_origins[jj].first;
1253 // apply constraint matrix from the left
1254 Number temp = 0.0;
1255 for (unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
1256 temp += c_pool.val[k] *
1257 phi->begin_dof_values()[comp * dofs_per_component +
1258 c_pool.col[k]][v];
1259
1260 // apply constraint matrix from the right
1262 [v][j + comp * c_pools[v].row_lid_to_gid.size()] +=
1263 temp * c_pool.val[c_pool.inverse_lookup_origins[jj].second];
1264 }
1265 }
1266 }
1267
1268 template <typename VectorType>
1269 inline void
1270 distribute_local_to_global(std::vector<VectorType *> &diagonal_global)
1271 {
1272 // STEP 4: assembly results: add into global vector
1273 const unsigned int n_fe_components =
1274 phi->get_dof_info().start_components.back();
1275
1276 if (n_fe_components == 1)
1277 AssertDimension(diagonal_global.size(), n_components);
1278
1279 for (unsigned int v = 0; v < n_lanes_filled; ++v)
1280 // if we have a block vector with components with the same
1281 // DoFHandler, we need to loop over all components manually and
1282 // need to apply the correct shift
1283 for (unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
1284 for (unsigned int comp = 0;
1285 comp < (n_fe_components == 1 ?
1286 static_cast<unsigned int>(n_components) :
1287 1);
1288 ++comp)
1290 *diagonal_global[n_fe_components == 1 ? comp : 0],
1293 [v][j + comp * c_pools[v].row_lid_to_gid.size()]);
1294 }
1295
1296 bool
1298 {
1300 }
1301
1302 private:
1305
1306 unsigned int dofs_per_component;
1307 unsigned int dofs_per_cell;
1308 unsigned int n_components;
1309
1310 unsigned int i;
1311
1312 unsigned int n_lanes_filled;
1313
1314 std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
1315
1316 // local storage: buffer so that we access the global vector once
1317 // note: may be larger then dofs_per_cell in the presence of
1318 // constraints!
1319 std::array<std::vector<Number>, n_lanes> diagonals_local_constrained;
1320
1321 std::map<
1323 std::vector<std::tuple<unsigned int, unsigned int, Number>>>
1325
1326 // scratch array
1328
1330 };
1331
1332 template <bool is_face,
1333 int dim,
1334 typename Number,
1335 typename VectorizedArrayType>
1336 bool
1339 const std::pair<unsigned int, unsigned int> &range,
1340 const unsigned int dof_no,
1341 const unsigned int quad_no,
1342 const unsigned int first_selected_component,
1343 const unsigned int fe_degree,
1344 const unsigned int n_q_points_1d,
1345 const bool is_interior_face = true)
1346 {
1347 const unsigned int static_n_q_points =
1348 is_face ? Utilities::pow(n_q_points_1d, dim - 1) :
1349 Utilities::pow(n_q_points_1d, dim);
1350
1351 unsigned int active_fe_index = 0;
1352 if (!is_face)
1353 active_fe_index = matrix_free.get_cell_active_fe_index(range);
1354 else if (is_interior_face)
1355 active_fe_index = matrix_free.get_face_range_category(range).first;
1356 else
1357 active_fe_index = matrix_free.get_face_range_category(range).second;
1358
1359 const auto init_data = ::internal::
1360 extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
1361 matrix_free,
1362 dof_no,
1363 first_selected_component,
1364 quad_no,
1365 fe_degree,
1366 static_n_q_points,
1367 active_fe_index,
1368 numbers::invalid_unsigned_int /*active_quad_index*/,
1369 numbers::invalid_unsigned_int /*face_type*/);
1370
1371 return init_data.shape_info->dofs_per_component_on_cell == 0;
1372 }
1373
1374 } // namespace internal
1375
1376
1377 template <int dim,
1378 int fe_degree,
1379 int n_q_points_1d,
1380 int n_components,
1381 typename Number,
1382 typename QuadOperation>
1383 class CellAction
1384 {
1385 public:
1387 const EvaluationFlags::EvaluationFlags evaluation_flags,
1390 , m_evaluation_flags(evaluation_flags)
1392 {}
1393
1394 KOKKOS_FUNCTION void
1395 operator()(const unsigned int cell,
1398 const Number *,
1399 Number *dst) const
1400 {
1401 Portable::
1402 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>
1403 fe_eval(gpu_data, shared_data);
1404 m_quad_operation.set_matrix_free_data(*gpu_data);
1405 m_quad_operation.set_cell(cell);
1406 constexpr int dofs_per_cell = decltype(fe_eval)::tensor_dofs_per_cell;
1407 typename decltype(fe_eval)::value_type
1408 diagonal[dofs_per_cell / n_components] = {};
1409 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1410 {
1411 const auto c = i % n_components;
1412
1413 Kokkos::parallel_for(
1414 Kokkos::TeamThreadRange(shared_data->team_member,
1415 dofs_per_cell / n_components),
1416 [&](int j) {
1417 typename decltype(fe_eval)::value_type val = {};
1418
1419 if constexpr (n_components == 1)
1420 {
1421 val = (i == j) ? 1 : 0;
1422 }
1423 else
1424 {
1425 val[c] = (i / n_components == j) ? 1 : 0;
1426 }
1427
1428 fe_eval.submit_dof_value(val, j);
1429 });
1430
1431 shared_data->team_member.team_barrier();
1432
1433 Portable::internal::
1434 resolve_hanging_nodes<dim, fe_degree, false, Number>(
1435 shared_data->team_member,
1436 gpu_data->constraint_weights,
1437 gpu_data->constraint_mask(cell * n_components + c),
1438 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
1439
1440 fe_eval.evaluate(m_evaluation_flags);
1441 fe_eval.apply_for_each_quad_point(m_quad_operation);
1442 fe_eval.integrate(m_integration_flags);
1443
1444 Portable::internal::
1445 resolve_hanging_nodes<dim, fe_degree, true, Number>(
1446 shared_data->team_member,
1447 gpu_data->constraint_weights,
1448 gpu_data->constraint_mask(cell * n_components + c),
1449 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
1450
1451 Kokkos::single(Kokkos::PerTeam(shared_data->team_member), [&] {
1452 if constexpr (n_components == 1)
1453 diagonal[i] = fe_eval.get_dof_value(i);
1454 else
1455 diagonal[i / n_components][i % n_components] =
1456 fe_eval.get_dof_value(i / n_components)[i % n_components];
1457 });
1458
1459 shared_data->team_member.team_barrier();
1460 }
1461
1462 Kokkos::single(Kokkos::PerTeam(shared_data->team_member), [&] {
1463 for (unsigned int i = 0; i < dofs_per_cell / n_components; ++i)
1464 fe_eval.submit_dof_value(diagonal[i], i);
1465 });
1466
1467 shared_data->team_member.team_barrier();
1468
1469 // We need to do the same as distribute_local_to_global but without
1470 // constraints since we have already taken care of them earlier
1471 if (gpu_data->use_coloring)
1472 {
1473 Kokkos::parallel_for(
1474 Kokkos::TeamThreadRange(shared_data->team_member, dofs_per_cell),
1475 [&](const int &i) {
1476 dst[gpu_data->local_to_global(cell, i)] +=
1477 shared_data->values(i % (dofs_per_cell / n_components),
1478 i / (dofs_per_cell / n_components));
1479 });
1480 }
1481 else
1482 {
1483 Kokkos::parallel_for(
1484 Kokkos::TeamThreadRange(shared_data->team_member, dofs_per_cell),
1485 [&](const int &i) {
1486 Kokkos::atomic_add(
1487 &dst[gpu_data->local_to_global(cell, i)],
1488 shared_data->values(i % (dofs_per_cell / n_components),
1489 i / (dofs_per_cell / n_components)));
1490 });
1491 }
1492 };
1493
1494 static constexpr unsigned int n_local_dofs =
1495 ::Utilities::pow(fe_degree + 1, dim) * n_components;
1496 static constexpr unsigned int n_q_points =
1497 ::Utilities::pow(n_q_points_1d, dim);
1498
1499 private:
1503 };
1504
1505
1506 template <int dim,
1507 int fe_degree,
1508 int n_q_points_1d,
1509 int n_components,
1510 typename Number,
1511 typename MemorySpace,
1512 typename QuadOperation>
1513 void
1515 const Portable::MatrixFree<dim, Number> &matrix_free,
1518 EvaluationFlags::EvaluationFlags evaluation_flags,
1520 const unsigned int dof_no,
1521 const unsigned int quad_no,
1522 const unsigned int first_selected_component,
1523 const unsigned int first_vector_component)
1524 {
1526 Assert(quad_no == 0, ExcNotImplemented());
1527 Assert(first_selected_component == 0, ExcNotImplemented());
1528 Assert(first_vector_component == 0, ExcNotImplemented());
1529
1531
1532
1533 CellAction<dim,
1534 fe_degree,
1535 n_q_points_1d,
1536 n_components,
1537 Number,
1539 cell_action(quad_operation, evaluation_flags, integration_flags);
1541 matrix_free.cell_loop(cell_action, dummy, diagonal_global);
1542
1543 matrix_free.set_constrained_values(Number(1.), diagonal_global);
1544 }
1545
1546 template <int dim,
1547 int fe_degree,
1548 int n_q_points_1d,
1549 int n_components,
1550 typename Number,
1551 typename VectorizedArrayType,
1552 typename VectorType>
1553 void
1556 VectorType &diagonal_global,
1557 const std::function<void(FEEvaluation<dim,
1558 fe_degree,
1559 n_q_points_1d,
1560 n_components,
1561 Number,
1562 VectorizedArrayType> &)>
1564 const unsigned int dof_no,
1565 const unsigned int quad_no,
1566 const unsigned int first_selected_component,
1567 const unsigned int first_vector_component)
1568 {
1569 compute_diagonal<dim,
1570 fe_degree,
1571 n_q_points_1d,
1572 n_components,
1573 Number,
1574 VectorizedArrayType,
1575 VectorType>(matrix_free,
1578 {},
1579 {},
1580 dof_no,
1581 quad_no,
1582 first_selected_component,
1583 first_vector_component);
1584 }
1585
1586 template <typename CLASS,
1587 int dim,
1588 int fe_degree,
1589 int n_q_points_1d,
1590 int n_components,
1591 typename Number,
1592 typename VectorizedArrayType,
1593 typename VectorType>
1594 void
1597 VectorType &diagonal_global,
1598 void (CLASS::*cell_operation)(FEEvaluation<dim,
1599 fe_degree,
1600 n_q_points_1d,
1601 n_components,
1602 Number,
1603 VectorizedArrayType> &) const,
1605 const unsigned int dof_no,
1606 const unsigned int quad_no,
1607 const unsigned int first_selected_component,
1608 const unsigned int first_vector_component)
1609 {
1610 compute_diagonal<dim,
1611 fe_degree,
1612 n_q_points_1d,
1613 n_components,
1614 Number,
1615 VectorizedArrayType,
1616 VectorType>(
1617 matrix_free,
1619 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
1620 dof_no,
1621 quad_no,
1622 first_selected_component,
1623 first_vector_component);
1624 }
1625
1626 template <int dim,
1627 int fe_degree,
1628 int n_q_points_1d,
1629 int n_components,
1630 typename Number,
1631 typename VectorizedArrayType,
1632 typename VectorType>
1633 void
1636 VectorType &diagonal_global,
1637 const std::function<void(FEEvaluation<dim,
1638 fe_degree,
1639 n_q_points_1d,
1640 n_components,
1641 Number,
1642 VectorizedArrayType> &)>
1644 const std::function<void(FEFaceEvaluation<dim,
1645 fe_degree,
1646 n_q_points_1d,
1647 n_components,
1648 Number,
1649 VectorizedArrayType> &,
1650 FEFaceEvaluation<dim,
1651 fe_degree,
1652 n_q_points_1d,
1653 n_components,
1654 Number,
1655 VectorizedArrayType> &)>
1657 const std::function<void(FEFaceEvaluation<dim,
1658 fe_degree,
1659 n_q_points_1d,
1660 n_components,
1661 Number,
1662 VectorizedArrayType> &)>
1664 const unsigned int dof_no,
1665 const unsigned int quad_no,
1666 const unsigned int first_selected_component,
1667 const unsigned int first_vector_component)
1668 {
1669 std::vector<typename ::internal::BlockVectorSelector<
1670 VectorType,
1672 diagonal_global_components(n_components);
1673
1674 for (unsigned int d = 0; d < n_components; ++d)
1675 diagonal_global_components[d] = ::internal::
1677 get_vector_component(diagonal_global, d + first_vector_component);
1678
1679 const auto &dof_info = matrix_free.get_dof_info(dof_no);
1680
1681 if (dof_info.start_components.back() == 1)
1682 for (unsigned int comp = 0; comp < n_components; ++comp)
1683 {
1685 ExcMessage("The finite element underlying this FEEvaluation "
1686 "object is scalar, but you requested " +
1687 std::to_string(n_components) +
1688 " components via the template argument in "
1689 "FEEvaluation. In that case, you must pass an "
1690 "std::vector<VectorType> or a BlockVector to " +
1691 "read_dof_values and distribute_local_to_global."));
1693 *diagonal_global_components[comp], matrix_free, dof_info);
1694 }
1695 else
1696 {
1698 *diagonal_global_components[0], matrix_free, dof_info);
1699 }
1700
1701 using FEEvalType = FEEvaluation<dim,
1702 fe_degree,
1703 n_q_points_1d,
1704 n_components,
1705 Number,
1706 VectorizedArrayType>;
1707
1708 using FEFaceEvalType = FEFaceEvaluation<dim,
1709 fe_degree,
1710 n_q_points_1d,
1711 n_components,
1712 Number,
1713 VectorizedArrayType>;
1714
1715 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1716 data_cell;
1717
1718 data_cell.dof_numbers = {dof_no};
1719 data_cell.quad_numbers = {quad_no};
1720 data_cell.n_components = {n_components};
1721 data_cell.first_selected_components = {first_selected_component};
1722 data_cell.batch_type = {0};
1723
1724 data_cell.op_create =
1725 [&](const std::pair<unsigned int, unsigned int> &range) {
1726 std::vector<
1727 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
1728 phi;
1729
1730 if (!internal::is_fe_nothing<false>(matrix_free,
1731 range,
1732 dof_no,
1733 quad_no,
1734 first_selected_component,
1735 fe_degree,
1736 n_q_points_1d))
1737 phi.emplace_back(std::make_unique<FEEvalType>(
1738 matrix_free, range, dof_no, quad_no, first_selected_component));
1739
1740 return phi;
1741 };
1742
1743 data_cell.op_reinit = [](auto &phi, const unsigned batch) {
1744 if (phi.size() == 1)
1745 static_cast<FEEvalType &>(*phi[0]).reinit(batch);
1746 };
1747
1748 if (cell_operation)
1749 data_cell.op_compute = [&](auto &phi) {
1750 cell_operation(static_cast<FEEvalType &>(*phi[0]));
1751 };
1752
1753 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1754 data_face;
1755
1756 data_face.dof_numbers = {dof_no, dof_no};
1757 data_face.quad_numbers = {quad_no, quad_no};
1758 data_face.n_components = {n_components, n_components};
1759 data_face.first_selected_components = {first_selected_component,
1760 first_selected_component};
1761 data_face.batch_type = {1, 2};
1762
1763 data_face.op_create =
1764 [&](const std::pair<unsigned int, unsigned int> &range) {
1765 std::vector<
1766 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1767 phi;
1768
1769 if (!internal::is_fe_nothing<true>(matrix_free,
1770 range,
1771 dof_no,
1772 quad_no,
1773 first_selected_component,
1774 fe_degree,
1775 n_q_points_1d,
1776 true) &&
1777 !internal::is_fe_nothing<true>(matrix_free,
1778 range,
1779 dof_no,
1780 quad_no,
1781 first_selected_component,
1782 fe_degree,
1783 n_q_points_1d,
1784 false))
1785 {
1786 phi.emplace_back(
1787 std::make_unique<FEFaceEvalType>(matrix_free,
1788 range,
1789 true,
1790 dof_no,
1791 quad_no,
1792 first_selected_component));
1793 phi.emplace_back(
1794 std::make_unique<FEFaceEvalType>(matrix_free,
1795 range,
1796 false,
1797 dof_no,
1798 quad_no,
1799 first_selected_component));
1800 }
1801
1802 return phi;
1803 };
1804
1805 data_face.op_reinit = [](auto &phi, const unsigned batch) {
1806 if (phi.size() == 2)
1807 {
1808 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
1809 static_cast<FEFaceEvalType &>(*phi[1]).reinit(batch);
1810 }
1811 };
1812
1813 if (face_operation)
1814 data_face.op_compute = [&](auto &phi) {
1815 face_operation(static_cast<FEFaceEvalType &>(*phi[0]),
1816 static_cast<FEFaceEvalType &>(*phi[1]));
1817 };
1818
1819 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1821
1822 data_boundary.dof_numbers = {dof_no};
1823 data_boundary.quad_numbers = {quad_no};
1824 data_boundary.n_components = {n_components};
1825 data_boundary.first_selected_components = {first_selected_component};
1826 data_boundary.batch_type = {1};
1827
1829 .op_create = [&](const std::pair<unsigned int, unsigned int> &range) {
1830 std::vector<
1831 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1832 phi;
1833
1834 if (!internal::is_fe_nothing<true>(matrix_free,
1835 range,
1836 dof_no,
1837 quad_no,
1838 first_selected_component,
1839 fe_degree,
1840 n_q_points_1d,
1841 true))
1842 phi.emplace_back(std::make_unique<FEFaceEvalType>(
1843 matrix_free, range, true, dof_no, quad_no, first_selected_component));
1844
1845 return phi;
1846 };
1847
1848 data_boundary.op_reinit = [](auto &phi, const unsigned batch) {
1849 if (phi.size() == 1)
1850 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
1851 };
1852
1854 data_boundary.op_compute = [&](auto &phi) {
1855 boundary_operation(static_cast<FEFaceEvalType &>(*phi[0]));
1856 };
1857
1858 internal::compute_diagonal(matrix_free,
1859 data_cell,
1860 data_face,
1864 }
1865
1866 namespace internal
1867 {
1868 template <int dim,
1869 typename Number,
1870 typename VectorizedArrayType,
1871 typename VectorType,
1872 typename VectorType2>
1873 void
1876 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1877 &data_cell,
1878 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1879 &data_face,
1880 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1882 VectorType &diagonal_global,
1883 std::vector<VectorType2 *> &diagonal_global_components)
1884 {
1885 // TODO: can we remove diagonal_global_components as argument?
1886
1887 int dummy = 0;
1888
1889 using Helper =
1890 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, false>;
1891
1892 using HelperFace =
1893 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, true>;
1894
1899
1900 const auto batch_operation =
1901 [&](auto &data,
1902 auto &scratch_data,
1903 const std::pair<unsigned int, unsigned int> &range) {
1904 if (!data.op_compute)
1905 return; // nothing to do
1906
1907 auto phi = data.op_create(range);
1908
1909 const unsigned int n_blocks = phi.size();
1910
1911 auto &helpers = scratch_data.get();
1912 helpers.resize(n_blocks);
1913
1914 for (unsigned int b = 0; b < n_blocks; ++b)
1915 helpers[b].initialize(*phi[b], matrix_free, data.n_components[b]);
1916
1917 for (unsigned int batch = range.first; batch < range.second; ++batch)
1918 {
1919 data.op_reinit(phi, batch);
1920
1921 for (unsigned int b = 0; b < n_blocks; ++b)
1922 helpers[b].reinit(batch);
1923
1924 if (n_blocks > 1)
1925 {
1926 Assert(std::all_of(helpers.begin(),
1927 helpers.end(),
1928 [](const auto &helper) {
1929 return helper.has_simple_constraints();
1930 }),
1932 }
1933
1934 for (unsigned int b = 0; b < n_blocks; ++b)
1935 {
1936 for (unsigned int i = 0;
1937 i < phi[b]->get_shape_info().dofs_per_component_on_cell *
1938 data.n_components[b];
1939 ++i)
1940 {
1941 for (unsigned int bb = 0; bb < n_blocks; ++bb)
1942 if (b == bb)
1943 helpers[bb].prepare_basis_vector(i);
1944 else
1945 helpers[bb].zero_basis_vector();
1946
1947 data.op_compute(phi);
1948 helpers[b].submit();
1949 }
1950
1951 helpers[b].distribute_local_to_global(
1953 }
1954 }
1955 };
1956
1957 const auto cell_operation_wrapped =
1958 [&](const auto &, auto &, const auto &, const auto range) {
1959 batch_operation(data_cell, scratch_data, range);
1960 };
1961
1962 const auto face_operation_wrapped =
1963 [&](const auto &, auto &, const auto &, const auto range) {
1965 };
1966
1967 const auto boundary_operation_wrapped =
1968 [&](const auto &, auto &, const auto &, const auto range) {
1970 };
1971
1972 if (data_face.op_compute || data_boundary.op_compute)
1973 matrix_free.template loop<VectorType, int>(cell_operation_wrapped,
1977 dummy,
1978 false);
1979 else
1982 dummy,
1983 false);
1984 }
1985 } // namespace internal
1986
1987 template <typename CLASS,
1988 int dim,
1989 int fe_degree,
1990 int n_q_points_1d,
1991 int n_components,
1992 typename Number,
1993 typename VectorizedArrayType,
1994 typename VectorType>
1995 void
1998 VectorType &diagonal_global,
1999 void (CLASS::*cell_operation)(FEEvaluation<dim,
2000 fe_degree,
2001 n_q_points_1d,
2002 n_components,
2003 Number,
2004 VectorizedArrayType> &) const,
2006 fe_degree,
2007 n_q_points_1d,
2008 n_components,
2009 Number,
2010 VectorizedArrayType> &,
2011 FEFaceEvaluation<dim,
2012 fe_degree,
2013 n_q_points_1d,
2014 n_components,
2015 Number,
2016 VectorizedArrayType> &)
2017 const,
2019 fe_degree,
2020 n_q_points_1d,
2021 n_components,
2022 Number,
2023 VectorizedArrayType> &)
2024 const,
2026 const unsigned int dof_no,
2027 const unsigned int quad_no,
2028 const unsigned int first_selected_component,
2029 const unsigned int first_vector_component)
2030 {
2031 compute_diagonal<dim,
2032 fe_degree,
2033 n_q_points_1d,
2034 n_components,
2035 Number,
2036 VectorizedArrayType,
2037 VectorType>(
2038 matrix_free,
2040 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
2041 [&](auto &phi_m, auto &phi_p) {
2043 },
2044 [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
2045 dof_no,
2046 quad_no,
2047 first_selected_component,
2048 first_vector_component);
2049 }
2050
2051 namespace internal
2052 {
2057 template <
2058 typename MatrixType,
2059 typename Number,
2060 std::enable_if_t<std::is_same_v<
2061 std::remove_const_t<
2062 std::remove_reference_t<typename MatrixType::value_type>>,
2063 std::remove_const_t<std::remove_reference_t<Number>>>> * = nullptr>
2066 const MatrixType &,
2067 const AffineConstraints<Number> &constraints,
2069 {
2070 return constraints;
2071 }
2072
2078 template <
2079 typename MatrixType,
2080 typename Number,
2081 std::enable_if_t<!std::is_same_v<
2082 std::remove_const_t<
2083 std::remove_reference_t<typename MatrixType::value_type>>,
2084 std::remove_const_t<std::remove_reference_t<Number>>>> * = nullptr>
2087 const MatrixType &,
2088 const AffineConstraints<Number> &constraints,
2091 {
2093 std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
2094 new_constraints->copy_from(constraints);
2095
2096 return *new_constraints;
2097 }
2098 } // namespace internal
2099
2100 template <int dim,
2101 int fe_degree,
2102 int n_q_points_1d,
2103 int n_components,
2104 typename Number,
2105 typename VectorizedArrayType,
2106 typename MatrixType>
2107 void
2111 MatrixType &matrix,
2112 const std::function<void(FEEvaluation<dim,
2113 fe_degree,
2114 n_q_points_1d,
2115 n_components,
2116 Number,
2117 VectorizedArrayType> &)>
2119 const unsigned int dof_no,
2120 const unsigned int quad_no,
2121 const unsigned int first_selected_component)
2122 {
2123 compute_matrix<dim,
2124 fe_degree,
2125 n_q_points_1d,
2126 n_components,
2127 Number,
2128 VectorizedArrayType,
2129 MatrixType>(matrix_free,
2131 matrix,
2133 {},
2134 {},
2135 dof_no,
2136 quad_no,
2137 first_selected_component);
2138 }
2139
2140 template <typename CLASS,
2141 int dim,
2142 int fe_degree,
2143 int n_q_points_1d,
2144 int n_components,
2145 typename Number,
2146 typename VectorizedArrayType,
2147 typename MatrixType>
2148 void
2151 const AffineConstraints<Number> &constraints,
2152 MatrixType &matrix,
2153 void (CLASS::*cell_operation)(FEEvaluation<dim,
2154 fe_degree,
2155 n_q_points_1d,
2156 n_components,
2157 Number,
2158 VectorizedArrayType> &) const,
2160 const unsigned int dof_no,
2161 const unsigned int quad_no,
2162 const unsigned int first_selected_component)
2163 {
2164 compute_matrix<dim,
2165 fe_degree,
2166 n_q_points_1d,
2167 n_components,
2168 Number,
2169 VectorizedArrayType,
2170 MatrixType>(
2171 matrix_free,
2172 constraints,
2173 matrix,
2174 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
2175 dof_no,
2176 quad_no,
2177 first_selected_component);
2178 }
2179
2180 template <int dim,
2181 int fe_degree,
2182 int n_q_points_1d,
2183 int n_components,
2184 typename Number,
2185 typename VectorizedArrayType,
2186 typename MatrixType>
2187 void
2191 MatrixType &matrix,
2192 const std::function<void(FEEvaluation<dim,
2193 fe_degree,
2194 n_q_points_1d,
2195 n_components,
2196 Number,
2197 VectorizedArrayType> &)>
2199 const std::function<void(FEFaceEvaluation<dim,
2200 fe_degree,
2201 n_q_points_1d,
2202 n_components,
2203 Number,
2204 VectorizedArrayType> &,
2205 FEFaceEvaluation<dim,
2206 fe_degree,
2207 n_q_points_1d,
2208 n_components,
2209 Number,
2210 VectorizedArrayType> &)>
2212 const std::function<void(FEFaceEvaluation<dim,
2213 fe_degree,
2214 n_q_points_1d,
2215 n_components,
2216 Number,
2217 VectorizedArrayType> &)>
2219 const unsigned int dof_no,
2220 const unsigned int quad_no,
2221 const unsigned int first_selected_component)
2222 {
2223 using FEEvalType = FEEvaluation<dim,
2224 fe_degree,
2225 n_q_points_1d,
2226 n_components,
2227 Number,
2228 VectorizedArrayType>;
2229
2230 using FEFaceEvalType = FEFaceEvaluation<dim,
2231 fe_degree,
2232 n_q_points_1d,
2233 n_components,
2234 Number,
2235 VectorizedArrayType>;
2236
2237 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2238 data_cell;
2239
2240 data_cell.dof_numbers = {dof_no};
2241 data_cell.quad_numbers = {quad_no};
2242 data_cell.n_components = {n_components};
2243 data_cell.first_selected_components = {first_selected_component};
2244 data_cell.batch_type = {0};
2245
2246 data_cell.op_create =
2247 [&](const std::pair<unsigned int, unsigned int> &range) {
2248 std::vector<
2249 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
2250 phi;
2251
2252 if (!internal::is_fe_nothing<false>(matrix_free,
2253 range,
2254 dof_no,
2255 quad_no,
2256 first_selected_component,
2257 fe_degree,
2258 n_q_points_1d))
2259 phi.emplace_back(std::make_unique<FEEvalType>(
2260 matrix_free, range, dof_no, quad_no, first_selected_component));
2261
2262 return phi;
2263 };
2264
2265 data_cell.op_reinit = [](auto &phi, const unsigned batch) {
2266 if (phi.size() == 1)
2267 static_cast<FEEvalType &>(*phi[0]).reinit(batch);
2268 };
2269
2270 if (cell_operation)
2271 data_cell.op_compute = [&](auto &phi) {
2272 cell_operation(static_cast<FEEvalType &>(*phi[0]));
2273 };
2274
2275 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2276 data_face;
2277
2278 data_face.dof_numbers = {dof_no, dof_no};
2279 data_face.quad_numbers = {quad_no, quad_no};
2280 data_face.n_components = {n_components, n_components};
2281 data_face.first_selected_components = {first_selected_component,
2282 first_selected_component};
2283 data_face.batch_type = {1, 2};
2284
2285 data_face.op_create =
2286 [&](const std::pair<unsigned int, unsigned int> &range) {
2287 std::vector<
2288 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2289 phi;
2290
2291 if (!internal::is_fe_nothing<true>(matrix_free,
2292 range,
2293 dof_no,
2294 quad_no,
2295 first_selected_component,
2296 fe_degree,
2297 n_q_points_1d,
2298 true) &&
2299 !internal::is_fe_nothing<true>(matrix_free,
2300 range,
2301 dof_no,
2302 quad_no,
2303 first_selected_component,
2304 fe_degree,
2305 n_q_points_1d,
2306 false))
2307 {
2308 phi.emplace_back(
2309 std::make_unique<FEFaceEvalType>(matrix_free,
2310 range,
2311 true,
2312 dof_no,
2313 quad_no,
2314 first_selected_component));
2315 phi.emplace_back(
2316 std::make_unique<FEFaceEvalType>(matrix_free,
2317 range,
2318 false,
2319 dof_no,
2320 quad_no,
2321 first_selected_component));
2322 }
2323
2324 return phi;
2325 };
2326
2327 data_face.op_reinit = [](auto &phi, const unsigned batch) {
2328 if (phi.size() == 2)
2329 {
2330 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
2331 static_cast<FEFaceEvalType &>(*phi[1]).reinit(batch);
2332 }
2333 };
2334
2335 if (face_operation)
2336 data_face.op_compute = [&](auto &phi) {
2337 face_operation(static_cast<FEFaceEvalType &>(*phi[0]),
2338 static_cast<FEFaceEvalType &>(*phi[1]));
2339 };
2340
2341 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2343
2344 data_boundary.dof_numbers = {dof_no};
2345 data_boundary.quad_numbers = {quad_no};
2346 data_boundary.n_components = {n_components};
2347 data_boundary.first_selected_components = {first_selected_component};
2348 data_boundary.batch_type = {1};
2349
2351 .op_create = [&](const std::pair<unsigned int, unsigned int> &range) {
2352 std::vector<
2353 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2354 phi;
2355
2356 if (!internal::is_fe_nothing<true>(matrix_free,
2357 range,
2358 dof_no,
2359 quad_no,
2360 first_selected_component,
2361 fe_degree,
2362 n_q_points_1d,
2363 true))
2364 phi.emplace_back(std::make_unique<FEFaceEvalType>(
2365 matrix_free, range, true, dof_no, quad_no, first_selected_component));
2366
2367 return phi;
2368 };
2369
2370 data_boundary.op_reinit = [](auto &phi, const unsigned batch) {
2371 if (phi.size() == 1)
2372 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
2373 };
2374
2376 data_boundary.op_compute = [&](auto &phi) {
2377 boundary_operation(static_cast<FEFaceEvalType &>(*phi[0]));
2378 };
2379
2380 internal::compute_matrix(
2381 matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix);
2382 }
2383
2384 namespace internal
2385 {
2386 template <int dim,
2387 typename Number,
2388 typename VectorizedArrayType,
2389 typename MatrixType>
2390 void
2394 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2395 &data_cell,
2396 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2397 &data_face,
2398 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2400 MatrixType &matrix)
2401 {
2402 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
2405 internal::create_new_affine_constraints_if_needed(
2407
2408 const auto batch_operation =
2409 [&matrix_free, &constraints, &matrix](
2410 auto &data, const std::pair<unsigned int, unsigned int> &range) {
2411 if (!data.op_compute)
2412 return; // nothing to do
2413
2414 auto phi = data.op_create(range);
2415
2416 const unsigned int n_blocks = phi.size();
2417
2418 if (n_blocks == 0)
2419 return;
2420
2421 Table<1, unsigned int> dofs_per_cell(n_blocks);
2422
2423 Table<1, std::vector<types::global_dof_index>> dof_indices(n_blocks);
2425 n_blocks, VectorizedArrayType::size());
2426 Table<1, std::vector<unsigned int>> lexicographic_numbering(n_blocks);
2427 Table<2,
2428 std::array<FullMatrix<typename MatrixType::value_type>,
2429 VectorizedArrayType::size()>>
2430 matrices(n_blocks, n_blocks);
2431
2432 for (unsigned int b = 0; b < n_blocks; ++b)
2433 {
2434 const auto &fe = matrix_free.get_dof_handler(data.dof_numbers[b])
2435 .get_fe(phi[b]->get_active_fe_index());
2436
2437 const auto component_base =
2438 matrix_free.get_dof_info(data.dof_numbers[b])
2439 .component_to_base_index[data.first_selected_components[b]];
2440 const auto component_in_base =
2441 data.first_selected_components[b] -
2442 matrix_free.get_dof_info(data.dof_numbers[b])
2443 .start_components[component_base];
2444
2445 const auto &shape_info = matrix_free.get_shape_info(
2446 data.dof_numbers[b],
2447 data.quad_numbers[b],
2449 phi[b]->get_active_fe_index(),
2450 phi[b]->get_active_quadrature_index());
2451
2452 dofs_per_cell[b] =
2453 shape_info.dofs_per_component_on_cell * data.n_components[b];
2454
2455 dof_indices[b].resize(fe.n_dofs_per_cell());
2456
2457 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2458 dof_indices_mf[b][v].resize(dofs_per_cell[b]);
2459
2460 lexicographic_numbering[b].insert(
2461 lexicographic_numbering[b].begin(),
2462 shape_info.lexicographic_numbering.begin() +
2463 component_in_base * shape_info.dofs_per_component_on_cell,
2464 shape_info.lexicographic_numbering.begin() +
2465 (component_in_base + data.n_components[b]) *
2466 shape_info.dofs_per_component_on_cell);
2467 }
2468
2469 for (unsigned int bj = 0; bj < n_blocks; ++bj)
2470 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2471 std::fill_n(matrices[bi][bj].begin(),
2472 VectorizedArrayType::size(),
2474 dofs_per_cell[bi], dofs_per_cell[bj]));
2475
2476 for (auto batch = range.first; batch < range.second; ++batch)
2477 {
2478 data.op_reinit(phi, batch);
2479
2480 const unsigned int n_filled_lanes =
2481 data.is_face ?
2482 matrix_free.n_active_entries_per_face_batch(batch) :
2483 matrix_free.n_active_entries_per_cell_batch(batch);
2484
2485 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2486 for (unsigned int b = 0; b < n_blocks; ++b)
2487 {
2488 unsigned int const cell_index =
2489 (data.batch_type[b] == 0) ?
2490 (batch * VectorizedArrayType::size() + v) :
2491 ((data.batch_type[b] == 1) ?
2492 matrix_free.get_face_info(batch).cells_interior[v] :
2493 matrix_free.get_face_info(batch).cells_exterior[v]);
2494
2495 const auto cell_iterator = matrix_free.get_cell_iterator(
2496 cell_index / VectorizedArrayType::size(),
2497 cell_index % VectorizedArrayType::size(),
2498 data.dof_numbers[b]);
2499
2500 if (matrix_free.get_mg_level() !=
2502 cell_iterator->get_mg_dof_indices(dof_indices[b]);
2503 else
2504 cell_iterator->get_dof_indices(dof_indices[b]);
2505
2506 for (unsigned int j = 0; j < dofs_per_cell[b]; ++j)
2507 dof_indices_mf[b][v][j] =
2508 dof_indices[b][lexicographic_numbering[b][j]];
2509 }
2510
2511 for (unsigned int bj = 0; bj < n_blocks; ++bj)
2512 {
2513 for (unsigned int j = 0; j < dofs_per_cell[bj]; ++j)
2514 {
2515 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2516 for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2517 phi[bi]->begin_dof_values()[i] =
2518 (bj == bi) ? static_cast<Number>(i == j) : 0.0;
2519
2520 data.op_compute(phi);
2521
2522 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2523 for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2524 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2525 matrices[bi][bj][v](i, j) =
2526 phi[bi]->begin_dof_values()[i][v];
2527 }
2528
2529 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2530 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2531 if (bi == bj)
2532 // specialization for blocks on the diagonal
2533 // to writing into diagonal elements of the
2534 // matrix if the corresponding degree of freedom
2535 // is constrained, see also the documentation
2536 // of AffineConstraints::distribute_local_to_global()
2537 constraints.distribute_local_to_global(
2538 matrices[bi][bi][v], dof_indices_mf[bi][v], matrix);
2539 else
2540 constraints.distribute_local_to_global(
2541 matrices[bi][bj][v],
2542 dof_indices_mf[bi][v],
2543 dof_indices_mf[bj][v],
2544 matrix);
2545 }
2546 }
2547 };
2548
2549 const auto cell_operation_wrapped =
2550 [&](const auto &, auto &, const auto &, const auto range) {
2552 };
2553
2554 const auto face_operation_wrapped =
2555 [&](const auto &, auto &, const auto &, const auto range) {
2557 };
2558
2559 const auto boundary_operation_wrapped =
2560 [&](const auto &, auto &, const auto &, const auto range) {
2562 };
2563
2564 if (data_face.op_compute || data_boundary.op_compute)
2565 {
2566 matrix_free.template loop<MatrixType, MatrixType>(
2570 matrix,
2571 matrix);
2572 }
2573 else
2574 matrix_free.template cell_loop<MatrixType, MatrixType>(
2575 cell_operation_wrapped, matrix, matrix);
2576
2577 matrix.compress(VectorOperation::add);
2578 }
2579 } // namespace internal
2580
2581 template <typename CLASS,
2582 int dim,
2583 int fe_degree,
2584 int n_q_points_1d,
2585 int n_components,
2586 typename Number,
2587 typename VectorizedArrayType,
2588 typename MatrixType>
2589 void
2592 const AffineConstraints<Number> &constraints,
2593 MatrixType &matrix,
2594 void (CLASS::*cell_operation)(FEEvaluation<dim,
2595 fe_degree,
2596 n_q_points_1d,
2597 n_components,
2598 Number,
2599 VectorizedArrayType> &) const,
2601 fe_degree,
2602 n_q_points_1d,
2603 n_components,
2604 Number,
2605 VectorizedArrayType> &,
2606 FEFaceEvaluation<dim,
2607 fe_degree,
2608 n_q_points_1d,
2609 n_components,
2610 Number,
2611 VectorizedArrayType> &)
2612 const,
2614 fe_degree,
2615 n_q_points_1d,
2616 n_components,
2617 Number,
2618 VectorizedArrayType> &)
2619 const,
2621 const unsigned int dof_no,
2622 const unsigned int quad_no,
2623 const unsigned int first_selected_component)
2624 {
2625 compute_matrix<dim,
2626 fe_degree,
2627 n_q_points_1d,
2628 n_components,
2629 Number,
2630 VectorizedArrayType,
2631 MatrixType>(
2632 matrix_free,
2633 constraints,
2634 matrix,
2635 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
2636 [&](auto &phi_m, auto &phi_p) {
2638 },
2639 [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
2640 dof_no,
2641 quad_no,
2642 first_selected_component);
2643 }
2644
2645#endif // DOXYGEN
2646
2647} // namespace MatrixFreeTools
2648
2650
2651
2652#endif
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) const
Definition tools.h:580
ObserverPointer< const MatrixFree< dim, Number, VectorizedArrayType > > matrix_free
Definition tools.h:647
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) const
Definition tools.h:546
void reinit(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const AdditionalData &additional_data=AdditionalData())
Definition tools.h:517
std::vector< unsigned int > n_components
Definition tools.h:47
std::vector< unsigned int > quad_numbers
Definition tools.h:46
std::vector< unsigned int > batch_type
Definition tools.h:49
std::vector< unsigned int > dof_numbers
Definition tools.h:45
std::function< void(std::vector< std::unique_ptr< FEEvalType > > &, const unsigned int)> op_reinit
Definition tools.h:57
std::function< void(std::vector< std::unique_ptr< FEEvalType > > &)> op_compute
Definition tools.h:59
std::function< std::vector< std::unique_ptr< FEEvalType > >(const std::pair< unsigned int, unsigned int > &)> op_create
Definition tools.h:54
std::vector< unsigned int > first_selected_components
Definition tools.h:48
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number, MemorySpaceType > &vec) const
void set_constrained_values(const Number val, VectorType &dst) const
void cell_loop(const Functor &func, const VectorType &src, VectorType &dst) const
constexpr void clear()
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
unsigned int level
Definition grid_out.cc:4632
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
void cell_action(IteratorType cell, DoFInfoBox< dim, DOFINFO > &dof_info, 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, const LoopControl &loop_control)
Definition loop.h:316
std::vector< index_type > data
Definition mpi.cc:740
std::size_t size
Definition mpi.cc:739
EvaluationFlags
The EvaluationFlags enum.
@ matrix
Contents is actually a matrix.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:49
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::ComputeMatrixScratchData< dim, VectorizedArrayType, false > &data_cell, const internal::ComputeMatrixScratchData< dim, VectorizedArrayType, true > &data_face, const internal::ComputeMatrixScratchData< dim, VectorizedArrayType, true > &data_boundary, VectorType &diagonal_global, std::vector< VectorType2 * > &diagonal_global_components)
void compute_matrix(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const AffineConstraints< Number > &constraints, const internal::ComputeMatrixScratchData< dim, VectorizedArrayType, false > &cell_operation, const internal::ComputeMatrixScratchData< dim, VectorizedArrayType, true > &face_operation, const internal::ComputeMatrixScratchData< dim, VectorizedArrayType, true > &boundary_operation, MatrixType &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 > &)> &cell_operation, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int first_vector_component=0)
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 > &)> &cell_operation, 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)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
std::uint8_t compressed_constraint_kind
Definition dof_info.h:86
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)
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)