Skip to content

Commit 3a3d284

Browse files
aabrown100-gitaabrown100-git
committed
Adding/moving templated Eigen matrix/tensor function implementations to header; using fixed size Matrix and Tensors, instead of just Tensors. Compiles and passes integration tests/unit tests that use NeoHookean
1 parent 30cb13e commit 3a3d284

File tree

4 files changed

+199
-95
lines changed

4 files changed

+199
-95
lines changed

‎Code/Source/svFSI/mat_fun.cpp‎

Lines changed: 0 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -829,23 +829,6 @@ ten_dyad_prod(const Array<double>& A, const Array<double>& B, const int nd)
829829
return C;
830830
}
831831

832-
833-
Eigen::Tensor<double, 4> ten_dyad_prod_eigen(const Eigen::Tensor<double, 2>& A, const Eigen::Tensor<double, 2>& B, const int nsd) {
834-
835-
Eigen::Tensor<double, 4> C(nsd, nsd, nsd, nsd);
836-
837-
for (int i = 0; i < nsd; ++i) {
838-
for (int j = 0; j < nsd; ++j) {
839-
for (int k = 0; k < nsd; ++k) {
840-
for (int l = 0; l < nsd; ++l) {
841-
C(i, j, k, l) = A(i, j) * B(k, l);
842-
}
843-
}
844-
}
845-
}
846-
847-
return C;
848-
}
849832
/// @brief Create a 4th order order symmetric identity tensor
850833
//
851834
Tensor4<double>

‎Code/Source/svFSI/mat_fun.h‎

Lines changed: 73 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@
3333
#include <eigen3/Eigen/Core>
3434
#include <eigen3/Eigen/Dense>
3535
#include <unsupported/Eigen/CXX11/Tensor>
36+
#include <stdexcept>
3637

3738
#include "Array.h"
3839
#include "Tensor4.h"
@@ -46,6 +47,32 @@
4647
/// \todo [TODO:DaveP] this should just be a namespace?
4748
//
4849
namespace mat_fun {
50+
// Helper function to convert Array<double> to Eigen::Matrix
51+
template <typename MatrixType>
52+
MatrixType toEigenMatrix(const Array<double>& src) {
53+
MatrixType mat;
54+
for (int i = 0; i < mat.rows(); ++i)
55+
for (int j = 0; j < mat.cols(); ++j)
56+
mat(i, j) = src(i, j);
57+
return mat;
58+
}
59+
60+
// Helper function to convert Eigen::Matrix to Array<double>
61+
template <typename MatrixType>
62+
void toArray(const MatrixType& mat, Array<double>& dest) {
63+
for (int i = 0; i < mat.rows(); ++i)
64+
for (int j = 0; j < mat.cols(); ++j)
65+
dest(i, j) = mat(i, j);
66+
}
67+
68+
// Helper function to convert a higher-dimensional array like Dm
69+
template <typename MatrixType>
70+
void copyDm(const MatrixType& mat, Array<double>& dest, int rows, int cols) {
71+
for (int i = 0; i < rows; ++i)
72+
for (int j = 0; j < cols; ++j)
73+
dest(i, j) = mat(i, j);
74+
}
75+
4976
double mat_ddot(const Array<double>& A, const Array<double>& B, const int nd);
5077
double mat_det(const Array<double>& A, const int nd);
5178
Array<double> mat_dev(const Array<double>& A, const int nd);
@@ -72,11 +99,56 @@ namespace mat_fun {
7299
Tensor4<double> ten_ddot_3424(const Tensor4<double>& A, const Tensor4<double>& B, const int nd);
73100

74101
Tensor4<double> ten_dyad_prod(const Array<double>& A, const Array<double>& B, const int nd);
75-
Eigen::Tensor<double, 4> ten_dyad_prod_eigen(const Eigen::Tensor<double, 2>& A, const Eigen::Tensor<double, 2>& B, const int nsd);
102+
103+
template <int nsd>
104+
Eigen::TensorFixedSize<double, Eigen::Sizes<nsd, nsd, nsd, nsd>>
105+
ten_dyad_prod_eigen(const Eigen::Matrix<double, nsd, nsd>& A, const Eigen::Matrix<double, nsd, nsd>& B) {
106+
// Create a fixed-size tensor to store the result
107+
Eigen::TensorFixedSize<double, Eigen::Sizes<nsd, nsd, nsd, nsd>> C;
108+
109+
// Perform the dyadic product (outer product) and assign values to the tensor
110+
for (int i = 0; i < nsd; ++i) {
111+
for (int j = 0; j < nsd; ++j) {
112+
for (int k = 0; k < nsd; ++k) {
113+
for (int l = 0; l < nsd; ++l) {
114+
C(i, j, k, l) = A(i, j) * B(k, l);
115+
}
116+
}
117+
}
118+
}
119+
120+
return C;
121+
}
122+
76123
Tensor4<double> ten_ids(const int nd);
77124
Array<double> ten_mddot(const Tensor4<double>& A, const Array<double>& B, const int nd);
78125

79126
Tensor4<double> ten_symm_prod(const Array<double>& A, const Array<double>& B, const int nd);
127+
128+
/// @brief Create a 4th order tensor from symmetric outer product of two matrices.
129+
///
130+
/// Reproduces 'FUNCTION TEN_SYMMPROD(A, B, nd) RESULT(C)'.
131+
//
132+
template <int nsd>
133+
Eigen::TensorFixedSize<double, Eigen::Sizes<nsd, nsd, nsd, nsd>>
134+
ten_symm_prod_eigen(const Eigen::Matrix<double, nsd, nsd>& A, const Eigen::Matrix<double, nsd, nsd>& B) {
135+
// Create a fixed-size tensor to store the result
136+
Eigen::TensorFixedSize<double, Eigen::Sizes<nsd, nsd, nsd, nsd>> C;
137+
138+
// Perform the symmetric product and assign values to the tensor
139+
for (int i = 0; i < nsd; ++i) {
140+
for (int j = 0; j < nsd; ++j) {
141+
for (int k = 0; k < nsd; ++k) {
142+
for (int l = 0; l < nsd; ++l) {
143+
C(i, j, k, l) = 0.5 * (A(i, k) * B(j, l) + A(i, l) * B(j, k));
144+
}
145+
}
146+
}
147+
}
148+
149+
return C;
150+
}
151+
80152
Tensor4<double> ten_transpose(const Tensor4<double>& A, const int nd);
81153

82154
Array<double> transpose(const Array<double>& A);

‎Code/Source/svFSI/mat_models.cpp‎

Lines changed: 126 additions & 73 deletions
Original file line numberDiff line numberDiff line change
@@ -119,6 +119,59 @@ void cc_to_voigt(const int nsd, const Tensor4<double>& CC, Array<double>& Dm)
119119
}
120120
}
121121

