Skip to content

Commit

Permalink
Modified implementation of cubic-CLE and ortho-CLE materials to hew m…
Browse files Browse the repository at this point in the history 
…ore closely to theoretical formulation
  • Loading branch information
@gateshian
gateshian committed Oct 9, 2024
1 parent ce496e4 commit 01b432a
Show file tree
Hide file tree
Showing 3 changed files with 66 additions and 64 deletions.
4 changes: 2 additions & 2 deletions Documentation/FEBio_User_Manual.lyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63371,7 +63371,7 @@ where

Here,

\begin_inset Formula ${\mathbf{E}\ensuremath{}}$
\begin_inset Formula $\mathbf{E}$
\end_inset

is the Lagrangian strain tensor and
Expand Down Expand Up @@ -75751,7 +75751,7 @@ Number of Taylor series terms
\begin_inset Text

\begin_layout Plain Layout
Shear modulus fir fiber-matrix interaction
Shear modulus for fiber-matrix interaction
\end_layout

\end_inset
Expand Down
63 changes: 32 additions & 31 deletions FEBioMech/FECubicCLE.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,9 +102,9 @@ mat3ds FECubicCLE::Stress(FEMaterialPoint& mp)

int i,j;
vec3d a0; // texture direction in reference configuration
mat3ds A[3]; // texture tensor in current configuration
double K[3]; // Ka
mat3ds BmI = pt.LeftCauchyGreen() - mat3dd(1); // B - I
mat3ds A0[3]; // texture tensor in reference configuration
double AE[3]; // E:A
mat3ds E = pt.Strain();

mat3d F = pt.m_F;
double J = pt.m_J;
Expand All @@ -115,23 +115,23 @@ mat3ds FECubicCLE::Stress(FEMaterialPoint& mp)
for (i=0; i<3; i++) { // Perform sum over all three texture directions
// Copy the texture direction in the reference configuration to a0
a0.x = Q[0][i]; a0.y = Q[1][i]; a0.z = Q[2][i];
A[i] = dyad(F*a0);
K[i] = 0.5*(A[i].tr() - 1);
A0[i] = dyad(a0);
AE[i] = a0*(E*a0);
}

lam[0][0] = (K[0] >= 0) ? m_lp1 : m_lm1;
lam[1][1] = (K[1] >= 0) ? m_lp1 : m_lm1;
lam[2][2] = (K[2] >= 0) ? m_lp1 : m_lm1;
lam[0][0] = (AE[0] >= 0) ? m_lp1 : m_lm1;
lam[1][1] = (AE[1] >= 0) ? m_lp1 : m_lm1;
lam[2][2] = (AE[2] >= 0) ? m_lp1 : m_lm1;

mat3ds s; s.zero();

for (i=0; i<3; ++i) {
s += (A[i]*BmI).sym()*(mu[i]);
s += (A0[i]*E).sym()*(mu[i]);
for (j=0; j<3; ++j) {
s += A[j]*(lam[i][j]*K[i]);
s += A0[j]*(lam[i][j]*AE[i]);
}
}
s /= J;
s = (F*s*F.transpose()).sym()/J;

return s;
}
Expand All @@ -154,10 +154,10 @@ tens4ds FECubicCLE::Tangent(FEMaterialPoint& mp)

int i,j;
vec3d a0; // texture direction in reference configuration
mat3ds A[3]; // texture tensor in current configuration
double K[3]; // Ka
mat3ds B = pt.LeftCauchyGreen();
mat3ds A[3]; // texture tensor in current configuration
double AE[3]; // E:A
mat3ds E = pt.Strain();

mat3d F = pt.m_F;
double J = pt.m_J;

Expand All @@ -168,16 +168,17 @@ tens4ds FECubicCLE::Tangent(FEMaterialPoint& mp)
// Copy the texture direction in the reference configuration to a0
a0.x = Q[0][i]; a0.y = Q[1][i]; a0.z = Q[2][i];
A[i] = dyad(F*a0);
K[i] = 0.5*(A[i].tr() - 1);
AE[i] = a0*(E*a0);
}

