Reference documentation for deal.II version GIT 5983d193e2 2023-05-27 16:50: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  {
287  static constexpr std::size_t width = VectorizedArrayTrait::width();
288 
290  std::pair<std::bitset<width>,
291  std::pair<Table<2, Number>, Table<2, Number>>>;
292 
294  : eps(std::sqrt(std::numeric_limits<ScalarNumber>::epsilon()))
295  {}
296 
297  bool
298  operator()(const MatrixPairType &left, const MatrixPairType &right) const
299  {
300  const auto &M_0 = left.second.first;
301  const auto &K_0 = left.second.second;
302  const auto &M_1 = right.second.first;
303  const auto &K_1 = right.second.second;
304 
305  std::bitset<width> mask;
306 
307  for (unsigned int v = 0; v < width; ++v)
308  mask[v] = left.first[v] && right.first[v];
309 
310  const FloatingPointComparator<Number> comparator(
311  eps, false /*use relative tolerance*/, mask);
312 
313  if (comparator(M_0, M_1))
314  return true;
315  else if (comparator(M_1, M_0))
316  return false;
317  else if (comparator(K_0, K_1))
318  return true;
319  else
320  return false;
321  }
322 
323  private:
325  };
326  } // namespace TensorProductMatrixSymmetricSum
327 } // namespace internal
328 
329 
330 
354 template <int dim, typename Number, int n_rows_1d = -1>
356 {
357  using MatrixPairType = std::pair<Table<2, Number>, Table<2, Number>>;
358 
359  using MatrixPairTypeWithMask = std::pair<
360  std::bitset<::internal::VectorizedArrayTrait<Number>::width()>,
362 
363 public:
368  {
373  const bool precompute_inverse_diagonal = true);
374 
379 
384  };
385 
390  const AdditionalData &additional_data = AdditionalData());
391 
396  void
397  reserve(const unsigned int size);
398 
404  template <typename T>
405  void
406  insert(const unsigned int index, const T &Ms, const T &Ks);
407 
412  void
414 
418  void
419  apply_inverse(const unsigned int index,
420  const ArrayView<Number> & dst_in,
421  const ArrayView<const Number> &src_in,
422  AlignedVector<Number> & tmp_array) const;
423 
427  std::size_t
429 
438  std::size_t
439  storage_size() const;
440 
441 private:
445  const bool compress_matrices;
446 
451 
456  std::vector<MatrixPairType> mass_and_derivative_matrices;
457 
462  std::map<
464  unsigned int,
467 
475  std::vector<unsigned int> indices;
476 
481 
486 
491 
496 
501 
505  std::vector<unsigned int> vector_ptr;
506 
510  std::vector<unsigned int> matrix_ptr;
511 
515  std::vector<unsigned int> vector_n_rows_1d;
516 };
517 
518 
519 /*----------------------- Inline functions ----------------------------------*/
520 
521 #ifndef DOXYGEN
522 
523 namespace internal
524 {
526  {
535  template <typename Number>
536  void
537  spectral_assembly(const Number * mass_matrix,
538  const Number * derivative_matrix,
539  const unsigned int n_rows,
540  const unsigned int n_cols,
541  Number * eigenvalues,
542  Number * eigenvectors)
543  {
544  Assert(n_rows == n_cols, ExcNotImplemented());
545 
546  std::vector<bool> constrained_dofs(n_rows, false);
547 
548  for (unsigned int i = 0; i < n_rows; ++i)
549  {
550  if (mass_matrix[i + i * n_rows] == 0.0)
551  {
552  Assert(derivative_matrix[i + i * n_rows] == 0.0,
553  ExcInternalError());
554 
555  for (unsigned int j = 0; j < n_rows; ++j)
556  {
557  Assert(derivative_matrix[i + j * n_rows] == 0,
558  ExcInternalError());
559  Assert(derivative_matrix[j + i * n_rows] == 0,
560  ExcInternalError());
561  }
562 
563  constrained_dofs[i] = true;
564  }
565  }
566 
567  const auto transpose_fill_nm = [&constrained_dofs](Number * out,
568  const Number * in,
569  const unsigned int n,
570  const unsigned int m) {
571  for (unsigned int mm = 0, c = 0; mm < m; ++mm)
572  for (unsigned int nn = 0; nn < n; ++nn, ++c)
573  out[mm + nn * m] =
574  (mm == nn && constrained_dofs[mm]) ? Number(1.0) : in[c];
575  };
576 
577  std::vector<::Vector<Number>> eigenvecs(n_rows);
578  LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
579  LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);
580 
581  transpose_fill_nm(&(mass_copy(0, 0)), mass_matrix, n_rows, n_cols);
582  transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
583 
584  deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy,
585  eigenvecs);
586  AssertDimension(eigenvecs.size(), n_rows);
587  for (unsigned int i = 0, c = 0; i < n_rows; ++i)
588  for (unsigned int j = 0; j < n_cols; ++j, ++c)
589  if (constrained_dofs[i] == false)
590  eigenvectors[c] = eigenvecs[j][i];
591 
592  for (unsigned int i = 0; i < n_rows; ++i, ++eigenvalues)
593  *eigenvalues = deriv_copy.eigenvalue(i).real();
594  }
595 
596 
597 
598  template <std::size_t dim, typename Number>
599  inline void
600  setup(const std::array<Table<2, Number>, dim> &mass_matrix,
601  const std::array<Table<2, Number>, dim> &derivative_matrix,
602  std::array<Table<2, Number>, dim> & eigenvectors,
603  std::array<AlignedVector<Number>, dim> & eigenvalues)
604  {
605  const unsigned int n_rows_1d = mass_matrix[0].n_cols();
606  (void)n_rows_1d;
607 
608  for (unsigned int dir = 0; dir < dim; ++dir)
609  {
610  AssertDimension(n_rows_1d, mass_matrix[dir].n_cols());
611  AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
612  AssertDimension(mass_matrix[dir].n_rows(),
613  derivative_matrix[dir].n_rows());
614  AssertDimension(mass_matrix[dir].n_rows(),
615  derivative_matrix[dir].n_cols());
616 
617  eigenvectors[dir].reinit(mass_matrix[dir].n_cols(),
618  mass_matrix[dir].n_rows());
619  eigenvalues[dir].resize(mass_matrix[dir].n_cols());
620  internal::TensorProductMatrixSymmetricSum::spectral_assembly<Number>(
621  &(mass_matrix[dir](0, 0)),
622  &(derivative_matrix[dir](0, 0)),
623  mass_matrix[dir].n_rows(),
624  mass_matrix[dir].n_cols(),
625  eigenvalues[dir].begin(),
626  &(eigenvectors[dir](0, 0)));
627  }
628  }
629 
630 
631 
632  template <std::size_t dim, typename Number, std::size_t n_lanes>
633  inline void
634  setup(
635  const std::array<Table<2, VectorizedArray<Number, n_lanes>>, dim>
636  &mass_matrix,
637  const std::array<Table<2, VectorizedArray<Number, n_lanes>>, dim>
638  &derivative_matrix,
641  &eigenvalues)
642  {
643  const unsigned int n_rows_1d = mass_matrix[0].n_cols();
644  constexpr unsigned int macro_size =
646  const std::size_t nm_flat_size_max = n_rows_1d * n_rows_1d * macro_size;
647  const std::size_t n_flat_size_max = n_rows_1d * macro_size;
648 
649  std::vector<Number> mass_matrix_flat;
650  std::vector<Number> deriv_matrix_flat;
651  std::vector<Number> eigenvalues_flat;
652  std::vector<Number> eigenvectors_flat;
653  mass_matrix_flat.resize(nm_flat_size_max);
654  deriv_matrix_flat.resize(nm_flat_size_max);
655  eigenvalues_flat.resize(n_flat_size_max);
656  eigenvectors_flat.resize(nm_flat_size_max);
657  std::array<unsigned int, macro_size> offsets_nm;
658  std::array<unsigned int, macro_size> offsets_n;
659  for (unsigned int dir = 0; dir < dim; ++dir)
660  {
661  AssertDimension(n_rows_1d, mass_matrix[dir].n_cols());
662  AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
663  AssertDimension(mass_matrix[dir].n_rows(),
664  derivative_matrix[dir].n_rows());
665  AssertDimension(mass_matrix[dir].n_rows(),
666  derivative_matrix[dir].n_cols());
667 
668  const unsigned int n_rows = mass_matrix[dir].n_rows();
669  const unsigned int n_cols = mass_matrix[dir].n_cols();
670  const unsigned int nm = n_rows * n_cols;
671  for (unsigned int vv = 0; vv < macro_size; ++vv)
672  offsets_nm[vv] = nm * vv;
673 
675  nm,
676  &(mass_matrix[dir](0, 0)),
677  offsets_nm.cbegin(),
678  mass_matrix_flat.data());
680  nm,
681  &(derivative_matrix[dir](0, 0)),
682  offsets_nm.cbegin(),
683  deriv_matrix_flat.data());
684 
685  const Number *mass_cbegin = mass_matrix_flat.data();
686  const Number *deriv_cbegin = deriv_matrix_flat.data();
687  Number * eigenvec_begin = eigenvectors_flat.data();
688  Number * eigenval_begin = eigenvalues_flat.data();
689  for (unsigned int lane = 0; lane < macro_size; ++lane)
690  internal::TensorProductMatrixSymmetricSum::spectral_assembly<
691  Number>(mass_cbegin + nm * lane,
692  deriv_cbegin + nm * lane,
693  n_rows,
694  n_cols,
695  eigenval_begin + n_rows * lane,
696  eigenvec_begin + nm * lane);
697 
698  eigenvalues[dir].resize(n_rows);
699  eigenvectors[dir].reinit(n_rows, n_cols);
700  for (unsigned int vv = 0; vv < macro_size; ++vv)
701  offsets_n[vv] = n_rows * vv;
703  eigenvalues_flat.data(),
704  offsets_n.cbegin(),
705  eigenvalues[dir].begin());
707  eigenvectors_flat.data(),
708  offsets_nm.cbegin(),
709  &(eigenvectors[dir](0, 0)));
710  }
711  }
712 
713 
714 
715  template <std::size_t dim, typename Number>
716  inline std::array<Table<2, Number>, dim>
717  convert(const std::array<Table<2, Number>, dim> &mass_matrix)
718  {
719  return mass_matrix;
720  }
721 
722 
723 
724  template <std::size_t dim, typename Number>
725  inline std::array<Table<2, Number>, dim>
726  convert(const std::array<FullMatrix<Number>, dim> &mass_matrix)
727  {
728  std::array<Table<2, Number>, dim> mass_copy;
729 
730  std::transform(mass_matrix.cbegin(),
731  mass_matrix.cend(),
732  mass_copy.begin(),
733  [](const FullMatrix<Number> &m) -> Table<2, Number> {
734  return m;
735  });
736 
737  return mass_copy;
738  }
739 
740 
741 
742  template <std::size_t dim, typename Number>
743  inline std::array<Table<2, Number>, dim>
744  convert(const Table<2, Number> &matrix)
745  {
746  std::array<Table<2, Number>, dim> matrices;
747 
748  std::fill(matrices.begin(), matrices.end(), matrix);
749 
750  return matrices;
751  }
752 
753 
754 
755  template <int n_rows_1d_templated, std::size_t dim, typename Number>
756  void
757  vmult(Number * dst,
758  const Number * src,
759  AlignedVector<Number> & tmp,
760  const unsigned int n_rows_1d_non_templated,
761  const std::array<const Number *, dim> &mass_matrix,
762  const std::array<const Number *, dim> &derivative_matrix)
763  {
764  const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
765  n_rows_1d_non_templated :
766  n_rows_1d_templated;
767  const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
768 
769  tmp.resize_fast(n * 2);
770  Number *t = tmp.begin();
771 
773  dim,
774  n_rows_1d_templated,
775  n_rows_1d_templated,
776  Number>
777  eval({}, {}, {}, n_rows_1d, n_rows_1d);
778 
779  if (dim == 1)
780  {
781  const Number *A = derivative_matrix[0];
782  eval.template apply<0, false, false>(A, src, dst);
783  }
784 
785  else if (dim == 2)
786  {
787  const Number *A0 = derivative_matrix[0];
788  const Number *M0 = mass_matrix[0];
789  const Number *A1 = derivative_matrix[1];
790  const Number *M1 = mass_matrix[1];
791  eval.template apply<0, false, false>(M0, src, t);
792  eval.template apply<1, false, false>(A1, t, dst);
793  eval.template apply<0, false, false>(A0, src, t);
794  eval.template apply<1, false, true>(M1, t, dst);
795  }
796 
797  else if (dim == 3)
798  {
799  const Number *A0 = derivative_matrix[0];
800  const Number *M0 = mass_matrix[0];
801  const Number *A1 = derivative_matrix[1];
802  const Number *M1 = mass_matrix[1];
803  const Number *A2 = derivative_matrix[2];
804  const Number *M2 = mass_matrix[2];
805  eval.template apply<0, false, false>(M0, src, t + n);
806  eval.template apply<1, false, false>(M1, t + n, t);
807  eval.template apply<2, false, false>(A2, t, dst);
808  eval.template apply<1, false, false>(A1, t + n, t);
809  eval.template apply<0, false, false>(A0, src, t + n);
810  eval.template apply<1, false, true>(M1, t + n, t);
811  eval.template apply<2, false, true>(M2, t, dst);
812  }
813 
814  else
815  AssertThrow(false, ExcNotImplemented());
816  }
817 
818 
819 
820  template <int n_rows_1d_templated, std::size_t dim, typename Number>
821  void
822  apply_inverse(Number * dst,
823  const Number * src,
825  const unsigned int n_rows_1d_non_templated,
826  const std::array<const Number *, dim> &eigenvectors,
827  const std::array<const Number *, dim> &eigenvalues,
828  const Number *inverted_eigenvalues = nullptr)
829  {
830  const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
831  n_rows_1d_non_templated :
832  n_rows_1d_templated;
833  const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
834 
835  tmp.resize_fast(n);
836  Number *t = tmp.begin();
837 
839  dim,
840  n_rows_1d_templated,
841  n_rows_1d_templated,
842  Number>
843  eval({}, {}, {}, n_rows_1d, n_rows_1d);
844 
845  // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
846  // --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
847  // --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
848  // while the eigenvectors are stored column-wise in S, i.e.
849  // rows correspond to dofs whereas columns to eigenvalue indices!
850  if (dim == 1)
851  {
852  const Number *S = eigenvectors[0];
853  eval.template apply<0, true, false>(S, src, t);
854 
855  for (unsigned int i = 0; i < n_rows_1d; ++i)
856  if (inverted_eigenvalues)
857  t[i] *= inverted_eigenvalues[i];
858  else
859  t[i] /= eigenvalues[0][i];
860 
861  eval.template apply<0, false, false>(S, t, dst);
862  }
863 
864  else if (dim == 2)
865  {
866  const Number *S0 = eigenvectors[0];
867  const Number *S1 = eigenvectors[1];
868  eval.template apply<0, true, false>(S0, src, t);
869  eval.template apply<1, true, false>(S1, t, dst);
870 
871  for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
872  for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
873  if (inverted_eigenvalues)
874  dst[c] *= inverted_eigenvalues[c];
875  else
876  dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
877 
878  eval.template apply<0, false, false>(S0, dst, t);
879  eval.template apply<1, false, false>(S1, t, dst);
880  }
881 
882  else if (dim == 3)
883  {
884  const Number *S0 = eigenvectors[0];
885  const Number *S1 = eigenvectors[1];
886  const Number *S2 = eigenvectors[2];
887  eval.template apply<0, true, false>(S0, src, t);
888  eval.template apply<1, true, false>(S1, t, dst);
889  eval.template apply<2, true, false>(S2, dst, t);
890 
891  for (unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
892  for (unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
893  for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
894  if (inverted_eigenvalues)
895  t[c] *= inverted_eigenvalues[c];
896  else
897  t[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] +
898  eigenvalues[0][i0]);
899 
900  eval.template apply<0, false, false>(S0, t, dst);
901  eval.template apply<1, false, false>(S1, dst, t);
902  eval.template apply<2, false, false>(S2, t, dst);
903  }
904 
905  else
906  Assert(false, ExcNotImplemented());
907  }
908 
909 
910 
911  template <int n_rows_1d_templated, std::size_t dim, typename Number>
912  void
913  select_vmult(Number * dst,
914  const Number * src,
915  AlignedVector<Number> & tmp,
916  const unsigned int n_rows_1d,
917  const std::array<const Number *, dim> &mass_matrix,
918  const std::array<const Number *, dim> &derivative_matrix);
919 
920 
921 
922  template <int n_rows_1d_templated, std::size_t dim, typename Number>
923  void
924  select_apply_inverse(Number * dst,
925  const Number * src,
926  AlignedVector<Number> & tmp,
927  const unsigned int n_rows_1d,
928  const std::array<const Number *, dim> &eigenvectors,
929  const std::array<const Number *, dim> &eigenvalues,
930  const Number *inverted_eigenvalues = nullptr);
931  } // namespace TensorProductMatrixSymmetricSum
932 } // namespace internal
933 
934 
935 template <int dim, typename Number, int n_rows_1d>
936 inline unsigned int
938 {
939  unsigned int m = mass_matrix[0].n_rows();
940  for (unsigned int d = 1; d < dim; ++d)
941  m *= mass_matrix[d].n_rows();
942  return m;
943 }
944 
945 
946 
947 template <int dim, typename Number, int n_rows_1d>
948 inline unsigned int
950 {
951  unsigned int n = mass_matrix[0].n_cols();
952  for (unsigned int d = 1; d < dim; ++d)
953  n *= mass_matrix[d].n_cols();
954  return n;
955 }
956 
957 
958 
959 template <int dim, typename Number, int n_rows_1d>
960 inline void
962  const ArrayView<Number> & dst_view,
963  const ArrayView<const Number> &src_view) const
964 {
965  std::lock_guard<std::mutex> lock(this->mutex);
966  this->vmult(dst_view, src_view, this->tmp_array);
967 }
968 
969 
970 
971 template <int dim, typename Number, int n_rows_1d>
972 inline void
974  const ArrayView<Number> & dst_view,
975  const ArrayView<const Number> &src_view,
976  AlignedVector<Number> & tmp_array) const
977 {
978  AssertDimension(dst_view.size(), this->m());
979  AssertDimension(src_view.size(), this->n());
980 
981  Number * dst = dst_view.begin();
982  const Number *src = src_view.begin();
983 
984  std::array<const Number *, dim> mass_matrix, derivative_matrix;
985 
986  for (unsigned int d = 0; d < dim; ++d)
987  {
988  mass_matrix[d] = &this->mass_matrix[d](0, 0);
989  derivative_matrix[d] = &this->derivative_matrix[d](0, 0);
990  }
991 
992  const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
993 
994  if (n_rows_1d != -1)
995  internal::TensorProductMatrixSymmetricSum::vmult<
996  n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
997  src,
998  tmp_array,
999  n_rows_1d_non_templated,
1000  mass_matrix,
1001  derivative_matrix);
1002  else
1003  internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
1004  dst,
1005  src,
1006  tmp_array,
1007  n_rows_1d_non_templated,
1008  mass_matrix,
1009  derivative_matrix);
1010 }
1011 
1012 
1013 
1014 template <int dim, typename Number, int n_rows_1d>
1015 inline void
1017  const ArrayView<Number> & dst_view,
1018  const ArrayView<const Number> &src_view) const
1019 {
1020  std::lock_guard<std::mutex> lock(this->mutex);
1021  this->apply_inverse(dst_view, src_view, this->tmp_array);
1022 }
1023 
1024 
1025 
1026 template <int dim, typename Number, int n_rows_1d>
1027 inline void
1029  const ArrayView<Number> & dst_view,
1030  const ArrayView<const Number> &src_view,
1031  AlignedVector<Number> & tmp_array) const
1032 {
1033  AssertDimension(dst_view.size(), this->n());
1034  AssertDimension(src_view.size(), this->m());
1035 
1036  Number * dst = dst_view.begin();
1037  const Number *src = src_view.begin();
1038 
1039  std::array<const Number *, dim> eigenvectors, eigenvalues;
1040 
1041  for (unsigned int d = 0; d < dim; ++d)
1042  {
1043  eigenvectors[d] = &this->eigenvectors[d](0, 0);
1044  eigenvalues[d] = this->eigenvalues[d].data();
1045  }
1046 
1047  const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
1048 
1049  if (n_rows_1d != -1)
1050  internal::TensorProductMatrixSymmetricSum::apply_inverse<
1051  n_rows_1d == -1 ? 0 : n_rows_1d>(
1052  dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1053  else
1054  internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1055  dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1056 }
1057 
1058 
1059 
1060 template <int dim, typename Number, int n_rows_1d>
1061 std::size_t
1063  const
1064 {
1066  MemoryConsumption::memory_consumption(derivative_matrix) +
1070 }
1071 
1072 
1073 
1074 template <int dim, typename Number, int n_rows_1d>
1075 template <typename T>
1078  const T &derivative_matrix)
1079 {
1080  reinit(mass_matrix, derivative_matrix);
1081 }
1082 
1083 
1084 
1085 template <int dim, typename Number, int n_rows_1d>
1086 template <typename T>
1087 inline void
1089  const T &mass_matrix,
1090  const T &derivative_matrix)
1091 {
1092  this->mass_matrix =
1093  internal::TensorProductMatrixSymmetricSum::convert<dim>(mass_matrix);
1094  this->derivative_matrix =
1095  internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1096 
1097  internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1098  this->derivative_matrix,
1099  this->eigenvectors,
1100  this->eigenvalues);
1101 }
1102 
1103 
1104 
1105 template <int dim, typename Number, int n_rows_1d>
1107  AdditionalData::AdditionalData(const bool compress_matrices,
1108  const bool precompute_inverse_diagonal)
1109  : compress_matrices(compress_matrices)
1110  , precompute_inverse_diagonal(precompute_inverse_diagonal)
1111 {}
1112 
1113 
1114 
1115 template <int dim, typename Number, int n_rows_1d>
1118  const AdditionalData &additional_data)
1119  : compress_matrices(additional_data.compress_matrices)
1120  , precompute_inverse_diagonal(additional_data.precompute_inverse_diagonal)
1121 {}
1122 
1123 
1124 
1125 template <int dim, typename Number, int n_rows_1d>
1126 void
1128  const unsigned int size)
1129 {
1130  if (compress_matrices == false)
1131  mass_and_derivative_matrices.resize(size * dim);
1132  else
1133  indices.assign(size * dim, numbers::invalid_unsigned_int);
1134 }
1135 
1136 
1137 
1138 template <int dim, typename Number, int n_rows_1d>
1139 template <typename T>
1140 void
1142  const unsigned int index,
1143  const T & Ms_in,
1144  const T & Ks_in)
1145 {
1146  const auto Ms =
1147  internal::TensorProductMatrixSymmetricSum::convert<dim>(Ms_in);
1148  const auto Ks =
1149  internal::TensorProductMatrixSymmetricSum::convert<dim>(Ks_in);
1150 
1151  for (unsigned int d = 0; d < dim; ++d)
1152  {
1153  if (compress_matrices == false)
1154  {
1155  const MatrixPairType matrix(Ms[d], Ks[d]);
1156  mass_and_derivative_matrices[index * dim + d] = matrix;
1157  }
1158  else
1159  {
1160  using VectorizedArrayTrait =
1162 
1163  std::bitset<VectorizedArrayTrait::width()> mask;
1164 
1165  for (unsigned int v = 0; v < VectorizedArrayTrait::width(); ++v)
1166  {
1167  typename VectorizedArrayTrait::value_type a = 0.0;
1168 
1169  for (unsigned int i = 0; i < Ms[d].size(0); ++i)
1170  for (unsigned int j = 0; j < Ms[d].size(1); ++j)
1171  {
1172  a += std::abs(VectorizedArrayTrait::get(Ms[d][i][j], v));
1173  a += std::abs(VectorizedArrayTrait::get(Ks[d][i][j], v));
1174  }
1175 
1176  mask[v] = (a != 0.0);
1177  }
1178 
1179  const MatrixPairTypeWithMask matrix{mask, {Ms[d], Ks[d]}};
1180 
1181  const auto ptr = cache.find(matrix);
1182 
1183  if (ptr != cache.end())
1184  {
1185  const auto ptr_index = ptr->second;
1186  indices[index * dim + d] = ptr_index;
1187 
1188  if ([&]() {
1189  for (unsigned int v = 0; v < VectorizedArrayTrait::width();
1190  ++v)
1191  if ((mask[v] == true) && (ptr->first.first[v] == false))
1192  return false;
1193 
1194  return true;
1195  }())
1196  {
1197  // nothing to do
1198  }
1199  else
1200  {
1201  auto mask_new = ptr->first.first;
1202  auto Ms_new = ptr->first.second.first;
1203  auto Ks_new = ptr->first.second.second;
1204 
1205  for (unsigned int v = 0; v < VectorizedArrayTrait::width();
1206  ++v)
1207  if (mask_new[v] == false && mask[v] == true)
1208  {
1209  mask_new[v] = true;
1210 
1211  for (unsigned int i = 0; i < Ms_new.size(0); ++i)
1212  for (unsigned int j = 0; j < Ms_new.size(1); ++j)
1213  {
1214  VectorizedArrayTrait::get(Ms_new[i][j], v) =
1215  VectorizedArrayTrait::get(Ms[d][i][j], v);
1216  VectorizedArrayTrait::get(Ks_new[i][j], v) =
1217  VectorizedArrayTrait::get(Ks[d][i][j], v);
1218  }
1219  }
1220 
1221  cache.erase(ptr);
1222 
1223  const MatrixPairTypeWithMask entry_new{mask_new,
1224  {Ms_new, Ks_new}};
1225 
1226  const auto ptr_ = cache.find(entry_new);
1227  AssertThrow(ptr_ == cache.end(), ExcNotImplemented());
1228 
1229  cache[entry_new] = ptr_index;
1230  }
1231  }
1232  else
1233  {
1234  const auto size = cache.size();
1235  indices[index * dim + d] = size;
1236  cache[matrix] = size;
1237  }
1238  }
1239  }
1240 }
1241 
1242 
1243 
1244 template <int dim, typename Number, int n_rows_1d>
1245 void
1247 {
1248  const auto store = [&](const unsigned int index,
1249  const MatrixPairType &M_and_K) {
1250  std::array<Table<2, Number>, 1> mass_matrix;
1251  mass_matrix[0] = M_and_K.first;
1252 
1253  std::array<Table<2, Number>, 1> derivative_matrix;
1254  derivative_matrix[0] = M_and_K.second;
1255 
1256  std::array<Table<2, Number>, 1> eigenvectors;
1257  std::array<AlignedVector<Number>, 1> eigenvalues;
1258 
1259  internal::TensorProductMatrixSymmetricSum::setup(mass_matrix,
1260  derivative_matrix,
1261  eigenvectors,
1262  eigenvalues);
1263 
1264  for (unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1265  i < mass_matrix[0].n_rows();
1266  ++i, ++v)
1267  {
1268  for (unsigned int j = 0; j < mass_matrix[0].n_cols(); ++j, ++m)
1269  {
1270  this->mass_matrices[m] = mass_matrix[0][i][j];
1271  this->derivative_matrices[m] = derivative_matrix[0][i][j];
1272  this->eigenvectors[m] = eigenvectors[0][i][j];
1273  }
1274 
1275  this->eigenvalues[v] = eigenvalues[0][i];
1276  }
1277  };
1278 
1279  if (compress_matrices == false)
1280  {
1281  // case 1) no compression requested
1282 
1283  AssertDimension(cache.size(), 0);
1284  AssertDimension(indices.size(), 0);
1285 
1286  this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1287  this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1288 
1289  for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1290  {
1291  const auto &M = mass_and_derivative_matrices[i].first;
1292 
1293  this->vector_ptr[i + 1] = M.n_rows();
1294  this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1295  }
1296 
1297  for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1298  {
1299  this->vector_ptr[i + 1] += this->vector_ptr[i];
1300  this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1301  }
1302 
1303  this->mass_matrices.resize_fast(matrix_ptr.back());
1304  this->derivative_matrices.resize_fast(matrix_ptr.back());
1305  this->eigenvectors.resize_fast(matrix_ptr.back());
1306  this->eigenvalues.resize_fast(vector_ptr.back());
1307 
1308  for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1309  store(i, mass_and_derivative_matrices[i]);
1310 
1311  mass_and_derivative_matrices.clear();
1312  }
1313  else if (cache.size() == indices.size())
1314  {
1315  // case 2) compression requested but none possible
1316 
1317  this->vector_ptr.resize(cache.size() + 1);
1318  this->matrix_ptr.resize(cache.size() + 1);
1319 
1320  std::map<unsigned int, MatrixPairType> inverted_cache;
1321 
1322  for (const auto &i : cache)
1323  inverted_cache[i.second] = i.first.second;
1324 
1325  for (unsigned int i = 0; i < indices.size(); ++i)
1326  {
1327  const auto &M = inverted_cache[indices[i]].first;
1328 
1329  this->vector_ptr[i + 1] = M.n_rows();
1330  this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1331  }
1332 
1333  for (unsigned int i = 0; i < cache.size(); ++i)
1334  {
1335  this->vector_ptr[i + 1] += this->vector_ptr[i];
1336  this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1337  }
1338 
1339  this->mass_matrices.resize_fast(matrix_ptr.back());
1340  this->derivative_matrices.resize_fast(matrix_ptr.back());
1341  this->eigenvectors.resize_fast(matrix_ptr.back());
1342  this->eigenvalues.resize_fast(vector_ptr.back());
1343 
1344  for (unsigned int i = 0; i < indices.size(); ++i)
1345  store(i, inverted_cache[indices[i]]);
1346 
1347  indices.clear();
1348  cache.clear();
1349  }
1350  else
1351  {
1352  // case 3) compress
1353 
1354  this->vector_ptr.resize(cache.size() + 1);
1355  this->matrix_ptr.resize(cache.size() + 1);
1356 
1357  for (const auto &i : cache)
1358  {
1359  const auto &M = i.first.second.first;
1360 
1361  this->vector_ptr[i.second + 1] = M.n_rows();
1362  this->matrix_ptr[i.second + 1] = M.n_rows() * M.n_cols();
1363  }
1364 
1365  for (unsigned int i = 0; i < cache.size(); ++i)
1366  {
1367  this->vector_ptr[i + 1] += this->vector_ptr[i];
1368  this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1369  }
1370 
1371  this->mass_matrices.resize_fast(matrix_ptr.back());
1372  this->derivative_matrices.resize_fast(matrix_ptr.back());
1373  this->eigenvectors.resize_fast(matrix_ptr.back());
1374  this->eigenvalues.resize_fast(vector_ptr.back());
1375 
1376  for (const auto &i : cache)
1377  store(i.second, i.first.second);
1378 
1379  cache.clear();
1380  }
1381 
1382  if (precompute_inverse_diagonal)
1383  {
1384  if (dim == 1)
1385  {
1386  // 1D case: simply invert 1D eigenvalues
1387  for (unsigned int i = 0; i < this->eigenvalues.size(); ++i)
1388  this->eigenvalues[i] = Number(1.0) / this->eigenvalues[i];
1389  std::swap(this->inverted_eigenvalues, eigenvalues);
1390  }
1391  else
1392  {
1393  // 2D and 3D case: we have 2 or 3 1d eigenvalues so that we
1394  // need to combine these
1395 
1396  // step 1) if eigenvalues/eigenvectors are compressed, we
1397  // need to compress the diagonal (the combination of ev
1398  // indices) as well. This is an optional step.
1399  std::vector<unsigned int> indices_ev;
1400 
1401  if (indices.size() > 0)
1402  {
1403  // 1a) create cache (ev indics -> diag index)
1404  const unsigned int n_cells = indices.size() / dim;
1405  std::map<std::array<unsigned int, dim>, unsigned int> cache_ev;
1406  std::vector<unsigned int> cache_ev_idx(n_cells);
1407 
1408  for (unsigned int i = 0, c = 0; i < n_cells; ++i)
1409  {
1410  std::array<unsigned int, dim> id;
1411 
1412  for (unsigned int d = 0; d < dim; ++d, ++c)
1413  id[d] = indices[c];
1414 
1415  const auto id_ptr = cache_ev.find(id);
1416 
1417  if (id_ptr == cache_ev.end())
1418  {
1419  const auto size = cache_ev.size();
1420  cache_ev_idx[i] = size;
1421  cache_ev[id] = size;
1422  }
1423  else
1424  {
1425  cache_ev_idx[i] = id_ptr->second;
1426  }
1427  }
1428 
1429  // 1b) store diagonal indices for each cell
1430  std::vector<unsigned int> new_indices;
1431  new_indices.reserve(indices.size() / dim * (dim + 1));
1432 
1433  for (unsigned int i = 0, c = 0; i < n_cells; ++i)
1434  {
1435  for (unsigned int d = 0; d < dim; ++d, ++c)
1436  new_indices.push_back(indices[c]);
1437  new_indices.push_back(cache_ev_idx[i]);
1438  }
1439 
1440  // 1c) transpose cache (diag index -> ev indices)
1441  indices_ev.resize(cache_ev.size() * dim);
1442  for (const auto &entry : cache_ev)
1443  for (unsigned int d = 0; d < dim; ++d)
1444  indices_ev[entry.second * dim + d] = entry.first[d];
1445 
1446  std::swap(this->indices, new_indices);
1447  }
1448 
1449  // step 2) allocate memory and set pointers
1450  const unsigned int n_diag =
1451  ((indices_ev.size() > 0) ? indices_ev.size() :
1452  (matrix_ptr.size() - 1)) /
1453  dim;
1454 
1455  std::vector<unsigned int> new_vector_ptr(n_diag + 1, 0);
1456  std::vector<unsigned int> new_vector_n_rows_1d(n_diag, 0);
1457 
1458  for (unsigned int i = 0; i < n_diag; ++i)
1459  {
1460  const unsigned int c = (indices_ev.size() > 0) ?
1461  indices_ev[dim * i + 0] :
1462  (dim * i + 0);
1463 
1464  const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c];
1465 
1466  new_vector_n_rows_1d[i] = n_rows;
1467  new_vector_ptr[i + 1] = Utilities::pow(n_rows, dim);
1468  }
1469 
1470  for (unsigned int i = 0; i < n_diag; ++i)
1471  new_vector_ptr[i + 1] += new_vector_ptr[i];
1472 
1473  this->inverted_eigenvalues.resize(new_vector_ptr.back());
1474 
1475  // step 3) loop over all unique diagonal entries and invert
1476  for (unsigned int i = 0; i < n_diag; ++i)
1477  {
1478  std::array<Number *, dim> evs;
1479 
1480  for (unsigned int d = 0; d < dim; ++d)
1481  evs[d] =
1482  &this
1483  ->eigenvalues[this->vector_ptr[(indices_ev.size() > 0) ?
1484  indices_ev[dim * i + d] :
1485  (dim * i + d)]];
1486 
1487  const unsigned int mm = new_vector_n_rows_1d[i];
1488  if (dim == 2)
1489  {
1490  for (unsigned int i1 = 0, c = 0; i1 < mm; ++i1)
1491  for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1492  this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1493  Number(1.0) / (evs[1][i1] + evs[0][i0]);
1494  }
1495  else
1496  {
1497  for (unsigned int i2 = 0, c = 0; i2 < mm; ++i2)
1498  for (unsigned int i1 = 0; i1 < mm; ++i1)
1499  for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1500  this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1501  Number(1.0) / (evs[2][i2] + evs[1][i1] + evs[0][i0]);
1502  }
1503  }
1504 
1505  // step 4) clean up
1506  std::swap(this->vector_ptr, new_vector_ptr);
1507  std::swap(this->vector_n_rows_1d, new_vector_n_rows_1d);
1508  }
1509 
1510  this->eigenvalues.clear();
1511  }
1512 }
1513 
1514 
1515 
1516 template <int dim, typename Number, int n_rows_1d>
1517 void
1519  apply_inverse(const unsigned int index,
1520  const ArrayView<Number> & dst_in,
1521  const ArrayView<const Number> &src_in,
1522  AlignedVector<Number> & tmp_array) const
1523 {
1524  Number * dst = dst_in.begin();
1525  const Number *src = src_in.begin();
1526 
1527  if (this->eigenvalues.empty() == false)
1528  {
1529  std::array<const Number *, dim> eigenvectors;
1530  std::array<const Number *, dim> eigenvalues;
1531  unsigned int n_rows_1d_non_templated = 0;
1532 
1533  for (unsigned int d = 0; d < dim; ++d)
1534  {
1535  const unsigned int translated_index =
1536  (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d);
1537 
1538  eigenvectors[d] =
1539  this->eigenvectors.data() + matrix_ptr[translated_index];
1540  eigenvalues[d] =
1541  this->eigenvalues.data() + vector_ptr[translated_index];
1542  n_rows_1d_non_templated =
1543  vector_ptr[translated_index + 1] - vector_ptr[translated_index];
1544  }
1545 
1546  if (n_rows_1d != -1)
1547  internal::TensorProductMatrixSymmetricSum::apply_inverse<
1548  n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
1549  src,
1550  tmp_array,
1551  n_rows_1d_non_templated,
1552  eigenvectors,
1553  eigenvalues);
1554  else
1555  internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1556  dst,
1557  src,
1558  tmp_array,
1559  n_rows_1d_non_templated,
1560  eigenvectors,
1561  eigenvalues);
1562  }
1563  else
1564  {
1565  std::array<const Number *, dim> eigenvectors;
1566  const Number * inverted_eigenvalues = nullptr;
1567  unsigned int n_rows_1d_non_templated = 0;
1568 
1569  for (unsigned int d = 0; d < dim; ++d)
1570  {
1571  const unsigned int translated_index =
1572  (indices.size() > 0) ?
1573  indices[((dim == 1) ? 1 : (dim + 1)) * index + d] :
1574  (dim * index + d);
1575 
1576  eigenvectors[d] =
1577  this->eigenvectors.data() + matrix_ptr[translated_index];
1578  }
1579 
1580  {
1581  const unsigned int translated_index =
1582  ((indices.size() > 0) && (dim != 1)) ?
1583  indices[(dim + 1) * index + dim] :
1584  index;
1585 
1586  inverted_eigenvalues =
1587  this->inverted_eigenvalues.data() + vector_ptr[translated_index];
1588  n_rows_1d_non_templated =
1589  (dim == 1) ?
1590  (vector_ptr[translated_index + 1] - vector_ptr[translated_index]) :
1591  vector_n_rows_1d[translated_index];
1592  }
1593 
1594  if (n_rows_1d != -1)
1595  internal::TensorProductMatrixSymmetricSum::apply_inverse<
1596  n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
1597  src,
1598  tmp_array,
1599  n_rows_1d_non_templated,
1600  eigenvectors,
1601  {},
1602  inverted_eigenvalues);
1603  else
1604  internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1605  dst,
1606  src,
1607  tmp_array,
1608  n_rows_1d_non_templated,
1609  eigenvectors,
1610  {},
1611  inverted_eigenvalues);
1612  }
1613 }
1614 
1615 
1616 
1617 template <int dim, typename Number, int n_rows_1d>
1618 std::size_t
1620  memory_consumption() const
1621 {
1622  return MemoryConsumption::memory_consumption(indices) +
1623  MemoryConsumption::memory_consumption(mass_matrices) +
1624  MemoryConsumption::memory_consumption(derivative_matrices) +
1629 }
1630 
1631 
1632 
1633 template <int dim, typename Number, int n_rows_1d>
1634 std::size_t
1636  storage_size() const
1637 {
1638  if (matrix_ptr.size() == 0)
1639  return 0; // if not initialized
1640 
1641  return matrix_ptr.size() - 1;
1642 }
1643 
1644 
1645 
1646 #endif
1647 
1649 
1650 #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
std::pair< std::bitset<::internal::VectorizedArrayTrait< Number >::width()>, MatrixPairType > MatrixPairTypeWithMask
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
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
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:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
@ 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)
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
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)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:447
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13826
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 Function &function, const unsigned int grainsize)
Definition: parallel.h:148
::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< 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)