122+
template <int nsd>
123+
void cc_to_voigt_eigen(const Eigen::TensorFixedSize<double, Eigen::Sizes<nsd, nsd, nsd, nsd>>& CC, Eigen::Matrix<double, 2*nsd, 2*nsd>& Dm)
124+
{
125+
if (nsd == 3) {
126+
Dm(0,0) = CC(0,0,0,0);
127+
Dm(0,1) = CC(0,0,1,1);
128+
Dm(0,2) = CC(0,0,2,2);
129+
Dm(0,3) = CC(0,0,0,1);
130+
Dm(0,4) = CC(0,0,1,2);
131+
Dm(0,5) = CC(0,0,2,0);
132+
133+
Dm(1,1) = CC(1,1,1,1);
134+
Dm(1,2) = CC(1,1,2,2);
135+
Dm(1,3) = CC(1,1,0,1);
136+
Dm(1,4) = CC(1,1,1,2);
137+
Dm(1,5) = CC(1,1,2,0);
138+
139+
Dm(2,2) = CC(2,2,2,2);
140+
Dm(2,3) = CC(2,2,0,1);
141+
Dm(2,4) = CC(2,2,1,2);
142+
Dm(2,5) = CC(2,2,2,0);
143+
144+
Dm(3,3) = CC(0,1,0,1);
145+
Dm(3,4) = CC(0,1,1,2);
146+
Dm(3,5) = CC(0,1,2,0);
147+
148+
Dm(4,4) = CC(1,2,1,2);
149+
Dm(4,5) = CC(1,2,2,0);
150+
151+
Dm(5,5) = CC(2,0,2,0);
152+
153+
for (int i = 1; i < 6; i++) {
154+
for (int j = 0; j <= i-1; j++) {
155+
Dm(i,j) = Dm(j,i);
156+
}
157+
}
158+
159+
} else if (nsd == 2) {
160+
Dm(0,0) = CC(0,0,0,0);
161+
Dm(0,1) = CC(0,0,1,1);
162+
Dm(0,2) = CC(0,0,0,1);
163+
164+
Dm(1,1) = CC(1,1,1,1);
165+
Dm(1,2) = CC(1,1,0,1);
166+
167+
Dm(2,2) = CC(0,1,0,1);
168+
169+
Dm(1,0) = Dm(0,1);
170+
Dm(2,0) = Dm(0,2);
171+
Dm(2,1) = Dm(1,2);
172+
}
173+
}
174+
122175
void voigt_to_cc(const int nsd, const Array<double>& Dm, Tensor4<double>& CC)
123176
{
124177
if (nsd == 3) {
@@ -187,67 +240,7 @@ void get_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibStrsT
187240
}
188241
}
189242

