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