33#include "eigen3/Eigen/Core"
34#include "eigen3/Eigen/Dense"
35#include "eigen3/unsupported/Eigen/CXX11/Tensor"
52 using Matrix = Eigen::Matrix<double, nsd, nsd>;
55 using Tensor = Eigen::TensorFixedSize<double, Eigen::Sizes<nsd, nsd, nsd, nsd>>;
58 template <
typename MatrixType>
59 MatrixType convert_to_eigen_matrix(
const Array<double>& src) {
61 for (
int i = 0; i < mat.rows(); ++i)
62 for (
int j = 0; j < mat.cols(); ++j)
63 mat(i, j) = src(i, j);
68 template <
typename MatrixType>
69 void convert_to_array(
const MatrixType& mat, Array<double>& dest) {
70 for (
int i = 0; i < mat.rows(); ++i)
71 for (
int j = 0; j < mat.cols(); ++j)
72 dest(i, j) = mat(i, j);
76 template <
typename MatrixType>
77 void copy_Dm(
const MatrixType& mat, Array<double>& dest,
int rows,
int cols) {
78 for (
int i = 0; i < rows; ++i)
79 for (
int j = 0; j < cols; ++j)
80 dest(i, j) = mat(i, j);
84 Eigen::Matrix<double, nsd, 1> cross_product(
const Eigen::Matrix<double, nsd, 1>& u,
const Eigen::Matrix<double, nsd, 1>& v) {
85 if constexpr (nsd == 2) {
86 return Eigen::Matrix<double, 2, 1>(v(1), - v(0));
88 else if constexpr (nsd == 3) {
92 throw std::runtime_error(
"[cross_product] Invalid number of spatial dimensions '" + std::to_string(nsd) +
"'. Valid dimensions are 2 or 3.");
96 double mat_ddot(
const Array<double>& A,
const Array<double>& B,
const int nd);
99 double double_dot_product(
const Matrix<nsd>& A,
const Matrix<nsd>& B) {
100 return A.cwiseProduct(B).sum();
103 double mat_det(
const Array<double>& A,
const int nd);
104 Array<double> mat_dev(
const Array<double>& A,
const int nd);
108 Array<double> mat_id(
const int nsd);
109 Array<double>
mat_inv(
const Array<double>& A,
const int nd,
bool debug =
false);
110 Array<double>
mat_inv_ge(
const Array<double>& A,
const int nd,
bool debug =
false);
111 Array<double>
mat_inv_lp(
const Array<double>& A,
const int nd);
114 Array<double>
mat_mul(
const Array<double>& A,
const Array<double>& B);
115 void mat_mul(
const Array<double>& A,
const Array<double>& B, Array<double>& result);
116 void mat_mul6x3(
const Array<double>& A,
const Array<double>& B, Array<double>& C);
118 Array<double>
mat_symm(
const Array<double>& A,
const int nd);
121 double mat_trace(
const Array<double>& A,
const int nd);
134 double_dot_product(
const Tensor<nsd>& A,
const std::array<int, 2>& dimsA,
135 const Tensor<nsd>& B,
const std::array<int, 2>& dimsB) {
138 Eigen::array<Eigen::IndexPair<int>, 2> contractionDims = {
139 Eigen::IndexPair<int>(dimsA[0], dimsB[0]),
140 Eigen::IndexPair<int>(dimsA[1], dimsB[1])
144 return A.contract(B, contractionDims);
167 for (
int i = 0; i < nsd; ++i) {
168 for (
int j = 0; j < nsd; ++j) {
169 for (
int k = 0; k < nsd; ++k) {
170 for (
int l = 0; l < nsd; ++l) {
171 C(i,j,k,l) = A(i,j) * B(k,l);
199 for (
int i = 0; i < nsd; ++i) {
200 for (
int j = 0; j < nsd; ++j) {
225 for (
int i = 0; i < nsd; ++i) {
226 for (
int j = 0; j < nsd; ++j) {
227 for (
int k = 0; k < nsd; ++k) {
228 for (
int l = 0; l < nsd; ++l) {
229 C(i,j,k,l) = 0.5 * (A(i,k) * B(j,l) + A(i,l) * B(j,k));
258 for (
int i = 0; i < nsd; ++i) {
259 for (
int j = 0; j < nsd; ++j) {
260 for (
int k = 0; k < nsd; ++k) {
261 for (
int l = 0; l < nsd; ++l) {
262 B(i,j,k,l) = A(k,l,i,j);
271 Array<double>
transpose(
const Array<double>& A);
The Tensor4 template class implements a simple interface to 4th order tensors.
Definition Tensor4.h:45
The Vector template class is used for storing int and double data.
Definition Vector.h:50
The classes defined here duplicate the data structures in the Fortran MATFUN module defined in MATFUN...
Definition mat_fun.cpp:41
Tensor< nsd > symmetric_dyadic_product(const Matrix< nsd > &A, const Matrix< nsd > &B)
Create a 4th order tensor from symmetric outer product of two matrices: C_ijkl = 0....
Definition mat_fun.h:219
void ten_init(const int nd)
Initialize tensor index pointer.
Definition mat_fun.cpp:787
Tensor4< double > ten_symm_prod(const Array< double > &A, const Array< double > &B, const int nd)
Create a 4th order tensor from symmetric outer product of two matrices.
Definition mat_fun.cpp:879
Array< double > ten_mddot(const Tensor4< double > &A, const Array< double > &B, const int nd)
Double dot product of a 4th order tensor and a 2nd order tensor.
Definition mat_fun.cpp:849
Array< double > mat_symm(const Array< double > &A, const int nd)
Symmetric part of a matrix, S = (A + A.T)/2.
Definition mat_fun.cpp:582
Tensor4< double > ten_ids(const int nd)
Create a 4th order order symmetric identity tensor.
Definition mat_fun.cpp:830
Array< double > transpose(const Array< double > &A)
Reproduces Fortran TRANSPOSE.
Definition mat_fun.cpp:915
Vector< double > mat_mul(const Array< double > &A, const Vector< double > &v)
Multiply a matrix by a vector.
Definition mat_fun.cpp:493
double mat_trace(const Array< double > &A, const int nd)
Trace of second order matrix of rank nd.
Definition mat_fun.cpp:614
Tensor4< double > ten_ddot(const Tensor4< double > &A, const Tensor4< double > &B, const int nd)
Double dot product of 2 4th order tensors T_ijkl = A_ijmn * B_klmn.
Definition mat_fun.cpp:653
Tensor4< double > ten_dyad_prod(const Array< double > &A, const Array< double > &B, const int nd)
Create a 4th order tensor from outer product of two matrices.
Definition mat_fun.cpp:811
Array< double > mat_inv_lp(const Array< double > &A, const int nd)
This function computes inverse of a square matrix using Lapack functions (DGETRF + DGETRI)
Definition mat_fun.cpp:423
Array< double > mat_inv(const Array< double > &A, const int nd, bool debug)
This function computes inverse of a square matrix.
Definition mat_fun.cpp:136
Array< double > mat_inv_ge(const Array< double > &Ain, const int n, bool debug)
This function computes inverse of a square matrix using Gauss Elimination method.
Definition mat_fun.cpp:193
Array< double > mat_symm_prod(const Vector< double > &u, const Vector< double > &v, const int nd)
Create a matrix from symmetric product of two vectors.
Definition mat_fun.cpp:599
Array< double > mat_dyad_prod(const Vector< double > &u, const Vector< double > &v, const int nd)
Create a matrix from outer product of two vectors.
Definition mat_fun.cpp:108
Tensor< nsd > fourth_order_identity()
Create a 4th order identity tensor: I_ijkl = 0.5 * (δ_ik * δ_jl + δ_il * δ_jk)
Definition mat_fun.h:193
double mat_ddot(const Array< double > &A, const Array< double > &B, const int nd)
Double dot product of 2 square matrices.
Definition mat_fun.cpp:47
Tensor4< double > ten_ddot_2412(const Tensor4< double > &A, const Tensor4< double > &B, const int nd)
T_ijkl = A_imjn * B_mnkl.
Definition mat_fun.cpp:698
Tensor4< double > ten_asym_prod12(const Array< double > &A, const Array< double > &B, const int nd)
Create a 4th order tensor from antisymmetric outer product of two matrices.
Definition mat_fun.cpp:631
Tensor< nsd > dyadic_product(const Matrix< nsd > &A, const Matrix< nsd > &B)
Compute the dyadic product of two 2nd order tensors A and B, C_ijkl = A_ij * B_kl.
Definition mat_fun.h:162