190-
/**
191-
* @brief Helper function to handle different number of spatial dimensions.
192-
*
193-
*/
194-
void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Array<double>& F, const int nfd,
195-
const Array<double> fl, const double ya, Array<double>& S, Array<double>& Dm, double& Ja)
196-
{
197-
// Number of spatial dimensions
198-
int nsd = com_mod.nsd;
199-
200-
if (nsd == 2) {
201-
Eigen::TensorFixedSize<double,Eigen::Sizes<2,2>> F_2D;
202-
F_2D.setValues( {{F(0,0), F(0,1)},{F(1,0), F(1,1)} });
203243

204-
Eigen::TensorFixedSize<double,Eigen::Sizes<2,2>> S_2D;
205-
S_2D.setValues( {{S(0,0), S(0,1)},{S(1,0), S(1,1)} });
206-
207-
Eigen::TensorFixedSize<double,Eigen::Sizes<4,4>> Dm_2D;
208-
Dm_2D.setValues( {{Dm(0,0), Dm(0,1), Dm(0,2), Dm(0,3)},
209-
{Dm(1,0), Dm(1,1), Dm(1,2), Dm(1,3)},
210-
{Dm(2,0), Dm(2,1), Dm(2,2), Dm(2,3)},
211-
{Dm(3,0), Dm(3,1), Dm(3,2), Dm(3,3)} });
212-
213-
_get_pk2cc<2>(com_mod, cep_mod, lDmn, F_2D, nfd, fl, ya, S_2D, Dm_2D, Ja);
214-
215-
S(0,0) = S_2D(0,0); S(0,1) = S_2D(0,1);
216-
S(1,0) = S_2D(1,0); S(1,1) = S_2D(1,1);
217-
218-
Dm(0,0) = Dm_2D(0,0); Dm(0,1) = Dm_2D(0,1);
219-
Dm(1,0) = Dm_2D(1,0); Dm(1,1) = Dm_2D(1,1);
220-
}
221-
else if (nsd == 3) {
222-
Eigen::TensorFixedSize<double,Eigen::Sizes<3,3>>F_3D;
223-
F_3D.setValues( {{F(0,0), F(0,1), F(0,2)},
224-
{F(1,0), F(1,1), F(1,2)},
225-
{F(2,0), F(2,1), F(2,2)} });
226-
227-
Eigen::TensorFixedSize<double,Eigen::Sizes<3,3>> S_3D;
228-
S_3D.setValues( {{S(0,0), S(0,1), S(0,2)},
229-
{S(1,0), S(1,1), S(1,2)},
230-
{S(2,0), S(2,1), S(2,2)} });
231-
Eigen::TensorFixedSize<double,Eigen::Sizes<6,6>>Dm_3D;
232-
Dm_3D.setValues( {{Dm(0,0), Dm(0,1), Dm(0,2), Dm(0,3), Dm(0,4), Dm(0,5)},
233-
{Dm(1,0), Dm(1,1), Dm(1,2), Dm(1,3), Dm(1,4), Dm(1,5)},
234-
{Dm(2,0), Dm(2,1), Dm(2,2), Dm(2,3), Dm(2,4), Dm(2,5)},
235-
{Dm(3,0), Dm(3,1), Dm(3,2), Dm(3,3), Dm(3,4), Dm(3,5)},
236-
{Dm(4,0), Dm(4,1), Dm(4,2), Dm(4,3), Dm(4,4), Dm(4,5)},
237-
{Dm(5,0), Dm(5,1), Dm(5,2), Dm(5,3), Dm(5,4), Dm(5,5)} });
238-
239-
_get_pk2cc<3>(com_mod, cep_mod, lDmn, F_3D, nfd, fl, ya, S_3D, Dm_3D, Ja);
240-
241-
S(0,0) = S_3D(0,0); S(0,1) = S_3D(0,1); S(0,2) = S_3D(0,2);
242-
S(1,0) = S_3D(1,0); S(1,1) = S_3D(1,1); S(1,2) = S_3D(1,2);
243-
S(2,0) = S_3D(2,0); S(2,1) = S_3D(2,1); S(2,2) = S_3D(2,2);
244-
245-
Dm(0,0) = Dm_3D(0,0); Dm(0,1) = Dm_3D(0,1); Dm(0,2) = Dm_3D(0,2);
246-
Dm(1,0) = Dm_3D(1,0); Dm(1,1) = Dm_3D(1,1); Dm(1,2) = Dm_3D(1,2);
247-
Dm(2,0) = Dm_3D(2,0); Dm(2,1) = Dm_3D(2,1); Dm(2,2) = Dm_3D(2,2);
248-
}
249-
250-
}
251244

