Skip to content

Fiber Stretch and Stretch Rate Computations #550

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
kko27 merged 15 commits into from
May 25, 2026
Merged

Fiber Stretch and Stretch Rate Computations #550

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
e1bec2f
Added comments for fib_strech for readability
kko27 May 14, 2026
1b7e827
added stretch rate and minor fixes
kko27 May 14, 2026
ae9d115
fiber stretch test case added
kko27 May 15, 2026
23617fd
changed fib stretch rate calc to FD method
kko27 May 18, 2026
a626e66
added reference solution for test case
kko27 May 19, 2026
41c8262
added documentation to test case
kko27 May 19, 2026
ea04798
minor formatting
kko27 May 19, 2026
a70e028
latex formatting
kko27 May 19, 2026
cc47640
explicitly expose compute_fib_stretch function
kko27 May 20, 2026
2c9164a
renamed I4f to fibStretch
kko27 May 21, 2026
8a5f5ce
Refactor deformation gradient in post.cpp
kko27 May 21, 2026
0817896
code cleanup suggested by PR comments
kko27 May 21, 2026
66bdb9b
simplified function calls and arguments
kko27 May 22, 2026
4726d54
addresses copilot comments on dt failure and inconsistent gnn call
kko27 May 22, 2026
b0cdc87
Addresses element domain check and dummy variable for gnn call comments
kko27 May 22, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Code/Source/solver/cep_ion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,7 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt

