Reference documentation for deal.II version 9.6.0
\(\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
22#include <deal.II/matrix_free/fe_evaluation.h>
23#include <deal.II/matrix_free/matrix_free.h>
24#include <deal.II/matrix_free/vector_access_internal.h>
25
26
27DEAL_II_NAMESPACE_OPEN
28
33namespace MatrixFreeTools
34{
35 namespace internal
36 {
37 template <int dim, typename Number, bool is_face_>
38 class ComputeMatrixScratchData
39 {
40 public:
41 using FEEvalType = FEEvaluationData<dim, Number, is_face_>;
42
43 std::vector<unsigned int> dof_numbers;
44 std::vector<unsigned int> quad_numbers;
45 std::vector<unsigned int> first_selected_components;
46 std::vector<unsigned int> batch_type;
47 static const bool is_face = is_face_;
48
49 std::function<std::vector<std::unique_ptr<FEEvalType>>(
50 const std::pair<unsigned int, unsigned int> &)>
51 op_create;
52 std::function<void(std::vector<std::unique_ptr<FEEvalType>> &,
53 const unsigned int)>
54 op_reinit;
55 std::function<void(std::vector<std::unique_ptr<FEEvalType>> &)>
56 op_compute;
57 };
58 } // namespace internal
59
65 template <int dim, typename AdditionalData>
66 void
67 categorize_by_boundary_ids(const Triangulation<dim> &tria,
68 AdditionalData &additional_data);
69
70
71
80 template <int dim,
81 int fe_degree,
82 int n_q_points_1d,
83 int n_components,
84 typename Number,
85 typename VectorizedArrayType,
86 typename VectorType>
87 void
88 compute_diagonal(
89 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
90 VectorType &diagonal_global,
91 const std::function<void(FEEvaluation<dim,
92 fe_degree,
93 n_q_points_1d,
94 n_components,
95 Number,
96 VectorizedArrayType> &)>
97 &cell_operation,
98 const unsigned int dof_no = 0,
99 const unsigned int quad_no = 0,
100 const unsigned int first_selected_component = 0);
101
102
103
107 template <typename CLASS,
108 int dim,
109 int fe_degree,
110 int n_q_points_1d,
111 int n_components,
112 typename Number,
113 typename VectorizedArrayType,
114 typename VectorType>
115 void
116 compute_diagonal(
117 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
118 VectorType &diagonal_global,
119 void (CLASS::*cell_operation)(FEEvaluation<dim,
120 fe_degree,
121 n_q_points_1d,
122 n_components,
123 Number,
124 VectorizedArrayType> &) const,
125 const CLASS *owning_class,
126 const unsigned int dof_no = 0,
127 const unsigned int quad_no = 0,
128 const unsigned int first_selected_component = 0);
129
130
131
142 template <int dim,
143 int fe_degree,
144 int n_q_points_1d,
145 int n_components,
146 typename Number,
147 typename VectorizedArrayType,
148 typename VectorType>
149 void
150 compute_diagonal(
151 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
152 VectorType &diagonal_global,
153 const std::function<void(FEEvaluation<dim,
154 fe_degree,
155 n_q_points_1d,
156 n_components,
157 Number,
158 VectorizedArrayType> &)>
159 &cell_operation,
160 const std::function<void(FEFaceEvaluation<dim,
161 fe_degree,
162 n_q_points_1d,
163 n_components,
164 Number,
165 VectorizedArrayType> &,
166 FEFaceEvaluation<dim,
167 fe_degree,
168 n_q_points_1d,
169 n_components,
170 Number,
171 VectorizedArrayType> &)>
172 &face_operation,
173 const std::function<void(FEFaceEvaluation<dim,
174 fe_degree,
175 n_q_points_1d,
176 n_components,
177 Number,
178 VectorizedArrayType> &)>
179 &boundary_operation,
180 const unsigned int dof_no = 0,
181 const unsigned int quad_no = 0,
182 const unsigned int first_selected_component = 0);
183
184
185
189 template <typename CLASS,
190 int dim,
191 int fe_degree,
192 int n_q_points_1d,
193 int n_components,
194 typename Number,
195 typename VectorizedArrayType,
196 typename VectorType>
197 void
198 compute_diagonal(
199 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
200 VectorType &diagonal_global,
201 void (CLASS::*cell_operation)(FEEvaluation<dim,
202 fe_degree,
203 n_q_points_1d,
204 n_components,
205 Number,
206 VectorizedArrayType> &) const,
207 void (CLASS::*face_operation)(FEFaceEvaluation<dim,
208 fe_degree,
209 n_q_points_1d,
210 n_components,
211 Number,
212 VectorizedArrayType> &,
213 FEFaceEvaluation<dim,
214 fe_degree,
215 n_q_points_1d,
216 n_components,
217 Number,
218 VectorizedArrayType> &)
219 const,
220 void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
221 fe_degree,
222 n_q_points_1d,
223 n_components,
224 Number,
225 VectorizedArrayType> &)
226 const,
227 const CLASS *owning_class,
228 const unsigned int dof_no = 0,
229 const unsigned int quad_no = 0,
230 const unsigned int first_selected_component = 0);
231
232
233
242 template <int dim,
243 int fe_degree,
244 int n_q_points_1d,
245 int n_components,
246 typename Number,
247 typename VectorizedArrayType,
248 typename MatrixType>
249 void
250 compute_matrix(
251 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
252 const AffineConstraints<Number> &constraints,
253 MatrixType &matrix,
254 const std::function<void(FEEvaluation<dim,
255 fe_degree,
256 n_q_points_1d,
257 n_components,
258 Number,
259 VectorizedArrayType> &)>
260 &cell_operation,
261 const unsigned int dof_no = 0,
262 const unsigned int quad_no = 0,
263 const unsigned int first_selected_component = 0);
264
265
266
270 template <typename CLASS,
271 int dim,
272 int fe_degree,
273 int n_q_points_1d,
274 int n_components,
275 typename Number,
276 typename VectorizedArrayType,
277 typename MatrixType>
278 void
279 compute_matrix(
280 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
281 const AffineConstraints<Number> &constraints,
282 MatrixType &matrix,
283 void (CLASS::*cell_operation)(FEEvaluation<dim,
284 fe_degree,
285 n_q_points_1d,
286 n_components,
287 Number,
288 VectorizedArrayType> &) const,
289 const CLASS *owning_class,
290 const unsigned int dof_no = 0,
291 const unsigned int quad_no = 0,
292 const unsigned int first_selected_component = 0);
293
294
295 namespace internal
296 {
301 template <int dim,
302 typename Number,
303 typename VectorizedArrayType,
304 typename MatrixType>
305 void
306 compute_matrix(
307 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
308 const AffineConstraints<Number> &constraints,
309 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
310 &cell_operation,
311 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
312 &face_operation,
313 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
314 &boundary_operation,
315 MatrixType &matrix);
316 } // namespace internal
317
318
319
329 template <int dim,
330 int fe_degree,
331 int n_q_points_1d,
332 int n_components,
333 typename Number,
334 typename VectorizedArrayType,
335 typename MatrixType>
336 void
337 compute_matrix(
338 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
339 const AffineConstraints<Number> &constraints,
340 MatrixType &matrix,
341 const std::function<void(FEEvaluation<dim,
342 fe_degree,
343 n_q_points_1d,
344 n_components,
345 Number,
346 VectorizedArrayType> &)>
347 &cell_operation,
348 const std::function<void(FEFaceEvaluation<dim,
349 fe_degree,
350 n_q_points_1d,
351 n_components,
352 Number,
353 VectorizedArrayType> &,
354 FEFaceEvaluation<dim,
355 fe_degree,
356 n_q_points_1d,
357 n_components,
358 Number,
359 VectorizedArrayType> &)>
360 &face_operation,
361 const std::function<void(FEFaceEvaluation<dim,
362 fe_degree,
363 n_q_points_1d,
364 n_components,
365 Number,
366 VectorizedArrayType> &)>
367 &boundary_operation,
368 const unsigned int dof_no = 0,
369 const unsigned int quad_no = 0,
370 const unsigned int first_selected_component = 0);
371
372
373
377 template <typename CLASS,
378 int dim,
379 int fe_degree,
380 int n_q_points_1d,
381 int n_components,
382 typename Number,
383 typename VectorizedArrayType,
384 typename MatrixType>
385 void
386 compute_matrix(
387 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
388 const AffineConstraints<Number> &constraints,
389 MatrixType &matrix,
390 void (CLASS::*cell_operation)(FEEvaluation<dim,
391 fe_degree,
392 n_q_points_1d,
393 n_components,
394 Number,
395 VectorizedArrayType> &) const,
396 void (CLASS::*face_operation)(FEFaceEvaluation<dim,
397 fe_degree,
398 n_q_points_1d,
399 n_components,
400 Number,
401 VectorizedArrayType> &,
402 FEFaceEvaluation<dim,
403 fe_degree,
404 n_q_points_1d,
405 n_components,
406 Number,
407 VectorizedArrayType> &)
408 const,
409 void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
410 fe_degree,
411 n_q_points_1d,
412 n_components,
413 Number,
414 VectorizedArrayType> &)
415 const,
416 const CLASS *owning_class,
417 const unsigned int dof_no = 0,
418 const unsigned int quad_no = 0,
419 const unsigned int first_selected_component = 0);
420
421
422
431 template <int dim,
432 typename Number,
433 typename VectorizedArrayType = VectorizedArray<Number>>
434 class ElementActivationAndDeactivationMatrixFree
435 {
436 public:
441 struct AdditionalData
442 {
446 AdditionalData(const unsigned int dof_index = 0)
447 : dof_index(dof_index)
448 {}
449
453 unsigned int dof_index;
454 };
455
465 void
466 reinit(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
467 const AdditionalData &additional_data = AdditionalData())
468 {
469 this->matrix_free = &matrix_free;
470
471 std::vector<unsigned int> valid_fe_indices;
472
473 const auto &fe_collection =
474 matrix_free.get_dof_handler(additional_data.dof_index)
475 .get_fe_collection();
476
477 for (unsigned int i = 0; i < fe_collection.size(); ++i)
478 if (fe_collection[i].n_dofs_per_cell() > 0)
479 valid_fe_indices.push_back(i);
480
481 // TODO: relax this so that arbitrary number of valid
482 // FEs can be accepted
483 AssertDimension(valid_fe_indices.size(), 1);
484
485 fe_index_valid = *valid_fe_indices.begin();
486 }
487
493 template <typename VectorTypeOut, typename VectorTypeIn>
494 void
495 cell_loop(const std::function<void(
496 const MatrixFree<dim, Number, VectorizedArrayType> &,
497 VectorTypeOut &,
498 const VectorTypeIn &,
499 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
500 VectorTypeOut &dst,
501 const VectorTypeIn &src,
502 const bool zero_dst_vector = false) const
503 {
504 const auto ebd_cell_operation = [&](const auto &matrix_free,
505 auto &dst,
506 const auto &src,
507 const auto &range) {
508 const auto category = matrix_free.get_cell_range_category(range);
509
510 if (category != fe_index_valid)
511 return;
512
513 cell_operation(matrix_free, dst, src, range);
514 };
515
516 matrix_free->template cell_loop<VectorTypeOut, VectorTypeIn>(
517 ebd_cell_operation, dst, src, zero_dst_vector);
518 }
519
527 template <typename VectorTypeOut, typename VectorTypeIn>
528 void
529 loop(const std::function<
530 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
531 VectorTypeOut &,
532 const VectorTypeIn &,
533 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
534 const std::function<
535 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
536 VectorTypeOut &,
537 const VectorTypeIn &,
538 const std::pair<unsigned int, unsigned int> &)> &face_operation,
539 const std::function<
540 void(const MatrixFree<dim, Number, VectorizedArrayType> &,
541 VectorTypeOut &,
542 const VectorTypeIn &,
543 const std::pair<unsigned int, unsigned int> &,
544 const bool)> &boundary_operation,
545 VectorTypeOut &dst,
546 const VectorTypeIn &src,
547 const bool zero_dst_vector = false) const
548 {
549 const auto ebd_cell_operation = [&](const auto &matrix_free,
550 auto &dst,
551 const auto &src,
552 const auto &range) {
553 const auto category = matrix_free.get_cell_range_category(range);
554
555 if (category != fe_index_valid)
556 return;
557
558 cell_operation(matrix_free, dst, src, range);
559 };
560
561 const auto ebd_internal_or_boundary_face_operation =
562 [&](const auto &matrix_free,
563 auto &dst,
564 const auto &src,
565 const auto &range) {
566 const auto category = matrix_free.get_face_range_category(range);
567
568 const unsigned int type =
569 static_cast<unsigned int>(category.first == fe_index_valid) +
570 static_cast<unsigned int>(category.second == fe_index_valid);
571
572 if (type == 0)
573 return; // deactivated face -> nothing to do
574
575 if (type == 1) // boundary face
576 boundary_operation(
577 matrix_free, dst, src, range, category.first == fe_index_valid);
578 else if (type == 2) // internal face
579 face_operation(matrix_free, dst, src, range);
580 };
581
582 matrix_free->template loop<VectorTypeOut, VectorTypeIn>(
583 ebd_cell_operation,
584 ebd_internal_or_boundary_face_operation,
585 ebd_internal_or_boundary_face_operation,
586 dst,
587 src,
588 zero_dst_vector);
589 }
590
591 private:
595 SmartPointer<const MatrixFree<dim, Number, VectorizedArrayType>>
596 matrix_free;
597
601 unsigned int fe_index_valid;
602 };
603
604 // implementations
605
606#ifndef DOXYGEN
607
608 template <int dim, typename AdditionalData>
609 void
610 categorize_by_boundary_ids(const Triangulation<dim> &tria,
611 AdditionalData &additional_data)
612 {
613 // ... determine if we are on an active or a multigrid level
614 const unsigned int level = additional_data.mg_level;
615 const bool is_mg = (level != numbers::invalid_unsigned_int);
616
617 // ... create empty list for the category of each cell
618 if (is_mg)
619 additional_data.cell_vectorization_category.assign(
620 std::distance(tria.begin(level), tria.end(level)), 0);
621 else
622 additional_data.cell_vectorization_category.assign(tria.n_active_cells(),
623 0);
624
625 // ... set up scaling factor
626 std::vector<unsigned int> factors(GeometryInfo<dim>::faces_per_cell);
627
628 auto bids = tria.get_boundary_ids();
629 std::sort(bids.begin(), bids.end());
630
631 {
632 unsigned int n_bids = bids.size() + 1;
633 int offset = 1;
634 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
635 i++, offset = offset * n_bids)
636 factors[i] = offset;
637 }
638
639 const auto to_category = [&](const auto &cell) {
640 unsigned int c_num = 0;
641 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
642 {
643 const auto face = cell->face(i);
644 if (face->at_boundary() && !cell->has_periodic_neighbor(i))
645 c_num +=
646 factors[i] * (1 + std::distance(bids.begin(),
647 std::find(bids.begin(),
648 bids.end(),
649 face->boundary_id())));
650 }
651 return c_num;
652 };
653
654 if (!is_mg)
655 {
656 for (auto cell = tria.begin_active(); cell != tria.end(); ++cell)
657 {
658 if (cell->is_locally_owned())
659 additional_data
660 .cell_vectorization_category[cell->active_cell_index()] =
661 to_category(cell);
662 }
663 }
664 else
665 {
666 for (auto cell = tria.begin(level); cell != tria.end(level); ++cell)
667 {
668 if (cell->is_locally_owned_on_level())
669 additional_data.cell_vectorization_category[cell->index()] =
670 to_category(cell);
671 }
672 }
673
674 // ... finalize set up of matrix_free
675 additional_data.hold_all_faces_to_owned_cells = true;
676 additional_data.cell_vectorization_categories_strict = true;
677 additional_data.mapping_update_flags_faces_by_cells =
678 additional_data.mapping_update_flags_inner_faces |
679 additional_data.mapping_update_flags_boundary_faces;
680 }
681
682 namespace internal
683 {
684 template <typename Number>
685 struct LocalCSR
686 {
687 std::vector<unsigned int> row_lid_to_gid;
688 std::vector<unsigned int> row;
689 std::vector<unsigned int> col;
690 std::vector<Number> val;
691
692 std::vector<unsigned int> inverse_lookup_rows;
693 std::vector<std::pair<unsigned int, unsigned int>> inverse_lookup_origins;
694 };
695
696 template <typename FEEvaluationType, bool is_face>
697 class ComputeDiagonalHelper
698 {
699 public:
700 using Number = typename FEEvaluationType::number_type;
701 using VectorizedArrayType = typename FEEvaluationType::NumberType;
702
703 static const unsigned int dim = FEEvaluationType::dimension;
704 static const unsigned int n_components = FEEvaluationType::n_components;
705 static const unsigned int n_lanes = VectorizedArrayType::size();
706
707 ComputeDiagonalHelper()
708 : phi(nullptr)
709 , dofs_per_component(0)
710 {}
711
712 ComputeDiagonalHelper(const ComputeDiagonalHelper &)
713 : phi(nullptr)
714 , dofs_per_component(0)
715 {}
716
717 void
718 initialize(FEEvaluationType &phi)
719 {
720 // if we are in hp mode and the number of unknowns changed, we must
721 // clear the map of entries
722 if (dofs_per_component != phi.dofs_per_component)
723 {
724 locally_relevant_constraints_hn_map.clear();
725 dofs_per_component = phi.dofs_per_component;
726 }
727 this->phi = &phi;
728 }
729
730 void
731 reinit(const unsigned int cell)
732 {
733 this->phi->reinit(cell);
734
735 // STEP 1: get relevant information from FEEvaluation
736 const auto &dof_info = phi->get_dof_info();
737 const unsigned int n_fe_components = dof_info.start_components.back();
738 const auto &matrix_free = phi->get_matrix_free();
739
740 // if we have a block vector with components with the same DoFHandler,
741 // each component is described with same set of constraints and
742 // we consider the shift in components only during access of the global
743 // vector
744 const unsigned int first_selected_component =
745 n_fe_components == 1 ? 0 : phi->get_first_selected_component();
746
747 this->n_lanes_filled =
748 is_face ? matrix_free.n_active_entries_per_face_batch(cell) :
749 matrix_free.n_active_entries_per_cell_batch(cell);
750
751 // STEP 2: setup CSR storage of transposed locally-relevant
752 // constraint matrix
753
754 // (constrained local index, global index of dof
755 // constraints, weight)
756 std::vector<std::tuple<unsigned int, unsigned int, Number>>
757 locally_relevant_constraints, locally_relevant_constraints_tmp;
758 locally_relevant_constraints.reserve(phi->dofs_per_cell);
759 std::vector<unsigned int> constraint_position;
760 std::vector<unsigned char> is_constrained_hn;
761
762 const std::array<unsigned int, n_lanes> &cells =
763 this->phi->get_cell_ids();
764
765 for (unsigned int v = 0; v < n_lanes_filled; ++v)
766 {
767 Assert(cells[v] != numbers::invalid_unsigned_int,
768 ExcInternalError());
769
770 const unsigned int *dof_indices;
771 unsigned int index_indicators, next_index_indicators;
772
773 const unsigned int start =
774 cells[v] * n_fe_components + first_selected_component;
775 dof_indices =
776 dof_info.dof_indices.data() + dof_info.row_starts[start].first;
777 index_indicators = dof_info.row_starts[start].second;
778 next_index_indicators = dof_info.row_starts[start + 1].second;
779
780 // STEP 2a: setup locally-relevant constraint matrix in a
781 // coordinate list (COO)
782 locally_relevant_constraints.clear();
783
784 if (n_components == 1 || n_fe_components == 1)
785 {
786 unsigned int ind_local = 0;
787 for (; index_indicators != next_index_indicators;
788 ++index_indicators, ++ind_local)
789 {
790 const std::pair<unsigned short, unsigned short> indicator =
791 dof_info.constraint_indicator[index_indicators];
792
793 for (unsigned int j = 0; j < indicator.first;
794 ++j, ++ind_local)
795 locally_relevant_constraints.emplace_back(ind_local,
796 dof_indices[j],
797 1.0);
798
799 dof_indices += indicator.first;
800
801 const Number *data_val =
802 matrix_free.constraint_pool_begin(indicator.second);
803 const Number *end_pool =
804 matrix_free.constraint_pool_end(indicator.second);
805
806 for (; data_val != end_pool; ++data_val, ++dof_indices)
807 locally_relevant_constraints.emplace_back(ind_local,
808 *dof_indices,
809 *data_val);
810 }
811
812 AssertIndexRange(ind_local, dofs_per_component + 1);
813
814 for (; ind_local < dofs_per_component;
815 ++dof_indices, ++ind_local)
816 locally_relevant_constraints.emplace_back(ind_local,
817 *dof_indices,
818 1.0);
819 }
820 else
821 {
822 // case with vector-valued finite elements where all
823 // components are included in one single vector. Assumption:
824 // first come all entries to the first component, then all
825 // entries to the second one, and so on. This is ensured by
826 // the way MatrixFree reads out the indices.
827 for (unsigned int comp = 0; comp < n_components; ++comp)
828 {
829 unsigned int ind_local = 0;
830
831 // check whether there is any constraint on the current
832 // cell
833 for (; index_indicators != next_index_indicators;
834 ++index_indicators, ++ind_local)
835 {
836 const std::pair<unsigned short, unsigned short>
837 indicator =
838 dof_info.constraint_indicator[index_indicators];
839
840 // run through values up to next constraint
841 for (unsigned int j = 0; j < indicator.first;
842 ++j, ++ind_local)
843 locally_relevant_constraints.emplace_back(
844 comp * dofs_per_component + ind_local,
845 dof_indices[j],
846 1.0);
847 dof_indices += indicator.first;
848
849 const Number *data_val =
850 matrix_free.constraint_pool_begin(indicator.second);
851 const Number *end_pool =
852 matrix_free.constraint_pool_end(indicator.second);
853
854 for (; data_val != end_pool; ++data_val, ++dof_indices)
855 locally_relevant_constraints.emplace_back(
856 comp * dofs_per_component + ind_local,
857 *dof_indices,
858 *data_val);
859 }
860
861 AssertIndexRange(ind_local, dofs_per_component + 1);
862
863 // get the dof values past the last constraint
864 for (; ind_local < dofs_per_component;
865 ++dof_indices, ++ind_local)
866 locally_relevant_constraints.emplace_back(
867 comp * dofs_per_component + ind_local,
868 *dof_indices,
869 1.0);
870
871 if (comp + 1 < n_components)
872 next_index_indicators =
873 dof_info.row_starts[start + comp + 2].second;
874 }
875 }
876
877 // we only need partial sortedness for the algorithm below in that
878 // all entries for a particular row must be adjacent. this is
879 // ensured by the way we fill the field, but check it again
880 for (unsigned int i = 1; i < locally_relevant_constraints.size();
881 ++i)
882 Assert(std::get<0>(locally_relevant_constraints[i]) >=
883 std::get<0>(locally_relevant_constraints[i - 1]),
884 ExcInternalError());
885
886 // STEP 2c: apply hanging-node constraints
887 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
888 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
889 dof_info.hanging_node_constraint_masks_comp
890 [phi->get_active_fe_index()][first_selected_component])
891 {
892 const auto mask =
893 dof_info.hanging_node_constraint_masks[cells[v]];
894
895 // cell has hanging nodes
896 if (mask != ::internal::MatrixFreeFunctions::
897 unconstrained_compressed_constraint_kind)
898 {
899 // check if hanging node internpolation matrix has been set
900 // up
901 if (locally_relevant_constraints_hn_map.find(mask) ==
902 locally_relevant_constraints_hn_map.end())
903 fill_constraint_type_into_map(mask);
904
905 const auto &locally_relevant_constraints_hn =
906 locally_relevant_constraints_hn_map[mask];
907
908 locally_relevant_constraints_tmp.clear();
909 if (locally_relevant_constraints_tmp.capacity() <
910 locally_relevant_constraints.size())
911 locally_relevant_constraints_tmp.reserve(
912 locally_relevant_constraints.size() +
913 locally_relevant_constraints_hn.size());
914
915 // combine with other constraints: to avoid binary
916 // searches, we first build a list of where constraints
917 // are pointing to, and then merge the two lists
918 constraint_position.assign(phi->dofs_per_cell,
919 numbers::invalid_unsigned_int);
920 for (auto &a : locally_relevant_constraints)
921 if (constraint_position[std::get<0>(a)] ==
922 numbers::invalid_unsigned_int)
923 constraint_position[std::get<0>(a)] =
924 std::distance(locally_relevant_constraints.data(),
925 &a);
926 is_constrained_hn.assign(phi->dofs_per_cell, false);
927 for (auto &hn : locally_relevant_constraints_hn)
928 is_constrained_hn[std::get<0>(hn)] = 1;
929
930 // not constrained from hanging nodes
931 for (const auto &a : locally_relevant_constraints)
932 if (is_constrained_hn[std::get<0>(a)] == 0)
933 locally_relevant_constraints_tmp.push_back(a);
934
935 // dof is constrained by hanging nodes: build transitive
936 // closure
937 for (const auto &hn : locally_relevant_constraints_hn)
938 if (constraint_position[std::get<1>(hn)] !=
939 numbers::invalid_unsigned_int)
940 {
941 AssertIndexRange(constraint_position[std::get<1>(hn)],
942 locally_relevant_constraints.size());
943 auto other = locally_relevant_constraints.begin() +
944 constraint_position[std::get<1>(hn)];
945 AssertDimension(std::get<0>(*other), std::get<1>(hn));
946
947 for (; other != locally_relevant_constraints.end() &&
948 std::get<0>(*other) == std::get<1>(hn);
949 ++other)
950 locally_relevant_constraints_tmp.emplace_back(
951 std::get<0>(hn),
952 std::get<1>(*other),
953 std::get<2>(hn) * std::get<2>(*other));
954 }
955
956 std::swap(locally_relevant_constraints,
957 locally_relevant_constraints_tmp);
958 }
959 }
960
961 // STEP 2d: transpose COO
962 std::sort(locally_relevant_constraints.begin(),
963 locally_relevant_constraints.end(),
964 [](const auto &a, const auto &b) {
965 if (std::get<1>(a) < std::get<1>(b))
966 return true;
967 return (std::get<1>(a) == std::get<1>(b)) &&
968 (std::get<0>(a) < std::get<0>(b));
969 });
970
971 // STEP 2e: translate COO to CRS
972 auto &c_pool = c_pools[v];
973 {
974 c_pool.row_lid_to_gid.clear();
975 c_pool.row.clear();
976 c_pool.row.push_back(0);
977 c_pool.col.clear();
978 c_pool.val.clear();
979
980 if (locally_relevant_constraints.size() > 0)
981 c_pool.row_lid_to_gid.emplace_back(
982 std::get<1>(locally_relevant_constraints.front()));
983 for (const auto &j : locally_relevant_constraints)
984 {
985 if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
986 {
987 c_pool.row_lid_to_gid.push_back(std::get<1>(j));
988 c_pool.row.push_back(c_pool.val.size());
989 }
990
991 c_pool.col.emplace_back(std::get<0>(j));
992 c_pool.val.emplace_back(std::get<2>(j));
993 }
994
995 if (c_pool.val.size() > 0)
996 c_pool.row.push_back(c_pool.val.size());
997
998 c_pool.inverse_lookup_rows.clear();
999 c_pool.inverse_lookup_rows.resize(1 + phi->dofs_per_cell);
1000 for (const unsigned int i : c_pool.col)
1001 c_pool.inverse_lookup_rows[1 + i]++;
1002 // transform to offsets
1003 std::partial_sum(c_pool.inverse_lookup_rows.begin(),
1004 c_pool.inverse_lookup_rows.end(),
1005 c_pool.inverse_lookup_rows.begin());
1006 AssertDimension(c_pool.inverse_lookup_rows.back(),
1007 c_pool.col.size());
1008
1009 c_pool.inverse_lookup_origins.resize(c_pool.col.size());
1010 std::vector<unsigned int> inverse_lookup_count(
1011 phi->dofs_per_cell);
1012 for (unsigned int row = 0; row < c_pool.row.size() - 1; ++row)
1013 for (unsigned int col = c_pool.row[row];
1014 col < c_pool.row[row + 1];
1015 ++col)
1016 {
1017 const unsigned int index = c_pool.col[col];
1018 c_pool.inverse_lookup_origins
1019 [c_pool.inverse_lookup_rows[index] +
1020 inverse_lookup_count[index]] = std::make_pair(row, col);
1021 ++inverse_lookup_count[index];
1022 }
1023 }
1024 }
1025
1026 // STEP 3: compute element matrix A_e, apply
1027 // locally-relevant constraints C_e^T * A_e * C_e, and get the
1028 // the diagonal entry
1029 // (C_e^T * A_e * C_e)(i,i)
1030 // or
1031 // C_e^T(i,:) * A_e * C_e(:,i).
1032 //
1033 // Since, we compute the element matrix column-by-column and as a
1034 // result never actually have the full element matrix, we actually
1035 // perform following steps:
1036 // 1) loop over all columns of the element matrix
1037 // a) compute column i
1038 // b) compute for each j (rows of C_e^T):
1039 // (C_e^T(j,:) * A_e(:,i)) * C_e(i,j)
1040 // or
1041 // (C_e^T(j,:) * A_e(:,i)) * C_e^T(j,i)
1042 // This gives a contribution the j-th entry of the
1043 // locally-relevant diagonal and comprises the multiplication
1044 // by the locally-relevant constraint matrix from the left and
1045 // the right. There is no contribution to the j-th vector
1046 // entry if the j-th row of C_e^T is empty or C_e^T(j,i) is
1047 // zero.
1048
1049 // set size locally-relevant diagonal
1050 for (unsigned int v = 0; v < n_lanes_filled; ++v)
1051 diagonals_local_constrained[v].assign(
1052 c_pools[v].row_lid_to_gid.size() *
1053 (n_fe_components == 1 ? n_components : 1),
1054 Number(0.0));
1055
1056 // check if fast path can be taken via FEEvaluation
1057 bool use_fast_path = true;
1058
1059 for (unsigned int v = 0; v < n_lanes_filled; ++v)
1060 {
1061 auto &c_pool = c_pools[v];
1062
1063 for (unsigned int i = 0; i < c_pool.row.size() - 1; ++i)
1064 {
1065 if ((c_pool.row[i + 1] - c_pool.row[i]) > 1)
1066 {
1067 use_fast_path = false;
1068 break;
1069 }
1070 else if (((c_pool.row[i + 1] - c_pool.row[i]) == 1) &&
1071 (c_pool.val[c_pool.row[i]] != 1.0))
1072 {
1073 use_fast_path = false;
1074 break;
1075 }
1076 }
1077
1078 if (use_fast_path == false)
1079 break;
1080 }
1081
1082 if (use_fast_path)
1083 {
1084 temp_values.resize(phi->dofs_per_cell);
1085 return;
1086 }
1087 else
1088 {
1089 temp_values.clear();
1090 }
1091 }
1092
1093 void
1094 fill_constraint_type_into_map(
1095 const ::internal::MatrixFreeFunctions::compressed_constraint_kind
1096 mask)
1097 {
1098 auto &constraints_hn = locally_relevant_constraints_hn_map[mask];
1099
1100 // assume that we constrain one face, i.e., (fe_degree + 1)^(dim-1)
1101 // unknowns - we might have more or less entries, but this is a good
1102 // first guess
1103 const unsigned int degree =
1104 phi->get_shape_info().data.front().fe_degree;
1105 constraints_hn.reserve(Utilities::pow(degree + 1, dim - 1));
1106
1107 // 1) collect hanging-node constraints for cell assuming
1108 // scalar finite element
1109 values_dofs.resize(dofs_per_component);
1110 std::array<
1111 ::internal::MatrixFreeFunctions::compressed_constraint_kind,
1112 VectorizedArrayType::size()>
1113 constraint_mask;
1114 constraint_mask.fill(::internal::MatrixFreeFunctions::
1115 unconstrained_compressed_constraint_kind);
1116 constraint_mask[0] = mask;
1117
1118 for (unsigned int i = 0; i < dofs_per_component; ++i)
1119 {
1120 for (unsigned int j = 0; j < dofs_per_component; ++j)
1121 values_dofs[j] = VectorizedArrayType();
1122 values_dofs[i] = Number(1);
1123
1124 ::internal::FEEvaluationHangingNodesFactory<
1125 dim,
1126 Number,
1127 VectorizedArrayType>::apply(1,
1128 degree,
1129 phi->get_shape_info(),
1130 false,
1131 constraint_mask,
1132 values_dofs.data());
1133
1134 const Number tolerance =
1135 std::max(Number(1e-12),
1136 std::numeric_limits<Number>::epsilon() * 16);
1137 for (unsigned int j = 0; j < dofs_per_component; ++j)
1138 if (std::abs(values_dofs[j][0]) > tolerance &&
1139 (j != i ||
1140 std::abs(values_dofs[j][0] - Number(1)) > tolerance))
1141 constraints_hn.emplace_back(j, i, values_dofs[j][0]);
1142 }
1143
1144 // 2) extend for multiple components
1145 const unsigned int n_hn_constraints = constraints_hn.size();
1146 constraints_hn.resize(n_hn_constraints * n_components);
1147
1148 for (unsigned int c = 1; c < n_components; ++c)
1149 for (unsigned int i = 0; i < n_hn_constraints; ++i)
1150 constraints_hn[c * n_hn_constraints + i] = std::make_tuple(
1151 std::get<0>(constraints_hn[i]) + c * dofs_per_component,
1152 std::get<1>(constraints_hn[i]) + c * dofs_per_component,
1153 std::get<2>(constraints_hn[i]));
1154 }
1155
1156 void
1157 prepare_basis_vector(const unsigned int i)
1158 {
1159 this->i = i;
1160
1161 // compute i-th column of element stiffness matrix:
1162 // this could be simply performed as done at the moment with
1163 // matrix-free operator evaluation applied to a ith-basis vector
1164 VectorizedArrayType *dof_values = phi->begin_dof_values();
1165 for (const unsigned int j : phi->dof_indices())
1166 dof_values[j] = VectorizedArrayType();
1167 dof_values[i] = Number(1);
1168 }
1169
1170 void
1171 zero_basis_vector()
1172 {
1173 VectorizedArrayType *dof_values = phi->begin_dof_values();
1174 for (const unsigned int j : phi->dof_indices())
1175 dof_values[j] = VectorizedArrayType();
1176 }
1177
1178 void
1179 submit()
1180 {
1181 if (!temp_values.empty())
1182 {
1183 temp_values[i] = phi->begin_dof_values()[i];
1184 return;
1185 }
1186
1187 // if we have a block vector with components with the same DoFHandler,
1188 // we need to figure out which component and which DoF within the
1189 // component are we currently considering
1190 const unsigned int n_fe_components =
1191 phi->get_dof_info().start_components.back();
1192 const unsigned int comp =
1193 n_fe_components == 1 ? i / dofs_per_component : 0;
1194 const unsigned int i_comp =
1195 n_fe_components == 1 ? (i % dofs_per_component) : i;
1196
1197 // apply local constraint matrix from left and from right:
1198 // loop over all rows of transposed constrained matrix
1199 for (unsigned int v = 0; v < n_lanes_filled; ++v)
1200 {
1201 const auto &c_pool = c_pools[v];
1202
1203 for (unsigned int jj = c_pool.inverse_lookup_rows[i_comp];
1204 jj < c_pool.inverse_lookup_rows[i_comp + 1];
1205 ++jj)
1206 {
1207 const unsigned int j = c_pool.inverse_lookup_origins[jj].first;
1208 // apply constraint matrix from the left
1209 Number temp = 0.0;
1210 for (unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
1211 temp += c_pool.val[k] *
1212 phi->begin_dof_values()[comp * dofs_per_component +
1213 c_pool.col[k]][v];
1214
1215 // apply constraint matrix from the right
1216 diagonals_local_constrained
1217 [v][j + comp * c_pools[v].row_lid_to_gid.size()] +=
1218 temp * c_pool.val[c_pool.inverse_lookup_origins[jj].second];
1219 }
1220 }
1221 }
1222
1223 template <typename VectorType>
1224 inline void
1225 distribute_local_to_global(
1226 std::array<VectorType *, n_components> &diagonal_global)
1227 {
1228 if (!temp_values.empty())
1229 {
1230 for (unsigned int j = 0; j < temp_values.size(); ++j)
1231 phi->begin_dof_values()[j] = temp_values[j];
1232
1233 phi->distribute_local_to_global(diagonal_global);
1234
1235 return;
1236 }
1237
1238 // STEP 4: assembly results: add into global vector
1239 const unsigned int n_fe_components =
1240 phi->get_dof_info().start_components.back();
1241
1242 for (unsigned int v = 0; v < n_lanes_filled; ++v)
1243 // if we have a block vector with components with the same
1244 // DoFHandler, we need to loop over all components manually and
1245 // need to apply the correct shift
1246 for (unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
1247 for (unsigned int comp = 0;
1248 comp < (n_fe_components == 1 ?
1249 static_cast<unsigned int>(n_components) :
1250 1);
1251 ++comp)
1252 ::internal::vector_access_add(
1253 *diagonal_global[n_fe_components == 1 ? comp : 0],
1254 c_pools[v].row_lid_to_gid[j],
1255 diagonals_local_constrained
1256 [v][j + comp * c_pools[v].row_lid_to_gid.size()]);
1257 }
1258
1259 bool
1260 use_fast_path() const
1261 {
1262 return !temp_values.empty();
1263 }
1264
1265 private:
1266 FEEvaluationType *phi;
1267
1268 unsigned int dofs_per_component;
1269
1270 unsigned int i;
1271
1272 unsigned int n_lanes_filled;
1273
1274 std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
1275
1276 // local storage: buffer so that we access the global vector once
1277 // note: may be larger then dofs_per_cell in the presence of
1278 // constraints!
1279 std::array<std::vector<Number>, n_lanes> diagonals_local_constrained;
1280
1281 std::map<
1282 ::internal::MatrixFreeFunctions::compressed_constraint_kind,
1283 std::vector<std::tuple<unsigned int, unsigned int, Number>>>
1284 locally_relevant_constraints_hn_map;
1285
1286 // scratch array
1287 AlignedVector<VectorizedArrayType> temp_values;
1288 AlignedVector<VectorizedArrayType> values_dofs;
1289 };
1290
1291 template <bool is_face,
1292 int dim,
1293 typename Number,
1294 typename VectorizedArrayType>
1295 bool
1296 is_fe_nothing(
1297 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1298 const std::pair<unsigned int, unsigned int> &range,
1299 const unsigned int dof_no,
1300 const unsigned int quad_no,
1301 const unsigned int first_selected_component,
1302 const unsigned int fe_degree,
1303 const unsigned int n_q_points_1d,
1304 const bool is_interior_face = true)
1305 {
1306 const unsigned int static_n_q_points =
1307 is_face ? Utilities::pow(n_q_points_1d, dim - 1) :
1308 Utilities::pow(n_q_points_1d, dim);
1309
1310 unsigned int active_fe_index = 0;
1311 if (!is_face)
1312 active_fe_index = matrix_free.get_cell_active_fe_index(range);
1313 else if (is_interior_face)
1314 active_fe_index = matrix_free.get_face_range_category(range).first;
1315 else
1316 active_fe_index = matrix_free.get_face_range_category(range).second;
1317
1318 const auto init_data = ::internal::
1319 extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
1320 matrix_free,
1321 dof_no,
1322 first_selected_component,
1323 quad_no,
1324 fe_degree,
1325 static_n_q_points,
1326 active_fe_index,
1327 numbers::invalid_unsigned_int /*active_quad_index*/,
1328 numbers::invalid_unsigned_int /*face_type*/);
1329
1330 return init_data.shape_info->dofs_per_component_on_cell == 0;
1331 }
1332
1333 } // namespace internal
1334
1335 template <int dim,
1336 int fe_degree,
1337 int n_q_points_1d,
1338 int n_components,
1339 typename Number,
1340 typename VectorizedArrayType,
1341 typename VectorType>
1342 void
1343 compute_diagonal(
1344 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1345 VectorType &diagonal_global,
1346 const std::function<void(FEEvaluation<dim,
1347 fe_degree,
1348 n_q_points_1d,
1349 n_components,
1350 Number,
1351 VectorizedArrayType> &)>
1352 &cell_operation,
1353 const unsigned int dof_no,
1354 const unsigned int quad_no,
1355 const unsigned int first_selected_component)
1356 {
1357 compute_diagonal<dim,
1358 fe_degree,
1359 n_q_points_1d,
1360 n_components,
1361 Number,
1362 VectorizedArrayType,
1363 VectorType>(matrix_free,
1364 diagonal_global,
1365 cell_operation,
1366 {},
1367 {},
1368 dof_no,
1369 quad_no,
1370 first_selected_component);
1371 }
1372
1373 template <typename CLASS,
1374 int dim,
1375 int fe_degree,
1376 int n_q_points_1d,
1377 int n_components,
1378 typename Number,
1379 typename VectorizedArrayType,
1380 typename VectorType>
1381 void
1382 compute_diagonal(
1383 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1384 VectorType &diagonal_global,
1385 void (CLASS::*cell_operation)(FEEvaluation<dim,
1386 fe_degree,
1387 n_q_points_1d,
1388 n_components,
1389 Number,
1390 VectorizedArrayType> &) const,
1391 const CLASS *owning_class,
1392 const unsigned int dof_no,
1393 const unsigned int quad_no,
1394 const unsigned int first_selected_component)
1395 {
1396 compute_diagonal<dim,
1397 fe_degree,
1398 n_q_points_1d,
1399 n_components,
1400 Number,
1401 VectorizedArrayType,
1402 VectorType>(
1403 matrix_free,
1404 diagonal_global,
1405 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
1406 dof_no,
1407 quad_no,
1408 first_selected_component);
1409 }
1410
1411 template <int dim,
1412 int fe_degree,
1413 int n_q_points_1d,
1414 int n_components,
1415 typename Number,
1416 typename VectorizedArrayType,
1417 typename VectorType>
1418 void
1419 compute_diagonal(
1420 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1421 VectorType &diagonal_global,
1422 const std::function<void(FEEvaluation<dim,
1423 fe_degree,
1424 n_q_points_1d,
1425 n_components,
1426 Number,
1427 VectorizedArrayType> &)>
1428 &cell_operation,
1429 const std::function<void(FEFaceEvaluation<dim,
1430 fe_degree,
1431 n_q_points_1d,
1432 n_components,
1433 Number,
1434 VectorizedArrayType> &,
1435 FEFaceEvaluation<dim,
1436 fe_degree,
1437 n_q_points_1d,
1438 n_components,
1439 Number,
1440 VectorizedArrayType> &)>
1441 &face_operation,
1442 const std::function<void(FEFaceEvaluation<dim,
1443 fe_degree,
1444 n_q_points_1d,
1445 n_components,
1446 Number,
1447 VectorizedArrayType> &)>
1448 &boundary_operation,
1449 const unsigned int dof_no,
1450 const unsigned int quad_no,
1451 const unsigned int first_selected_component)
1452 {
1453 int dummy = 0;
1454
1455 std::array<typename ::internal::BlockVectorSelector<
1456 VectorType,
1457 IsBlockVector<VectorType>::value>::BaseVectorType *,
1458 n_components>
1459 diagonal_global_components;
1460
1461 for (unsigned int d = 0; d < n_components; ++d)
1462 diagonal_global_components[d] = ::internal::
1463 BlockVectorSelector<VectorType, IsBlockVector<VectorType>::value>::
1464 get_vector_component(diagonal_global, d + first_selected_component);
1465
1466 const auto &dof_info = matrix_free.get_dof_info(dof_no);
1467
1468 if (dof_info.start_components.back() == 1)
1469 for (unsigned int comp = 0; comp < n_components; ++comp)
1470 {
1471 Assert(diagonal_global_components[comp] != nullptr,
1472 ExcMessage("The finite element underlying this FEEvaluation "
1473 "object is scalar, but you requested " +
1474 std::to_string(n_components) +
1475 " components via the template argument in "
1476 "FEEvaluation. In that case, you must pass an "
1477 "std::vector<VectorType> or a BlockVector to " +
1478 "read_dof_values and distribute_local_to_global."));
1479 ::internal::check_vector_compatibility(
1480 *diagonal_global_components[comp], matrix_free, dof_info);
1481 }
1482 else
1483 {
1484 ::internal::check_vector_compatibility(
1485 *diagonal_global_components[0], matrix_free, dof_info);
1486 }
1487
1488 using Helper =
1489 internal::ComputeDiagonalHelper<FEEvaluation<dim,
1490 fe_degree,
1491 n_q_points_1d,
1492 n_components,
1493 Number,
1494 VectorizedArrayType>,
1495 false>;
1496
1497 using HelperFace =
1498 internal::ComputeDiagonalHelper<FEFaceEvaluation<dim,
1499 fe_degree,
1500 n_q_points_1d,
1501 n_components,
1502 Number,
1503 VectorizedArrayType>,
1504 true>;
1505
1506 Threads::ThreadLocalStorage<Helper> scratch_data;
1507 Threads::ThreadLocalStorage<HelperFace> scratch_data_m;
1508 Threads::ThreadLocalStorage<HelperFace> scratch_data_p;
1509
1510 const auto cell_operation_wrapped =
1511 [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1512 VectorType &,
1513 const int &,
1514 const std::pair<unsigned int, unsigned int> &range) {
1515 if (!cell_operation)
1516 return;
1517
1518 // shortcut for FE_Nothing cells
1519 if (internal::is_fe_nothing<false>(matrix_free,
1520 range,
1521 dof_no,
1522 quad_no,
1523 first_selected_component,
1524 fe_degree,
1525 n_q_points_1d))
1526 return;
1527
1528 Helper &helper = scratch_data.get();
1529
1530 FEEvaluation<dim,
1531 fe_degree,
1532 n_q_points_1d,
1533 n_components,
1534 Number,
1535 VectorizedArrayType>
1536 phi(matrix_free, range, dof_no, quad_no, first_selected_component);
1537 helper.initialize(phi);
1538
1539 for (unsigned int cell = range.first; cell < range.second; ++cell)
1540 {
1541 helper.reinit(cell);
1542
1543 for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
1544 {
1545 helper.prepare_basis_vector(i);
1546 cell_operation(phi);
1547 helper.submit();
1548 }
1549
1550 helper.distribute_local_to_global(diagonal_global_components);
1551 }
1552 };
1553
1554 const auto face_operation_wrapped =
1555 [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1556 VectorType &,
1557 const int &,
1558 const std::pair<unsigned int, unsigned int> &range) {
1559 if (!face_operation)
1560 return;
1561
1562 // shortcut for faces containing purely FE_Nothing
1563 if (internal::is_fe_nothing<true>(matrix_free,
1564 range,
1565 dof_no,
1566 quad_no,
1567 first_selected_component,
1568 fe_degree,
1569 n_q_points_1d,
1570 true) ||
1571 internal::is_fe_nothing<true>(matrix_free,
1572 range,
1573 dof_no,
1574 quad_no,
1575 first_selected_component,
1576 fe_degree,
1577 n_q_points_1d,
1578 false))
1579 return;
1580
1581 HelperFace &helper_m = scratch_data_m.get();
1582 HelperFace &helper_p = scratch_data_p.get();
1583
1584 FEFaceEvaluation<dim,
1585 fe_degree,
1586 n_q_points_1d,
1587 n_components,
1588 Number,
1589 VectorizedArrayType>
1590 phi_m(matrix_free,
1591 range,
1592 true,
1593 dof_no,
1594 quad_no,
1595 first_selected_component);
1596 FEFaceEvaluation<dim,
1597 fe_degree,
1598 n_q_points_1d,
1599 n_components,
1600 Number,
1601 VectorizedArrayType>
1602 phi_p(matrix_free,
1603 range,
1604 false,
1605 dof_no,
1606 quad_no,
1607 first_selected_component);
1608
1609 helper_m.initialize(phi_m);
1610 helper_p.initialize(phi_p);
1611
1612 for (unsigned int face = range.first; face < range.second; ++face)
1613 {
1614 helper_m.reinit(face);
1615 helper_p.reinit(face);
1616
1617 // make check only if both adjacent cells have DoFs
1618 Assert(helper_m.use_fast_path() && helper_p.use_fast_path(),
1619 ExcNotImplemented());
1620
1621 for (unsigned int i = 0; i < phi_m.dofs_per_cell; ++i)
1622 {
1623 helper_m.prepare_basis_vector(i);
1624 helper_p.zero_basis_vector();
1625 face_operation(phi_m, phi_p);
1626 helper_m.submit();
1627 }
1628
1629 helper_m.distribute_local_to_global(diagonal_global_components);
1630
1631 for (unsigned int i = 0; i < phi_p.dofs_per_cell; ++i)
1632 {
1633 helper_m.zero_basis_vector();
1634 helper_p.prepare_basis_vector(i);
1635 face_operation(phi_m, phi_p);
1636 helper_p.submit();
1637 }
1638
1639 helper_p.distribute_local_to_global(diagonal_global_components);
1640 }
1641 };
1642
1643 const auto boundary_operation_wrapped =
1644 [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1645 VectorType &,
1646 const int &,
1647 const std::pair<unsigned int, unsigned int> &range) {
1648 if (!boundary_operation)
1649 return;
1650
1651 // shortcut for faces containing purely FE_Nothing
1652 if (internal::is_fe_nothing<true>(matrix_free,
1653 range,
1654 dof_no,
1655 quad_no,
1656 first_selected_component,
1657 fe_degree,
1658 n_q_points_1d,
1659 true))
1660 return;
1661
1662 HelperFace &helper = scratch_data_m.get();
1663
1664 FEFaceEvaluation<dim,
1665 fe_degree,
1666 n_q_points_1d,
1667 n_components,
1668 Number,
1669 VectorizedArrayType>
1670 phi(matrix_free,
1671 range,
1672 true,
1673 dof_no,
1674 quad_no,
1675 first_selected_component);
1676 helper.initialize(phi);
1677
1678 for (unsigned int face = range.first; face < range.second; ++face)
1679 {
1680 helper.reinit(face);
1681
1682 for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
1683 {
1684 helper.prepare_basis_vector(i);
1685 boundary_operation(phi);
1686 helper.submit();
1687 }
1688
1689 helper.distribute_local_to_global(diagonal_global_components);
1690 }
1691 };
1692
1693 if (face_operation || boundary_operation)
1694 matrix_free.template loop<VectorType, int>(cell_operation_wrapped,
1695 face_operation_wrapped,
1696 boundary_operation_wrapped,
1697 diagonal_global,
1698 dummy,
1699 false);
1700 else
1701 matrix_free.template cell_loop<VectorType, int>(cell_operation_wrapped,
1702 diagonal_global,
1703 dummy,
1704 false);
1705 }
1706
1707 template <typename CLASS,
1708 int dim,
1709 int fe_degree,
1710 int n_q_points_1d,
1711 int n_components,
1712 typename Number,
1713 typename VectorizedArrayType,
1714 typename VectorType>
1715 void
1716 compute_diagonal(
1717 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1718 VectorType &diagonal_global,
1719 void (CLASS::*cell_operation)(FEEvaluation<dim,
1720 fe_degree,
1721 n_q_points_1d,
1722 n_components,
1723 Number,
1724 VectorizedArrayType> &) const,
1725 void (CLASS::*face_operation)(FEFaceEvaluation<dim,
1726 fe_degree,
1727 n_q_points_1d,
1728 n_components,
1729 Number,
1730 VectorizedArrayType> &,
1731 FEFaceEvaluation<dim,
1732 fe_degree,
1733 n_q_points_1d,
1734 n_components,
1735 Number,
1736 VectorizedArrayType> &)
1737 const,
1738 void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
1739 fe_degree,
1740 n_q_points_1d,
1741 n_components,
1742 Number,
1743 VectorizedArrayType> &)
1744 const,
1745 const CLASS *owning_class,
1746 const unsigned int dof_no,
1747 const unsigned int quad_no,
1748 const unsigned int first_selected_component)
1749 {
1750 compute_diagonal<dim,
1751 fe_degree,
1752 n_q_points_1d,
1753 n_components,
1754 Number,
1755 VectorizedArrayType,
1756 VectorType>(
1757 matrix_free,
1758 diagonal_global,
1759 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
1760 [&](auto &phi_m, auto &phi_p) {
1761 (owning_class->*face_operation)(phi_m, phi_p);
1762 },
1763 [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
1764 dof_no,
1765 quad_no,
1766 first_selected_component);
1767 }
1768
1769 namespace internal
1770 {
1775 template <
1776 typename MatrixType,
1777 typename Number,
1778 std::enable_if_t<std::is_same_v<
1779 std::remove_const_t<
1780 std::remove_reference_t<typename MatrixType::value_type>>,
1781 std::remove_const_t<std::remove_reference_t<Number>>>> * = nullptr>
1782 const AffineConstraints<typename MatrixType::value_type> &
1783 create_new_affine_constraints_if_needed(
1784 const MatrixType &,
1785 const AffineConstraints<Number> &constraints,
1786 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>> &)
1787 {
1788 return constraints;
1789 }
1790
1796 template <
1797 typename MatrixType,
1798 typename Number,
1799 std::enable_if_t<!std::is_same_v<
1800 std::remove_const_t<
1801 std::remove_reference_t<typename MatrixType::value_type>>,
1802 std::remove_const_t<std::remove_reference_t<Number>>>> * = nullptr>
1803 const AffineConstraints<typename MatrixType::value_type> &
1804 create_new_affine_constraints_if_needed(
1805 const MatrixType &,
1806 const AffineConstraints<Number> &constraints,
1807 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
1808 &new_constraints)
1809 {
1810 new_constraints =
1811 std::make_unique<AffineConstraints<typename MatrixType::value_type>>();
1812 new_constraints->copy_from(constraints);
1813
1814 return *new_constraints;
1815 }
1816 } // namespace internal
1817
1818 template <int dim,
1819 int fe_degree,
1820 int n_q_points_1d,
1821 int n_components,
1822 typename Number,
1823 typename VectorizedArrayType,
1824 typename MatrixType>
1825 void
1826 compute_matrix(
1827 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1828 const AffineConstraints<Number> &constraints_in,
1829 MatrixType &matrix,
1830 const std::function<void(FEEvaluation<dim,
1831 fe_degree,
1832 n_q_points_1d,
1833 n_components,
1834 Number,
1835 VectorizedArrayType> &)>
1836 &cell_operation,
1837 const unsigned int dof_no,
1838 const unsigned int quad_no,
1839 const unsigned int first_selected_component)
1840 {
1841 compute_matrix<dim,
1842 fe_degree,
1843 n_q_points_1d,
1844 n_components,
1845 Number,
1846 VectorizedArrayType,
1847 MatrixType>(matrix_free,
1848 constraints_in,
1849 matrix,
1850 cell_operation,
1851 {},
1852 {},
1853 dof_no,
1854 quad_no,
1855 first_selected_component);
1856 }
1857
1858 template <typename CLASS,
1859 int dim,
1860 int fe_degree,
1861 int n_q_points_1d,
1862 int n_components,
1863 typename Number,
1864 typename VectorizedArrayType,
1865 typename MatrixType>
1866 void
1867 compute_matrix(
1868 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1869 const AffineConstraints<Number> &constraints,
1870 MatrixType &matrix,
1871 void (CLASS::*cell_operation)(FEEvaluation<dim,
1872 fe_degree,
1873 n_q_points_1d,
1874 n_components,
1875 Number,
1876 VectorizedArrayType> &) const,
1877 const CLASS *owning_class,
1878 const unsigned int dof_no,
1879 const unsigned int quad_no,
1880 const unsigned int first_selected_component)
1881 {
1882 compute_matrix<dim,
1883 fe_degree,
1884 n_q_points_1d,
1885 n_components,
1886 Number,
1887 VectorizedArrayType,
1888 MatrixType>(
1889 matrix_free,
1890 constraints,
1891 matrix,
1892 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
1893 dof_no,
1894 quad_no,
1895 first_selected_component);
1896 }
1897
1898 template <int dim,
1899 int fe_degree,
1900 int n_q_points_1d,
1901 int n_components,
1902 typename Number,
1903 typename VectorizedArrayType,
1904 typename MatrixType>
1905 void
1906 compute_matrix(
1907 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
1908 const AffineConstraints<Number> &constraints_in,
1909 MatrixType &matrix,
1910 const std::function<void(FEEvaluation<dim,
1911 fe_degree,
1912 n_q_points_1d,
1913 n_components,
1914 Number,
1915 VectorizedArrayType> &)>
1916 &cell_operation,
1917 const std::function<void(FEFaceEvaluation<dim,
1918 fe_degree,
1919 n_q_points_1d,
1920 n_components,
1921 Number,
1922 VectorizedArrayType> &,
1923 FEFaceEvaluation<dim,
1924 fe_degree,
1925 n_q_points_1d,
1926 n_components,
1927 Number,
1928 VectorizedArrayType> &)>
1929 &face_operation,
1930 const std::function<void(FEFaceEvaluation<dim,
1931 fe_degree,
1932 n_q_points_1d,
1933 n_components,
1934 Number,
1935 VectorizedArrayType> &)>
1936 &boundary_operation,
1937 const unsigned int dof_no,
1938 const unsigned int quad_no,
1939 const unsigned int first_selected_component)
1940 {
1941 using FEEvalType = FEEvaluation<dim,
1942 fe_degree,
1943 n_q_points_1d,
1944 n_components,
1945 Number,
1946 VectorizedArrayType>;
1947
1948 using FEFaceEvalType = FEFaceEvaluation<dim,
1949 fe_degree,
1950 n_q_points_1d,
1951 n_components,
1952 Number,
1953 VectorizedArrayType>;
1954
1955 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
1956 data_cell;
1957
1958 data_cell.dof_numbers = {dof_no};
1959 data_cell.quad_numbers = {quad_no};
1960 data_cell.first_selected_components = {first_selected_component};
1961 data_cell.batch_type = {0};
1962
1963 data_cell.op_create =
1964 [&](const std::pair<unsigned int, unsigned int> &range) {
1965 std::vector<
1966 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
1967 phi;
1968
1969 if (!internal::is_fe_nothing<false>(matrix_free,
1970 range,
1971 dof_no,
1972 quad_no,
1973 first_selected_component,
1974 fe_degree,
1975 n_q_points_1d))
1976 phi.emplace_back(std::make_unique<FEEvalType>(
1977 matrix_free, range, dof_no, quad_no, first_selected_component));
1978
1979 return phi;
1980 };
1981
1982 data_cell.op_reinit = [](auto &phi, const unsigned batch) {
1983 static_cast<FEEvalType &>(*phi[0]).reinit(batch);
1984 };
1985
1986 if (cell_operation)
1987 data_cell.op_compute = [&](auto &phi) {
1988 cell_operation(static_cast<FEEvalType &>(*phi[0]));
1989 };
1990
1991 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
1992 data_face;
1993
1994 data_face.dof_numbers = {dof_no, dof_no};
1995 data_face.quad_numbers = {quad_no, quad_no};
1996 data_face.first_selected_components = {first_selected_component,
1997 first_selected_component};
1998 data_face.batch_type = {1, 2};
1999
2000 data_face.op_create =
2001 [&](const std::pair<unsigned int, unsigned int> &range) {
2002 std::vector<
2003 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2004 phi;
2005
2006 if (!internal::is_fe_nothing<true>(matrix_free,
2007 range,
2008 dof_no,
2009 quad_no,
2010 first_selected_component,
2011 fe_degree,
2012 n_q_points_1d,
2013 true) &&
2014 !internal::is_fe_nothing<true>(matrix_free,
2015 range,
2016 dof_no,
2017 quad_no,
2018 first_selected_component,
2019 fe_degree,
2020 n_q_points_1d,
2021 false))
2022 {
2023 phi.emplace_back(
2024 std::make_unique<FEFaceEvalType>(matrix_free,
2025 range,
2026 true,
2027 dof_no,
2028 quad_no,
2029 first_selected_component));
2030 phi.emplace_back(
2031 std::make_unique<FEFaceEvalType>(matrix_free,
2032 range,
2033 false,
2034 dof_no,
2035 quad_no,
2036 first_selected_component));
2037 }
2038
2039 return phi;
2040 };
2041
2042 data_face.op_reinit = [](auto &phi, const unsigned batch) {
2043 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
2044 static_cast<FEFaceEvalType &>(*phi[1]).reinit(batch);
2045 };
2046
2047 if (face_operation)
2048 data_face.op_compute = [&](auto &phi) {
2049 face_operation(static_cast<FEFaceEvalType &>(*phi[0]),
2050 static_cast<FEFaceEvalType &>(*phi[1]));
2051 };
2052
2053 internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2054 data_boundary;
2055
2056 data_boundary.dof_numbers = {dof_no};
2057 data_boundary.quad_numbers = {quad_no};
2058 data_boundary.first_selected_components = {first_selected_component};
2059 data_boundary.batch_type = {1};
2060
2061 data_boundary
2062 .op_create = [&](const std::pair<unsigned int, unsigned int> &range) {
2063 std::vector<
2064 std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
2065 phi;
2066
2067 if (!internal::is_fe_nothing<true>(matrix_free,
2068 range,
2069 dof_no,
2070 quad_no,
2071 first_selected_component,
2072 fe_degree,
2073 n_q_points_1d,
2074 true))
2075 phi.emplace_back(std::make_unique<FEFaceEvalType>(
2076 matrix_free, range, true, dof_no, quad_no, first_selected_component));
2077
2078 return phi;
2079 };
2080
2081 data_boundary.op_reinit = [](auto &phi, const unsigned batch) {
2082 static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
2083 };
2084
2085 if (boundary_operation)
2086 data_boundary.op_compute = [&](auto &phi) {
2087 boundary_operation(static_cast<FEFaceEvalType &>(*phi[0]));
2088 };
2089
2090 internal::compute_matrix(
2091 matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix);
2092 }
2093
2094 namespace internal
2095 {
2096 template <int dim,
2097 typename Number,
2098 typename VectorizedArrayType,
2099 typename MatrixType>
2100 void
2101 compute_matrix(
2102 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
2103 const AffineConstraints<Number> &constraints_in,
2104 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
2105 &data_cell,
2106 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2107 &data_face,
2108 const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
2109 &data_boundary,
2110 MatrixType &matrix)
2111 {
2112 std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
2113 constraints_for_matrix;
2114 const AffineConstraints<typename MatrixType::value_type> &constraints =
2115 internal::create_new_affine_constraints_if_needed(
2116 matrix, constraints_in, constraints_for_matrix);
2117
2118 const auto batch_operation =
2119 [&matrix_free, &constraints, &matrix](
2120 auto &data, const std::pair<unsigned int, unsigned int> &range) {
2121 if (!data.op_compute)
2122 return; // nothing to do
2123
2124 auto phi = data.op_create(range);
2125
2126 const unsigned int n_blocks = phi.size();
2127
2128 if (n_blocks == 0)
2129 return;
2130
2131 Table<1, unsigned int> dofs_per_cell(n_blocks);
2132
2133 Table<1, std::vector<types::global_dof_index>> dof_indices(n_blocks);
2134 Table<2, std::vector<types::global_dof_index>> dof_indices_mf(
2135 n_blocks, VectorizedArrayType::size());
2136 Table<1, std::vector<unsigned int>> lexicographic_numbering(n_blocks);
2137 Table<2,
2138 std::array<FullMatrix<typename MatrixType::value_type>,
2139 VectorizedArrayType::size()>>
2140 matrices(n_blocks, n_blocks);
2141
2142 for (unsigned int b = 0; b < n_blocks; ++b)
2143 {
2144 const auto &shape_info = matrix_free.get_shape_info(
2145 data.dof_numbers[b],
2146 data.quad_numbers[b],
2147 data.first_selected_components[b],
2148 phi[b]->get_active_fe_index(),
2149 phi[b]->get_active_quadrature_index());
2150
2151 dofs_per_cell[b] =
2152 shape_info.dofs_per_component_on_cell * shape_info.n_components;
2153
2154 dof_indices[b].resize(dofs_per_cell[b]);
2155
2156 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2157 dof_indices_mf[b][v].resize(dofs_per_cell[b]);
2158
2159 lexicographic_numbering[b] = shape_info.lexicographic_numbering;
2160 }
2161
2162 for (unsigned int bj = 0; bj < n_blocks; ++bj)
2163 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2164 std::fill_n(matrices[bi][bj].begin(),
2165 VectorizedArrayType::size(),
2166 FullMatrix<typename MatrixType::value_type>(
2167 dofs_per_cell[bi], dofs_per_cell[bj]));
2168
2169 for (auto batch = range.first; batch < range.second; ++batch)
2170 {
2171 data.op_reinit(phi, batch);
2172
2173 const unsigned int n_filled_lanes =
2174 data.is_face ?
2175 matrix_free.n_active_entries_per_face_batch(batch) :
2176 matrix_free.n_active_entries_per_cell_batch(batch);
2177
2178 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2179 for (unsigned int b = 0; b < n_blocks; ++b)
2180 {
2181 unsigned int const cell_index =
2182 (data.batch_type[b] == 0) ?
2183 (batch * VectorizedArrayType::size() + v) :
2184 ((data.batch_type[b] == 1) ?
2185 matrix_free.get_face_info(batch).cells_interior[v] :
2186 matrix_free.get_face_info(batch).cells_exterior[v]);
2187
2188 const auto cell_iterator = matrix_free.get_cell_iterator(
2189 cell_index / VectorizedArrayType::size(),
2190 cell_index % VectorizedArrayType::size(),
2191 data.dof_numbers[b]);
2192
2193 if (matrix_free.get_mg_level() !=
2194 numbers::invalid_unsigned_int)
2195 cell_iterator->get_mg_dof_indices(dof_indices[b]);
2196 else
2197 cell_iterator->get_dof_indices(dof_indices[b]);
2198
2199 for (unsigned int j = 0; j < dof_indices[b].size(); ++j)
2200 dof_indices_mf[b][v][j] =
2201 dof_indices[b][lexicographic_numbering[b][j]];
2202 }
2203
2204 for (unsigned int bj = 0; bj < n_blocks; ++bj)
2205 {
2206 for (unsigned int j = 0; j < dofs_per_cell[bj]; ++j)
2207 {
2208 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2209 for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2210 phi[bi]->begin_dof_values()[i] =
2211 (bj == bi) ? static_cast<Number>(i == j) : 0.0;
2212
2213 data.op_compute(phi);
2214
2215 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2216 for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
2217 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2218 matrices[bi][bj][v](i, j) =
2219 phi[bi]->begin_dof_values()[i][v];
2220 }
2221
2222 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2223 for (unsigned int bi = 0; bi < n_blocks; ++bi)
2224 if (bi == bj)
2225 // specialization for blocks on the diagonal
2226 // to writing into diagonal elements of the
2227 // matrix if the corresponding degree of freedom
2228 // is constrained, see also the documentation
2229 // of AffineConstraints::distribute_local_to_global()
2230 constraints.distribute_local_to_global(
2231 matrices[bi][bi][v], dof_indices_mf[bi][v], matrix);
2232 else
2233 constraints.distribute_local_to_global(
2234 matrices[bi][bj][v],
2235 dof_indices_mf[bi][v],
2236 dof_indices_mf[bj][v],
2237 matrix);
2238 }
2239 }
2240 };
2241
2242 const auto cell_operation_wrapped =
2243 [&](const auto &, auto &, const auto &, const auto range) {
2244 batch_operation(data_cell, range);
2245 };
2246
2247 const auto face_operation_wrapped =
2248 [&](const auto &, auto &, const auto &, const auto range) {
2249 batch_operation(data_face, range);
2250 };
2251
2252 const auto boundary_operation_wrapped =
2253 [&](const auto &, auto &, const auto &, const auto range) {
2254 batch_operation(data_boundary, range);
2255 };
2256
2257 if (data_face.op_compute || data_boundary.op_compute)
2258 {
2259 matrix_free.template loop<MatrixType, MatrixType>(
2260 cell_operation_wrapped,
2261 face_operation_wrapped,
2262 boundary_operation_wrapped,
2263 matrix,
2264 matrix);
2265 }
2266 else
2267 matrix_free.template cell_loop<MatrixType, MatrixType>(
2268 cell_operation_wrapped, matrix, matrix);
2269
2270 matrix.compress(VectorOperation::add);
2271 }
2272 } // namespace internal
2273
2274 template <typename CLASS,
2275 int dim,
2276 int fe_degree,
2277 int n_q_points_1d,
2278 int n_components,
2279 typename Number,
2280 typename VectorizedArrayType,
2281 typename MatrixType>
2282 void
2283 compute_matrix(
2284 const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
2285 const AffineConstraints<Number> &constraints,
2286 MatrixType &matrix,
2287 void (CLASS::*cell_operation)(FEEvaluation<dim,
2288 fe_degree,
2289 n_q_points_1d,
2290 n_components,
2291 Number,
2292 VectorizedArrayType> &) const,
2293 void (CLASS::*face_operation)(FEFaceEvaluation<dim,
2294 fe_degree,
2295 n_q_points_1d,
2296 n_components,
2297 Number,
2298 VectorizedArrayType> &,
2299 FEFaceEvaluation<dim,
2300 fe_degree,
2301 n_q_points_1d,
2302 n_components,
2303 Number,
2304 VectorizedArrayType> &)
2305 const,
2306 void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
2307 fe_degree,
2308 n_q_points_1d,
2309 n_components,
2310 Number,
2311 VectorizedArrayType> &)
2312 const,
2313 const CLASS *owning_class,
2314 const unsigned int dof_no,
2315 const unsigned int quad_no,
2316 const unsigned int first_selected_component)
2317 {
2318 compute_matrix<dim,
2319 fe_degree,
2320 n_q_points_1d,
2321 n_components,
2322 Number,
2323 VectorizedArrayType,
2324 MatrixType>(
2325 matrix_free,
2326 constraints,
2327 matrix,
2328 [&](auto &phi) { (owning_class->*cell_operation)(phi); },
2329 [&](auto &phi_m, auto &phi_p) {
2330 (owning_class->*face_operation)(phi_m, phi_p);
2331 },
2332 [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
2333 dof_no,
2334 quad_no,
2335 first_selected_component);
2336 }
2337
2338#endif // DOXYGEN
2339
2340} // namespace MatrixFreeTools
2341
2342DEAL_II_NAMESPACE_CLOSE
2343
2344
2345#endif
AffineConstraints::distribute_local_to_global
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
Definition affine_constraints.h:2651
AffineConstraints::copy_from
void copy_from(const AffineConstraints< other_number > &other)
Definition affine_constraints.h:2778
MatrixFreeTools::ElementActivationAndDeactivationMatrixFree::loop
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:529
MatrixFreeTools::ElementActivationAndDeactivationMatrixFree::cell_loop
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:495
MatrixFreeTools::ElementActivationAndDeactivationMatrixFree::matrix_free
SmartPointer< const MatrixFree< dim, Number, VectorizedArrayType > > matrix_free
Definition tools.h:596
MatrixFreeTools::ElementActivationAndDeactivationMatrixFree::reinit
void reinit(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const AdditionalData &additional_data=AdditionalData())
Definition tools.h:466
MatrixFreeTools::internal::ComputeMatrixScratchData::quad_numbers
std::vector< unsigned int > quad_numbers
Definition tools.h:44
MatrixFreeTools::internal::ComputeMatrixScratchData::batch_type
std::vector< unsigned int > batch_type
Definition tools.h:46
MatrixFreeTools::internal::ComputeMatrixScratchData::dof_numbers
std::vector< unsigned int > dof_numbers
Definition tools.h:43
MatrixFreeTools::internal::ComputeMatrixScratchData::op_reinit
std::function< void(std::vector< std::unique_ptr< FEEvalType > > &, const unsigned int)> op_reinit
Definition tools.h:54
MatrixFreeTools::internal::ComputeMatrixScratchData::op_compute
std::function< void(std::vector< std::unique_ptr< FEEvalType > > &)> op_compute
Definition tools.h:56
MatrixFreeTools::internal::ComputeMatrixScratchData::op_create
std::function< std::vector< std::unique_ptr< FEEvalType > >(const std::pair< unsigned int, unsigned int > &)> op_create
Definition tools.h:51
MatrixFreeTools::internal::ComputeMatrixScratchData::first_selected_components
std::vector< unsigned int > first_selected_components
Definition tools.h:45
MatrixFree::get_cell_active_fe_index
unsigned int get_cell_active_fe_index(const std::pair< unsigned int, unsigned int > range) const
MatrixFree::n_active_entries_per_cell_batch
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
MatrixFree::get_mg_level
unsigned int get_mg_level() const
MatrixFree::get_dof_info
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
MatrixFree::constraint_pool_begin
const Number * constraint_pool_begin(const unsigned int pool_index) const
MatrixFree::get_dof_handler
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
MatrixFree::get_cell_iterator
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
MatrixFree::get_face_range_category
std::pair< unsigned int, unsigned int > get_face_range_category(const std::pair< unsigned int, unsigned int > face_batch_range) const
MatrixFree::constraint_pool_end
const Number * constraint_pool_end(const unsigned int pool_index) const
MatrixFree::n_active_entries_per_face_batch
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
MatrixFree::get_shape_info
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
Threads::ThreadLocalStorage
A class that provides a separate storage location on each thread that accesses the object.
Definition thread_local_storage.h:105
Triangulation::begin
cell_iterator begin(const unsigned int level=0) const
Triangulation::n_active_cells
unsigned int n_active_cells() const
Triangulation::end
cell_iterator end() const
Triangulation::get_boundary_ids
virtual std::vector< types::boundary_id > get_boundary_ids() const
Triangulation::begin_active
active_cell_iterator begin_active(const unsigned int level=0) const
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
fe_evaluation.h
level
unsigned int level
Definition grid_out.cc:4626
cell_index
unsigned int cell_index
Definition grid_tools_topology.cc:783
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1638
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1985
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2053
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
matrix_free.h
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:57
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition operators.h:49
MatrixFreeTools::internal::compute_matrix
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)
MatrixFreeTools::compute_matrix
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)
MatrixFreeTools::compute_diagonal
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)
MatrixFreeTools::categorize_by_boundary_ids
void categorize_by_boundary_ids(const Triangulation< dim > &tria, AdditionalData &additional_data)
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition trilinos_sparse_matrix.cc:50
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
internal::MatrixFreeFunctions::compressed_constraint_kind
std::uint8_t compressed_constraint_kind
Definition dof_info.h:86
internal::MatrixFreeFunctions::unconstrained_compressed_constraint_kind
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
Definition hanging_nodes_internal.h:84
internal::vector_access_add
void vector_access_add(VectorType &vec, const unsigned int entry, const typename VectorType::value_type &val)
Definition vector_access_internal.h:94
internal::check_vector_compatibility
void check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
Definition vector_access_internal.h:179
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition matrix_block.h:617
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
std
STL namespace.
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6943
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Definition vectorization.h:6885
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6927
vector_access_internal.h
Generated by doxygen 1.11.0