Dear all,
I think I caught a bug in FEBio source code. Specifically, it is the calculation of tangent modulus for Uncoupled Ogden phenomenological model.
In febiosource-2.2.6\FEBioMech\FEOgdenMaterial.cpp,
Error #1: line 411 - 412:
beta1 += m_c[l]/m_m[l]*(lamd[i][l] - (lamd[0][l] + lamd[1][l] + lamd[2][l])/3.0);
beta3 += m_c[l]/m_m[l]*(lamd[j][l] - (lamd[0][l] + lamd[1][l] + lamd[2][l])/3.0);
Should be the same with line 150-151 (I swapped the order of two lines though):
beta3 += m_c[l]/m_m[l]*(lamd[i][l] - (lamd[0][l] + lamd[1][l] + lamd[2][l])/3.0);
beta1 += m_c[l]/m_m[l]*(lamd[j][l] - (lamd[0][l] + lamd[1][l] + lamd[2][l])/3.0);
Error #2: line 416 - 417:
for (l=0; l<MAX_TERMS; ++l) g11 += m_c[l]*(lamd[i][l]/3.0 + (lamd[0][l]+lamd[1][l]+lamd[2][l])/9.0);
for (l=0; l<MAX_TERMS; ++l) g33 += m_c[l]*(lamd[j][l]/3.0 + (lamd[0][l]+lamd[1][l]+lamd[2][l])/9.0);
Should be:
for (l=0; l<MAX_TERMS; ++l) g33 += m_c[l]*(lamd[i][l]/3.0 + (lamd[0][l]+lamd[1][l]+lamd[2][l])/9.0);
for (l=0; l<MAX_TERMS; ++l) g11 += m_c[l]*(lamd[j][l]/3.0 + (lamd[0][l]+lamd[1][l]+lamd[2][l])/9.0);
The following are the rationale.
I think it is obvious that all code comes from the paper "Simo JC, Taylor RL. Quasi-incompressible finite elasticity in principal stretches. continuum basis and numerical algorithms. Comput Methods Appl Mech Eng 85: 273?310, 1991."
The problematic code is in implementing the equation (2.70). The paper was denoting the distinct eigenvalue as #3, and the double eigenvalues as #1 and #2. In the code at line 395 - 397, it is pretty clear that the code conforms with the paper; however at the problematic lines reported above, the #1 was confused as the distinct eigenvalue, whereas it should be the doubled eigenvalue. Vice versa for the #3.
PS: Sorry if this is the wrong place to post this - I assume FEBio does not have a github organization?
Shawn
I think I caught a bug in FEBio source code. Specifically, it is the calculation of tangent modulus for Uncoupled Ogden phenomenological model.
In febiosource-2.2.6\FEBioMech\FEOgdenMaterial.cpp,
Error #1: line 411 - 412:
beta1 += m_c[l]/m_m[l]*(lamd[i][l] - (lamd[0][l] + lamd[1][l] + lamd[2][l])/3.0);
beta3 += m_c[l]/m_m[l]*(lamd[j][l] - (lamd[0][l] + lamd[1][l] + lamd[2][l])/3.0);
Should be the same with line 150-151 (I swapped the order of two lines though):
beta3 += m_c[l]/m_m[l]*(lamd[i][l] - (lamd[0][l] + lamd[1][l] + lamd[2][l])/3.0);
beta1 += m_c[l]/m_m[l]*(lamd[j][l] - (lamd[0][l] + lamd[1][l] + lamd[2][l])/3.0);
Error #2: line 416 - 417:
for (l=0; l<MAX_TERMS; ++l) g11 += m_c[l]*(lamd[i][l]/3.0 + (lamd[0][l]+lamd[1][l]+lamd[2][l])/9.0);
for (l=0; l<MAX_TERMS; ++l) g33 += m_c[l]*(lamd[j][l]/3.0 + (lamd[0][l]+lamd[1][l]+lamd[2][l])/9.0);
Should be:
for (l=0; l<MAX_TERMS; ++l) g33 += m_c[l]*(lamd[i][l]/3.0 + (lamd[0][l]+lamd[1][l]+lamd[2][l])/9.0);
for (l=0; l<MAX_TERMS; ++l) g11 += m_c[l]*(lamd[j][l]/3.0 + (lamd[0][l]+lamd[1][l]+lamd[2][l])/9.0);
The following are the rationale.
I think it is obvious that all code comes from the paper "Simo JC, Taylor RL. Quasi-incompressible finite elasticity in principal stretches. continuum basis and numerical algorithms. Comput Methods Appl Mech Eng 85: 273?310, 1991."
The problematic code is in implementing the equation (2.70). The paper was denoting the distinct eigenvalue as #3, and the double eigenvalues as #1 and #2. In the code at line 395 - 397, it is pretty clear that the code conforms with the paper; however at the problematic lines reported above, the #1 was confused as the distinct eigenvalue, whereas it should be the doubled eigenvalue. Vice versa for the #3.
PS: Sorry if this is the wrong place to post this - I assume FEBio does not have a github organization?
Shawn
Comment