deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
tensor_product_matrix.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 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_tensor_product_matrix_h
16#define dealii_tensor_product_matrix_h
17
18
19#include <deal.II/base/config.h>
20
23#include <deal.II/base/mutex.h>
25
27
29
30#include <bitset>
31
32
34
35// Forward declarations
36#ifndef DOXYGEN
37template <typename>
38class Vector;
39template <typename>
40class FullMatrix;
41#endif
42
114template <int dim, typename Number, int n_rows_1d = -1>
116{
117public:
122 using value_type = Number;
123
128 static constexpr int n_rows_1d_static = n_rows_1d;
129
134
139 template <typename T>
141 const T &derivative_matrix);
142
160 template <typename T>
161 void
163
169 unsigned int
170 m() const;
171
177 unsigned int
178 n() const;
179
193 void
194 vmult(const ArrayView<Number> &dst, const ArrayView<const Number> &src) const;
195
201 void
203 const ArrayView<const Number> &src,
204 AlignedVector<Number> &tmp) const;
205
212 void
214 const ArrayView<const Number> &src) const;
215
219 std::size_t
221
222protected:
226 std::array<Table<2, Number>, dim> mass_matrix;
227
231 std::array<Table<2, Number>, dim> derivative_matrix;
232
237 std::array<AlignedVector<Number>, dim> eigenvalues;
238
243 std::array<Table<2, Number>, dim> eigenvectors;
244
245private:
250
255};
256
257
258
259namespace internal
260{
262 {
263 template <typename Number>
265 {
269 static constexpr std::size_t width = VectorizedArrayTrait::width();
270
272 std::pair<std::bitset<width>,
273 std::pair<Table<2, Number>, Table<2, Number>>>;
274
276 : eps(std::sqrt(std::numeric_limits<ScalarNumber>::epsilon()))
277 {}
278
279 bool
280 operator()(const MatrixPairType &left, const MatrixPairType &right) const
281 {
282 const auto &M_0 = left.second.first;
283 const auto &K_0 = left.second.second;
284 const auto &M_1 = right.second.first;
285 const auto &K_1 = right.second.second;
286
287 std::bitset<width> mask;
288
289 for (unsigned int v = 0; v < width; ++v)
290 mask[v] = left.first[v] && right.first[v];
291
293 eps, false /*use relative tolerance*/, mask);
294
295 if (comparator(M_0, M_1))
296 return true;
297 else if (comparator(M_1, M_0))
298 return false;
299 else if (comparator(K_0, K_1))
300 return true;
301 else
302 return false;
303 }
304
305 private:
307 };
308 } // namespace TensorProductMatrixSymmetricSum
309} // namespace internal
310
311
312
336template <int dim, typename Number, int n_rows_1d = -1>
338{
339 using MatrixPairType = std::pair<Table<2, Number>, Table<2, Number>>;
340
341 using MatrixPairTypeWithMask = std::pair<
342 std::bitset<::internal::VectorizedArrayTrait<Number>::width()>,
344
345public:
350 {
355 const bool precompute_inverse_diagonal = true);
356
361
366 };
367
372 const AdditionalData &additional_data = AdditionalData());
373
378 void
379 reserve(const unsigned int size);
380
386 template <typename T>
387 void
388 insert(const unsigned int index, const T &Ms, const T &Ks);
389
394 void
396
400 void
401 apply_inverse(const unsigned int index,
403 const ArrayView<const Number> &src_in) const;
404
408 std::size_t
410
419 std::size_t
421
422private:
427
432
437 std::vector<MatrixPairType> mass_and_derivative_matrices;
438
443 std::map<
445 unsigned int,
448
456 std::vector<unsigned int> indices;
457
462
467
472
477
482
486 std::vector<unsigned int> vector_ptr;
487
491 std::vector<unsigned int> matrix_ptr;
492
496 std::vector<unsigned int> vector_n_rows_1d;
497};
498
499
500/*----------------------- Inline functions ----------------------------------*/
501
502#ifndef DOXYGEN
503
504namespace internal
505{
507 {
516 template <typename Number>
517 void
518 spectral_assembly(const Number *mass_matrix,
519 const Number *derivative_matrix,
520 const unsigned int n_rows,
521 const unsigned int n_cols,
522 Number *eigenvalues,
523 Number *eigenvectors)
524 {
525 Assert(n_rows == n_cols, ExcNotImplemented());
526
527 std::vector<bool> constrained_dofs(n_rows, false);
528
529 for (unsigned int i = 0; i < n_rows; ++i)
530 {
531 if (mass_matrix[i + i * n_rows] == 0.0)
532 {
533 Assert(derivative_matrix[i + i * n_rows] == 0.0,
535
536 for (unsigned int j = 0; j < n_rows; ++j)
537 {
538 Assert(derivative_matrix[i + j * n_rows] == 0,
540 Assert(derivative_matrix[j + i * n_rows] == 0,
542 }
543
544 constrained_dofs[i] = true;
545 }
546 }
547
548 const auto transpose_fill_nm = [&constrained_dofs](Number *out,
549 const Number *in,
550 const unsigned int n,
551 const unsigned int m) {
552 for (unsigned int mm = 0, c = 0; mm < m; ++mm)
553 for (unsigned int nn = 0; nn < n; ++nn, ++c)
554 out[mm + nn * m] =
555 (mm == nn && constrained_dofs[mm]) ? Number(1.0) : in[c];
556 };
557
558 std::vector<::Vector<Number>> eigenvecs(n_rows);
559 LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
560 LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);
561
562 transpose_fill_nm(&(mass_copy(0, 0)), mass_matrix, n_rows, n_cols);
563 transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
564
565 deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy,
566 eigenvecs);
567 AssertDimension(eigenvecs.size(), n_rows);
568 for (unsigned int i = 0, c = 0; i < n_rows; ++i)
569 for (unsigned int j = 0; j < n_cols; ++j, ++c)
570 if (constrained_dofs[i] == false)
571 eigenvectors[c] = eigenvecs[j][i];
572
573 for (unsigned int i = 0; i < n_rows; ++i, ++eigenvalues)
574 *eigenvalues = deriv_copy.eigenvalue(i).real();
575 }
576
577
578
579 template <std::size_t dim, typename Number>
580 inline void
581 setup(const std::array<Table<2, Number>, dim> &mass_matrix,
582 const std::array<Table<2, Number>, dim> &derivative_matrix,
583 std::array<Table<2, Number>, dim> &eigenvectors,
584 std::array<AlignedVector<Number>, dim> &eigenvalues)
585 {
586 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
588
589 for (unsigned int dir = 0; dir < dim; ++dir)
590 {
591 AssertDimension(n_rows_1d, mass_matrix[dir].n_cols());
592 AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
593 AssertDimension(mass_matrix[dir].n_rows(),
594 derivative_matrix[dir].n_rows());
595 AssertDimension(mass_matrix[dir].n_rows(),
596 derivative_matrix[dir].n_cols());
597
598 eigenvectors[dir].reinit(mass_matrix[dir].n_cols(),
599 mass_matrix[dir].n_rows());
600 eigenvalues[dir].resize(mass_matrix[dir].n_cols());
601 internal::TensorProductMatrixSymmetricSum::spectral_assembly<Number>(
602 &(mass_matrix[dir](0, 0)),
603 &(derivative_matrix[dir](0, 0)),
604 mass_matrix[dir].n_rows(),
605 mass_matrix[dir].n_cols(),
606 eigenvalues[dir].begin(),
607 &(eigenvectors[dir](0, 0)));
608 }
609 }
610
611
612
613 template <std::size_t dim, typename Number, std::size_t n_lanes>
614 inline void
615 setup(
616 const std::array<Table<2, VectorizedArray<Number, n_lanes>>, dim>
617 &mass_matrix,
618 const std::array<Table<2, VectorizedArray<Number, n_lanes>>, dim>
619 &derivative_matrix,
623 {
624 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
625 constexpr unsigned int macro_size =
627 const std::size_t nm_flat_size_max = n_rows_1d * n_rows_1d * macro_size;
628 const std::size_t n_flat_size_max = n_rows_1d * macro_size;
629
630 std::vector<Number> mass_matrix_flat;
631 std::vector<Number> deriv_matrix_flat;
632 std::vector<Number> eigenvalues_flat;
633 std::vector<Number> eigenvectors_flat;
638 std::array<unsigned int, macro_size> offsets_nm;
639 std::array<unsigned int, macro_size> offsets_n;
640 for (unsigned int dir = 0; dir < dim; ++dir)
641 {
642 AssertDimension(n_rows_1d, mass_matrix[dir].n_cols());
643 AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
644 AssertDimension(mass_matrix[dir].n_rows(),
645 derivative_matrix[dir].n_rows());
646 AssertDimension(mass_matrix[dir].n_rows(),
647 derivative_matrix[dir].n_cols());
648
649 const unsigned int n_rows = mass_matrix[dir].n_rows();
650 const unsigned int n_cols = mass_matrix[dir].n_cols();
651 const unsigned int nm = n_rows * n_cols;
652 for (unsigned int vv = 0; vv < macro_size; ++vv)
653 offsets_nm[vv] = nm * vv;
654
656 false,
657 nm,
658 &(mass_matrix[dir](0, 0)),
659 offsets_nm.data(),
660 mass_matrix_flat.data());
662 false,
663 nm,
664 &(derivative_matrix[dir](0, 0)),
665 offsets_nm.data(),
666 deriv_matrix_flat.data());
667
668 const Number *mass_cbegin = mass_matrix_flat.data();
669 const Number *deriv_cbegin = deriv_matrix_flat.data();
670 Number *eigenvec_begin = eigenvectors_flat.data();
671 Number *eigenval_begin = eigenvalues_flat.data();
672 for (unsigned int lane = 0; lane < macro_size; ++lane)
673 internal::TensorProductMatrixSymmetricSum::spectral_assembly<
674 Number>(mass_cbegin + nm * lane,
675 deriv_cbegin + nm * lane,
676 n_rows,
677 n_cols,
678 eigenval_begin + n_rows * lane,
679 eigenvec_begin + nm * lane);
680
681 eigenvalues[dir].resize(n_rows);
682 eigenvectors[dir].reinit(n_rows, n_cols);
683 for (unsigned int vv = 0; vv < macro_size; ++vv)
684 offsets_n[vv] = n_rows * vv;
686 n_rows,
687 eigenvalues_flat.data(),
688 offsets_n.data(),
689 eigenvalues[dir].begin());
691 nm,
692 eigenvectors_flat.data(),
693 offsets_nm.data(),
694 &(eigenvectors[dir](0, 0)));
695 }
696 }
697
698
699
700 template <std::size_t dim, typename Number>
701 inline std::array<Table<2, Number>, dim>
702 convert(const std::array<Table<2, Number>, dim> &mass_matrix)
703 {
704 return mass_matrix;
705 }
706
707
708
709 template <std::size_t dim, typename Number>
710 inline std::array<Table<2, Number>, dim>
711 convert(const std::array<FullMatrix<Number>, dim> &mass_matrix)
712 {
713 std::array<Table<2, Number>, dim> mass_copy;
714
715 std::transform(mass_matrix.cbegin(),
716 mass_matrix.cend(),
717 mass_copy.begin(),
718 [](const FullMatrix<Number> &m) -> Table<2, Number> {
719 return m;
720 });
721
722 return mass_copy;
723 }
724
725
726
727 template <std::size_t dim, typename Number>
728 inline std::array<Table<2, Number>, dim>
729 convert(const Table<2, Number> &matrix)
730 {
731 std::array<Table<2, Number>, dim> matrices;
732
733 std::fill(matrices.begin(), matrices.end(), matrix);
734
735 return matrices;
736 }
737
738
739
740 template <int n_rows_1d_templated, std::size_t dim, typename Number>
741 void
742 vmult(Number *dst,
743 const Number *src,
745 const unsigned int n_rows_1d_non_templated,
746 const std::array<const Number *, dim> &mass_matrix,
747 const std::array<const Number *, dim> &derivative_matrix)
748 {
749 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
752 const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
753
754 tmp.resize_fast(n * 2);
755 Number *t = tmp.begin();
756
758 dim,
761 Number>
762 eval({}, {}, {}, n_rows_1d, n_rows_1d);
763
764 if (dim == 1)
765 {
766 const Number *A = derivative_matrix[0];
767 eval.template apply<0, false, false>(A, src, dst);
768 }
769
770 else if (dim == 2)
771 {
772 const Number *A0 = derivative_matrix[0];
773 const Number *M0 = mass_matrix[0];
774 const Number *A1 = derivative_matrix[1];
775 const Number *M1 = mass_matrix[1];
776 eval.template apply<0, false, false>(M0, src, t);
777 eval.template apply<1, false, false>(A1, t, dst);
778 eval.template apply<0, false, false>(A0, src, t);
779 eval.template apply<1, false, true>(M1, t, dst);
780 }
781
782 else if (dim == 3)
783 {
784 const Number *A0 = derivative_matrix[0];
785 const Number *M0 = mass_matrix[0];
786 const Number *A1 = derivative_matrix[1];
787 const Number *M1 = mass_matrix[1];
788 const Number *A2 = derivative_matrix[2];
789 const Number *M2 = mass_matrix[2];
790 eval.template apply<0, false, false>(M0, src, t + n);
791 eval.template apply<1, false, false>(M1, t + n, t);
792 eval.template apply<2, false, false>(A2, t, dst);
793 eval.template apply<1, false, false>(A1, t + n, t);
794 eval.template apply<0, false, false>(A0, src, t + n);
795 eval.template apply<1, false, true>(M1, t + n, t);
796 eval.template apply<2, false, true>(M2, t, dst);
797 }
798
799 else
801 }
802
803
804
805 template <int n_rows_1d_templated, std::size_t dim, typename Number>
806 void
807 apply_inverse(Number *dst,
808 const Number *src,
809 const unsigned int n_rows_1d_non_templated,
810 const std::array<const Number *, dim> &eigenvectors,
811 const std::array<const Number *, dim> &eigenvalues,
812 const Number *inverted_eigenvalues = nullptr)
813 {
814 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
817
819 dim,
822 Number>
823 eval({}, {}, {}, n_rows_1d, n_rows_1d);
824
825 // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
826 // --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
827 // --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
828 // while the eigenvectors are stored column-wise in S, i.e.
829 // rows correspond to dofs whereas columns to eigenvalue indices!
830 if (dim == 1)
831 {
832 const Number *S = eigenvectors[0];
833 eval.template apply<0, true, false>(S, src, dst);
834
835 for (unsigned int i = 0; i < n_rows_1d; ++i)
836 if (inverted_eigenvalues)
837 dst[i] *= inverted_eigenvalues[i];
838 else
839 dst[i] /= eigenvalues[0][i];
840
841 eval.template apply<0, false, false>(S, dst, dst);
842 }
843
844 else if (dim == 2)
845 {
846 const Number *S0 = eigenvectors[0];
847 const Number *S1 = eigenvectors[1];
848 eval.template apply<0, true, false>(S0, src, dst);
849 eval.template apply<1, true, false>(S1, dst, dst);
850
851 for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
852 for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
853 if (inverted_eigenvalues)
854 dst[c] *= inverted_eigenvalues[c];
855 else
856 dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
857
858 eval.template apply<1, false, false>(S1, dst, dst);
859 eval.template apply<0, false, false>(S0, dst, dst);
860 }
861
862 else if (dim == 3)
863 {
864 const Number *S0 = eigenvectors[0];
865 const Number *S1 = eigenvectors[1];
866 const Number *S2 = eigenvectors[2];
867 eval.template apply<0, true, false>(S0, src, dst);
868 eval.template apply<1, true, false>(S1, dst, dst);
869 eval.template apply<2, true, false>(S2, dst, dst);
870
871 for (unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
872 for (unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
873 for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
874 if (inverted_eigenvalues)
875 dst[c] *= inverted_eigenvalues[c];
876 else
877 dst[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] +
878 eigenvalues[0][i0]);
879
880 eval.template apply<2, false, false>(S2, dst, dst);
881 eval.template apply<1, false, false>(S1, dst, dst);
882 eval.template apply<0, false, false>(S0, dst, dst);
883 }
884
885 else
887 }
888
889
890
891 template <int n_rows_1d_templated, std::size_t dim, typename Number>
892 void
893 select_vmult(Number *dst,
894 const Number *src,
896 const unsigned int n_rows_1d,
897 const std::array<const Number *, dim> &mass_matrix,
898 const std::array<const Number *, dim> &derivative_matrix);
899
900
901
902 template <int n_rows_1d_templated, std::size_t dim, typename Number>
903 void
904 select_apply_inverse(Number *dst,
905 const Number *src,
906 const unsigned int n_rows_1d,
907 const std::array<const Number *, dim> &eigenvectors,
908 const std::array<const Number *, dim> &eigenvalues,
909 const Number *inverted_eigenvalues = nullptr);
910 } // namespace TensorProductMatrixSymmetricSum
911} // namespace internal
912
913
914template <int dim, typename Number, int n_rows_1d>
915inline unsigned int
917{
918 unsigned int m = mass_matrix[0].n_rows();
919 for (unsigned int d = 1; d < dim; ++d)
920 m *= mass_matrix[d].n_rows();
921 return m;
922}
923
924
925
926template <int dim, typename Number, int n_rows_1d>
927inline unsigned int
929{
930 unsigned int n = mass_matrix[0].n_cols();
931 for (unsigned int d = 1; d < dim; ++d)
932 n *= mass_matrix[d].n_cols();
933 return n;
934}
935
936
937
938template <int dim, typename Number, int n_rows_1d>
939inline void
943{
944 std::lock_guard<std::mutex> lock(this->mutex);
945 this->vmult(dst_view, src_view, this->tmp_array);
946}
947
948
949
950template <int dim, typename Number, int n_rows_1d>
951inline void
955 AlignedVector<Number> &tmp_array) const
956{
957 AssertDimension(dst_view.size(), this->m());
958 AssertDimension(src_view.size(), this->n());
959
960 Number *dst = dst_view.begin();
961 const Number *src = src_view.begin();
962
963 std::array<const Number *, dim> mass_matrix, derivative_matrix;
964
965 for (unsigned int d = 0; d < dim; ++d)
966 {
967 mass_matrix[d] = &this->mass_matrix[d](0, 0);
968 derivative_matrix[d] = &this->derivative_matrix[d](0, 0);
969 }
970
971 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
972
973 if (n_rows_1d != -1)
974 internal::TensorProductMatrixSymmetricSum::vmult<
975 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
976 src,
977 tmp_array,
980 derivative_matrix);
981 else
982 internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
983 dst,
984 src,
985 tmp_array,
987 mass_matrix,
988 derivative_matrix);
989}
990
991
992
993template <int dim, typename Number, int n_rows_1d>
994inline void
998{
999 AssertDimension(dst_view.size(), this->n());
1000 AssertDimension(src_view.size(), this->m());
1001
1002 Number *dst = dst_view.begin();
1003 const Number *src = src_view.begin();
1004
1005 std::array<const Number *, dim> eigenvectors, eigenvalues;
1006
1007 for (unsigned int d = 0; d < dim; ++d)
1008 {
1009 eigenvectors[d] = &this->eigenvectors[d](0, 0);
1010 eigenvalues[d] = this->eigenvalues[d].data();
1011 }
1012
1013 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
1014
1015 if (n_rows_1d != -1)
1016 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1017 n_rows_1d == -1 ? 0 : n_rows_1d>(
1019 else
1020 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1022}
1023
1024
1025
1026template <int dim, typename Number, int n_rows_1d>
1027std::size_t
1029 const
1030{
1031 return MemoryConsumption::memory_consumption(mass_matrix) +
1032 MemoryConsumption::memory_consumption(derivative_matrix) +
1036}
1037
1038
1039
1040template <int dim, typename Number, int n_rows_1d>
1041template <typename T>
1043 TensorProductMatrixSymmetricSum(const T &mass_matrix,
1044 const T &derivative_matrix)
1045{
1046 reinit(mass_matrix, derivative_matrix);
1047}
1048
1049
1050
1051template <int dim, typename Number, int n_rows_1d>
1052template <typename T>
1053inline void
1055 const T &mass_matrix,
1056 const T &derivative_matrix)
1057{
1058 this->mass_matrix =
1059 internal::TensorProductMatrixSymmetricSum::convert<dim>(mass_matrix);
1060 this->derivative_matrix =
1061 internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1062
1063 internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1064 this->derivative_matrix,
1065 this->eigenvectors,
1066 this->eigenvalues);
1067}
1068
1069
1070
1071template <int dim, typename Number, int n_rows_1d>
1073 AdditionalData::AdditionalData(const bool compress_matrices,
1074 const bool precompute_inverse_diagonal)
1075 : compress_matrices(compress_matrices)
1076 , precompute_inverse_diagonal(precompute_inverse_diagonal)
1077{}
1078
1079
1080
1081template <int dim, typename Number, int n_rows_1d>
1084 const AdditionalData &additional_data)
1085 : compress_matrices(additional_data.compress_matrices)
1086 , precompute_inverse_diagonal(additional_data.precompute_inverse_diagonal)
1087{}
1088
1089
1090
1091template <int dim, typename Number, int n_rows_1d>
1092void
1094 const unsigned int size)
1095{
1096 if (compress_matrices == false)
1097 mass_and_derivative_matrices.resize(size * dim);
1098 else
1099 indices.assign(size * dim, numbers::invalid_unsigned_int);
1100}
1101
1102
1103
1104template <int dim, typename Number, int n_rows_1d>
1105template <typename T>
1106void
1108 const unsigned int index,
1109 const T &Ms_in,
1110 const T &Ks_in)
1111{
1112 const auto Ms =
1113 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ms_in);
1114 const auto Ks =
1115 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ks_in);
1116
1117 for (unsigned int d = 0; d < dim; ++d)
1118 {
1119 if (compress_matrices == false)
1120 {
1121 const MatrixPairType matrix(Ms[d], Ks[d]);
1122 mass_and_derivative_matrices[index * dim + d] = matrix;
1123 }
1124 else
1125 {
1126 using VectorizedArrayTrait =
1128
1129 std::bitset<VectorizedArrayTrait::width()> mask;
1130
1131 for (unsigned int v = 0; v < VectorizedArrayTrait::width(); ++v)
1132 {
1133 typename VectorizedArrayTrait::value_type a = 0.0;
1134
1135 for (unsigned int i = 0; i < Ms[d].size(0); ++i)
1136 for (unsigned int j = 0; j < Ms[d].size(1); ++j)
1137 {
1138 a += std::abs(VectorizedArrayTrait::get(Ms[d][i][j], v));
1139 a += std::abs(VectorizedArrayTrait::get(Ks[d][i][j], v));
1140 }
1141
1142 mask[v] = (a != 0.0);
1143 }
1144
1145 const MatrixPairTypeWithMask matrix{mask, {Ms[d], Ks[d]}};
1146
1147 const auto ptr = cache.find(matrix);
1148
1149 if (ptr != cache.end())
1150 {
1151 const auto ptr_index = ptr->second;
1152 indices[index * dim + d] = ptr_index;
1153
1154 if ([&]() {
1155 for (unsigned int v = 0; v < VectorizedArrayTrait::width();
1156 ++v)
1157 if ((mask[v] == true) && (ptr->first.first[v] == false))
1158 return false;
1159
1160 return true;
1161 }())
1162 {
1163 // nothing to do
1164 }
1165 else
1166 {
1167 auto mask_new = ptr->first.first;
1168 auto Ms_new = ptr->first.second.first;
1169 auto Ks_new = ptr->first.second.second;
1170
1171 for (unsigned int v = 0; v < VectorizedArrayTrait::width();
1172 ++v)
1173 if (mask_new[v] == false && mask[v] == true)
1174 {
1175 mask_new[v] = true;
1176
1177 for (unsigned int i = 0; i < Ms_new.size(0); ++i)
1178 for (unsigned int j = 0; j < Ms_new.size(1); ++j)
1179 {
1180 VectorizedArrayTrait::get(Ms_new[i][j], v) =
1181 VectorizedArrayTrait::get(Ms[d][i][j], v);
1182 VectorizedArrayTrait::get(Ks_new[i][j], v) =
1183 VectorizedArrayTrait::get(Ks[d][i][j], v);
1184 }
1185 }
1186
1187 cache.erase(ptr);
1188
1189 const MatrixPairTypeWithMask entry_new{mask_new,
1190 {Ms_new, Ks_new}};
1191
1192 const auto ptr_ = cache.find(entry_new);
1193 AssertThrow(ptr_ == cache.end(), ExcNotImplemented());
1194
1195 cache[entry_new] = ptr_index;
1196 }
1197 }
1198 else
1199 {
1200 const auto size = cache.size();
1201 indices[index * dim + d] = size;
1202 cache[matrix] = size;
1203 }
1204 }
1205 }
1206}
1207
1208
1209
1210template <int dim, typename Number, int n_rows_1d>
1211void
1213{
1214 const auto store = [&](const unsigned int index,
1215 const MatrixPairType &M_and_K) {
1216 std::array<Table<2, Number>, 1> mass_matrix;
1217 mass_matrix[0] = M_and_K.first;
1218
1219 std::array<Table<2, Number>, 1> derivative_matrix;
1220 derivative_matrix[0] = M_and_K.second;
1221
1222 std::array<Table<2, Number>, 1> eigenvectors;
1223 std::array<AlignedVector<Number>, 1> eigenvalues;
1224
1225 internal::TensorProductMatrixSymmetricSum::setup(mass_matrix,
1226 derivative_matrix,
1228 eigenvalues);
1229
1230 for (unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1231 i < mass_matrix[0].n_rows();
1232 ++i, ++v)
1233 {
1234 for (unsigned int j = 0; j < mass_matrix[0].n_cols(); ++j, ++m)
1235 {
1236 this->mass_matrices[m] = mass_matrix[0][i][j];
1237 this->derivative_matrices[m] = derivative_matrix[0][i][j];
1238 this->eigenvectors[m] = eigenvectors[0][i][j];
1239 }
1240
1241 this->eigenvalues[v] = eigenvalues[0][i];
1242 }
1243 };
1244
1245 if (compress_matrices == false)
1246 {
1247 // case 1) no compression requested
1248
1249 AssertDimension(cache.size(), 0);
1250 AssertDimension(indices.size(), 0);
1251
1252 this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1253 this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1254
1255 for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1256 {
1257 const auto &M = mass_and_derivative_matrices[i].first;
1258
1259 this->vector_ptr[i + 1] = M.n_rows();
1260 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1261 }
1262
1263 for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1264 {
1265 this->vector_ptr[i + 1] += this->vector_ptr[i];
1266 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1267 }
1268
1269 this->mass_matrices.resize_fast(matrix_ptr.back());
1270 this->derivative_matrices.resize_fast(matrix_ptr.back());
1271 this->eigenvectors.resize_fast(matrix_ptr.back());
1272 this->eigenvalues.resize_fast(vector_ptr.back());
1273
1274 for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1275 store(i, mass_and_derivative_matrices[i]);
1276
1277 mass_and_derivative_matrices.clear();
1278 }
1279 else if (cache.size() == indices.size())
1280 {
1281 // case 2) compression requested but none possible
1282
1283 this->vector_ptr.resize(cache.size() + 1);
1284 this->matrix_ptr.resize(cache.size() + 1);
1285
1286 std::map<unsigned int, MatrixPairType> inverted_cache;
1287
1288 for (const auto &i : cache)
1290
1291 for (unsigned int i = 0; i < indices.size(); ++i)
1292 {
1293 const auto &M = inverted_cache[indices[i]].first;
1294
1295 this->vector_ptr[i + 1] = M.n_rows();
1296 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1297 }
1298
1299 for (unsigned int i = 0; i < cache.size(); ++i)
1300 {
1301 this->vector_ptr[i + 1] += this->vector_ptr[i];
1302 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1303 }
1304
1305 this->mass_matrices.resize_fast(matrix_ptr.back());
1306 this->derivative_matrices.resize_fast(matrix_ptr.back());
1307 this->eigenvectors.resize_fast(matrix_ptr.back());
1308 this->eigenvalues.resize_fast(vector_ptr.back());
1309
1310 for (unsigned int i = 0; i < indices.size(); ++i)
1311 store(i, inverted_cache[indices[i]]);
1312
1313 indices.clear();
1314 cache.clear();
1315 }
1316 else
1317 {
1318 // case 3) compress
1319
1320 this->vector_ptr.resize(cache.size() + 1);
1321 this->matrix_ptr.resize(cache.size() + 1);
1322
1323 for (const auto &i : cache)
1324 {
1325 const auto &M = i.first.second.first;
1326
1327 this->vector_ptr[i.second + 1] = M.n_rows();
1328 this->matrix_ptr[i.second + 1] = M.n_rows() * M.n_cols();
1329 }
1330
1331 for (unsigned int i = 0; i < cache.size(); ++i)
1332 {
1333 this->vector_ptr[i + 1] += this->vector_ptr[i];
1334 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1335 }
1336
1337 this->mass_matrices.resize_fast(matrix_ptr.back());
1338 this->derivative_matrices.resize_fast(matrix_ptr.back());
1339 this->eigenvectors.resize_fast(matrix_ptr.back());
1340 this->eigenvalues.resize_fast(vector_ptr.back());
1341
1342 for (const auto &i : cache)
1343 store(i.second, i.first.second);
1344
1345 cache.clear();
1346 }
1347
1348 if (precompute_inverse_diagonal)
1349 {
1350 if (dim == 1)
1351 {
1352 // 1D case: simply invert 1D eigenvalues
1353 for (unsigned int i = 0; i < this->eigenvalues.size(); ++i)
1354 this->eigenvalues[i] = Number(1.0) / this->eigenvalues[i];
1355 std::swap(this->inverted_eigenvalues, eigenvalues);
1356 }
1357 else
1358 {
1359 // 2D and 3D case: we have 2 or 3 1d eigenvalues so that we
1360 // need to combine these
1361
1362 // step 1) if eigenvalues/eigenvectors are compressed, we
1363 // need to compress the diagonal (the combination of ev
1364 // indices) as well. This is an optional step.
1365 std::vector<unsigned int> indices_ev;
1366
1367 if (indices.size() > 0)
1368 {
1369 // 1a) create cache (ev indics -> diag index)
1370 const unsigned int n_cells = indices.size() / dim;
1371 std::map<std::array<unsigned int, dim>, unsigned int> cache_ev;
1372 std::vector<unsigned int> cache_ev_idx(n_cells);
1373
1374 for (unsigned int i = 0, c = 0; i < n_cells; ++i)
1375 {
1376 std::array<unsigned int, dim> id;
1377
1378 for (unsigned int d = 0; d < dim; ++d, ++c)
1379 id[d] = indices[c];
1380
1381 const auto id_ptr = cache_ev.find(id);
1382
1383 if (id_ptr == cache_ev.end())
1384 {
1385 const auto size = cache_ev.size();
1386 cache_ev_idx[i] = size;
1387 cache_ev[id] = size;
1388 }
1389 else
1390 {
1391 cache_ev_idx[i] = id_ptr->second;
1392 }
1393 }
1394
1395 // 1b) store diagonal indices for each cell
1396 std::vector<unsigned int> new_indices;
1397 new_indices.reserve(indices.size() / dim * (dim + 1));
1398
1399 for (unsigned int i = 0, c = 0; i < n_cells; ++i)
1400 {
1401 for (unsigned int d = 0; d < dim; ++d, ++c)
1402 new_indices.push_back(indices[c]);
1403 new_indices.push_back(cache_ev_idx[i]);
1404 }
1405
1406 // 1c) transpose cache (diag index -> ev indices)
1407 indices_ev.resize(cache_ev.size() * dim);
1408 for (const auto &entry : cache_ev)
1409 for (unsigned int d = 0; d < dim; ++d)
1410 indices_ev[entry.second * dim + d] = entry.first[d];
1411
1412 std::swap(this->indices, new_indices);
1413 }
1414
1415 // step 2) allocate memory and set pointers
1416 const unsigned int n_diag =
1417 ((indices_ev.size() > 0) ? indices_ev.size() :
1418 (matrix_ptr.size() - 1)) /
1419 dim;
1420
1421 std::vector<unsigned int> new_vector_ptr(n_diag + 1, 0);
1422 std::vector<unsigned int> new_vector_n_rows_1d(n_diag, 0);
1423
1424 for (unsigned int i = 0; i < n_diag; ++i)
1425 {
1426 const unsigned int c = (indices_ev.size() > 0) ?
1427 indices_ev[dim * i + 0] :
1428 (dim * i + 0);
1429
1430 const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c];
1431
1432 new_vector_n_rows_1d[i] = n_rows;
1433 new_vector_ptr[i + 1] = Utilities::pow(n_rows, dim);
1434 }
1435
1436 for (unsigned int i = 0; i < n_diag; ++i)
1437 new_vector_ptr[i + 1] += new_vector_ptr[i];
1438
1439 this->inverted_eigenvalues.resize(new_vector_ptr.back());
1440
1441 // step 3) loop over all unique diagonal entries and invert
1442 for (unsigned int i = 0; i < n_diag; ++i)
1443 {
1444 std::array<Number *, dim> evs;
1445
1446 for (unsigned int d = 0; d < dim; ++d)
1447 evs[d] =
1448 &this
1449 ->eigenvalues[this->vector_ptr[(indices_ev.size() > 0) ?
1450 indices_ev[dim * i + d] :
1451 (dim * i + d)]];
1452
1453 const unsigned int mm = new_vector_n_rows_1d[i];
1454 if (dim == 2)
1455 {
1456 for (unsigned int i1 = 0, c = 0; i1 < mm; ++i1)
1457 for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1458 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1459 Number(1.0) / (evs[1][i1] + evs[0][i0]);
1460 }
1461 else
1462 {
1463 for (unsigned int i2 = 0, c = 0; i2 < mm; ++i2)
1464 for (unsigned int i1 = 0; i1 < mm; ++i1)
1465 for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1466 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1467 Number(1.0) / (evs[2][i2] + evs[1][i1] + evs[0][i0]);
1468 }
1469 }
1470
1471 // step 4) clean up
1472 std::swap(this->vector_ptr, new_vector_ptr);
1473 std::swap(this->vector_n_rows_1d, new_vector_n_rows_1d);
1474 }
1475
1476 this->eigenvalues.clear();
1477 }
1478}
1479
1480
1481
1482template <int dim, typename Number, int n_rows_1d>
1483void
1485 apply_inverse(const unsigned int index,
1487 const ArrayView<const Number> &src_in) const
1488{
1489 Number *dst = dst_in.begin();
1490 const Number *src = src_in.begin();
1491
1492 if (this->eigenvalues.empty() == false)
1493 {
1494 std::array<const Number *, dim> eigenvectors;
1495 std::array<const Number *, dim> eigenvalues;
1496 unsigned int n_rows_1d_non_templated = 0;
1497
1498 for (unsigned int d = 0; d < dim; ++d)
1499 {
1500 const unsigned int translated_index =
1501 (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d);
1502
1503 eigenvectors[d] =
1504 this->eigenvectors.data() + matrix_ptr[translated_index];
1505 eigenvalues[d] =
1506 this->eigenvalues.data() + vector_ptr[translated_index];
1508 vector_ptr[translated_index + 1] - vector_ptr[translated_index];
1509 }
1510
1511 if (n_rows_1d != -1)
1512 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1513 n_rows_1d == -1 ? 0 : n_rows_1d>(
1515 else
1516 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1518 }
1519 else
1520 {
1521 std::array<const Number *, dim> eigenvectors;
1522 const Number *inverted_eigenvalues = nullptr;
1523 unsigned int n_rows_1d_non_templated = 0;
1524
1525 for (unsigned int d = 0; d < dim; ++d)
1526 {
1527 const unsigned int translated_index =
1528 (indices.size() > 0) ?
1529 indices[((dim == 1) ? 1 : (dim + 1)) * index + d] :
1530 (dim * index + d);
1531
1532 eigenvectors[d] =
1533 this->eigenvectors.data() + matrix_ptr[translated_index];
1534 }
1535
1536 {
1537 const unsigned int translated_index =
1538 ((indices.size() > 0) && (dim != 1)) ?
1539 indices[(dim + 1) * index + dim] :
1540 index;
1541
1542 inverted_eigenvalues =
1543 this->inverted_eigenvalues.data() + vector_ptr[translated_index];
1545 (dim == 1) ?
1546 (vector_ptr[translated_index + 1] - vector_ptr[translated_index]) :
1547 vector_n_rows_1d[translated_index];
1548 }
1549
1550 if (n_rows_1d != -1)
1551 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1552 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
1553 src,
1556 {},
1557 inverted_eigenvalues);
1558 else
1559 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1560 dst,
1561 src,
1564 {},
1565 inverted_eigenvalues);
1566 }
1567}
1568
1569
1570
1571template <int dim, typename Number, int n_rows_1d>
1572std::size_t
1574 memory_consumption() const
1575{
1578 MemoryConsumption::memory_consumption(derivative_matrices) +
1583}
1584
1585
1586
1587template <int dim, typename Number, int n_rows_1d>
1588std::size_t
1590 storage_size() const
1591{
1592 if (matrix_ptr.empty())
1593 return 0; // if not initialized
1594
1595 return matrix_ptr.size() - 1;
1596}
1597
1598
1599
1600#endif
1601
1603
1604#endif
void apply_inverse(const unsigned int index, const ArrayView< Number > &dst_in, const ArrayView< const Number > &src_in) const
std::vector< MatrixPairType > mass_and_derivative_matrices
std::pair< std::bitset<::internal::VectorizedArrayTrait< Number >::width()>, MatrixPairType > MatrixPairTypeWithMask
std::pair< Table< 2, Number >, Table< 2, Number > > MatrixPairType
void reserve(const unsigned int size)
void insert(const unsigned int index, const T &Ms, const T &Ks)
std::map< MatrixPairTypeWithMask, unsigned int, internal::TensorProductMatrixSymmetricSum::MatrixPairComparator< Number > > cache
TensorProductMatrixSymmetricSumCollection(const AdditionalData &additional_data=AdditionalData())
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src, AlignedVector< Number > &tmp) const
std::array< Table< 2, Number >, dim > eigenvectors
std::array< Table< 2, Number >, dim > derivative_matrix
void reinit(const T &mass_matrix, const T &derivative_matrix)
void apply_inverse(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
std::size_t memory_consumption() const
std::array< Table< 2, Number >, dim > mass_matrix
std::array< AlignedVector< Number >, dim > eigenvalues
TensorProductMatrixSymmetricSum(const T &mass_matrix, const T &derivative_matrix)
constexpr void clear()
friend class Tensor
Definition tensor.h:865
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
std::size_t size
Definition mpi.cc:739
@ matrix
Contents is actually a matrix.
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
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:967
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14889
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
AdditionalData(const bool compress_matrices=true, const bool precompute_inverse_diagonal=true)
bool operator()(const MatrixPairType &left, const MatrixPairType &right) const
std::pair< std::bitset< width >, std::pair< Table< 2, Number >, Table< 2, Number > > > MatrixPairType
static constexpr std::size_t width()
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)