252245

253246
/**
@@ -269,8 +262,8 @@ else if (nsd == 3) {
269262
* @return None, but modifies S, Dm, and Ja in place.
270263
*/
271264
template<size_t nsd>
272-
void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Eigen::TensorFixedSize<double, Eigen::Sizes<nsd,nsd>>& F, const int nfd,
273-
const Array<double>& fl, const double ya, Eigen::TensorFixedSize<double, Eigen::Sizes<nsd,nsd>>& S, Eigen::TensorFixedSize<double, Eigen::Sizes<2*nsd,2*nsd>>& Dm, double& Ja)
265+
void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Eigen::Matrix<double, nsd, nsd>& F, const int nfd,
266+
const Eigen::Matrix<double, nsd, Eigen::Dynamic> fl, const double ya, Eigen::Matrix<double, nsd, nsd>& S, Eigen::Matrix<double, 2*nsd, 2*nsd>& Dm, double& Ja)
274267
{
275268
using namespace consts;
276269
using namespace mat_fun;
@@ -290,7 +283,6 @@ void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDm
290283
// ustruct flag
291284
bool ustruct = (lDmn.phys == EquationType::phys_ustruct);
292285

293-
294286
S.setZero();
295287
Dm.setZero();
296288

@@ -313,7 +305,12 @@ void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDm
313305
auto Fe = F;
314306
auto Fa = Eigen::Matrix<double, nsd, nsd>::Identity();
315307
auto Fai = Fa;
316-
308+
309+
// if (cep_mod.cem.aStrain) {
310+
// actv_strain(com_mod, cep_mod, ya, nfd, fl, Fa);
311+
// Fai = Fa.inverse();
312+
// Fe = F * Fai;
313+
// }
317314

318315
Ja = Fa.determinant();
319316
double J = Fe.determinant();
@@ -342,27 +339,24 @@ void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDm
342339

343340
// Now, compute isochoric and total stress, elasticity tensors
344341
//
345-
Eigen::Tensor<double, 4> CC(nsd,nsd,nsd,nsd);
342+
Eigen::TensorFixedSize<double, Eigen::Sizes<nsd,nsd,nsd,nsd>> CC;
346343

347344
switch (stM.isoType) {
348345

349346
// NeoHookean model
350347
case ConstitutiveModelType::stIso_nHook: {
351348
double g1 = 2.0 * stM.C10;
352-
auto Sb = g1*Idm;
349+
Eigen::Matrix<double, nsd, nsd> Sb = g1*Idm;
353350

354351
// Fiber reinforcement/active stress
355-
//Sb += Tfa * mat_dyad_prod(fl.col(0), fl.col(0), nsd);
352+
Sb += Tfa * (fl.col(0) * fl.col(0).transpose());
356353

357354
double r1 = g1 * Inv1 / nd;
358355
S = J2d*Sb - r1*Ci;
359356

360-
//CC = ten_dyad_prod_eigen(Ci, S, nsd);
361-
CC = (-2.0/nd) * ( ten_dyad_prod_eigen(Ci, S, nsd) + ten_dyad_prod_eigen(S, Ci, nsd));
362-
//CC = (-2.0/nd) * (Eigen::kroneckerProduct(Ci, S.transpose()) + Eigen::kroneckerProduct(S , Ci.transpose()));
357+
CC = (-2.0/nd) * ( ten_dyad_prod_eigen<nsd>(Ci, S) + ten_dyad_prod_eigen<nsd>(S, Ci));
363358
S += p*J*Ci;
364-
//CC += 2.0*(r1 - p*J) * ten_symm_prod(Ci, Ci, nsd) + (pl*J - 2.0*r1/nd) * ten_dyad_prod(Ci, Ci, nsd);
365-
//CC += 2.0*(r1 - p*J) * (Ci * Ci.transpose()).symmetricView() + (pl*J - 2.0*r1/nd) * (Ci * Ci.transpose());
359+
CC += 2.0*(r1 - p*J) * ten_symm_prod_eigen<nsd>(Ci, Ci) + (pl*J - 2.0*r1/nd) * ten_dyad_prod_eigen<nsd>(Ci, Ci);
366360

367361
} break;
368362

@@ -371,9 +365,68 @@ void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDm
371365
}
372366

373367
// Convert to Voigt Notation
374-
//cc_to_voigt(nsd, CC, Dm);
368+
cc_to_voigt_eigen<nsd>(CC, Dm);
375369
}
376370

