Reference documentation for deal.II version GIT f0f8c7fe18 2023-03-21 21:25: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\}}\)
tensor_product_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2021 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_tensor_product_matrix_h
17 #define dealii_tensor_product_matrix_h
18 
19 
20 #include <deal.II/base/config.h>
21 
24 #include <deal.II/base/mutex.h>
26 
28 
30 
32 
33 // Forward declarations
34 #ifndef DOXYGEN
35 template <typename>
36 class Vector;
37 template <typename>
38 class FullMatrix;
39 #endif
40 
115 template <int dim, typename Number, int n_rows_1d = -1>
117 {
118 public:
123  using value_type = Number;
124 
129  static constexpr int n_rows_1d_static = n_rows_1d;
130 
135 
140  template <typename T>
142  const T &derivative_matrix);
143 
161  template <typename T>
162  void
164 
170  unsigned int
171  m() const;
172 
178  unsigned int
179  n() const;
180 
194  void
195  vmult(const ArrayView<Number> &dst, const ArrayView<const Number> &src) const;
196 
202  void
204  const ArrayView<const Number> &src,
205  AlignedVector<Number> & tmp) const;
206 
220  void
222  const ArrayView<const Number> &src) const;
223 
229  void
231  const ArrayView<const Number> &src,
232  AlignedVector<Number> & tmp) const;
233 
237  std::size_t
239 
240 protected:
244  std::array<Table<2, Number>, dim> mass_matrix;
245 
249  std::array<Table<2, Number>, dim> derivative_matrix;
250 
255  std::array<AlignedVector<Number>, dim> eigenvalues;
256 
261  std::array<Table<2, Number>, dim> eigenvectors;
262 
263 private:
268 
273 };
274 
275 
276 
277 namespace internal
278 {
280  {
281  template <typename Number>
283  {
284  using MatrixPairType = std::pair<Table<2, Number>, Table<2, Number>>;
288  static constexpr std::size_t width = VectorizedArrayTrait::width;
289 
291  : eps(std::sqrt(std::numeric_limits<ScalarNumber>::epsilon()))
292  {}
293 
294  bool
295  operator()(const MatrixPairType &left, const MatrixPairType &right) const
296  {
297  const auto &M_0 = left.first;
298  const auto &K_0 = left.second;
299  const auto &M_1 = right.first;
300  const auto &K_1 = right.second;
301 
302  std::bitset<width> mask;
303 
304  for (unsigned int v = 0; v < width; ++v)
305  {
306  ScalarNumber a = 0.0;
307  ScalarNumber b = 0.0;
308 
309  for (unsigned int i = 0; i < M_0.size(0); ++i)
310  for (unsigned int j = 0; j < M_0.size(1); ++j)
311  {
312  a += std::abs(VectorizedArrayTrait::get(M_0[i][j], v));
313  a += std::abs(VectorizedArrayTrait::get(K_0[i][j], v));
314  b += std::abs(VectorizedArrayTrait::get(M_1[i][j], v));
315  b += std::abs(VectorizedArrayTrait::get(K_1[i][j], v));
316  }
317 
318  mask[v] = (a != 0.0) && (b != 0.0);
319  }
320 
321  const FloatingPointComparator<Number> comparator(
322  eps, false /*use relative tolerance*/, mask);
323 
324  if (comparator(M_0, M_1))
325  return true;
326  else if (comparator(M_1, M_0))
327  return false;
328  else if (comparator(K_0, K_1))
329  return true;
330  else
331  return false;
332  }
333 
334  private:
336  };
337  } // namespace TensorProductMatrixSymmetricSum
338 } // namespace internal
339 
340 
341 
365 template <int dim, typename Number, int n_rows_1d = -1>
367 {
368  using MatrixPairType = std::pair<Table<2, Number>, Table<2, Number>>;
369 
370 public:
375  {
379  AdditionalData(const bool compress_matrices = true);
380 
385  };
386 
391  const AdditionalData &additional_data = AdditionalData());
392 
397  void
398  reserve(const unsigned int size);
399 
405  template <typename T>
406  void
407  insert(const unsigned int index, const T &Ms, const T &Ks);
408 
413  void
415 
419  void
420  apply_inverse(const unsigned int index,
421  const ArrayView<Number> & dst_in,
422  const ArrayView<const Number> &src_in,
423  AlignedVector<Number> & tmp_array) const;
424 
428  std::size_t
430 
439  std::size_t
440  storage_size() const;
441 
442 private:
446  const bool compress_matrices;
447 
452  std::vector<MatrixPairType> mass_and_derivative_matrices;
453 
458  std::map<
460  unsigned int,
463 
471  std::vector<unsigned int> indices;
472 
477 
482 
487 
492 
496  std::vector<unsigned int> vector_ptr;
497 
501  std::vector<unsigned int> matrix_ptr;
502 };
503 
504 
505 /*----------------------- Inline functions ----------------------------------*/
506 
507 #ifndef DOXYGEN
508 
509 namespace internal
510 {
512  {
521  template <typename Number>
522  void
523  spectral_assembly(const Number * mass_matrix,
524  const Number * derivative_matrix,
525  const unsigned int n_rows,
526  const unsigned int n_cols,
527  Number * eigenvalues,
528  Number * eigenvectors)
529  {
530  Assert(n_rows == n_cols, ExcNotImplemented());
531 
532  std::vector<bool> constrained_dofs(n_rows, false);
533 
534  for (unsigned int i = 0; i < n_rows; ++i)
535  {
536  if (mass_matrix[i + i * n_rows] == 0.0)
537  {
538  Assert(derivative_matrix[i + i * n_rows] == 0.0,
539  ExcInternalError());
540 
541  for (unsigned int j = 0; j < n_rows; ++j)
542  {
543  Assert(derivative_matrix[i + j * n_rows] == 0,
544  ExcInternalError());
545  Assert(derivative_matrix[j + i * n_rows] == 0,
546  ExcInternalError());
547  }
548 
549  constrained_dofs[i] = true;
550  }
551  }
552 
553  const auto transpose_fill_nm = [&constrained_dofs](Number * out,
554  const Number * in,
555  const unsigned int n,
556  const unsigned int m) {
557  for (unsigned int mm = 0, c = 0; mm < m; ++mm)
558  for (unsigned int nn = 0; nn < n; ++nn, ++c)
559  out[mm + nn * m] =
560  (mm == nn && constrained_dofs[mm]) ? Number(1.0) : in[c];
561  };
562 
563  std::vector<::Vector<Number>> eigenvecs(n_rows);
564  LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
565  LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);
566 
567  transpose_fill_nm(&(mass_copy(0, 0)), mass_matrix, n_rows, n_cols);
568  transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
569 
570  deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy,
571  eigenvecs);
572  AssertDimension(eigenvecs.size(), n_rows);
573  for (unsigned int i = 0, c = 0; i < n_rows; ++i)
574  for (unsigned int j = 0; j < n_cols; ++j, ++c)
575  if (constrained_dofs[i] == false)
576  eigenvectors[c] = eigenvecs[j][i];
577 
578  for (unsigned int i = 0; i < n_rows; ++i, ++eigenvalues)
579  *eigenvalues = deriv_copy.eigenvalue(i).real();
580  }
581 
582 
583 
584  template <std::size_t dim, typename Number>
585  inline void
586  setup(const std::array<Table<2, Number>, dim> &mass_matrix,
587  const std::array<Table<2, Number>, dim> &derivative_matrix,
588  std::array<Table<2, Number>, dim> & eigenvectors,
589  std::array<AlignedVector<Number>, dim> & eigenvalues)
590  {
591  const unsigned int n_rows_1d = mass_matrix[0].n_cols();
592  (void)n_rows_1d;
593 
594  for (unsigned int dir = 0; dir < dim; ++dir)
595  {
596  AssertDimension(n_rows_1d, mass_matrix[dir].n_cols());
597  AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
598  AssertDimension(mass_matrix[dir].n_rows(),
599  derivative_matrix[dir].n_rows());
600  AssertDimension(mass_matrix[dir].n_rows(),
601  derivative_matrix[dir].n_cols());
602 
603  eigenvectors[dir].reinit(mass_matrix[dir].n_cols(),
604  mass_matrix[dir].n_rows());
605  eigenvalues[dir].resize(mass_matrix[dir].n_cols());
606  internal::TensorProductMatrixSymmetricSum::spectral_assembly<Number>(
607  &(mass_matrix[dir](0, 0)),
608  &(derivative_matrix[dir](0, 0)),
609  mass_matrix[dir].n_rows(),
610  mass_matrix[dir].n_cols(),
611  eigenvalues[dir].begin(),
612  &(eigenvectors[dir](0, 0)));
613  }
614  }
615 
616 
617 
618  template <std::size_t dim, typename Number, std::size_t n_lanes>
619  inline void
620  setup(
621  const std::array<Table<2, VectorizedArray<Number, n_lanes>>, dim>
622  &mass_matrix,
623  const std::array<Table<2, VectorizedArray<Number, n_lanes>>, dim>
624  &derivative_matrix,
627  &eigenvalues)
628  {
629  const unsigned int n_rows_1d = mass_matrix[0].n_cols();
630  constexpr unsigned int macro_size =
632  const std::size_t nm_flat_size_max = n_rows_1d * n_rows_1d * macro_size;
633  const std::size_t n_flat_size_max = n_rows_1d * macro_size;
634 
635  std::vector<Number> mass_matrix_flat;
636  std::vector<Number> deriv_matrix_flat;
637  std::vector<Number> eigenvalues_flat;
638  std::vector<Number> eigenvectors_flat;
639  mass_matrix_flat.resize(nm_flat_size_max);
640  deriv_matrix_flat.resize(nm_flat_size_max);
641  eigenvalues_flat.resize(n_flat_size_max);
642  eigenvectors_flat.resize(nm_flat_size_max);
643  std::array<unsigned int, macro_size> offsets_nm;
644  std::array<unsigned int, macro_size> offsets_n;
645  for (unsigned int dir = 0; dir < dim; ++dir)
646  {
647  AssertDimension(n_rows_1d, mass_matrix[dir].n_cols());
648  AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
649  AssertDimension(mass_matrix[dir].n_rows(),
650  derivative_matrix[dir].n_rows());
651  AssertDimension(mass_matrix[dir].n_rows(),
652  derivative_matrix[dir].n_cols());
653 
654  const unsigned int n_rows = mass_matrix[dir].n_rows();
655  const unsigned int n_cols = mass_matrix[dir].n_cols();
656  const unsigned int nm = n_rows * n_cols;
657  for (unsigned int vv = 0; vv < macro_size; ++vv)
658  offsets_nm[vv] = nm * vv;
659 
661  nm,
662  &(mass_matrix[dir](0, 0)),
663  offsets_nm.cbegin(),
664  mass_matrix_flat.data());
666  nm,
667  &(derivative_matrix[dir](0, 0)),
668  offsets_nm.cbegin(),
669  deriv_matrix_flat.data());
670 
671  const Number *mass_cbegin = mass_matrix_flat.data();
672  const Number *deriv_cbegin = deriv_matrix_flat.data();
673  Number * eigenvec_begin = eigenvectors_flat.data();
674  Number * eigenval_begin = eigenvalues_flat.data();
675  for (unsigned int lane = 0; lane < macro_size; ++lane)
676  internal::TensorProductMatrixSymmetricSum::spectral_assembly<
677  Number>(mass_cbegin + nm * lane,
678  deriv_cbegin + nm * lane,
679  n_rows,
680  n_cols,
681  eigenval_begin + n_rows * lane,
682  eigenvec_begin + nm * lane);
683 
684  eigenvalues[dir].resize(n_rows);
685  eigenvectors[dir].reinit(n_rows, n_cols);
686  for (unsigned int vv = 0; vv < macro_size; ++vv)
687  offsets_n[vv] = n_rows * vv;
689  eigenvalues_flat.data(),
690  offsets_n.cbegin(),
691  eigenvalues[dir].begin());
693  eigenvectors_flat.data(),
694  offsets_nm.cbegin(),
695  &(eigenvectors[dir](0, 0)));
696  }
697  }
698 
699 
700 
701  template <std::size_t dim, typename Number>
702  inline std::array<Table<2, Number>, dim>
703  convert(const std::array<Table<2, Number>, dim> &mass_matrix)
704  {
705  return mass_matrix;
706  }
707 
708 
709 
710  template <std::size_t dim, typename Number>
711  inline std::array<Table<2, Number>, dim>
712  convert(const std::array<FullMatrix<Number>, dim> &mass_matrix)
713  {
714  std::array<Table<2, Number>, dim> mass_copy;
715 
716  std::transform(mass_matrix.cbegin(),
717  mass_matrix.cend(),
718  mass_copy.begin(),
719  [](const FullMatrix<Number> &m) -> Table<2, Number> {
720  return m;
721  });
722 
723  return mass_copy;
724  }
725 
726 
727 
728  template <std::size_t dim, typename Number>
729  inline std::array<Table<2, Number>, dim>
730  convert(const Table<2, Number> &matrix)
731  {
732  std::array<Table<2, Number>, dim> matrices;
733 
734  std::fill(matrices.begin(), matrices.end(), matrix);
735 
736  return matrices;
737  }
738 
739 
740 
741  template <int n_rows_1d_templated, std::size_t dim, typename Number>
742  void
743  vmult(Number * dst,
744  const Number * src,
745  AlignedVector<Number> & tmp,
746  const unsigned int n_rows_1d_non_templated,
747  const std::array<const Number *, dim> &mass_matrix,
748  const std::array<const Number *, dim> &derivative_matrix)
749  {
750  const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
751  n_rows_1d_non_templated :
752  n_rows_1d_templated;
753  const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
754 
755  tmp.resize_fast(n * 2);
756  Number *t = tmp.begin();
757 
759  dim,
760  n_rows_1d_templated,
761  n_rows_1d_templated,
762  Number>
763  eval({}, {}, {}, n_rows_1d, n_rows_1d);
764 
765  if (dim == 1)
766  {
767  const Number *A = derivative_matrix[0];
768  eval.template apply<0, false, false>(A, src, dst);
769  }
770 
771  else if (dim == 2)
772  {
773  const Number *A0 = derivative_matrix[0];
774  const Number *M0 = mass_matrix[0];
775  const Number *A1 = derivative_matrix[1];
776  const Number *M1 = mass_matrix[1];
777  eval.template apply<0, false, false>(M0, src, t);
778  eval.template apply<1, false, false>(A1, t, dst);
779  eval.template apply<0, false, false>(A0, src, t);
780  eval.template apply<1, false, true>(M1, t, dst);
781  }
782 
783  else if (dim == 3)
784  {
785  const Number *A0 = derivative_matrix[0];
786  const Number *M0 = mass_matrix[0];
787  const Number *A1 = derivative_matrix[1];
788  const Number *M1 = mass_matrix[1];
789  const Number *A2 = derivative_matrix[2];
790  const Number *M2 = mass_matrix[2];
791  eval.template apply<0, false, false>(M0, src, t + n);
792  eval.template apply<1, false, false>(M1, t + n, t);
793  eval.template apply<2, false, false>(A2, t, dst);
794  eval.template apply<1, false, false>(A1, t + n, t);
795  eval.template apply<0, false, false>(A0, src, t + n);
796  eval.template apply<1, false, true>(M1, t + n, t);
797  eval.template apply<2, false, true>(M2, t, dst);
798  }
799 
800  else
801  AssertThrow(false, ExcNotImplemented());
802  }
803 
804 
805 
806  template <int n_rows_1d_templated, std::size_t dim, typename Number>
807  void
808  apply_inverse(Number * dst,
809  const Number * src,
811  const unsigned int n_rows_1d_non_templated,
812  const std::array<const Number *, dim> &eigenvectors,
813  const std::array<const Number *, dim> &eigenvalues)
814  {
815  const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
816  n_rows_1d_non_templated :
817  n_rows_1d_templated;
818  const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
819 
820  tmp.resize_fast(n);
821  Number *t = tmp.begin();
822 
824  dim,
825  n_rows_1d_templated,
826  n_rows_1d_templated,
827  Number>
828  eval({}, {}, {}, n_rows_1d, n_rows_1d);
829 
830  // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
831  // --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
832  // --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
833  // while the eigenvectors are stored column-wise in S, i.e.
834  // rows correspond to dofs whereas columns to eigenvalue indices!
835  if (dim == 1)
836  {
837  const Number *S = eigenvectors[0];
838  eval.template apply<0, true, false>(S, src, t);
839  for (unsigned int i = 0; i < n_rows_1d; ++i)
840  t[i] /= eigenvalues[0][i];
841  eval.template apply<0, false, false>(S, t, 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, t);
849  eval.template apply<1, true, false>(S1, t, dst);
850  for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
851  for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
852  dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
853  eval.template apply<0, false, false>(S0, dst, t);
854  eval.template apply<1, false, false>(S1, t, dst);
855  }
856 
857  else if (dim == 3)
858  {
859  const Number *S0 = eigenvectors[0];
860  const Number *S1 = eigenvectors[1];
861  const Number *S2 = eigenvectors[2];
862  eval.template apply<0, true, false>(S0, src, t);
863  eval.template apply<1, true, false>(S1, t, dst);
864  eval.template apply<2, true, false>(S2, dst, t);
865  for (unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
866  for (unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
867  for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
868  t[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] +
869  eigenvalues[0][i0]);
870  eval.template apply<0, false, false>(S0, t, dst);
871  eval.template apply<1, false, false>(S1, dst, t);
872  eval.template apply<2, false, false>(S2, t, dst);
873  }
874 
875  else
876  Assert(false, ExcNotImplemented());
877  }
878 
879 
880 
881  template <int n_rows_1d_templated, std::size_t dim, typename Number>
882  void
883  select_vmult(Number * dst,
884  const Number * src,
885  AlignedVector<Number> & tmp,
886  const unsigned int n_rows_1d,
887  const std::array<const Number *, dim> &mass_matrix,
888  const std::array<const Number *, dim> &derivative_matrix);
889 
890 
891 
892  template <int n_rows_1d_templated, std::size_t dim, typename Number>
893  void
894  select_apply_inverse(Number * dst,
895  const Number * src,
896  AlignedVector<Number> & tmp,
897  const unsigned int n_rows_1d,
898  const std::array<const Number *, dim> &eigenvectors,
899  const std::array<const Number *, dim> &eigenvalues);
900  } // namespace TensorProductMatrixSymmetricSum
901 } // namespace internal
902 
903 
904 template <int dim, typename Number, int n_rows_1d>
905 inline unsigned int
907 {
908  unsigned int m = mass_matrix[0].n_rows();
909  for (unsigned int d = 1; d < dim; ++d)
910  m *= mass_matrix[d].n_rows();
911  return m;
912 }
913 
914 
915 
916 template <int dim, typename Number, int n_rows_1d>
917 inline unsigned int
919 {
920  unsigned int n = mass_matrix[0].n_cols();
921  for (unsigned int d = 1; d < dim; ++d)
922  n *= mass_matrix[d].n_cols();
923  return n;
924 }
925 
926 
927 
928 template <int dim, typename Number, int n_rows_1d>
929 inline void
931  const ArrayView<Number> & dst_view,
932  const ArrayView<const Number> &src_view) const
933 {
934  std::lock_guard<std::mutex> lock(this->mutex);
935  this->vmult(dst_view, src_view, this->tmp_array);
936 }
937 
938 
939 
940 template <int dim, typename Number, int n_rows_1d>
941 inline void
943  const ArrayView<Number> & dst_view,
944  const ArrayView<const Number> &src_view,
945  AlignedVector<Number> & tmp_array) const
946 {
947  AssertDimension(dst_view.size(), this->m());
948  AssertDimension(src_view.size(), this->n());
949 
950  Number * dst = dst_view.begin();
951  const Number *src = src_view.begin();
952 
953  std::array<const Number *, dim> mass_matrix, derivative_matrix;
954 
955  for (unsigned int d = 0; d < dim; ++d)
956  {
957  mass_matrix[d] = &this->mass_matrix[d](0, 0);
958  derivative_matrix[d] = &this->derivative_matrix[d](0, 0);
959  }
960 
961  const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
962 
963  if (n_rows_1d != -1)
964  internal::TensorProductMatrixSymmetricSum::vmult<
965  n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
966  src,
967  tmp_array,
968  n_rows_1d_non_templated,
969  mass_matrix,
970  derivative_matrix);
971  else
972  internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
973  dst,
974  src,
975  tmp_array,
976  n_rows_1d_non_templated,
977  mass_matrix,
978  derivative_matrix);
979 }
980 
981 
982 
983 template <int dim, typename Number, int n_rows_1d>
984 inline void
986  const ArrayView<Number> & dst_view,
987  const ArrayView<const Number> &src_view) const
988 {
989  std::lock_guard<std::mutex> lock(this->mutex);
990  this->apply_inverse(dst_view, src_view, this->tmp_array);
991 }
992 
993 
994 
995 template <int dim, typename Number, int n_rows_1d>
996 inline void
998  const ArrayView<Number> & dst_view,
999  const ArrayView<const Number> &src_view,
1000  AlignedVector<Number> & tmp_array) const
1001 {
1002  AssertDimension(dst_view.size(), this->n());
1003  AssertDimension(src_view.size(), this->m());
1004 
1005  Number * dst = dst_view.begin();
1006  const Number *src = src_view.begin();
1007 
1008  std::array<const Number *, dim> eigenvectors, eigenvalues;
1009 
1010  for (unsigned int d = 0; d < dim; ++d)
1011  {
1012  eigenvectors[d] = &this->eigenvectors[d](0, 0);
1013  eigenvalues[d] = this->eigenvalues[d].data();
1014  }
1015 
1016  const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
1017 
1018  if (n_rows_1d != -1)
1019  internal::TensorProductMatrixSymmetricSum::apply_inverse<
1020  n_rows_1d == -1 ? 0 : n_rows_1d>(
1021  dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1022  else
1023  internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1024  dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1025 }
1026 
1027 
1028 
1029 template <int dim, typename Number, int n_rows_1d>
1030 std::size_t
1032  const
1033 {
1035  MemoryConsumption::memory_consumption(derivative_matrix) +
1039 }
1040 
1041 
1042 
1043 template <int dim, typename Number, int n_rows_1d>
1044 template <typename T>
1047  const T &derivative_matrix)
1048 {
1049  reinit(mass_matrix, derivative_matrix);
1050 }
1051 
1052 
1053 
1054 template <int dim, typename Number, int n_rows_1d>
1055 template <typename T>
1056 inline void
1058  const T &mass_matrix,
1059  const T &derivative_matrix)
1060 {
1061  this->mass_matrix =
1062  internal::TensorProductMatrixSymmetricSum::convert<dim>(mass_matrix);
1063  this->derivative_matrix =
1064  internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1065 
1066  internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1067  this->derivative_matrix,
1068  this->eigenvectors,
1069  this->eigenvalues);
1070 }
1071 
1072 
1073 
1074 template <int dim, typename Number, int n_rows_1d>
1076  AdditionalData::AdditionalData(const bool compress_matrices)
1077  : compress_matrices(compress_matrices)
1078 {}
1079 
1080 
1081 
1082 template <int dim, typename Number, int n_rows_1d>
1085  const AdditionalData &additional_data)
1086  : compress_matrices(additional_data.compress_matrices)
1087 {}
1088 
1089 
1090 
1091 template <int dim, typename Number, int n_rows_1d>
1092 void
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 
1104 template <int dim, typename Number, int n_rows_1d>
1105 template <typename T>
1106 void
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  const MatrixPairType matrix(Ms[d], Ks[d]);
1120 
1121  if (compress_matrices == false)
1122  {
1123  mass_and_derivative_matrices[index * dim + d] = matrix;
1124  }
1125  else
1126  {
1127  const auto ptr = cache.find(matrix);
1128 
1129  if (ptr != cache.end())
1130  indices[index * dim + d] = ptr->second;
1131  else
1132  {
1133  const auto size = cache.size();
1134  indices[index * dim + d] = size;
1135  cache[matrix] = size;
1136  }
1137  }
1138  }
1139 }
1140 
1141 
1142 
1143 template <int dim, typename Number, int n_rows_1d>
1144 void
1146 {
1147  const auto store = [&](const unsigned int index,
1148  const MatrixPairType &M_and_K) {
1149  std::array<Table<2, Number>, 1> mass_matrix;
1150  mass_matrix[0] = M_and_K.first;
1151 
1152  std::array<Table<2, Number>, 1> derivative_matrix;
1153  derivative_matrix[0] = M_and_K.second;
1154 
1155  std::array<Table<2, Number>, 1> eigenvectors;
1156  std::array<AlignedVector<Number>, 1> eigenvalues;
1157 
1158  internal::TensorProductMatrixSymmetricSum::setup(mass_matrix,
1159  derivative_matrix,
1160  eigenvectors,
1161  eigenvalues);
1162 
1163  for (unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1164  i < mass_matrix[0].n_rows();
1165  ++i, ++v)
1166  {
1167  for (unsigned int j = 0; j < mass_matrix[0].n_cols(); ++j, ++m)
1168  {
1169  this->mass_matrices[m] = mass_matrix[0][i][j];
1170  this->derivative_matrices[m] = derivative_matrix[0][i][j];
1171  this->eigenvectors[m] = eigenvectors[0][i][j];
1172  }
1173 
1174  this->eigenvalues[v] = eigenvalues[0][i];
1175  }
1176  };
1177 
1178  if (compress_matrices == false)
1179  {
1180  // case 1) no compression requested
1181 
1182  AssertDimension(cache.size(), 0);
1183  AssertDimension(indices.size(), 0);
1184 
1185  this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1186  this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1187 
1188  for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1189  {
1190  const auto &M = mass_and_derivative_matrices[i].first;
1191 
1192  this->vector_ptr[i + 1] = M.n_rows();
1193  this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1194  }
1195 
1196  for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1197  {
1198  this->vector_ptr[i + 1] += this->vector_ptr[i];
1199  this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1200  }
1201 
1202  this->mass_matrices.resize_fast(matrix_ptr.back());
1203  this->derivative_matrices.resize_fast(matrix_ptr.back());
1204  this->eigenvectors.resize_fast(matrix_ptr.back());
1205  this->eigenvalues.resize_fast(vector_ptr.back());
1206 
1207  for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1208  store(i, mass_and_derivative_matrices[i]);
1209 
1210  mass_and_derivative_matrices.clear();
1211  }
1212  else if (cache.size() == indices.size())
1213  {
1214  // case 2) compression requested but none possible
1215 
1216  this->vector_ptr.resize(cache.size() + 1);
1217  this->matrix_ptr.resize(cache.size() + 1);
1218 
1219  std::map<unsigned int, MatrixPairType> inverted_cache;
1220 
1221  for (const auto &i : cache)
1222  inverted_cache[i.second] = i.first;
1223 
1224  for (unsigned int i = 0; i < indices.size(); ++i)
1225  {
1226  const auto &M = inverted_cache[indices[i]].first;
1227 
1228  this->vector_ptr[i + 1] = M.n_rows();
1229  this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1230  }
1231 
1232  for (unsigned int i = 0; i < cache.size(); ++i)
1233  {
1234  this->vector_ptr[i + 1] += this->vector_ptr[i];
1235  this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1236  }
1237 
1238  this->mass_matrices.resize_fast(matrix_ptr.back());
1239  this->derivative_matrices.resize_fast(matrix_ptr.back());
1240  this->eigenvectors.resize_fast(matrix_ptr.back());
1241  this->eigenvalues.resize_fast(vector_ptr.back());
1242 
1243  for (unsigned int i = 0; i < indices.size(); ++i)
1244  store(i, inverted_cache[indices[i]]);
1245 
1246  indices.clear();
1247  cache.clear();
1248  }
1249  else
1250  {
1251  // case 2) compression requested but none possible
1252 
1253  this->vector_ptr.resize(cache.size() + 1);
1254  this->matrix_ptr.resize(cache.size() + 1);
1255 
1256  for (const auto &i : cache)
1257  {
1258  const auto &M = i.first.first;
1259 
1260  this->vector_ptr[i.second + 1] = M.n_rows();
1261  this->matrix_ptr[i.second + 1] = M.n_rows() * M.n_cols();
1262  }
1263 
1264  for (unsigned int i = 0; i < cache.size(); ++i)
1265  {
1266  this->vector_ptr[i + 1] += this->vector_ptr[i];
1267  this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1268  }
1269 
1270  this->mass_matrices.resize_fast(matrix_ptr.back());
1271  this->derivative_matrices.resize_fast(matrix_ptr.back());
1272  this->eigenvectors.resize_fast(matrix_ptr.back());
1273  this->eigenvalues.resize_fast(vector_ptr.back());
1274 
1275  for (const auto &i : cache)
1276  store(i.second, i.first);
1277 
1278  cache.clear();
1279  }
1280 }
1281 
1282 
1283 
1284 template <int dim, typename Number, int n_rows_1d>
1285 void
1287  apply_inverse(const unsigned int index,
1288  const ArrayView<Number> & dst_in,
1289  const ArrayView<const Number> &src_in,
1290  AlignedVector<Number> & tmp_array) const
1291 {
1292  Number * dst = dst_in.begin();
1293  const Number *src = src_in.begin();
1294 
1295  std::array<const Number *, dim> eigenvectors, eigenvalues;
1296  unsigned int n_rows_1d_non_templated = 0;
1297 
1298  for (unsigned int d = 0; d < dim; ++d)
1299  {
1300  const unsigned int translated_index =
1301  (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d);
1302 
1303  eigenvectors[d] =
1304  this->eigenvectors.data() + matrix_ptr[translated_index];
1305  eigenvalues[d] = this->eigenvalues.data() + vector_ptr[translated_index];
1306  n_rows_1d_non_templated =
1307  vector_ptr[translated_index + 1] - vector_ptr[translated_index];
1308  }
1309 
1310  if (n_rows_1d != -1)
1311  internal::TensorProductMatrixSymmetricSum::apply_inverse<
1312  n_rows_1d == -1 ? 0 : n_rows_1d>(
1313  dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1314  else
1315  internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1316  dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1317 }
1318 
1319 
1320 
1321 template <int dim, typename Number, int n_rows_1d>
1322 std::size_t
1324  memory_consumption() const
1325 {
1326  return MemoryConsumption::memory_consumption(indices) +
1327  MemoryConsumption::memory_consumption(mass_matrices) +
1328  MemoryConsumption::memory_consumption(derivative_matrices) +
1333 }
1334 
1335 
1336 
1337 template <int dim, typename Number, int n_rows_1d>
1338 std::size_t
1340  storage_size() const
1341 {
1342  if (matrix_ptr.size() == 0)
1343  return 0; // if not initialized
1344 
1345  return matrix_ptr.size() - 1;
1346 }
1347 
1348 
1349 
1350 #endif
1351 
1353 
1354 #endif
void resize_fast(const size_type new_size)
iterator begin()
iterator begin() const
Definition: array_view.h:594
std::size_t size() const
Definition: array_view.h:576
std::complex< typename numbers::NumberTraits< number >::real_type > eigenvalue(const size_type i) const
void compute_generalized_eigenvalues_symmetric(LAPACKFullMatrix< number > &B, const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, std::vector< Vector< number >> &eigenvectors, const types::blas_int itype=1)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
std::vector< MatrixPairType > mass_and_derivative_matrices
void apply_inverse(const unsigned int index, const ArrayView< Number > &dst_in, const ArrayView< const Number > &src_in, AlignedVector< Number > &tmp_array) const
std::pair< Table< 2, Number >, Table< 2, Number > > MatrixPairType
std::map< MatrixPairType, unsigned int, internal::TensorProductMatrixSymmetricSum::MatrixPairComparator< Number > > cache
void reserve(const unsigned int size)
void insert(const unsigned int index, const T &Ms, const T &Ks)
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
void apply_inverse(const ArrayView< Number > &dst, const ArrayView< const Number > &src, AlignedVector< Number > &tmp) const
TensorProductMatrixSymmetricSum(const T &mass_matrix, const T &derivative_matrix)
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
@ matrix
Contents is actually a matrix.
@ eigenvalues
Eigenvalue vector is filled.
static const char A
static const char T
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: l2.h:58
std::enable_if_t< std::is_fundamental< T >::value, 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)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
static const unsigned int invalid_unsigned_int
Definition: types.h:213
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:133
AdditionalData(const bool compress_matrices=true)
std::pair< Table< 2, Number >, Table< 2, Number > > MatrixPairType
bool operator()(const MatrixPairType &left, const MatrixPairType &right) const
static constexpr std::size_t width
static Number & get(Number &value, unsigned int c)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
void vectorized_load_and_transpose(const unsigned int n_entries, const Number *in, const unsigned int *offsets, VectorizedArray< Number, width > *out)
void vectorized_transpose_and_store(const bool add_into, const unsigned int n_entries, const VectorizedArray< Number, width > *in, const unsigned int *offsets, Number *out)