deal.II version GIT relicensing-3301-g5b03c77952 2025-05-11 13:40:00+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#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, dof_no);
1356 else if (is_interior_face)
1357 active_fe_index =
1358 matrix_free.get_face_range_category(range, dof_no).first;
1359 else
1360 active_fe_index =
1361 matrix_free.get_face_range_category(range, dof_no).second;
1362
1363 const auto init_data = ::internal::
1364 extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
1365 matrix_free,
1366 dof_no,
1367 first_selected_component,
1368 quad_no,
1369 fe_degree,
1370 static_n_q_points,
1371 active_fe_index,
1372 numbers::invalid_unsigned_int /*active_quad_index*/,
1373 numbers::invalid_unsigned_int /*face_type*/);
1374
1375 return init_data.shape_info->dofs_per_component_on_cell == 0;
1376 }
1377
1378
1379
1383 template <int dim,
1384 int fe_degree,
1385 int n_q_points_1d,
1386 int n_components,
1387 typename Number,
1388 typename QuadOperation>
1389 class ComputeDiagonalCellAction
1390 {
1391 public:
1392 ComputeDiagonalCellAction(
1393 const QuadOperation &quad_operation,
1394 const EvaluationFlags::EvaluationFlags evaluation_flags,
1395 const EvaluationFlags::EvaluationFlags integration_flags)
1396 : m_quad_operation(quad_operation)
1397 , m_evaluation_flags(evaluation_flags)
1398 , m_integration_flags(integration_flags)
1399 {}
1400
1401 KOKKOS_FUNCTION void
1402 operator()(const typename Portable::MatrixFree<dim, Number>::Data *data,
1405 {
1406 Portable::
1407 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>
1408 fe_eval(data);
1410 *gpu_data = data->precomputed_data;
1411 const int cell = data->cell_index;
1412
1413 constexpr int dofs_per_cell = decltype(fe_eval)::tensor_dofs_per_cell;
1414 typename decltype(fe_eval)::value_type
1415 diagonal[dofs_per_cell / n_components] = {};
1416 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1417 {
1418 const auto c = i % n_components;
1419
1420 Kokkos::parallel_for(
1421 Kokkos::TeamThreadRange(data->team_member,
1422 dofs_per_cell / n_components),
1423 [&](int j) {
1424 typename decltype(fe_eval)::value_type val = {};
1425
1426 if constexpr (n_components == 1)
1427 {
1428 val = (i == j) ? 1 : 0;
1429 }
1430 else
1431 {
1432 val[c] = (i / n_components == j) ? 1 : 0;
1433 }
1434
1435 fe_eval.submit_dof_value(val, j);
1436 });
1437
1438 data->team_member.team_barrier();
1439
1440 Portable::internal::
1441 resolve_hanging_nodes<dim, fe_degree, false, Number>(
1442 data->team_member,
1443 gpu_data->constraint_weights,
1444 gpu_data->constraint_mask(cell * n_components + c),
1445 Kokkos::subview(data->shared_data->values, Kokkos::ALL, c));
1446
1447 fe_eval.evaluate(m_evaluation_flags);
1448 fe_eval.apply_for_each_quad_point(m_quad_operation);
1449 fe_eval.integrate(m_integration_flags);
1450
1451 Portable::internal::
1452 resolve_hanging_nodes<dim, fe_degree, true, Number>(
1453 data->team_member,
1454 gpu_data->constraint_weights,
1455 gpu_data->constraint_mask(cell * n_components + c),
1456 Kokkos::subview(data->shared_data->values, Kokkos::ALL, c));
1457
1458 Kokkos::single(Kokkos::PerTeam(data->team_member), [&] {
1459 if constexpr (n_components == 1)
1460 diagonal[i] = fe_eval.get_dof_value(i);
1461 else
1462 diagonal[i / n_components][i % n_components] =
1463 fe_eval.get_dof_value(i / n_components)[i % n_components];
1464 });
1465
1466 data->team_member.team_barrier();
1467 }
1468
1469 Kokkos::single(Kokkos::PerTeam(data->team_member), [&] {
1470 for (unsigned int i = 0; i < dofs_per_cell / n_components; ++i)
1471 fe_eval.submit_dof_value(diagonal[i], i);
1472 });
1473
1474 data->team_member.team_barrier();
1475
1476 // We need to do the same as distribute_local_to_global but without
1477 // constraints since we have already taken care of them earlier
1478 if (gpu_data->use_coloring)
1479 {
1480 Kokkos::parallel_for(
1481 Kokkos::TeamThreadRange(data->team_member, dofs_per_cell),
1482 [&](const int &i) {
1483 dst[gpu_data->local_to_global(i, cell)] +=
1484 data->shared_data->values(i % (dofs_per_cell / n_components),
1485 i / (dofs_per_cell / n_components));
1486 });
1487 }
1488 else
1489 {
1490 Kokkos::parallel_for(
1491 Kokkos::TeamThreadRange(data->team_member, dofs_per_cell),
1492 [&](const int &i) {
1493 Kokkos::atomic_add(&dst[gpu_data->local_to_global(i, cell)],
1494 data->shared_data->values(
1495 i % (dofs_per_cell / n_components),
1496 i / (dofs_per_cell / n_components)));
1497 });
1498 }
1499 };
1500
1501 static constexpr unsigned int n_local_dofs =
1502 ::Utilities::pow(fe_degree + 1, dim) * n_components;
1503 static constexpr unsigned int n_q_points =
1504 ::Utilities::pow(n_q_points_1d, dim);
1505
1506 private:
1507 const QuadOperation m_quad_operation;
1508 const EvaluationFlags::EvaluationFlags m_evaluation_flags;
1509 const EvaluationFlags::EvaluationFlags m_integration_flags;
1510 };
1511
1512 } // namespace internal
1513
1514
1515
1516 template <int dim,
1517 int fe_degree,
1518 int n_q_points_1d,
1519 int n_components,
1520 typename Number,
1521 typename MemorySpace,
1522 typename QuadOperation>
1523 void
1525 const Portable::MatrixFree<dim, Number> &matrix_free,
1527 const QuadOperation &quad_operation,
1528 EvaluationFlags::EvaluationFlags evaluation_flags,
1529 EvaluationFlags::EvaluationFlags integration_flags,
1530 const unsigned int dof_no,
1531 const unsigned int quad_no,
1532 const unsigned int first_selected_component,
1533 const unsigned int first_vector_component)
1534 {
1535 Assert(dof_no == 0, ExcNotImplemented());
1536 Assert(quad_no == 0, ExcNotImplemented());
1537 Assert(first_selected_component == 0, ExcNotImplemented());
1538 Assert(first_vector_component == 0, ExcNotImplemented());
1539
1540 matrix_free.initialize_dof_vector(diagonal_global);
1541
1542
1543 internal::ComputeDiagonalCellAction<dim,
1544 fe_degree,
1545 n_q_points_1d,
1546 n_components,
1547 Number,
1548 QuadOperation>
1549 cell_action(quad_operation, evaluation_flags, integration_flags);
1551 matrix_free.cell_loop(cell_action, dummy, diagonal_global);
1552
1553 matrix_free.set_constrained_values(Number(1.), diagonal_global);
1554 }
1555
1556 template <int dim,
1557 int fe_degree,
1558 int n_q_points_1d,
1559 int n_components,
1560 typename Number,
1561 typename VectorizedArrayType,
1562 typename VectorType>
1563 void
1566 VectorType &diagonal_global,
1567 const std::function<void(FEEvaluation<dim,
1568 fe_degree,
1569 n_q_points_1d,
1570 n_components,
1571 Number,
1572 VectorizedArrayType> &)>
1573 &cell_operation,
1574 const unsigned int dof_no,
1575 const unsigned int quad_no,
1576 const unsigned int first_selected_component,
1577 const unsigned int first_vector_component)
1578 {
1579 compute_diagonal<dim,
1580 fe_degree,
1581 n_q_points_1d,
1582 n_components,
1583 Number,
1584 VectorizedArrayType,
1585 VectorType>(matrix_free,
1586 diagonal_global,
1587 cell_operation,
1588 {},
1589 {},
1590 dof_no,
1591 quad_no,
1592 first_selected_component,
1593 first_vector_component);
1594 }
1595
1596 template <typename CLASS,
1597 int dim,
1598 int fe_degree,
1599 int n_q_points_1d,
1600 int n_components,
1601 typename Number,
1602 typename VectorizedArrayType,
1603 typename VectorType>
1604 void
1607 VectorType &diagonal_global,
1608 void (CLASS::*cell_operation)(FEEvaluation<dim,
1609 fe_degree,
1610 n_q_points_1d,
1611 n_components,
1612 Number,
1613 VectorizedArrayType> &) const,
1614 const CLASS *owning_class,
1615 const unsigned int dof_no,
1616 const unsigned int quad_no,
1617 const unsigned int first_selected_component,
1618 const unsigned int first_vector_component)
1619 {
1620 compute_diagonal<dim,
1621 fe_degree,
1622 n_q_points_1d,
1623 n_components,
1624 Number,
1625 VectorizedArrayType,
1626 VectorType>(
1627 matrix_free,
1628 diagonal_global,
1629 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
1630 dof_no,
1631 quad_no,
1632 first_selected_component,
1633 first_vector_component);
1634 }
1635
1636 template <int dim,
1637 int fe_degree,
1638 int n_q_points_1d,
1639 int n_components,
1640 typename Number,
1641 typename VectorizedArrayType,
1642 typename VectorType>
1643 void
1646 VectorType &diagonal_global,
1647 const std::function<void(FEEvaluation<dim,
1648 fe_degree,
1649 n_q_points_1d,
1650 n_components,
1651 Number,
1652 VectorizedArrayType> &)>
1653 &cell_operation,
1654 const std::function<void(FEFaceEvaluation<dim,
1655 fe_degree,
1656 n_q_points_1d,
1657 n_components,
1658 Number,
1659 VectorizedArrayType> &,
1660 FEFaceEvaluation<dim,
1661 fe_degree,
1662 n_q_points_1d,
1663 n_components,
1664 Number,
1665 VectorizedArrayType> &)>
1666 &face_operation,
1667 const std::function<void(FEFaceEvaluation<dim,
1668 fe_degree,
1669 n_q_points_1d,
1670 n_components,
1671 Number,
1672 VectorizedArrayType> &)>
1673 &boundary_operation,
1674 const unsigned int dof_no,
1675 const unsigned int quad_no,
1676 const unsigned int first_selected_component,
1677 const unsigned int first_vector_component)
1678 {
1679 std::vector<typename ::internal::BlockVectorSelector<
1680 VectorType,
1681 IsBlockVector<VectorType>::value>::BaseVectorType *>
1682 diagonal_global_components(n_components);
1683
1684 for (unsigned int d = 0; d < n_components; ++d)
1685 diagonal_global_components[d] = ::internal::
1686 BlockVectorSelector<VectorType, IsBlockVector<VectorType>::value>::
1687 get_vector_component(diagonal_global, d + first_vector_component);
1688
1689 const auto &dof_info = matrix_free.get_dof_info(dof_no);
1690
1691 if (dof_info.start_components.back() == 1)
1692 for (unsigned int comp = 0; comp < n_components; ++comp)
1693 {
1694 Assert(diagonal_global_components[comp] != nullptr,
1695 ExcMessage("The finite element underlying this FEEvaluation "
1696 "object is scalar, but you requested " +
1697 std::to_string(n_components) +
1698 " components via the template argument in "
1699 "FEEvaluation. In that case, you must pass an "
1700 "std::vector<VectorType> or a BlockVector to " +
1701 "read_dof_values and distribute_local_to_global."));
1703 *diagonal_global_components[comp], matrix_free, dof_info);
1704 }
1705 else
1706 {
1708 *diagonal_global_components[0], matrix_free, dof_info);
1709 }
1710
1711 using FEEvalType = FEEvaluation<dim,
1712 fe_degree,
1713 n_q_points_1d,
1714 n_components,
1715 Number,
1716 VectorizedArrayType>;
1717
1718 using FEFaceEvalType = FEFaceEvaluation<dim,
1719 fe_degree,
1720 n_q_points_1d,
1721 n_components,
1722 Number,
1723 VectorizedArrayType>;
1724
1725 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1726 data_cell;
1727
1728 data_cell.dof_numbers = {dof_no};
1729 data_cell.quad_numbers = {quad_no};
1730 data_cell.n_components = {n_components};
1731 data_cell.first_selected_components = {first_selected_component};
1732 data_cell.batch_type = {0};
1733
1734 data_cell.op_create =
1735 [&](const std::pair<unsigned int, unsigned int> &range) {
1736 std::vector<
1737 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
1738 phi;
1739
1740 if (!internal::is_fe_nothing<false>(matrix_free,
1741 range,
1742 dof_no,
1743 quad_no,
1744 first_selected_component,
1745 fe_degree,
1746 n_q_points_1d))
1747 phi.emplace_back(std::make_unique<FEEvalType>(
1748 matrix_free, range, dof_no, quad_no, first_selected_component));
1749
1750 return phi;
1751 };
1752
1753 data_cell.op_reinit = [](auto &phi, const unsigned batch) {
1754 if (phi.size() == 1)
1755 static_cast<FEEvalType &>(*phi[0]).reinit(batch);
1756 };
1757
1758 if (cell_operation)
1759 data_cell.op_compute = [&](auto &phi) {
1760 cell_operation(static_cast<FEEvalType &>(*phi[0]));
1761 };
1762
1763 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1764 data_face;
1765
1766 data_face.dof_numbers = {dof_no, dof_no};
1767 data_face.quad_numbers = {quad_no, quad_no};
1768 data_face.n_components = {n_components, n_components};
1769 data_face.first_selected_components = {first_selected_component,
1770 first_selected_component};
1771 data_face.batch_type = {1, 2};
1772
1773 data_face.op_create =
1774 [&](const std::pair<unsigned int, unsigned int> &range) {
1775 std::vector<
1776 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1777 phi;
1778
1779 if (!internal::is_fe_nothing<true>(matrix_free,
1780 range,
1781 dof_no,
1782 quad_no,
1783 first_selected_component,
1784 fe_degree,
1785 n_q_points_1d,
1786 true) &&
1787 !internal::is_fe_nothing<true>(matrix_free,
1788 range,
1789 dof_no,
1790 quad_no,
1791 first_selected_component,
1792 fe_degree,
1793 n_q_points_1d,
1794 false))
1795 {
1796 phi.emplace_back(
1797 std::make_unique<FEFaceEvalType>(matrix_free,
1798 range,
1799 true,
1800 dof_no,
1801 quad_no,
1802 first_selected_component));
1803 phi.emplace_back(
1804 std::make_unique<FEFaceEvalType>(matrix_free,
1805 range,
1806 false,
1807 dof_no,
1808 quad_no,
1809 first_selected_component));
1810 }
1811
1812 return phi;
1813 };
1814
1815 data_face.op_reinit = [](auto &phi, const unsigned batch) {
1816 if (phi.size() == 2)
1817 {
1818 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
1819 static_cast<FEFaceEvalType &>(*phi[1]).reinit(batch);
1820 }
1821 };
1822
1823 if (face_operation)
1824 data_face.op_compute = [&](auto &phi) {
1825 face_operation(static_cast<FEFaceEvalType &>(*phi[0]),
1826 static_cast<FEFaceEvalType &>(*phi[1]));
1827 };
1828
1829 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1830 data_boundary;
1831
1832 data_boundary.dof_numbers = {dof_no};
1833 data_boundary.quad_numbers = {quad_no};
1834 data_boundary.n_components = {n_components};
1835 data_boundary.first_selected_components = {first_selected_component};
1836 data_boundary.batch_type = {1};
1837
1838 data_boundary
1839 .op_create = [&](const std::pair<unsigned int, unsigned int> &range) {
1840 std::vector<
1841 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
1842 phi;
1843
1844 if (!internal::is_fe_nothing<true>(matrix_free,
1845 range,
1846 dof_no,
1847 quad_no,
1848 first_selected_component,
1849 fe_degree,
1850 n_q_points_1d,
1851 true))
1852 phi.emplace_back(std::make_unique<FEFaceEvalType>(
1853 matrix_free, range, true, dof_no, quad_no, first_selected_component));
1854
1855 return phi;
1856 };
1857
1858 data_boundary.op_reinit = [](auto &phi, const unsigned batch) {
1859 if (phi.size() == 1)
1860 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
1861 };
1862
1863 if (boundary_operation)
1864 data_boundary.op_compute = [&](auto &phi) {
1865 boundary_operation(static_cast<FEFaceEvalType &>(*phi[0]));
1866 };
1867
1868 internal::compute_diagonal(matrix_free,
1869 data_cell,
1870 data_face,
1871 data_boundary,
1872 diagonal_global,
1873 diagonal_global_components);
1874 }
1875
1876 namespace internal
1877 {
1878 template <int dim,
1879 typename Number,
1880 typename VectorizedArrayType,
1881 typename VectorType,
1882 typename VectorType2>
1883 void
1886 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1887 &data_cell,
1888 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1889 &data_face,
1890 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1891 &data_boundary,
1892 VectorType &diagonal_global,
1893 std::vector<VectorType2 *> &diagonal_global_components)
1894 {
1895 // TODO: can we remove diagonal_global_components as argument?
1896
1897 int dummy = 0;
1898
1899 using Helper =
1900 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, false>;
1901
1902 using HelperFace =
1903 internal::ComputeDiagonalHelper<dim, VectorizedArrayType, true>;
1904
1907 scratch_data_internal;
1909
1910 const auto batch_operation =
1911 [&](auto &data,
1912 auto &scratch_data,
1913 const std::pair<unsigned int, unsigned int> &range) {
1914 if (!data.op_compute)
1915 return; // nothing to do
1916
1917 auto phi = data.op_create(range);
1918
1919 const unsigned int n_blocks = phi.size();
1920
1921 auto &helpers = scratch_data.get();
1922 helpers.resize(n_blocks);
1923
1924 for (unsigned int b = 0; b < n_blocks; ++b)
1925 helpers[b].initialize(*phi[b], matrix_free, data.n_components[b]);
1926
1927 for (unsigned int batch = range.first; batch < range.second; ++batch)
1928 {
1929 data.op_reinit(phi, batch);
1930
1931 for (unsigned int b = 0; b < n_blocks; ++b)
1932 helpers[b].reinit(batch);
1933
1934 if (n_blocks > 1)
1935 {
1936 Assert(std::all_of(helpers.begin(),
1937 helpers.end(),
1938 [](const auto &helper) {
1939 return helper.has_simple_constraints();
1940 }),
1942 }
1943
1944 for (unsigned int b = 0; b < n_blocks; ++b)
1945 {
1946 for (unsigned int i = 0;
1947 i < phi[b]->get_shape_info().dofs_per_component_on_cell *
1948 data.n_components[b];
1949 ++i)
1950 {
1951 for (unsigned int bb = 0; bb < n_blocks; ++bb)
1952 if (b == bb)
1953 helpers[bb].prepare_basis_vector(i);
1954 else
1955 helpers[bb].zero_basis_vector();
1956
1957 data.op_compute(phi);
1958 helpers[b].submit();
1959 }
1960
1961 helpers[b].distribute_local_to_global(
1962 diagonal_global_components);
1963 }
1964 }
1965 };
1966
1967 const auto cell_operation_wrapped =
1968 [&](const auto &, auto &, const auto &, const auto range) {
1969 batch_operation(data_cell, scratch_data, range);
1970 };
1971
1972 const auto face_operation_wrapped =
1973 [&](const auto &, auto &, const auto &, const auto range) {
1974 batch_operation(data_face, scratch_data_internal, range);
1975 };
1976
1977 const auto boundary_operation_wrapped =
1978 [&](const auto &, auto &, const auto &, const auto range) {
1979 batch_operation(data_boundary, scratch_data_bc, range);
1980 };
1981
1982 if (data_face.op_compute || data_boundary.op_compute)
1983 matrix_free.template loop<VectorType, int>(cell_operation_wrapped,
1984 face_operation_wrapped,
1985 boundary_operation_wrapped,
1986 diagonal_global,
1987 dummy,
1988 false);
1989 else
1990 matrix_free.template cell_loop<VectorType, int>(cell_operation_wrapped,
1991 diagonal_global,
1992 dummy,
1993 false);
1994 }
1995 } // namespace internal
1996
1997 template <typename CLASS,
1998 int dim,
1999 int fe_degree,
2000 int n_q_points_1d,
2001 int n_components,
2002 typename Number,
2003 typename VectorizedArrayType,
2004 typename VectorType>
2005 void
2008 VectorType &diagonal_global,
2009 void (CLASS::*cell_operation)(FEEvaluation<dim,
2010 fe_degree,
2011 n_q_points_1d,
2012 n_components,
2013 Number,
2014 VectorizedArrayType> &) const,
2015 void (CLASS::*face_operation)(FEFaceEvaluation<dim,
2016 fe_degree,
2017 n_q_points_1d,
2018 n_components,
2019 Number,
2020 VectorizedArrayType> &,
2021 FEFaceEvaluation<dim,
2022 fe_degree,
2023 n_q_points_1d,
2024 n_components,
2025 Number,
2026 VectorizedArrayType> &)
2027 const,
2028 void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
2029 fe_degree,
2030 n_q_points_1d,
2031 n_components,
2032 Number,
2033 VectorizedArrayType> &)
2034 const,
2035 const CLASS *owning_class,
2036 const unsigned int dof_no,
2037 const unsigned int quad_no,
2038 const unsigned int first_selected_component,
2039 const unsigned int first_vector_component)
2040 {
2041 compute_diagonal<dim,
2042 fe_degree,
2043 n_q_points_1d,
2044 n_components,
2045 Number,
2046 VectorizedArrayType,
2047 VectorType>(
2048 matrix_free,
2049 diagonal_global,
2050 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
2051 [&](auto &phi_m, auto &phi_p) {
2052 (owning_class->*face_operation)(phi_m, phi_p);
2053 },
2054 [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
2055 dof_no,
2056 quad_no,
2057 first_selected_component,
2058 first_vector_component);
2059 }
2060
2061 namespace internal
2062 {
2067 template <
2068 typename MatrixType,
2069 typename Number,
2070 std::enable_if_t<std::is_same_v<
2071 std::remove_const_t<
2072 std::remove_reference_t<typename MatrixType::value_type>>,
2073 std::remove_const_t<std::remove_reference_t<Number>>>> * = nullptr>
2075 create_new_affine_constraints_if_needed(
2076 const MatrixType &,
2077 const AffineConstraints<Number> &constraints,
2079 {
2080 return constraints;
2081 }
2082
2088 template <
2089 typename MatrixType,
2090 typename Number,
2091 std::enable_if_t<!std::is_same_v<
2092 std::remove_const_t<
2093 std::remove_reference_t<typename MatrixType::value_type>>,
2094 std::remove_const_t<std::remove_reference_t<Number>>>> * = nullptr>
2096 create_new_affine_constraints_if_needed(
2097 const MatrixType &,
2098 const AffineConstraints<Number> &constraints,
2100 &new_constraints)
2101 {
2102 new_constraints =
2103 std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
2104 new_constraints->copy_from(constraints);
2105
2106 return *new_constraints;
2107 }
2108 } // namespace internal
2109
2110 template <int dim,
2111 int fe_degree,
2112 int n_q_points_1d,
2113 int n_components,
2114 typename Number,
2115 typename VectorizedArrayType,
2116 typename MatrixType>
2117 void
2120 const AffineConstraints<Number> &constraints_in,
2121 MatrixType &matrix,
2122 const std::function<void(FEEvaluation<dim,
2123 fe_degree,
2124 n_q_points_1d,
2125 n_components,
2126 Number,
2127 VectorizedArrayType> &)>
2128 &cell_operation,
2129 const unsigned int dof_no,
2130 const unsigned int quad_no,
2131 const unsigned int first_selected_component)
2132 {
2133 compute_matrix<dim,
2134 fe_degree,
2135 n_q_points_1d,
2136 n_components,
2137 Number,
2138 VectorizedArrayType,
2139 MatrixType>(matrix_free,
2140 constraints_in,
2141 matrix,
2142 cell_operation,
2143 {},
2144 {},
2145 dof_no,
2146 quad_no,
2147 first_selected_component);
2148 }
2149
2150 template <typename CLASS,
2151 int dim,
2152 int fe_degree,
2153 int n_q_points_1d,
2154 int n_components,
2155 typename Number,
2156 typename VectorizedArrayType,
2157 typename MatrixType>
2158 void
2161 const AffineConstraints<Number> &constraints,
2162 MatrixType &matrix,
2163 void (CLASS::*cell_operation)(FEEvaluation<dim,
2164 fe_degree,
2165 n_q_points_1d,
2166 n_components,
2167 Number,
2168 VectorizedArrayType> &) const,
2169 const CLASS *owning_class,
2170 const unsigned int dof_no,
2171 const unsigned int quad_no,
2172 const unsigned int first_selected_component)
2173 {
2174 compute_matrix<dim,
2175 fe_degree,
2176 n_q_points_1d,
2177 n_components,
2178 Number,
2179 VectorizedArrayType,
2180 MatrixType>(
2181 matrix_free,
2182 constraints,
2183 matrix,
2184 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
2185 dof_no,
2186 quad_no,
2187 first_selected_component);
2188 }
2189
2190 template <int dim,
2191 int fe_degree,
2192 int n_q_points_1d,
2193 int n_components,
2194 typename Number,
2195 typename VectorizedArrayType,
2196 typename MatrixType>
2197 void
2200 const AffineConstraints<Number> &constraints_in,
2201 MatrixType &matrix,
2202 const std::function<void(FEEvaluation<dim,
2203 fe_degree,
2204 n_q_points_1d,
2205 n_components,
2206 Number,
2207 VectorizedArrayType> &)>
2208 &cell_operation,
2209 const std::function<void(FEFaceEvaluation<dim,
2210 fe_degree,
2211 n_q_points_1d,
2212 n_components,
2213 Number,
2214 VectorizedArrayType> &,
2215 FEFaceEvaluation<dim,
2216 fe_degree,
2217 n_q_points_1d,
2218 n_components,
2219 Number,
2220 VectorizedArrayType> &)>
2221 &face_operation,
2222 const std::function<void(FEFaceEvaluation<dim,
2223 fe_degree,
2224 n_q_points_1d,
2225 n_components,
2226 Number,
2227 VectorizedArrayType> &)>
2228 &boundary_operation,
2229 const unsigned int dof_no,
2230 const unsigned int quad_no,
2231 const unsigned int first_selected_component)
2232 {
2233 using FEEvalType = FEEvaluation<dim,
2234 fe_degree,
2235 n_q_points_1d,
2236 n_components,
2237 Number,
2238 VectorizedArrayType>;
2239
2240 using FEFaceEvalType = FEFaceEvaluation<dim,
2241 fe_degree,
2242 n_q_points_1d,
2243 n_components,
2244 Number,
2245 VectorizedArrayType>;
2246
2247 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2248 data_cell;
2249
2250 data_cell.dof_numbers = {dof_no};
2251 data_cell.quad_numbers = {quad_no};
2252 data_cell.n_components = {n_components};
2253 data_cell.first_selected_components = {first_selected_component};
2254 data_cell.batch_type = {0};
2255
2256 data_cell.op_create =
2257 [&](const std::pair<unsigned int, unsigned int> &range) {
2258 std::vector<
2259 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
2260 phi;
2261
2262 if (!internal::is_fe_nothing<false>(matrix_free,
2263 range,
2264 dof_no,
2265 quad_no,
2266 first_selected_component,
2267 fe_degree,
2268 n_q_points_1d))
2269 phi.emplace_back(std::make_unique<FEEvalType>(
2270 matrix_free, range, dof_no, quad_no, first_selected_component));
2271
2272 return phi;
2273 };
2274
2275 data_cell.op_reinit = [](auto &phi, const unsigned batch) {
2276 if (phi.size() == 1)
2277 static_cast<FEEvalType &>(*phi[0]).reinit(batch);
2278 };
2279
2280 if (cell_operation)
2281 data_cell.op_compute = [&](auto &phi) {
2282 cell_operation(static_cast<FEEvalType &>(*phi[0]));
2283 };
2284
2285 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2286 data_face;
2287
2288 data_face.dof_numbers = {dof_no, dof_no};
2289 data_face.quad_numbers = {quad_no, quad_no};
2290 data_face.n_components = {n_components, n_components};
2291 data_face.first_selected_components = {first_selected_component,
2292 first_selected_component};
2293 data_face.batch_type = {1, 2};
2294
2295 data_face.op_create =
2296 [&](const std::pair<unsigned int, unsigned int> &range) {
2297 std::vector<
2298 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2299 phi;
2300
2301 if (!internal::is_fe_nothing<true>(matrix_free,
2302 range,
2303 dof_no,
2304 quad_no,
2305 first_selected_component,
2306 fe_degree,
2307 n_q_points_1d,
2308 true) &&
2309 !internal::is_fe_nothing<true>(matrix_free,
2310 range,
2311 dof_no,
2312 quad_no,
2313 first_selected_component,
2314 fe_degree,
2315 n_q_points_1d,
2316 false))
2317 {
2318 phi.emplace_back(
2319 std::make_unique<FEFaceEvalType>(matrix_free,
2320 range,
2321 true,
2322 dof_no,
2323 quad_no,
2324 first_selected_component));
2325 phi.emplace_back(
2326 std::make_unique<FEFaceEvalType>(matrix_free,
2327 range,
2328 false,
2329 dof_no,
2330 quad_no,
2331 first_selected_component));
2332 }
2333
2334 return phi;
2335 };
2336
2337 data_face.op_reinit = [](auto &phi, const unsigned batch) {
2338 if (phi.size() == 2)
2339 {
2340 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
2341 static_cast<FEFaceEvalType &>(*phi[1]).reinit(batch);
2342 }
2343 };
2344
2345 if (face_operation)
2346 data_face.op_compute = [&](auto &phi) {
2347 face_operation(static_cast<FEFaceEvalType &>(*phi[0]),
2348 static_cast<FEFaceEvalType &>(*phi[1]));
2349 };
2350
2351 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2352 data_boundary;
2353
2354 data_boundary.dof_numbers = {dof_no};
2355 data_boundary.quad_numbers = {quad_no};
2356 data_boundary.n_components = {n_components};
2357 data_boundary.first_selected_components = {first_selected_component};
2358 data_boundary.batch_type = {1};
2359
2360 data_boundary
2361 .op_create = [&](const std::pair<unsigned int, unsigned int> &range) {
2362 std::vector<
2363 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2364 phi;
2365
2366 if (!internal::is_fe_nothing<true>(matrix_free,
2367 range,
2368 dof_no,
2369 quad_no,
2370 first_selected_component,
2371 fe_degree,
2372 n_q_points_1d,
2373 true))
2374 phi.emplace_back(std::make_unique<FEFaceEvalType>(
2375 matrix_free, range, true, dof_no, quad_no, first_selected_component));
2376
2377 return phi;
2378 };
2379
2380 data_boundary.op_reinit = [](auto &phi, const unsigned batch) {
2381 if (phi.size() == 1)
2382 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
2383 };
2384
2385 if (boundary_operation)
2386 data_boundary.op_compute = [&](auto &phi) {
2387 boundary_operation(static_cast<FEFaceEvalType &>(*phi[0]));
2388 };
2389
2390 internal::compute_matrix(
2391 matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix);
2392 }
2393
2394 namespace internal
2395 {
2396 template <int dim,
2397 typename Number,
2398 typename VectorizedArrayType,
2399 typename MatrixType>
2400 void
2403 const AffineConstraints<Number> &constraints_in,
2404 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2405 &data_cell,
2406 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2407 &data_face,
2408 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2409 &data_boundary,
2410 MatrixType &matrix)
2411 {
2412 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
2413 constraints_for_matrix;
2415 internal::create_new_affine_constraints_if_needed(
2416 matrix, constraints_in, constraints_for_matrix);
2417
2418 const auto batch_operation =
2419 [&matrix_free, &constraints, &matrix](
2420 auto &data, const std::pair<unsigned int, unsigned int> &range) {
2421 if (!data.op_compute)
2422 return; // nothing to do
2423
2424 auto phi = data.op_create(range);
2425
2426 const unsigned int n_blocks = phi.size();
2427
2428 if (n_blocks == 0)
2429 return;
2430
2431 Table<1, unsigned int> dofs_per_cell(n_blocks);
2432
2433 Table<1, std::vector<types::global_dof_index>> dof_indices(n_blocks);
2435 n_blocks, VectorizedArrayType::size());
2436 Table<1, std::vector<unsigned int>> lexicographic_numbering(n_blocks);
2437 Table<2,
2438 std::array<FullMatrix<typename MatrixType::value_type>,
2439 VectorizedArrayType::size()>>
2440 matrices(n_blocks, n_blocks);
2441
2442 for (unsigned int b = 0; b < n_blocks; ++b)
2443 {
2444 const auto &fe = matrix_free.get_dof_handler(data.dof_numbers[b])
2445 .get_fe(phi[b]->get_active_fe_index());
2446
2447 const auto component_base =
2448 matrix_free.get_dof_info(data.dof_numbers[b])
2449 .component_to_base_index[data.first_selected_components[b]];
2450 const auto component_in_base =
2451 data.first_selected_components[b] -
2452 matrix_free.get_dof_info(data.dof_numbers[b])
2453 .start_components[component_base];
2454
2455 const auto &shape_info = matrix_free.get_shape_info(
2456 data.dof_numbers[b],
2457 data.quad_numbers[b],
2458 component_base,
2459 phi[b]->get_active_fe_index(),
2460 phi[b]->get_active_quadrature_index());
2461
2462 dofs_per_cell[b] =
2463 shape_info.dofs_per_component_on_cell * data.n_components[b];
2464
2465 dof_indices[b].resize(fe.n_dofs_per_cell());
2466
2467 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2468 dof_indices_mf[b][v].resize(dofs_per_cell[b]);
2469
2470 lexicographic_numbering[b].insert(
2471 lexicographic_numbering[b].begin(),
2472 shape_info.lexicographic_numbering.begin() +
2473 component_in_base * shape_info.dofs_per_component_on_cell,
2474 shape_info.lexicographic_numbering.begin() +
2475 (component_in_base + data.n_components[b]) *
2476 shape_info.dofs_per_component_on_cell);
2477 }
2478
2479 for (unsigned int bj = 0; bj < n_blocks; ++bj)
2480 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2481 std::fill_n(matrices[bi][bj].begin(),
2482 VectorizedArrayType::size(),
2484 dofs_per_cell[bi], dofs_per_cell[bj]));
2485
2486 for (auto batch = range.first; batch < range.second; ++batch)
2487 {
2488 data.op_reinit(phi, batch);
2489
2490 const unsigned int n_filled_lanes =
2491 data.is_face ?
2492 matrix_free.n_active_entries_per_face_batch(batch) :
2493 matrix_free.n_active_entries_per_cell_batch(batch);
2494
2495 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2496 for (unsigned int b = 0; b < n_blocks; ++b)
2497 {
2498 unsigned int const cell_index =
2499 (data.batch_type[b] == 0) ?
2500 (batch * VectorizedArrayType::size() + v) :
2501 ((data.batch_type[b] == 1) ?
2502 matrix_free.get_face_info(batch).cells_interior[v] :
2503 matrix_free.get_face_info(batch).cells_exterior[v]);
2504
2505 const auto cell_iterator = matrix_free.get_cell_iterator(
2506 cell_index / VectorizedArrayType::size(),
2507 cell_index % VectorizedArrayType::size(),
2508 data.dof_numbers[b]);
2509
2510 if (matrix_free.get_mg_level() !=
2512 cell_iterator->get_mg_dof_indices(dof_indices[b]);
2513 else
2514 cell_iterator->get_dof_indices(dof_indices[b]);
2515
2516 for (unsigned int j = 0; j < dofs_per_cell[b]; ++j)
2517 dof_indices_mf[b][v][j] =
2518 dof_indices[b][lexicographic_numbering[b][j]];
2519 }
2520
2521 for (unsigned int bj = 0; bj < n_blocks; ++bj)
2522 {
2523 for (unsigned int j = 0; j < dofs_per_cell[bj]; ++j)
2524 {
2525 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2526 for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2527 phi[bi]->begin_dof_values()[i] =
2528 (bj == bi) ? static_cast<Number>(i == j) : 0.0;
2529
2530 data.op_compute(phi);
2531
2532 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2533 for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2534 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2535 matrices[bi][bj][v](i, j) =
2536 phi[bi]->begin_dof_values()[i][v];
2537 }
2538
2539 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2540 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2541 if (bi == bj)
2542 // specialization for blocks on the diagonal
2543 // to writing into diagonal elements of the
2544 // matrix if the corresponding degree of freedom
2545 // is constrained, see also the documentation
2546 // of AffineConstraints::distribute_local_to_global()
2547 constraints.distribute_local_to_global(
2548 matrices[bi][bi][v], dof_indices_mf[bi][v], matrix);
2549 else
2550 constraints.distribute_local_to_global(
2551 matrices[bi][bj][v],
2552 dof_indices_mf[bi][v],
2553 dof_indices_mf[bj][v],
2554 matrix);
2555 }
2556 }
2557 };
2558
2559 const auto cell_operation_wrapped =
2560 [&](const auto &, auto &, const auto &, const auto range) {
2561 batch_operation(data_cell, range);
2562 };
2563
2564 const auto face_operation_wrapped =
2565 [&](const auto &, auto &, const auto &, const auto range) {
2566 batch_operation(data_face, range);
2567 };
2568
2569 const auto boundary_operation_wrapped =
2570 [&](const auto &, auto &, const auto &, const auto range) {
2571 batch_operation(data_boundary, range);
2572 };
2573
2574 if (data_face.op_compute || data_boundary.op_compute)
2575 {
2576 matrix_free.template loop<MatrixType, MatrixType>(
2577 cell_operation_wrapped,
2578 face_operation_wrapped,
2579 boundary_operation_wrapped,
2580 matrix,
2581 matrix);
2582 }
2583 else
2584 matrix_free.template cell_loop<MatrixType, MatrixType>(
2585 cell_operation_wrapped, matrix, matrix);
2586
2587 matrix.compress(VectorOperation::add);
2588 }
2589 } // namespace internal
2590
2591 template <typename CLASS,
2592 int dim,
2593 int fe_degree,
2594 int n_q_points_1d,
2595 int n_components,
2596 typename Number,
2597 typename VectorizedArrayType,
2598 typename MatrixType>
2599 void
2602 const AffineConstraints<Number> &constraints,
2603 MatrixType &matrix,
2604 void (CLASS::*cell_operation)(FEEvaluation<dim,
2605 fe_degree,
2606 n_q_points_1d,
2607 n_components,
2608 Number,
2609 VectorizedArrayType> &) const,
2610 void (CLASS::*face_operation)(FEFaceEvaluation<dim,
2611 fe_degree,
2612 n_q_points_1d,
2613 n_components,
2614 Number,
2615 VectorizedArrayType> &,
2616 FEFaceEvaluation<dim,
2617 fe_degree,
2618 n_q_points_1d,
2619 n_components,
2620 Number,
2621 VectorizedArrayType> &)
2622 const,
2623 void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
2624 fe_degree,
2625 n_q_points_1d,
2626 n_components,
2627 Number,
2628 VectorizedArrayType> &)
2629 const,
2630 const CLASS *owning_class,
2631 const unsigned int dof_no,
2632 const unsigned int quad_no,
2633 const unsigned int first_selected_component)
2634 {
2635 compute_matrix<dim,
2636 fe_degree,
2637 n_q_points_1d,
2638 n_components,
2639 Number,
2640 VectorizedArrayType,
2641 MatrixType>(
2642 matrix_free,
2643 constraints,
2644 matrix,
2645 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
2646 [&](auto &phi_m, auto &phi_p) {
2647 (owning_class->*face_operation)(phi_m, phi_p);
2648 },
2649 [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
2650 dof_no,
2651 quad_no,
2652 first_selected_component);
2653 }
2654
2655#endif // DOXYGEN
2656
2657} // namespace MatrixFreeTools
2658
2660
2661
2662#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 n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
unsigned int get_mg_level() const
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) 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
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range, const unsigned int dof_handler_index=numbers::invalid_unsigned_int) 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:238
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