Skip to content

Commit 30cb13e

Browse files
Marisa BazziMarisa Bazzi
Marisa Bazzi
authored and
Marisa Bazzi
committed
Change from Matrix to Tensor
1 parent dccdc64 commit 30cb13e

File tree

4 files changed

+75
-23
lines changed

4 files changed

+75
-23
lines changed

‎Code/Source/svFSI/mat_fun.cpp‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@
3636
*/
3737

3838
#include "consts.h"
39+
3940
#include "utils.h"
4041

4142
#include <math.h>
@@ -828,6 +829,23 @@ ten_dyad_prod(const Array<double>& A, const Array<double>& B, const int nd)
828829
return C;
829830
}
830831

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+
}
831849
/// @brief Create a 4th order order symmetric identity tensor
832850
//
833851
Tensor4<double>

‎Code/Source/svFSI/mat_fun.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,9 @@
3030

3131
#ifndef MAT_FUN_H
3232
#define MAT_FUN_H
33+
#include <eigen3/Eigen/Core>
34+
#include <eigen3/Eigen/Dense>
35+
#include <unsupported/Eigen/CXX11/Tensor>
3336

3437
#include "Array.h"
3538
#include "Tensor4.h"
@@ -69,6 +72,7 @@ namespace mat_fun {
6972
Tensor4<double> ten_ddot_3424(const Tensor4<double>& A, const Tensor4<double>& B, const int nd);
7073

7174
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);
7276
Tensor4<double> ten_ids(const int nd);
7377
Array<double> ten_mddot(const Tensor4<double>& A, const Array<double>& B, const int nd);
7478

‎Code/Source/svFSI/mat_models.cpp‎

Lines changed: 49 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -191,31 +191,60 @@ void get_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibStrsT
191191
* @brief Helper function to handle different number of spatial dimensions.
192192
*
193193
*/
194-
void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Eigen::MatrixXd& F, const int nfd,
195-
const Array<double> fl, const double ya, Eigen::MatrixXd& S, Eigen::MatrixXd& Dm, double& Ja)
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)
196196
{
197197
// Number of spatial dimensions
198-
int nsd = F.rows();
198+
int nsd = com_mod.nsd;
199199

200200
if (nsd == 2) {
201-
Eigen::Matrix2d F_2D = F.topLeftCorner(2,2);
202-
Eigen::Matrix2d S_2D = S.topLeftCorner(2,2);
203-
Eigen::Matrix2d Dm_2D = Dm.topLeftCorner(3,3);
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)} });
203+
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)} });
204212

205213
_get_pk2cc<2>(com_mod, cep_mod, lDmn, F_2D, nfd, fl, ya, S_2D, Dm_2D, Ja);
206214

207-
S.topLeftCorner(2,2) = S_2D;
208-
Dm.topLeftCorner(3,3) = Dm_2D;
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);
209220
}
210-
else if (nsd == 3) {
211-
Eigen::Matrix3d F_3D = F.topLeftCorner(3,3);
212-
Eigen::Matrix3d S_3D = S.topLeftCorner(3,3);
213-
Eigen::Matrix3d Dm_3D = Dm.topLeftCorner(6,6);
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)} });
214238

215239
_get_pk2cc<3>(com_mod, cep_mod, lDmn, F_3D, nfd, fl, ya, S_3D, Dm_3D, Ja);
216240

217-
S.topLeftCorner(3,3) = S_3D;
218-
Dm.topLeftCorner(6,6) = Dm_3D;
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);
219248
}
220249

221250
}
@@ -240,8 +269,8 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn
240269
* @return None, but modifies S, Dm, and Ja in place.
241270
*/
242271
template<size_t nsd>
243-
void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Eigen::Matrix<double, nsd, nsd>& F, const int nfd,
244-
const Array<double>& fl, const double ya, Eigen::Matrix<double, nsd, nsd>& S, Eigen::Matrix<double, nsd, nsd>& Dm, double& Ja)
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)
245274
{
246275
using namespace consts;
247276
using namespace mat_fun;
@@ -284,7 +313,7 @@ void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDm
284313
auto Fe = F;
285314
auto Fa = Eigen::Matrix<double, nsd, nsd>::Identity();
286315
auto Fai = Fa;
287-
316+
288317

289318
Ja = Fa.determinant();
290319
double J = Fe.determinant();
@@ -328,8 +357,9 @@ void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDm
328357
double r1 = g1 * Inv1 / nd;
329358
S = J2d*Sb - r1*Ci;
330359

331-
//CC = (-2.0/nd) * ( ten_dyad_prod(Ci, S, nsd) + ten_dyad_prod(S, Ci, nsd));
332-
CC = (-2.0/nd) * (Ci * S.transpose() + S * Ci.transpose());
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()));
333363
S += p*J*Ci;
334364
//CC += 2.0*(r1 - p*J) * ten_symm_prod(Ci, Ci, nsd) + (pl*J - 2.0*r1/nd) * ten_dyad_prod(Ci, Ci, nsd);
335365
//CC += 2.0*(r1 - p*J) * (Ci * Ci.transpose()).symmetricView() + (pl*J - 2.0*r1/nd) * (Ci * Ci.transpose());

‎Code/Source/svFSI/mat_models.h‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -53,12 +53,12 @@ void voigt_to_cc(const int nsd, const Array<double>& Dm, Tensor4<double>& CC);
5353

5454
void get_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibStrsType& Tfl, double& g);
5555

56-
void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Eigen::MatrixXd& F, const int nfd,
57-
const Array<double>& fl, const double ya, Eigen::MatrixXd& S, Eigen::MatrixXd& Dm, double& Ja);
56+
void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Array<double>& F, const int nfd,
57+
const Array<double>& fl, const double ya, Array<double>& S, Array<double>& Dm, double& Ja);
5858

5959
template<size_t nsd>
60-
void _get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Eigen::Matrix<double, nsd, nsd>& F, const int nfd,
61-
const Array<double>& fl, const double ya, Eigen::Matrix<double, nsd, nsd>& S, Eigen::Matrix<double, nsd, nsd>& Dm, double& Ja);
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);
6262

6363
void get_pk2cc_shlc(const ComMod& com_mod, const dmnType& lDmn, const int nfd, const Array<double>& fNa0,
6464
const Array<double>& gg_0, const Array<double>& gg_x, double& g33, Vector<double>& Sml, Array<double>& Dml);

0 commit comments

Comments
 (0)