371+
/**
372+
* @brief Get the 2nd Piola-Kirchhoff stress tensor and material elasticity tensor.
373+
*
374+
* This is a wrapper function for the templated function _get_pk2cc.
375+
*
376+
*/
377+
void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Array<double>& F, const int nfd,
378+
const Array<double>& fl, const double ya, Array<double>& S, Array<double>& Dm, double& Ja)
379+
{
380+
// Number of spatial dimensions
381+
int nsd = com_mod.nsd;
382+
383+
if (nsd == 2) {
384+
// Copy deformation gradient to Eigen matrix
385+
auto F_2D = mat_fun::toEigenMatrix<Eigen::Matrix2d>(F);
386+
387+
// Copy fiber directions to Eigen matrix
388+
Eigen::Matrix<double, 2, Eigen::Dynamic> fl_2D(2, nfd);
389+
for (int i = 0; i < nfd; i++) {
390+
fl_2D(0, i) = fl(0, i);
391+
fl_2D(1, i) = fl(1, i);
392+
}
393+
394+
// Initialize stress and elasticity tensors
395+
Eigen::Matrix2d S_2D = Eigen::Matrix2d::Zero();
396+
Eigen::Matrix4d Dm_2D = Eigen::Matrix4d::Zero();
397+
398+
// Call templated function
399+
_get_pk2cc<2>(com_mod, cep_mod, lDmn, F_2D, nfd, fl_2D, ya, S_2D, Dm_2D, Ja);
400+
401+
// Copy results back
402+
mat_fun::toArray(S_2D, S);
403+
mat_fun::copyDm(Dm_2D, Dm, 4, 4);
404+
405+
} else if (nsd == 3) {
406+
// Copy deformation gradient to Eigen matrix
407+
auto F_3D = mat_fun::toEigenMatrix<Eigen::Matrix3d>(F);
408+
409+
// Copy fiber directions to Eigen matrix
410+
Eigen::Matrix<double, 3, Eigen::Dynamic> fl_3D(3, nfd);
411+
for (int i = 0; i < nfd; i++) {
412+
fl_3D(0, i) = fl(0, i);
413+
fl_3D(1, i) = fl(1, i);
414+
fl_3D(2, i) = fl(2, i);
415+
}
416+
417+
// Initialize stress and elasticity tensors
418+
Eigen::Matrix3d S_3D = Eigen::Matrix3d::Zero();
419+
Eigen::Matrix<double, 6, 6> Dm_3D;
420+
Dm_3D.setZero();
421+
422+
// Call templated function
423+
_get_pk2cc<3>(com_mod, cep_mod, lDmn, F_3D, nfd, fl_3D, ya, S_3D, Dm_3D, Ja);
424+
425+
// Copy results back
426+
mat_fun::toArray(S_3D, S);
427+
mat_fun::copyDm(Dm_3D, Dm, 6, 6);
428+
}
429+
}
377430

378431
/// @brief Compute 2nd Piola-Kirchhoff stress and material stiffness tensors
379432
/// for compressible shell elements.

‎Code/Source/svFSI/mat_models.h‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -56,10 +56,6 @@ void get_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibStrsT
5656
void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Array<double>& F, const int nfd,
5757
const Array<double>& fl, const double ya, Array<double>& S, Array<double>& Dm, double& Ja);
5858

59-
template<size_t nsd>
60-
void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Eigen::TensorFixedSize<double, Eigen::Sizes<nsd,nsd>>& F, const int nfd,
61-
const Array<double>& fl, const double ya, Eigen::TensorFixedSize<double, Eigen::Sizes<nsd,nsd>>& S, Eigen::TensorFixedSize<double, Eigen::Sizes<2*nsd,2*nsd>>& Dm, double& Ja);
62-
6359
void get_pk2cc_shlc(const ComMod& com_mod, const dmnType& lDmn, const int nfd, const Array<double>& fNa0,
6460
const Array<double>& gg_0, const Array<double>& gg_x, double& g33, Vector<double>& Sml, Array<double>& Dml);
6561

0 commit comments

Comments
 (0)