(HOWTO) Insertion of a new fiber/contraction model

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts
  • HeikoStark
    Junior Member
    • Mar 2009
    • 26

    (HOWTO) Insertion of a new fiber/contraction model

    A little framework for insertion of new fiber/contraction models in FEBio. It implements in all materials with active_contraction the new models. Feel free to discuss.
    1. Get the new attribute values in FEBioImport.cpp (line ~760)

      Code:
      else if (tag == "active_contraction")
      {
      	const char* szlc = tag.AttributeValue("lc", true);
      	int lc = 0;
      	if (szlc) lc = atoi(szlc);
      	pm->m_fib.m_lcna = lc;
      	tag.value(pm->m_fib.m_ascl);
      	[B]pm->m_fib.m_model = 1;[/B]
      
      ...
      
      }
      else if [B](tag == "active_contraction_new")[/B]
      {
      	const char* szlc = tag.AttributeValue("lc", true);
      	int lc = 0;
      	if (szlc) lc = atoi(szlc);
      	pm->m_fib.m_lcna = lc;
      	tag.value(pm->m_fib.m_ascl);
      	[B]pm->m_fib.m_model = 2;[/B]
      
      ... [B]new tag values[/B] ...
      
      }
    2. Define the new parameters in FEMaterial.h:

      Code:
      class FEFiberMaterial : public FEMaterial
      {
      public:
      	FEFiberMaterial()
      	{
      		m_plc = 0;
      		m_lcna = -1;
      		m_ascl = 0;
      		[B]m_model = 0;	[/B]
      
      ...
      
      	double	m_ascl;		//!< activation scale factor
      [B]	int  	m_model;	//!< model[/B]
      // Model 1	
      	double	m_ca0;		//!< intracellular calcium concentration
      	double	m_beta;		//!< shape of peak isometric tension-sarcomere length relation
      	double	m_l0;		//!< unloaded length
      	double	m_refl;	 	//!< sarcomere length
      // Model 2
      	
      ... [B]new attributes[/B] ...
      
      // we need load curve data for active contraction
      	FELoadCurve* m_plc;	//!< pointer to current load curve values
      };
    3. Work with the new model in FEFibermaterial.cpp

      Code:
      mat3ds FEFiberMaterial::Stress(FEMaterialPoint &mp)
      {
      	FEElasticMaterialPoint& pt = *mp.ExtractData<FEElasticMaterialPoint>();
      
      ...
      	
        	mat3ds s;
      	
      	if [B](m_model == 1)[/B]
      	{	
      
      ...
      
      	}
      	else if [B](m_model == 2)[/B]
      	{
      
      ... [B]the new model[/B] ...
      	
      	}
      	else
      	{
      		s.zero ();
      	}
      	return s;
      }


    Cheers!
    Last edited by HeikoStark; 01-03-2011, 06:52 AM.
  • maas
    Lead Code Developer
    • Nov 2007
    • 3721

    #2
    Hi Heiko,

    Thanks for sharing this. I hope it was not too hard to dig through the source code to figure this out. I know we don't have that much information for developers yet.

    I do have a few remarks about your implementation. First, I recently made some changes to the FEBioImport file. You probably want to grab the newest version of the web (posted on dec 9 2010) to make sure your modifications still works. Second, In future versions of FEBio I envision all material implementations to have their own separate class. I know that for the implementation of fiber models at this point that might be a little too cumbersome, but this is most likely the road we will take. Please keep this in mind.

    I you have any questions about implementation details or suggestions on how the code should evolve, please feel free to share them in the FEBio developers forum.

    Thanks,

    Steve.
    Department of Bioengineering, University of Utah
    Scientific Computing and Imaging institute, University of Utah

    Comment

    • HeikoStark
      Junior Member
      • Mar 2009
      • 26

      #3
      Hi!

      Thanks for the newest version!

      I think my implementation is only a interim solution. The best way would be to split up the strain-energy function in a isotropic, anisotropic and volume part. And to split up the anisotropic part in a passive and active fiber contribution. With this it should be possible to combine different models.

      Isotropic
      Neo-Hookean
      Mooney-Rivlin
      Anisotropic (One or two fibers)
      Passive
      Different tissues
      ...
      Active
      Static or dynamic
      Fibertypes
      ...
      Volume

      W = WIso + WAni + WVol
      WAni = Wpassive + Wactive

      But that is much work.

      Cheers.
      Last edited by HeikoStark; 12-15-2010, 06:10 AM.

      Comment

      Working...
      X