if (msh.nFn != 0) {
Vector<double> sA(msh.nNo);
post::fib_strech(simulation, iEq, msh, solutions, sA);
post::fib_stretch(com_mod, iEq, msh, solutions.current.get_displacement(), sA);
for (int a = 0; a < msh.nNo; a++) {
int Ac = msh.gN(a);
I4f(Ac) = sA(a);
Expand Down
9 changes: 6 additions & 3 deletions Code/Source/solver/consts.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -342,9 +342,11 @@ enum class OutputNameType
outGrp_divV = 521,
outGrp_Visc = 522,
outGrp_fS = 523,
outGrp_C = 524,
outGrp_C = 524,
outGrp_I1 = 525,
outGrp_ionicState = 526,
outGrp_fibStretch = 527,
outGrp_fibStretchRate = 528,

out_velocity = 599,
out_pressure = 598,
Expand Down Expand Up @@ -373,7 +375,9 @@ enum class OutputNameType
out_viscosity = 575,
out_fibStrn = 574,
out_CGstrain = 573,
out_CGInv1 = 572
out_CGInv1 = 572,
out_fibStretch = 571,
out_fibStretchRate = 570
};

/// @brief Simulation output file types.
Expand Down Expand Up @@ -497,4 +501,3 @@ enum class LinearAlgebraType {
};

#endif

212 changes: 88 additions & 124 deletions Code/Source/solver/post.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

#include "post.h"

#include "FE/Common/FEException.h"
#include "all_fun.h"
#include "fluid.h"
#include "fs.h"
Expand All @@ -17,6 +18,35 @@

namespace post {

namespace {

Array<double> deformation_gradient(const Array<double>& Nx, const Array<double>& dl,
int nsd, int nNo, int eq_start)
{
auto F = mat_fun::mat_id(nsd);
for (int a = 0; a < nNo; a++) {
if (nsd == 3) {
F(0,0) = F(0,0) + Nx(0,a)*dl(eq_start,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(eq_start,a);
F(0,2) = F(0,2) + Nx(2,a)*dl(eq_start,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(eq_start+1,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(eq_start+1,a);
F(1,2) = F(1,2) + Nx(2,a)*dl(eq_start+1,a);
F(2,0) = F(2,0) + Nx(0,a)*dl(eq_start+2,a);
F(2,1) = F(2,1) + Nx(1,a)*dl(eq_start+2,a);
F(2,2) = F(2,2) + Nx(2,a)*dl(eq_start+2,a);
} else {
F(0,0) = F(0,0) + Nx(0,a)*dl(eq_start,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(eq_start,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(eq_start+1,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(eq_start+1,a);
}
}
return F;
}

}

void all_post(Simulation* simulation, Array<double>& res, const SolutionStates& solutions,
consts::OutputNameType outGrp, const int iEq)
{
Expand Down Expand Up @@ -443,17 +473,8 @@ void div_post(Simulation* simulation, const mshType& lM, Array<double>& res, con
vx(2,0) = vx(2,0) + Nx(0,a)*yl(k,a);
vx(2,1) = vx(2,1) + Nx(1,a)*yl(k,a);
vx(2,2) = vx(2,2) + Nx(2,a)*yl(k,a);

F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(0,2) = F(0,2) + Nx(2,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
F(1,2) = F(1,2) + Nx(2,a)*dl(j,a);
F(2,0) = F(2,0) + Nx(0,a)*dl(k,a);
F(2,1) = F(2,1) + Nx(1,a)*dl(k,a);
F(2,2) = F(2,2) + Nx(2,a)*dl(k,a);
}
F = deformation_gradient(Nx, dl, nsd, eNoN, i);

auto Fi = mat_fun::mat_inv(F,3);

Expand All @@ -468,12 +489,8 @@ void div_post(Simulation* simulation, const mshType& lM, Array<double>& res, con
vx(0,1) = vx(0,1) + Nx(1,a)*yl(i,a);
vx(1,0) = vx(1,0) + Nx(0,a)*yl(j,a);
vx(1,1) = vx(1,1) + Nx(1,a)*yl(j,a);

F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
}
F = deformation_gradient(Nx, dl, nsd, eNoN, i);

auto Fi = mat_fun::mat_inv(F,2);
VxFi(0) = vx(0,0)*Fi(0,0) + vx(0,1)*Fi(1,0);
Expand Down Expand Up @@ -562,31 +579,12 @@ void fib_algn_post(Simulation* simulation, const mshType& lM, Array<double>& res
double Jac = 0.0;

if (g == 0 || !lM.lShpF) {
auto Nx = lM.Nx.slice(g);
nn::gnn(eNoN, nsd, nsd, Nx, xl, Nx, Jac, F);
auto Nxi = lM.Nx.slice(g);
nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, F);
}

double w = lM.w(g)*Jac;
auto F = mat_fun::mat_id(nsd);

for (int a = 0; a < eNoN; a++) {
if (nsd == 3) {
F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(0,2) = F(0,2) + Nx(2,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
F(1,2) = F(1,2) + Nx(2,a)*dl(j,a);
F(2,0) = F(2,0) + Nx(0,a)*dl(k,a);
F(2,1) = F(2,1) + Nx(1,a)*dl(k,a);
F(2,2) = F(2,2) + Nx(2,a)*dl(k,a);
} else {
F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
}
}
auto F = deformation_gradient(Nx, dl, nsd, eNoN, i);
Comment thread
kko27 marked this conversation as resolved.
for (int iFn = 0; iFn < 2; iFn++) {
for (int i = 0; i < nsd; i++) {
auto fN_col = fN.col(iFn);
Expand Down Expand Up @@ -684,26 +682,7 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra

double w = lM.w(g) * Jac;
N = lM.N.col(g);
F = mat_fun::mat_id(nsd);

for (int a = 0; a < eNoN; a++) {
if (nsd == 3) {
F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(0,2) = F(0,2) + Nx(2,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
F(1,2) = F(1,2) + Nx(2,a)*dl(j,a);
F(2,0) = F(2,0) + Nx(0,a)*dl(k,a);
F(2,1) = F(2,1) + Nx(1,a)*dl(k,a);
F(2,2) = F(2,2) + Nx(2,a)*dl(k,a);
} else {
F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
}
}
F = deformation_gradient(Nx, dl, nsd, eNoN, i);

for (int iFn = 0; iFn < lM.nFn; iFn++) {
for (int i = 0; i < nsd; i++) {
Expand Down Expand Up @@ -740,45 +719,38 @@ void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Arra

}

/// @brief Compute fiber stretch based on 4th invariant: I_{4,f}
/// @brief Compute fiber stretch based on 4th invariant: λ = sqrt(I_{4,f})
//
void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector<double>& res)
void fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM,
const Array<double>& lD, Vector<double>& res)
{
using namespace consts;

auto& com_mod = simulation->com_mod;
auto& cm = com_mod.cm;
auto& cm_mod = simulation->cm_mod;
const auto& lD = solutions.old.get_displacement();
auto& eq = com_mod.eq[iEq];

int nsd = com_mod.nsd;
int tnNo = com_mod.tnNo;
int tDof = com_mod.tDof;

// [NOTE] Setting gobal variable 'dof'.
com_mod.dof = eq.dof;

int eNoN = lM.eNoN;
int i = eq.s;
int j = i + 1;
int k = j + 1;

Vector<double> sA(tnNo);
Vector<double> sF(tnNo);
Array<double> xl(nsd,eNoN);
Vector<double> sA(tnNo);
Vector<double> sF(tnNo);
Array<double> xl(nsd,eNoN);
Array<double> dl(tDof,eNoN);
Array<double> Nx(nsd,eNoN);
Vector<double> N(eNoN);
Array<double> Nx(nsd,eNoN);

for (int e = 0; e < lM.nEl; e++) {
int cDmn = all_fun::domain(com_mod, lM, iEq, e);
int cDmn = all_fun::domain(com_mod, lM, iEq, e);
auto cPhys = eq.dmn[cDmn].phys;
if (cPhys != EquationType::phys_struct && cPhys != EquationType::phys_ustruct) {
continue;
}

if (lM.eType == ElementType::NRB) {
//CALL NRBNNX(lM, e)
}

for (int a = 0; a < eNoN; a++) {
for (int a = 0; a < eNoN; a++) {
int Ac = lM.IEN(a,e);
xl.set_col(a, com_mod.x.col(Ac));
dl.set_col(a, lD.col(Ac));
Expand All @@ -788,40 +760,25 @@ void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const
double Jac = 0.0;
Array<double> F(nsd,nsd);
if (g == 0 || !lM.lShpF) {
auto Nx = lM.Nx.slice(g);
nn::gnn(eNoN, nsd, nsd, Nx, xl, Nx, Jac, F);
auto Nxi = lM.Nx.slice(g);
Array<double> dummy_ksix(nsd,nsd);
nn::gnn(eNoN, nsd, nsd, Nxi, xl, Nx, Jac, dummy_ksix);
}
Comment thread
kko27 marked this conversation as resolved.

double w = lM.w(g)*Jac;
auto N = lM.N.col(g);
F = mat_fun::mat_id(nsd);

for (int a = 0; a < eNoN; a++) {
if (nsd == 3) {
F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(0,2) = F(0,2) + Nx(2,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
F(1,2) = F(1,2) + Nx(2,a)*dl(j,a);
F(2,0) = F(2,0) + Nx(0,a)*dl(k,a);
F(2,1) = F(2,1) + Nx(1,a)*dl(k,a);
F(2,2) = F(2,2) + Nx(2,a)*dl(k,a);
} else {
F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
}
}
// Compute Deformation Gradient: F = I + grad(u)
F = deformation_gradient(Nx, dl, nsd, eNoN, i);

// Compute fiber stretch based on 4th invariant: I_{4,f} = F.fN.F.fN
auto fl = mat_fun::mat_mul(F, lM.fN.rows(0,nsd-1,e));
double I4f = utils::norm(fl);
double lambda = sqrt(utils::norm(fl));
Comment thread
michelebucelli marked this conversation as resolved.

for (int a = 0; a < eNoN; a++) {
// L2 projection from integration points to nodes
double w = lM.w(g)*Jac;
auto N = lM.N.col(g);
for (int a = 0; a < eNoN; a++) {
int Ac = lM.IEN(a,e);
sA(Ac) = sA(Ac) + w*N(a);
sF(Ac) = sF(Ac) + w*N(a)*I4f;
sF(Ac) = sF(Ac) + w*N(a)*lambda;
}
}
}
Expand All @@ -837,7 +794,33 @@ void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const
res(a) = res(a) + sF(Ac) / sA(Ac);
}
}
}

/// @brief Compute fiber stretch rate dλ/dt via backward finite difference.
//
void fib_stretch_rate(const ComMod& com_mod, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector<double>& res)
{
const double dt = com_mod.dt;
int nNo = lM.nNo;

if (dt <= 0.0) {
svmp::raise<svmp::FE::InvalidArgumentException>(
SVMP_HERE,
"[fib_stretch_rate] Expected com_mod.dt > 0, but got " + std::to_string(dt) + ".");
}

if (res.size() != nNo) {
svmp::raise<svmp::FE::InvalidArgumentException>(
SVMP_HERE,
"[fib_stretch_rate] Expected res size " + std::to_string(nNo) + ", but got " + std::to_string(res.size()) + ".");
}

Vector<double> lambda_old(nNo);

fib_stretch(com_mod, iEq, lM, solutions.current.get_displacement(), res);
fib_stretch(com_mod, iEq, lM, solutions.old.get_displacement(), lambda_old);

res = (res - lambda_old) / dt;
}

void post(Simulation* simulation, const mshType& lM, Array<double>& res, const SolutionStates& solutions,
Expand Down Expand Up @@ -1803,26 +1786,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array<double>
Je = Je + w;

auto Im = mat_fun::mat_id(nsd);
auto F = Im;

for (int a = 0; a < fs.eNoN; a++) {
if (nsd == 3) {
F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(0,2) = F(0,2) + Nx(2,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
F(1,2) = F(1,2) + Nx(2,a)*dl(j,a);
F(2,0) = F(2,0) + Nx(0,a)*dl(k,a);
F(2,1) = F(2,1) + Nx(1,a)*dl(k,a);
F(2,2) = F(2,2) + Nx(2,a)*dl(k,a);
} else {
F(0,0) = F(0,0) + Nx(0,a)*dl(i,a);
F(0,1) = F(0,1) + Nx(1,a)*dl(i,a);
F(1,0) = F(1,0) + Nx(0,a)*dl(j,a);
F(1,1) = F(1,1) + Nx(1,a)*dl(j,a);
}
}
auto F = deformation_gradient(Nx, dl, nsd, fs.eNoN, i);

double detF = mat_fun::mat_det(F, nsd);

Expand Down
5 changes: 3 additions & 2 deletions Code/Source/solver/post.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,9 @@ void fib_algn_post(Simulation* simulation, const mshType& lM, Array<double>& res

void fib_dir_post(Simulation* simulation, const mshType& lM, const int nFn, Array<double>& res, const SolutionStates& solutions, const int iEq);

void fib_strech(Simulation* simulation, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector<double>& res);
void fib_stretch(const ComMod& com_mod, const int iEq, const mshType& lM, const Array<double>& lD, Vector<double>& res);

void fib_stretch_rate(const ComMod& com_mod, const int iEq, const mshType& lM, const SolutionStates& solutions, Vector<double>& res);

void post(Simulation* simulation, const mshType& lM, Array<double>& res, const SolutionStates& solutions,
consts::OutputNameType outGrp, const int iEq);
Expand All @@ -38,4 +40,3 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array<double>
};

#endif

Loading
Loading