(HOWTO) Insertion of a new fiber/contraction model

[HeikoStark]

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)
   
   

   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:
   
   

   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
   
   

   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!

[maas]

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.

[HeikoStark]

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.