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