lam[0][0] = (K[0] >= 0) ? m_lp1 : m_lm1;
lam[1][1] = (K[1] >= 0) ? m_lp1 : m_lm1;
lam[2][2] = (K[2] >= 0) ? m_lp1 : m_lm1;
lam[0][0] = (AE[0] >= 0) ? m_lp1 : m_lm1;
lam[1][1] = (AE[1] >= 0) ? m_lp1 : m_lm1;
lam[2][2] = (AE[2] >= 0) ? m_lp1 : m_lm1;

tens4ds c;
c.zero();


mat3ds B = pt.LeftCauchyGreen();
for (i=0; i<3; ++i) {
c += dyad4s(A[i], B)*mu[i];
for (j=0; j<3; ++j) {
Expand Down Expand Up @@ -207,9 +208,9 @@ double FECubicCLE::StrainEnergyDensity(FEMaterialPoint& mp)

int i,j;
vec3d a0; // texture direction in reference configuration
mat3ds A0; // texture tensor in current configuration
double K[3], L[3]; // Ka
mat3ds E = (pt.RightCauchyGreen() - mat3dd(1))/2;
mat3ds A0; // texture tensor in reference configuration
double AE[3], AE2[3]; // Ka
mat3ds E = pt.Strain();

// get the local coordinate systems
mat3d Q = GetLocalCS(mp);
Expand All @@ -218,20 +219,20 @@ double FECubicCLE::StrainEnergyDensity(FEMaterialPoint& mp)
// Copy the texture direction in the reference configuration to a0
a0.x = Q[0][i]; a0.y = Q[1][i]; a0.z = Q[2][i];
A0 = dyad(a0);
K[i] = (A0*E).trace();
L[i] = (A0*E*E).trace();
AE[i] = a0*(E*a0);
AE2[i] = a0*(E*E*a0);
}

lam[0][0] = (K[0] >= 0) ? m_lp1 : m_lm1;
lam[1][1] = (K[1] >= 0) ? m_lp1 : m_lm1;
lam[2][2] = (K[2] >= 0) ? m_lp1 : m_lm1;
lam[0][0] = (AE[0] >= 0) ? m_lp1 : m_lm1;
lam[1][1] = (AE[1] >= 0) ? m_lp1 : m_lm1;
lam[2][2] = (AE[2] >= 0) ? m_lp1 : m_lm1;

double W = 0;

for (i=0; i<3; ++i) {
W += L[i]*mu[i];
W += AE2[i]*mu[i];
for (j=0; j<3; ++j) {
W += K[i]*K[j]*lam[i][j]/2;
W += AE[i]*AE[j]*lam[i][j]/2;
}
}

Expand Down
63 changes: 32 additions & 31 deletions FEBioMech/FEOrthotropicCLE.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,10 +108,10 @@ mat3ds FEOrthotropicCLE::Stress(FEMaterialPoint& mp)

int i,j;
vec3d a0; // texture direction in reference configuration
mat3ds A[3]; // texture tensor in current configuration
double K[3]; // Ka
mat3ds BmI = pt.LeftCauchyGreen() - mat3dd(1); // B - I
mat3ds A0[3]; // texture tensor in reference configuration
double AE[3]; // E:A
mat3ds E = pt.Strain();

mat3d F = pt.m_F;
double J = pt.m_J;

Expand All @@ -121,25 +121,25 @@ mat3ds FEOrthotropicCLE::Stress(FEMaterialPoint& mp)
for (i=0; i<3; i++) { // Perform sum over all three texture directions
// Copy the texture direction in the reference configuration to a0
a0.x = Q[0][i]; a0.y = Q[1][i]; a0.z = Q[2][i];
A[i] = dyad(F*a0);
K[i] = 0.5*(A[i].tr() - 1);
A0[i] = dyad(a0);
AE[i] = a0*(E*a0);
}

lam[0][0] = (K[0] >= 0) ? lp11 : lm11;
lam[1][1] = (K[1] >= 0) ? lp22 : lm22;
lam[2][2] = (K[2] >= 0) ? lp33 : lm33;
lam[0][0] = (AE[0] >= 0) ? lp11 : lm11;
lam[1][1] = (AE[1] >= 0) ? lp22 : lm22;
lam[2][2] = (AE[2] >= 0) ? lp33 : lm33;

mat3ds s;
s.zero();

for (i=0; i<3; ++i) {
s += (A[i]*BmI).sym()*(mu[i]);
s += (A0[i]*E).sym()*(mu[i]);
for (j=0; j<3; ++j) {
s += A[j]*(lam[i][j]*K[i]);
s += A0[j]*(lam[i][j]*AE[i]);
}
}
s /= J;
s = (F*s*F.transpose()).sym()/J;

return s;
}

