svMultiPhysics
Loading...
Searching...
No Matches
mat_fun.h
1/* Copyright (c) Stanford University, The Regents of the University of California, and others.
2 *
3 * All Rights Reserved.
4 *
5 * See Copyright-SimVascular.txt for additional details.
6 *
7 * Permission is hereby granted, free of charge, to any person obtaining
8 * a copy of this software and associated documentation files (the
9 * "Software"), to deal in the Software without restriction, including
10 * without limitation the rights to use, copy, modify, merge, publish,
11 * distribute, sublicense, and/or sell copies of the Software, and to
12 * permit persons to whom the Software is furnished to do so, subject
13 * to the following conditions:
14 *
15 * The above copyright notice and this permission notice shall be included
16 * in all copies or substantial portions of the Software.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
19 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
20 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
21 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
22 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
25 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
26 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
27 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30
31#ifndef MAT_FUN_H
32#define MAT_FUN_H
33#include "eigen3/Eigen/Core"
34#include "eigen3/Eigen/Dense"
35#include "eigen3/unsupported/Eigen/CXX11/Tensor"
36#include <stdexcept>
37
38#include "Array.h"
39#include "Tensor4.h"
40#include "Vector.h"
41
42/// @brief The classes defined here duplicate the data structures in the
43/// Fortran MATFUN module defined in MATFUN.f.
44///
45/// This module defines data structures for generally performed matrix and tensor operations.
46///
47/// \todo [TODO:DaveP] this should just be a namespace?
48//
49namespace mat_fun {
50 // Define templated type aliases for Eigen matrices and tensors for convenience
51 template<size_t nsd>
52 using Matrix = Eigen::Matrix<double, nsd, nsd>;
53
54 template<size_t nsd>
55 using Tensor = Eigen::TensorFixedSize<double, Eigen::Sizes<nsd, nsd, nsd, nsd>>;
56
57 // Function to convert Array<double> to Eigen::Matrix
58 template <typename MatrixType>
59 MatrixType convert_to_eigen_matrix(const Array<double>& src) {
60 MatrixType mat;
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);
64 return mat;
65 }
66
67 // Function to convert Eigen::Matrix to Array<double>
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);
73 }
74
75 // Function to convert a higher-dimensional array like Dm
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);
81 }
82
83 template <int nsd>
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));
87 }
88 else if constexpr (nsd == 3) {
89 return u.cross(v);
90 }
91 else {
92 throw std::runtime_error("[cross_product] Invalid number of spatial dimensions '" + std::to_string(nsd) + "'. Valid dimensions are 2 or 3.");
93 }
94 }
95
96 double mat_ddot(const Array<double>& A, const Array<double>& B, const int nd);
97
98 template <int nsd>
99 double double_dot_product(const Matrix<nsd>& A, const Matrix<nsd>& B) {
100 return A.cwiseProduct(B).sum();
101 }
102
103 double mat_det(const Array<double>& A, const int nd);
104 Array<double> mat_dev(const Array<double>& A, const int nd);
105
106 Array<double> mat_dyad_prod(const Vector<double>& u, const Vector<double>& v, const int nd);
107
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);
112
113 Vector<double> mat_mul(const Array<double>& A, const Vector<double>& v);
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);
117
118 Array<double> mat_symm(const Array<double>& A, const int nd);
119 Array<double> mat_symm_prod(const Vector<double>& u, const Vector<double>& v, const int nd);
120
121 double mat_trace(const Array<double>& A, const int nd);
122
123 Tensor4<double> ten_asym_prod12(const Array<double>& A, const Array<double>& B, const int nd);
124 Tensor4<double> ten_ddot(const Tensor4<double>& A, const Tensor4<double>& B, const int nd);
125 Tensor4<double> ten_ddot_2412(const Tensor4<double>& A, const Tensor4<double>& B, const int nd);
126 Tensor4<double> ten_ddot_3424(const Tensor4<double>& A, const Tensor4<double>& B, const int nd);
127
128 /**
129 * @brief Contracts two 4th order tensors A and B over two dimensions,
130 *
131 */
132 template <int nsd>
133 Tensor<nsd>
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) {
136
137 // Define the contraction dimensions
138 Eigen::array<Eigen::IndexPair<int>, 2> contractionDims = {
139 Eigen::IndexPair<int>(dimsA[0], dimsB[0]), // Contract A's dimsA[0] with B's dimsB[0]
140 Eigen::IndexPair<int>(dimsA[1], dimsB[1]) // Contract A's dimsA[1] with B's dimsB[1]
141 };
142
143 // Return the double dot product
144 return A.contract(B, contractionDims);
145
146 // For some reason, in this case the Eigen::Tensor contract function is
147 // faster than a for loop implementation.
148 }
149
150 Tensor4<double> ten_dyad_prod(const Array<double>& A, const Array<double>& B, const int nd);
151
152 /**
153 * @brief Compute the dyadic product of two 2nd order tensors A and B, C_ijkl = A_ij * B_kl
154 *
155 * @tparam nsd, the number of spatial dimensions
156 * @param A, the first 2nd order tensor
157 * @param B, the second 2nd order tensor
158 * @return Tensor<nsd>
159 */
160 template <int nsd>
161 Tensor<nsd>
162 dyadic_product(const Matrix<nsd>& A, const Matrix<nsd>& B) {
163 // Initialize the result tensor
164 Tensor<nsd> C;
165
166 // Compute the dyadic product: C_ijkl = A_ij * B_kl
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);
172 }
173 }
174 }
175 }
176 // For some reason, in this case the Eigen::Tensor contract function is
177 // slower than the for loop implementation
178
179 return C;
180 }
181
182 Tensor4<double> ten_ids(const int nd);
183
184 /**
185 * @brief Create a 4th order identity tensor:
186 * I_ijkl = 0.5 * (δ_ik * δ_jl + δ_il * δ_jk)
187 *
188 * @tparam nsd, the number of spatial dimensions
189 * @return Tensor<nsd>
190 */
191 template <int nsd>
192 Tensor<nsd>
194 // Initialize as zero
195 Tensor<nsd> I;
196 I.setZero();
197
198 // Set only non-zero entries
199 for (int i = 0; i < nsd; ++i) {
200 for (int j = 0; j < nsd; ++j) {
201 I(i,j,i,j) += 0.5;
202 I(i,j,j,i) += 0.5;
203 }
204 }
205
206 return I;
207 }
208
209 Array<double> ten_mddot(const Tensor4<double>& A, const Array<double>& B, const int nd);
210
211 Tensor4<double> ten_symm_prod(const Array<double>& A, const Array<double>& B, const int nd);
212
213 /// @brief Create a 4th order tensor from symmetric outer product of two matrices: C_ijkl = 0.5 * (A_ik * B_jl + A_il * B_jk)
214 ///
215 /// Reproduces 'FUNCTION TEN_SYMMPROD(A, B, nd) RESULT(C)'.
216 //
217 template <int nsd>
218 Tensor<nsd>
219 symmetric_dyadic_product(const Matrix<nsd>& A, const Matrix<nsd>& B) {
220
221 // Initialize the result tensor
222 Tensor<nsd> C;
223
224 // Compute the symmetric product: C_ijkl = 0.5 * (A_ik * B_jl + A_il * B_jk)
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));
230 }
231 }
232 }
233 }
234 // For some reason, in this case the for loop implementation is faster
235 // than the Eigen::Tensor contract method
236
237 // Return the symmetric product
238 return C;
239 }
240
241 Tensor4<double> ten_transpose(const Tensor4<double>& A, const int nd);
242
243 /**
244 * @brief Performs a tensor transpose operation on a 4th order tensor A, B_ijkl = A_klij
245 *
246 * @tparam nsd, the number of spatial dimensions
247 * @param A, the input 4th order tensor
248 * @return Tensor<nsd>
249 */
250 template <int nsd>
251 Tensor<nsd>
252 transpose(const Tensor<nsd>& A) {
253
254 // Initialize the result tensor
255 Tensor<nsd> B;
256
257 // Permute the tensor indices to perform the transpose operation
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);
263 }
264 }
265 }
266 }
267
268 return B;
269 }
270
271 Array<double> transpose(const Array<double>& A);
272
273 void ten_init(const int nd);
274
275};
276
277#endif
278
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