Reference documentation for deal.II version 9.5.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// Copyright (C) 2020 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_matrix_free_tools_h
17#define dealii_matrix_free_tools_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/grid/tria.h>
22
26
27
29
35{
41 template <int dim, typename AdditionalData>
42 void
44 AdditionalData & additional_data);
45
54 template <int dim,
55 int fe_degree,
56 int n_q_points_1d,
57 int n_components,
58 typename Number,
59 typename VectorizedArrayType,
60 typename VectorType>
61 void
64 VectorType & diagonal_global,
65 const std::function<void(FEEvaluation<dim,
66 fe_degree,
67 n_q_points_1d,
68 n_components,
69 Number,
70 VectorizedArrayType> &)> &local_vmult,
71 const unsigned int dof_no = 0,
72 const unsigned int quad_no = 0,
73 const unsigned int first_selected_component = 0);
74
78 template <typename CLASS,
79 int dim,
80 int fe_degree,
81 int n_q_points_1d,
82 int n_components,
83 typename Number,
84 typename VectorizedArrayType,
85 typename VectorType>
86 void
89 VectorType & diagonal_global,
90 void (CLASS::*cell_operation)(FEEvaluation<dim,
91 fe_degree,
92 n_q_points_1d,
93 n_components,
94 Number,
95 VectorizedArrayType> &) const,
96 const CLASS * owning_class,
97 const unsigned int dof_no = 0,
98 const unsigned int quad_no = 0,
99 const unsigned int first_selected_component = 0);
100
101
110 template <int dim,
111 int fe_degree,
112 int n_q_points_1d,
113 int n_components,
114 typename Number,
115 typename VectorizedArrayType,
116 typename MatrixType>
117 void
120 const AffineConstraints<Number> & constraints,
121 MatrixType & matrix,
122 const std::function<void(FEEvaluation<dim,
123 fe_degree,
124 n_q_points_1d,
125 n_components,
126 Number,
127 VectorizedArrayType> &)> &local_vmult,
128 const unsigned int dof_no = 0,
129 const unsigned int quad_no = 0,
130 const unsigned int first_selected_component = 0);
131
135 template <typename CLASS,
136 int dim,
137 int fe_degree,
138 int n_q_points_1d,
139 int n_components,
140 typename Number,
141 typename VectorizedArrayType,
142 typename MatrixType>
143 void
146 const AffineConstraints<Number> & constraints,
147 MatrixType & matrix,
148 void (CLASS::*cell_operation)(FEEvaluation<dim,
149 fe_degree,
150 n_q_points_1d,
151 n_components,
152 Number,
153 VectorizedArrayType> &) const,
154 const CLASS * owning_class,
155 const unsigned int dof_no = 0,
156 const unsigned int quad_no = 0,
157 const unsigned int first_selected_component = 0);
158
159
160
169 template <int dim,
170 typename Number,
171 typename VectorizedArrayType = VectorizedArray<Number>>
173 {
174 public:
180 {
184 AdditionalData(const unsigned int dof_index = 0)
186 {}
187
191 unsigned int dof_index;
192 };
193
203 void
205 const AdditionalData &additional_data = AdditionalData())
206 {
207 this->matrix_free = &matrix_free;
208
209 std::vector<unsigned int> valid_fe_indices;
210
211 const auto &fe_collection =
212 matrix_free.get_dof_handler(additional_data.dof_index)
213 .get_fe_collection();
214
215 for (unsigned int i = 0; i < fe_collection.size(); ++i)
216 if (fe_collection[i].n_dofs_per_cell() > 0)
217 valid_fe_indices.push_back(i);
218
219 // TODO: relax this so that arbitrary number of valid
220 // FEs can be accepted
221 AssertDimension(valid_fe_indices.size(), 1);
222
223 fe_index_valid = *valid_fe_indices.begin();
224 }
225
231 template <typename VectorTypeOut, typename VectorTypeIn>
232 void
233 cell_loop(const std::function<void(
235 VectorTypeOut &,
236 const VectorTypeIn &,
237 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
238 VectorTypeOut & dst,
239 const VectorTypeIn & src,
240 const bool zero_dst_vector = false) const
241 {
242 const auto ebd_cell_operation = [&](const auto &matrix_free,
243 auto & dst,
244 const auto &src,
245 const auto &range) {
246 const auto category = matrix_free.get_cell_range_category(range);
247
248 if (category != fe_index_valid)
249 return;
250
251 cell_operation(matrix_free, dst, src, range);
252 };
253
254 matrix_free->template cell_loop<VectorTypeOut, VectorTypeIn>(
255 ebd_cell_operation, dst, src, zero_dst_vector);
256 }
257
265 template <typename VectorTypeOut, typename VectorTypeIn>
266 void
267 loop(const std::function<
269 VectorTypeOut &,
270 const VectorTypeIn &,
271 const std::pair<unsigned int, unsigned int> &)> &cell_operation,
272 const std::function<
274 VectorTypeOut &,
275 const VectorTypeIn &,
276 const std::pair<unsigned int, unsigned int> &)> &face_operation,
277 const std::function<
279 VectorTypeOut &,
280 const VectorTypeIn &,
281 const std::pair<unsigned int, unsigned int> &,
282 const bool)> &boundary_operation,
283 VectorTypeOut & dst,
284 const VectorTypeIn & src,
285 const bool zero_dst_vector = false) const
286 {
287 const auto ebd_cell_operation = [&](const auto &matrix_free,
288 auto & dst,
289 const auto &src,
290 const auto &range) {
291 const auto category = matrix_free.get_cell_range_category(range);
292
293 if (category != fe_index_valid)
294 return;
295
296 cell_operation(matrix_free, dst, src, range);
297 };
298
299 const auto ebd_internal_or_boundary_face_operation =
300 [&](const auto &matrix_free,
301 auto & dst,
302 const auto &src,
303 const auto &range) {
304 const auto category = matrix_free.get_face_range_category(range);
305
306 const unsigned int type =
307 static_cast<unsigned int>(category.first == fe_index_valid) +
308 static_cast<unsigned int>(category.second == fe_index_valid);
309
310 if (type == 0)
311 return; // deactivated face -> nothing to do
312
313 if (type == 1) // boundary face
314 boundary_operation(
315 matrix_free, dst, src, range, category.first == fe_index_valid);
316 else if (type == 2) // internal face
317 face_operation(matrix_free, dst, src, range);
318 };
319
320 matrix_free->template loop<VectorTypeOut, VectorTypeIn>(
321 ebd_cell_operation,
322 ebd_internal_or_boundary_face_operation,
323 ebd_internal_or_boundary_face_operation,
324 dst,
325 src,
326 zero_dst_vector);
327 }
328
329 private:
335
339 unsigned int fe_index_valid;
340 };
341
342 // implementations
343
344#ifndef DOXYGEN
345
346 template <int dim, typename AdditionalData>
347 void
349 AdditionalData & additional_data)
350 {
351 // ... determine if we are on an active or a multigrid level
352 const unsigned int level = additional_data.mg_level;
353 const bool is_mg = (level != numbers::invalid_unsigned_int);
354
355 // ... create empty list for the category of each cell
356 if (is_mg)
357 additional_data.cell_vectorization_category.assign(
358 std::distance(tria.begin(level), tria.end(level)), 0);
359 else
360 additional_data.cell_vectorization_category.assign(tria.n_active_cells(),
361 0);
362
363 // ... set up scaling factor
364 std::vector<unsigned int> factors(GeometryInfo<dim>::faces_per_cell);
365
366 auto bids = tria.get_boundary_ids();
367 std::sort(bids.begin(), bids.end());
368
369 {
370 unsigned int n_bids = bids.size() + 1;
371 int offset = 1;
372 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
373 i++, offset = offset * n_bids)
374 factors[i] = offset;
375 }
376
377 const auto to_category = [&](const auto &cell) {
378 unsigned int c_num = 0;
379 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
380 {
381 const auto face = cell->face(i);
382 if (face->at_boundary() && !cell->has_periodic_neighbor(i))
383 c_num +=
384 factors[i] * (1 + std::distance(bids.begin(),
385 std::find(bids.begin(),
386 bids.end(),
387 face->boundary_id())));
388 }
389 return c_num;
390 };
391
392 if (!is_mg)
393 {
394 for (auto cell = tria.begin_active(); cell != tria.end(); ++cell)
395 {
396 if (cell->is_locally_owned())
397 additional_data
398 .cell_vectorization_category[cell->active_cell_index()] =
399 to_category(cell);
400 }
401 }
402 else
403 {
404 for (auto cell = tria.begin(level); cell != tria.end(level); ++cell)
405 {
406 if (cell->is_locally_owned_on_level())
407 additional_data.cell_vectorization_category[cell->index()] =
408 to_category(cell);
409 }
410 }
411
412 // ... finalize set up of matrix_free
413 additional_data.hold_all_faces_to_owned_cells = true;
414 additional_data.cell_vectorization_categories_strict = true;
415 additional_data.mapping_update_flags_faces_by_cells =
416 additional_data.mapping_update_flags_inner_faces |
417 additional_data.mapping_update_flags_boundary_faces;
418 }
419
420 namespace internal
421 {
422 template <typename Number>
423 struct LocalCSR
424 {
425 std::vector<unsigned int> row_lid_to_gid;
426 std::vector<unsigned int> row;
427 std::vector<unsigned int> col;
428 std::vector<Number> val;
429
430 std::vector<unsigned int> inverse_lookup_rows;
431 std::vector<std::pair<unsigned int, unsigned int>> inverse_lookup_origins;
432 };
433
434 template <int dim,
435 int fe_degree,
436 int n_q_points_1d,
437 int n_components,
438 typename Number,
439 typename VectorizedArrayType>
440 class ComputeDiagonalHelper
441 {
442 public:
443 static const unsigned int n_lanes = VectorizedArrayType::size();
444
445 ComputeDiagonalHelper()
446 : phi(nullptr)
447 , dofs_per_component(0)
448 {}
449
450 ComputeDiagonalHelper(const ComputeDiagonalHelper &)
451 : phi(nullptr)
452 , dofs_per_component(0)
453 {}
454
455 void
456 initialize(FEEvaluation<dim,
457 fe_degree,
458 n_q_points_1d,
459 n_components,
460 Number,
461 VectorizedArrayType> &phi)
462 {
463 // if we are in hp mode and the number of unknowns changed, we must
464 // clear the map of entries
465 if (dofs_per_component != phi.dofs_per_component)
466 {
467 locally_relevant_constraints_hn_map.clear();
468 dofs_per_component = phi.dofs_per_component;
469 }
470 this->phi = &phi;
471 }
472
473 void
474 reinit(const unsigned int cell)
475 {
476 this->phi->reinit(cell);
477
478 // STEP 1: get relevant information from FEEvaluation
479 const auto & dof_info = phi->get_dof_info();
480 const unsigned int n_fe_components = dof_info.start_components.back();
481 const auto & matrix_free = phi->get_matrix_free();
482
483 // if we have a block vector with components with the same DoFHandler,
484 // each component is described with same set of constraints and
485 // we consider the shift in components only during access of the global
486 // vector
487 const unsigned int first_selected_component =
488 n_fe_components == 1 ? 0 : phi->get_first_selected_component();
489
490 const unsigned int n_lanes_filled =
491 matrix_free.n_active_entries_per_cell_batch(cell);
492
493 // STEP 2: setup CSR storage of transposed locally-relevant
494 // constraint matrix
495
496 // (constrained local index, global index of dof
497 // constraints, weight)
498 std::vector<std::tuple<unsigned int, unsigned int, Number>>
499 locally_relevant_constraints, locally_relevant_constraints_tmp;
500 locally_relevant_constraints.reserve(phi->dofs_per_cell);
501 std::vector<unsigned int> constraint_position;
502 std::vector<unsigned char> is_constrained_hn;
503
504 for (unsigned int v = 0; v < n_lanes_filled; ++v)
505 {
506 const unsigned int *dof_indices;
507 unsigned int index_indicators, next_index_indicators;
508
509 const unsigned int start =
510 (cell * n_lanes + v) * n_fe_components + first_selected_component;
511 dof_indices =
512 dof_info.dof_indices.data() + dof_info.row_starts[start].first;
513 index_indicators = dof_info.row_starts[start].second;
514 next_index_indicators = dof_info.row_starts[start + 1].second;
515
516 // STEP 2a: setup locally-relevant constraint matrix in a
517 // coordinate list (COO)
518 locally_relevant_constraints.clear();
519
520 if (n_components == 1 || n_fe_components == 1)
521 {
522 unsigned int ind_local = 0;
523 for (; index_indicators != next_index_indicators;
524 ++index_indicators, ++ind_local)
525 {
526 const std::pair<unsigned short, unsigned short> indicator =
527 dof_info.constraint_indicator[index_indicators];
528
529 for (unsigned int j = 0; j < indicator.first;
530 ++j, ++ind_local)
531 locally_relevant_constraints.emplace_back(ind_local,
532 dof_indices[j],
533 1.0);
534
535 dof_indices += indicator.first;
536
537 const Number *data_val =
538 matrix_free.constraint_pool_begin(indicator.second);
539 const Number *end_pool =
540 matrix_free.constraint_pool_end(indicator.second);
541
542 for (; data_val != end_pool; ++data_val, ++dof_indices)
543 locally_relevant_constraints.emplace_back(ind_local,
544 *dof_indices,
545 *data_val);
546 }
547
548 AssertIndexRange(ind_local, dofs_per_component + 1);
549
550 for (; ind_local < dofs_per_component;
551 ++dof_indices, ++ind_local)
552 locally_relevant_constraints.emplace_back(ind_local,
553 *dof_indices,
554 1.0);
555 }
556 else
557 {
558 // case with vector-valued finite elements where all
559 // components are included in one single vector. Assumption:
560 // first come all entries to the first component, then all
561 // entries to the second one, and so on. This is ensured by
562 // the way MatrixFree reads out the indices.
563 for (unsigned int comp = 0; comp < n_components; ++comp)
564 {
565 unsigned int ind_local = 0;
566
567 // check whether there is any constraint on the current
568 // cell
569 for (; index_indicators != next_index_indicators;
570 ++index_indicators, ++ind_local)
571 {
572 const std::pair<unsigned short, unsigned short>
573 indicator =
574 dof_info.constraint_indicator[index_indicators];
575
576 // run through values up to next constraint
577 for (unsigned int j = 0; j < indicator.first;
578 ++j, ++ind_local)
579 locally_relevant_constraints.emplace_back(
580 comp * dofs_per_component + ind_local,
581 dof_indices[j],
582 1.0);
583 dof_indices += indicator.first;
584
585 const Number *data_val =
586 matrix_free.constraint_pool_begin(indicator.second);
587 const Number *end_pool =
588 matrix_free.constraint_pool_end(indicator.second);
589
590 for (; data_val != end_pool; ++data_val, ++dof_indices)
591 locally_relevant_constraints.emplace_back(
592 comp * dofs_per_component + ind_local,
593 *dof_indices,
594 *data_val);
595 }
596
597 AssertIndexRange(ind_local, dofs_per_component + 1);
598
599 // get the dof values past the last constraint
600 for (; ind_local < dofs_per_component;
601 ++dof_indices, ++ind_local)
602 locally_relevant_constraints.emplace_back(
603 comp * dofs_per_component + ind_local,
604 *dof_indices,
605 1.0);
606
607 if (comp + 1 < n_components)
608 next_index_indicators =
609 dof_info.row_starts[start + comp + 2].second;
610 }
611 }
612
613 // we only need partial sortedness for the algorithm below in that
614 // all entries for a particular row must be adjacent. this is
615 // ensured by the way we fill the field, but check it again
616 for (unsigned int i = 1; i < locally_relevant_constraints.size();
617 ++i)
618 Assert(std::get<0>(locally_relevant_constraints[i]) >=
619 std::get<0>(locally_relevant_constraints[i - 1]),
621
622 // STEP 2c: apply hanging-node constraints
623 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
624 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
625 dof_info.hanging_node_constraint_masks_comp
626 [phi->get_active_fe_index()][first_selected_component])
627 {
628 const auto mask =
629 dof_info.hanging_node_constraint_masks[cell * n_lanes + v];
630
631 // cell has hanging nodes
632 if (mask != ::internal::MatrixFreeFunctions::
634 {
635 // check if hanging node internpolation matrix has been set
636 // up
637 if (locally_relevant_constraints_hn_map.find(mask) ==
638 locally_relevant_constraints_hn_map.end())
639 fill_constraint_type_into_map(mask);
640
641 const auto &locally_relevant_constraints_hn =
642 locally_relevant_constraints_hn_map[mask];
643
644 locally_relevant_constraints_tmp.clear();
645 if (locally_relevant_constraints_tmp.capacity() <
646 locally_relevant_constraints.size())
647 locally_relevant_constraints_tmp.reserve(
648 locally_relevant_constraints.size() +
649 locally_relevant_constraints_hn.size());
650
651 // combine with other constraints: to avoid binary
652 // searches, we first build a list of where constraints
653 // are pointing to, and then merge the two lists
654 constraint_position.assign(phi->dofs_per_cell,
656 for (auto &a : locally_relevant_constraints)
657 if (constraint_position[std::get<0>(a)] ==
658 numbers::invalid_unsigned_int)
659 constraint_position[std::get<0>(a)] =
660 std::distance(locally_relevant_constraints.data(),
661 &a);
662 is_constrained_hn.assign(phi->dofs_per_cell, false);
663 for (auto &hn : locally_relevant_constraints_hn)
664 is_constrained_hn[std::get<0>(hn)] = 1;
665
666 // not constrained from hanging nodes
667 for (const auto &a : locally_relevant_constraints)
668 if (is_constrained_hn[std::get<0>(a)] == 0)
669 locally_relevant_constraints_tmp.push_back(a);
670
671 // dof is constrained by hanging nodes: build transitive
672 // closure
673 for (const auto &hn : locally_relevant_constraints_hn)
674 if (constraint_position[std::get<1>(hn)] !=
675 numbers::invalid_unsigned_int)
676 {
677 AssertIndexRange(constraint_position[std::get<1>(hn)],
678 locally_relevant_constraints.size());
679 auto other = locally_relevant_constraints.begin() +
680 constraint_position[std::get<1>(hn)];
681 AssertDimension(std::get<0>(*other), std::get<1>(hn));
682
683 for (; other != locally_relevant_constraints.end() &&
684 std::get<0>(*other) == std::get<1>(hn);
685 ++other)
686 locally_relevant_constraints_tmp.emplace_back(
687 std::get<0>(hn),
688 std::get<1>(*other),
689 std::get<2>(hn) * std::get<2>(*other));
690 }
691
692 std::swap(locally_relevant_constraints,
693 locally_relevant_constraints_tmp);
694 }
695 }
696
697 // STEP 2d: transpose COO
698 std::sort(locally_relevant_constraints.begin(),
699 locally_relevant_constraints.end(),
700 [](const auto &a, const auto &b) {
701 if (std::get<1>(a) < std::get<1>(b))
702 return true;
703 return (std::get<1>(a) == std::get<1>(b)) &&
704 (std::get<0>(a) < std::get<0>(b));
705 });
706
707 // STEP 2e: translate COO to CRS
708 auto &c_pool = c_pools[v];
709 {
710 c_pool.row_lid_to_gid.clear();
711 c_pool.row.clear();
712 c_pool.row.push_back(0);
713 c_pool.col.clear();
714 c_pool.val.clear();
715
716 if (locally_relevant_constraints.size() > 0)
717 c_pool.row_lid_to_gid.emplace_back(
718 std::get<1>(locally_relevant_constraints.front()));
719 for (const auto &j : locally_relevant_constraints)
720 {
721 if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
722 {
723 c_pool.row_lid_to_gid.push_back(std::get<1>(j));
724 c_pool.row.push_back(c_pool.val.size());
725 }
726
727 c_pool.col.emplace_back(std::get<0>(j));
728 c_pool.val.emplace_back(std::get<2>(j));
729 }
730
731 if (c_pool.val.size() > 0)
732 c_pool.row.push_back(c_pool.val.size());
733
734 c_pool.inverse_lookup_rows.clear();
735 c_pool.inverse_lookup_rows.resize(1 + phi->dofs_per_cell);
736 for (const unsigned int i : c_pool.col)
737 c_pool.inverse_lookup_rows[1 + i]++;
738 // transform to offsets
739 std::partial_sum(c_pool.inverse_lookup_rows.begin(),
740 c_pool.inverse_lookup_rows.end(),
741 c_pool.inverse_lookup_rows.begin());
742 AssertDimension(c_pool.inverse_lookup_rows.back(),
743 c_pool.col.size());
744
745 c_pool.inverse_lookup_origins.resize(c_pool.col.size());
746 std::vector<unsigned int> inverse_lookup_count(
747 phi->dofs_per_cell);
748 for (unsigned int row = 0; row < c_pool.row.size() - 1; ++row)
749 for (unsigned int col = c_pool.row[row];
750 col < c_pool.row[row + 1];
751 ++col)
752 {
753 const unsigned int index = c_pool.col[col];
754 c_pool.inverse_lookup_origins
755 [c_pool.inverse_lookup_rows[index] +
756 inverse_lookup_count[index]] = std::make_pair(row, col);
757 ++inverse_lookup_count[index];
758 }
759 }
760 }
761
762 // STEP 3: compute element matrix A_e, apply
763 // locally-relevant constraints C_e^T * A_e * C_e, and get the
764 // the diagonal entry
765 // (C_e^T * A_e * C_e)(i,i)
766 // or
767 // C_e^T(i,:) * A_e * C_e(:,i).
768 //
769 // Since, we compute the element matrix column-by-column and as a
770 // result never actually have the full element matrix, we actually
771 // perform following steps:
772 // 1) loop over all columns of the element matrix
773 // a) compute column i
774 // b) compute for each j (rows of C_e^T):
775 // (C_e^T(j,:) * A_e(:,i)) * C_e(i,j)
776 // or
777 // (C_e^T(j,:) * A_e(:,i)) * C_e^T(j,i)
778 // This gives a contribution the j-th entry of the
779 // locally-relevant diagonal and comprises the multiplication
780 // by the locally-relevant constraint matrix from the left and
781 // the right. There is no contribution to the j-th vector
782 // entry if the j-th row of C_e^T is empty or C_e^T(j,i) is
783 // zero.
784
785 // set size locally-relevant diagonal
786 for (unsigned int v = 0; v < n_lanes_filled; ++v)
787 diagonals_local_constrained[v].assign(
788 c_pools[v].row_lid_to_gid.size() *
789 (n_fe_components == 1 ? n_components : 1),
790 Number(0.0));
791 }
792
793 void
794 fill_constraint_type_into_map(
795 const ::internal::MatrixFreeFunctions::compressed_constraint_kind
796 mask)
797 {
798 auto &constraints_hn = locally_relevant_constraints_hn_map[mask];
799
800 // assume that we constrain one face, i.e., (fe_degree + 1)^(dim-1)
801 // unknowns - we might have more or less entries, but this is a good
802 // first guess
803 const unsigned int degree =
804 phi->get_shape_info().data.front().fe_degree;
805 constraints_hn.reserve(Utilities::pow(degree + 1, dim - 1));
806
807 // 1) collect hanging-node constraints for cell assuming
808 // scalar finite element
809 values_dofs.resize(dofs_per_component);
810 std::array<
812 VectorizedArrayType::size()>
813 constraint_mask;
814 constraint_mask.fill(::internal::MatrixFreeFunctions::
816 constraint_mask[0] = mask;
817
818 for (unsigned int i = 0; i < dofs_per_component; ++i)
819 {
820 for (unsigned int j = 0; j < dofs_per_component; ++j)
821 values_dofs[j] = VectorizedArrayType();
822 values_dofs[i] = Number(1);
823
825 dim,
826 Number,
827 VectorizedArrayType>::apply(1,
828 degree,
829 phi->get_shape_info(),
830 false,
831 constraint_mask,
832 values_dofs.data());
833
834 const Number tolerance =
835 std::max(Number(1e-12),
836 std::numeric_limits<Number>::epsilon() * 16);
837 for (unsigned int j = 0; j < dofs_per_component; ++j)
838 if (std::abs(values_dofs[j][0]) > tolerance &&
839 (j != i ||
840 std::abs(values_dofs[j][0] - Number(1)) > tolerance))
841 constraints_hn.emplace_back(j, i, values_dofs[j][0]);
842 }
843
844 // 2) extend for multiple components
845 const unsigned int n_hn_constraints = constraints_hn.size();
846 constraints_hn.resize(n_hn_constraints * n_components);
847
848 for (unsigned int c = 1; c < n_components; ++c)
849 for (unsigned int i = 0; i < n_hn_constraints; ++i)
850 constraints_hn[c * n_hn_constraints + i] = std::make_tuple(
851 std::get<0>(constraints_hn[i]) + c * dofs_per_component,
852 std::get<1>(constraints_hn[i]) + c * dofs_per_component,
853 std::get<2>(constraints_hn[i]));
854 }
855
856 void
857 prepare_basis_vector(const unsigned int i)
858 {
859 this->i = i;
860
861 // compute i-th column of element stiffness matrix:
862 // this could be simply performed as done at the moment with
863 // matrix-free operator evaluation applied to a ith-basis vector
864 for (unsigned int j = 0; j < phi->dofs_per_cell; ++j)
865 phi->begin_dof_values()[j] = VectorizedArrayType();
866 phi->begin_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<
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>::value> * = 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<
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>::value> * = 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:267
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:233
SmartPointer< const MatrixFree< dim, Number, VectorizedArrayType > > matrix_free
Definition tools.h:334
void reinit(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const AdditionalData &additional_data=AdditionalData())
Definition tools.h:204
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
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
unsigned int level
Definition grid_out.cc:4618
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
@ 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:447
std::uint8_t compressed_constraint_kind
Definition dof_info.h:83
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:213
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::Triangulation< dim, spacedim > & tria