A bug found in FEBio source code (likely a typo)

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts
  • yw5aj
    Junior Member
    • Apr 2013
    • 13

    A bug found in FEBio source code (likely a typo)

    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
  • yw5aj
    Junior Member
    • Apr 2013
    • 13

    #2
    Hi,

    Would anyone be interested to discuss/fix this one? Thanks,

    Shawn

    Comment

    • ateshian
      Developer
      • Dec 2007
      • 1824

      #3
      Hi Shawn,

      Sorry for the delay in responding to this issue. Since the test suite problems that use Ogden materials were running fine, this was not treated as an urgent concern, but we did not forget it.

      I checked the validity of the tangent calculation in the existing code using the tangent diagnostic tool with this file: tangent_OG.feb (febio2 -d tangent_OG.feb). It confirmed that the tangent calculation was incorrect. I implemented the corrections you mentioned in your post and the results improved, but the tangent diagnostics still failed overall. This meant that there were additional errors to the ones you found.

      To fix the problem, it was simpler to rewrite the whole routine from scratch. I copied the algorithm for stress and tangent calculations from the compressible "Ogden unconstrained" material and modified them to apply to the "Ogden" uncoupled formulation. This fixed the problem (the tangent diagnostic passed successfully). I checked the problems in the test suite (co28.feb, co29.feb, co37.feb, ma01.feb, ma16.feb) and one of the problems showed significant improvements with the fixed algorithm (co37.feb, which went down from 478 to 359 equilibrium iterations), whereas the others ran just as well or only slightly better (thankfully...).

      This fix will be available in the next FEBio release. Thanks again for going through the code and checking it. This was very helpful.

      Best,

      Gerard

      Comment

      Working...
      X
      😀
      😂
      🥰
      😘
      🤢
      😎
      😞
      😡
      👍
      👎