15#ifndef dealii_tensor_product_matrix_h
16#define dealii_tensor_product_matrix_h
114template <
int dim,
typename Number,
int n_rows_1d = -1>
139 template <
typename T>
160 template <
typename T>
263 template <
typename Number>
272 std::pair<std::bitset<width>,
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;
287 std::bitset<width>
mask;
289 for (
unsigned int v = 0; v <
width; ++v)
290 mask[v] = left.first[v] && right.first[v];
336template <
int dim,
typename Number,
int n_rows_1d = -1>
342 std::bitset<::internal::VectorizedArrayTrait<Number>::width()>,
386 template <
typename T>
516 template <
typename Number>
519 const Number *derivative_matrix,
520 const unsigned int n_rows,
521 const unsigned int n_cols,
527 std::vector<bool> constrained_dofs(n_rows,
false);
529 for (
unsigned int i = 0; i < n_rows; ++i)
531 if (mass_matrix[i + i * n_rows] == 0.0)
533 Assert(derivative_matrix[i + i * n_rows] == 0.0,
536 for (
unsigned int j = 0;
j < n_rows; ++
j)
538 Assert(derivative_matrix[i +
j * n_rows] == 0,
540 Assert(derivative_matrix[
j + i * n_rows] == 0,
544 constrained_dofs[i] =
true;
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)
555 (
mm ==
nn && constrained_dofs[
mm]) ? Number(1.0) : in[c];
558 std::vector<::Vector<Number>>
eigenvecs(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)
573 for (
unsigned int i = 0; i < n_rows; ++i, ++
eigenvalues)
579 template <std::
size_t dim,
typename Number>
586 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
589 for (
unsigned int dir = 0; dir < dim; ++dir)
594 derivative_matrix[dir].n_rows());
596 derivative_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(),
613 template <std::
size_t dim,
typename Number, std::
size_t n_lanes>
624 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
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)
645 derivative_matrix[dir].n_rows());
647 derivative_matrix[dir].n_cols());
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;
658 &(mass_matrix[dir](0, 0)),
664 &(derivative_matrix[dir](0, 0)),
672 for (
unsigned int lane = 0; lane <
macro_size; ++lane)
673 internal::TensorProductMatrixSymmetricSum::spectral_assembly<
700 template <std::
size_t dim,
typename Number>
701 inline std::array<Table<2, Number>, dim>
709 template <std::
size_t dim,
typename Number>
710 inline std::array<Table<2, Number>, dim>
713 std::array<Table<2, Number>, dim>
mass_copy;
715 std::transform(mass_matrix.cbegin(),
727 template <std::
size_t dim,
typename Number>
728 inline std::array<Table<2, Number>, dim>
731 std::array<Table<2, Number>, dim> matrices;
733 std::fill(matrices.begin(), matrices.end(), matrix);
740 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
746 const std::array<const Number *, dim> &mass_matrix,
747 const std::array<const Number *, dim> &derivative_matrix)
752 const unsigned int n = Utilities::fixed_power<dim>(
n_rows_1d);
754 tmp.resize_fast(n * 2);
755 Number *t = tmp.begin();
766 const Number *A = derivative_matrix[0];
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];
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];
805 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
807 apply_inverse(Number *dst,
811 const std::array<const Number *, dim> &
eigenvalues,
812 const Number *inverted_eigenvalues =
nullptr)
835 for (
unsigned int i = 0; i <
n_rows_1d; ++i)
836 if (inverted_eigenvalues)
837 dst[i] *= inverted_eigenvalues[i];
853 if (inverted_eigenvalues)
854 dst[c] *= inverted_eigenvalues[c];
874 if (inverted_eigenvalues)
875 dst[c] *= inverted_eigenvalues[c];
891 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
897 const std::array<const Number *, dim> &mass_matrix,
898 const std::array<const Number *, dim> &derivative_matrix);
902 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
908 const std::array<const Number *, dim> &
eigenvalues,
909 const Number *inverted_eigenvalues =
nullptr);
914template <
int dim,
typename Number,
int n_rows_1d>
919 for (
unsigned int d = 1;
d < dim; ++
d)
920 m *= mass_matrix[d].n_rows();
926template <
int dim,
typename Number,
int n_rows_1d>
931 for (
unsigned int d = 1;
d < dim; ++
d)
932 n *= mass_matrix[d].n_cols();
938template <
int dim,
typename Number,
int n_rows_1d>
944 std::lock_guard<std::mutex> lock(this->mutex);
950template <
int dim,
typename Number,
int n_rows_1d>
961 const Number *src =
src_view.begin();
963 std::array<const Number *, dim>
mass_matrix, derivative_matrix;
965 for (
unsigned int d = 0;
d < dim; ++
d)
968 derivative_matrix[
d] = &this->derivative_matrix[
d](0, 0);
974 internal::TensorProductMatrixSymmetricSum::vmult<
982 internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
993template <
int dim,
typename Number,
int n_rows_1d>
1003 const Number *src =
src_view.begin();
1007 for (
unsigned int d = 0;
d < dim; ++
d)
1016 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1020 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1026template <
int dim,
typename Number,
int n_rows_1d>
1040template <
int dim,
typename Number,
int n_rows_1d>
1041template <
typename T>
1044 const T &derivative_matrix)
1046 reinit(mass_matrix, derivative_matrix);
1051template <
int dim,
typename Number,
int n_rows_1d>
1052template <
typename T>
1055 const T &mass_matrix,
1056 const T &derivative_matrix)
1059 internal::TensorProductMatrixSymmetricSum::convert<dim>(mass_matrix);
1060 this->derivative_matrix =
1061 internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1063 internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1064 this->derivative_matrix,
1071template <
int dim,
typename Number,
int n_rows_1d>
1074 const bool precompute_inverse_diagonal)
1075 : compress_matrices(compress_matrices)
1076 , precompute_inverse_diagonal(precompute_inverse_diagonal)
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)
1091template <
int dim,
typename Number,
int n_rows_1d>
1094 const unsigned int size)
1096 if (compress_matrices ==
false)
1097 mass_and_derivative_matrices.resize(
size * dim);
1104template <
int dim,
typename Number,
int n_rows_1d>
1105template <
typename T>
1108 const unsigned int index,
1113 internal::TensorProductMatrixSymmetricSum::convert<dim>(
Ms_in);
1115 internal::TensorProductMatrixSymmetricSum::convert<dim>(
Ks_in);
1117 for (
unsigned int d = 0;
d < dim; ++
d)
1119 if (compress_matrices ==
false)
1122 mass_and_derivative_matrices[
index * dim +
d] =
matrix;
1126 using VectorizedArrayTrait =
1129 std::bitset<VectorizedArrayTrait::width()>
mask;
1131 for (
unsigned int v = 0; v < VectorizedArrayTrait::width(); ++v)
1133 typename VectorizedArrayTrait::value_type a = 0.0;
1135 for (
unsigned int i = 0; i <
Ms[
d].size(0); ++i)
1136 for (
unsigned int j = 0;
j <
Ms[
d].size(1); ++
j)
1138 a +=
std::abs(VectorizedArrayTrait::get(
Ms[d][i][
j], v));
1139 a +=
std::abs(VectorizedArrayTrait::get(
Ks[d][i][
j], v));
1142 mask[v] = (a != 0.0);
1147 const auto ptr = cache.find(matrix);
1149 if (ptr != cache.end())
1155 for (
unsigned int v = 0; v < VectorizedArrayTrait::width();
1157 if ((mask[v] ==
true) && (ptr->first.first[v] ==
false))
1168 auto Ms_new = ptr->first.second.first;
1169 auto Ks_new = ptr->first.second.second;
1171 for (
unsigned int v = 0; v < VectorizedArrayTrait::width();
1173 if (
mask_new[v] ==
false && mask[v] ==
true)
1177 for (
unsigned int i = 0; i <
Ms_new.size(0); ++i)
1178 for (
unsigned int j = 0;
j <
Ms_new.size(1); ++
j)
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);
1200 const auto size = cache.size();
1210template <
int dim,
typename Number,
int n_rows_1d>
1214 const auto store = [&](
const unsigned int index,
1215 const MatrixPairType &
M_and_K) {
1219 std::array<Table<2, Number>, 1> derivative_matrix;
1220 derivative_matrix[0] =
M_and_K.second;
1225 internal::TensorProductMatrixSymmetricSum::setup(mass_matrix,
1230 for (
unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1237 this->derivative_matrices[m] = derivative_matrix[0][i][
j];
1245 if (compress_matrices ==
false)
1252 this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1253 this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1255 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1257 const auto &
M = mass_and_derivative_matrices[i].first;
1259 this->vector_ptr[i + 1] =
M.n_rows();
1260 this->matrix_ptr[i + 1] =
M.n_rows() *
M.n_cols();
1263 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1265 this->vector_ptr[i + 1] += this->vector_ptr[i];
1266 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
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());
1274 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1275 store(i, mass_and_derivative_matrices[i]);
1277 mass_and_derivative_matrices.
clear();
1279 else if (cache.size() == indices.size())
1283 this->vector_ptr.resize(cache.size() + 1);
1284 this->matrix_ptr.resize(cache.size() + 1);
1288 for (
const auto &i : cache)
1291 for (
unsigned int i = 0; i < indices.size(); ++i)
1295 this->vector_ptr[i + 1] =
M.n_rows();
1296 this->matrix_ptr[i + 1] =
M.n_rows() *
M.n_cols();
1299 for (
unsigned int i = 0; i < cache.size(); ++i)
1301 this->vector_ptr[i + 1] += this->vector_ptr[i];
1302 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
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());
1310 for (
unsigned int i = 0; i < indices.size(); ++i)
1320 this->vector_ptr.resize(cache.size() + 1);
1321 this->matrix_ptr.resize(cache.size() + 1);
1323 for (
const auto &i : cache)
1325 const auto &
M = i.first.second.first;
1327 this->vector_ptr[i.second + 1] =
M.n_rows();
1328 this->matrix_ptr[i.second + 1] =
M.n_rows() *
M.n_cols();
1331 for (
unsigned int i = 0; i < cache.size(); ++i)
1333 this->vector_ptr[i + 1] += this->vector_ptr[i];
1334 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
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());
1342 for (
const auto &i : cache)
1348 if (precompute_inverse_diagonal)
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);
1367 if (indices.size() > 0)
1370 const unsigned int n_cells = indices.size() / dim;
1371 std::map<std::array<unsigned int, dim>,
unsigned int>
cache_ev;
1374 for (
unsigned int i = 0, c = 0; i <
n_cells; ++i)
1376 std::array<unsigned int, dim> id;
1378 for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1397 new_indices.reserve(indices.size() / dim * (dim + 1));
1399 for (
unsigned int i = 0, c = 0; i <
n_cells; ++i)
1401 for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1410 indices_ev[entry.second * dim + d] = entry.first[d];
1416 const unsigned int n_diag =
1418 (matrix_ptr.size() - 1)) /
1424 for (
unsigned int i = 0; i <
n_diag; ++i)
1426 const unsigned int c = (
indices_ev.size() > 0) ?
1430 const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c];
1436 for (
unsigned int i = 0; i <
n_diag; ++i)
1442 for (
unsigned int i = 0; i <
n_diag; ++i)
1444 std::array<Number *, dim>
evs;
1446 for (
unsigned int d = 0;
d < dim; ++
d)
1456 for (
unsigned int i1 = 0, c = 0;
i1 <
mm; ++
i1)
1457 for (
unsigned int i0 = 0;
i0 <
mm; ++
i0, ++c)
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)
1476 this->eigenvalues.clear();
1482template <
int dim,
typename Number,
int n_rows_1d>
1489 Number *dst =
dst_in.begin();
1490 const Number *src =
src_in.begin();
1492 if (this->eigenvalues.empty() ==
false)
1498 for (
unsigned int d = 0;
d < dim; ++
d)
1501 (indices.size() > 0) ? indices[dim * index + d] : (dim *
index +
d);
1512 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1516 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1522 const Number *inverted_eigenvalues =
nullptr;
1525 for (
unsigned int d = 0;
d < dim; ++
d)
1528 (indices.size() > 0) ?
1529 indices[((dim == 1) ? 1 : (dim + 1)) *
index +
d] :
1538 ((indices.size() > 0) && (dim != 1)) ?
1539 indices[(dim + 1) *
index + dim] :
1542 inverted_eigenvalues =
1551 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1557 inverted_eigenvalues);
1559 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1565 inverted_eigenvalues);
1571template <
int dim,
typename Number,
int n_rows_1d>
1587template <
int dim,
typename Number,
int n_rows_1d>
1592 if (matrix_ptr.empty())
1595 return matrix_ptr.size() - 1;
std::size_t storage_size() const
AlignedVector< Number > mass_matrices
const bool precompute_inverse_diagonal
std::size_t memory_consumption() const
std::vector< unsigned int > matrix_ptr
void apply_inverse(const unsigned int index, const ArrayView< Number > &dst_in, const ArrayView< const Number > &src_in) const
AlignedVector< Number > eigenvectors
std::vector< unsigned int > vector_ptr
const bool compress_matrices
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)
AlignedVector< Number > inverted_eigenvalues
std::vector< unsigned int > indices
std::vector< unsigned int > vector_n_rows_1d
std::map< MatrixPairTypeWithMask, unsigned int, internal::TensorProductMatrixSymmetricSum::MatrixPairComparator< Number > > cache
AlignedVector< Number > eigenvalues
AlignedVector< Number > derivative_matrices
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
AlignedVector< Number > tmp_array
static constexpr int n_rows_1d_static
std::size_t memory_consumption() const
TensorProductMatrixSymmetricSum()=default
std::array< Table< 2, Number >, dim > mass_matrix
std::array< AlignedVector< Number >, dim > eigenvalues
TensorProductMatrixSymmetricSum(const T &mass_matrix, const T &derivative_matrix)
static constexpr std::size_t size()
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
@ matrix
Contents is actually a matrix.
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
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)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
bool precompute_inverse_diagonal
AdditionalData(const bool compress_matrices=true, const bool precompute_inverse_diagonal=true)
typename VectorizedArrayTrait::value_type ScalarNumber
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
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)