Expand All @@ -162,9 +162,9 @@ tens4ds FEOrthotropicCLE::Tangent(FEMaterialPoint& mp)
int i,j;
vec3d a0; // texture direction in reference configuration
mat3ds A[3]; // texture tensor in current configuration
double K[3]; // Ka
mat3ds B = pt.LeftCauchyGreen();
double AE[3]; // E:A
mat3ds E = pt.Strain();

mat3d F = pt.m_F;
double J = pt.m_J;

Expand All @@ -175,16 +175,17 @@ tens4ds FEOrthotropicCLE::Tangent(FEMaterialPoint& mp)
// Copy the texture direction in the reference configuration to a0
a0.x = Q[0][i]; a0.y = Q[1][i]; a0.z = Q[2][i];
A[i] = dyad(F*a0);
K[i] = 0.5*(A[i].tr() - 1);
AE[i] = a0*(E*a0);
}

lam[0][0] = (K[0] >= 0) ? lp11 : lm11;
lam[1][1] = (K[1] >= 0) ? lp22 : lm22;
lam[2][2] = (K[2] >= 0) ? lp33 : lm33;
lam[0][0] = (AE[0] >= 0) ? lp11 : lm11;
lam[1][1] = (AE[1] >= 0) ? lp22 : lm22;
lam[2][2] = (AE[2] >= 0) ? lp33 : lm33;

tens4ds c;
c.zero();

mat3ds B = pt.LeftCauchyGreen();
for (i=0; i<3; ++i) {
c += dyad4s(A[i], B)*mu[i];
for (j=0; j<3; ++j) {
Expand Down Expand Up @@ -214,31 +215,31 @@ double FEOrthotropicCLE::StrainEnergyDensity(FEMaterialPoint& mp)

int i,j;
vec3d a0; // texture direction in reference configuration
mat3ds A0; // texture tensor in current configuration
double K[3], L[3]; // Ka
mat3ds E = (pt.RightCauchyGreen() - mat3dd(1))/2;
mat3ds A0; // texture tensor in reference configuration
double AE[3], AE2[3]; // Ka
mat3ds E = pt.Strain();

// get local coordinates
mat3d Q = GetLocalCS(mp);

for (i=0; i<3; i++) { // Perform sum over all three texture directions
// Copy the texture direction in the reference configuration to a0
a0.x = Q[0][i]; a0.y = Q[1][i]; a0.z = Q[2][i];
A0 = dyad(a0);
K[i] = (A0*E).trace();
L[i] = (A0*E*E).trace();
AE[i] = a0*(E*a0);
AE2[i] = a0*(E*E*a0);
}

lam[0][0] = (K[0] >= 0) ? lp11 : lm11;
lam[1][1] = (K[1] >= 0) ? lp22 : lm22;
lam[2][2] = (K[2] >= 0) ? lp33 : lm33;
lam[0][0] = (AE[0] >= 0) ? lp11 : lm11;
lam[1][1] = (AE[1] >= 0) ? lp22 : lm22;
lam[2][2] = (AE[2] >= 0) ? lp33 : lm33;

double W = 0;

for (i=0; i<3; ++i) {
W += L[i]*mu[i];
W += AE2[i]*mu[i];
for (j=0; j<3; ++j) {
W += K[i]*K[j]*lam[i][j]/2;
W += AE[i]*AE[j]*lam[i][j]/2;
}
}

Expand Down

0 comments on commit 01b432a

Please sign in to comment.