Example 5d4: Active Contraction based on a user-defined CellML function

This examples demonstrates how a user-specific function can be defined within CellML and used to simulate actice contraction. Here, a constant Calcium level is asumed. The user-specific function is used to enhance the 2nd Piola Kirchhoff stress tensor. The 'ipacti' file demonstrates how one chooses a CellMl-defined function, what type of stress this user-defined function resembles, and in which direction it should be added to. This example is based on a modified Mooney Rivlin constitutive law which can be enhanced by a passive and active stress response that depends on the fibrestretch. The example is twofold. First, we just elongate a unit cube in one direction by .1 (10%) and apply only a passive stress ($alpha=0). Then we keep the node values of the cube fixed and increase the activation level to a $alpha=0.2. This implies that we add to the 2nd PKST a passive and active stress component. This is explained in more detail in Rohrle o., Pullan A., "Three-dimensional finite element analysis of muscle forces during mastication", 2006 (submitted)

The comfile run by this example is as follows:

#  
#  EXAMPLE 5d4: Active Contraction based on a user-defined CellML function
#
#	This examples demonstrates how a user-specific function can be defined
#	within CellML and used to simulate actice contraction. Here, a constant
#	Calcium level is asumed. The user-specific function is used to enhance
#	the 2nd Piola Kirchhoff stress tensor. The 'ipacti' file demonstrates
#	how one chooses a CellMl-defined function, what type of stress this 
#	user-defined function resembles, and in which direction it should be
#	added to.
#
#	This example is based on a modified Mooney Rivlin constitutive law which
# 	can be enhanced by a passive and active stress response that depends on
#	the fibrestretch. The example is twofold. First, we just elongate a unit
#	cube in one direction by .1 (10%) and apply only a passive stress ($alpha=0).
#	Then we keep the node values of the cube fixed and increase the activation
#	level to a $alpha=0.2. This implies that we add to the 2nd PKST a passive
#	and active stress component.
#
#	This is explained in more detail in Rohrle o., Pullan A., "Three-dimensional 
#	finite element analysis of muscle forces during mastication", 2006 (submitted)
#


# If the example path is not set, default to current directory
if (! defined $example) {
    $example = "./";
}
# Drop off the trailing / in the example path
$chopped = chop $example;
if ($chopped ne "/") {
    $example .= $chopped;
}

# Set the output directory
$output = ".";
if (! -d $output) {
    mkdir $output,0700;
}


# Definition of basic variables

FEM define para ;r; $example/RunParameters
FEM define coordinate 3,1
FEM define base ;r; $example/Activation

# Definition of a unit cube
FEM define node ;r; $example/UnitCube
FEM define elem ;r; $example/UnitCube
FEM define fibr ;r; $example/UnitCube30
FEM define elem ;r; $example/UnitCube      fibre
FEM define fiel ;r; $example/UnitCube
FEM define elem ;r; $example/UnitCube      field

FEM group elem 1 as ALL_ELEMENTS

#
# 	Setup the mechanics problem, in particular that the constitutive law is referred

#	
FEM define equa ;r; $example/Activation
FEM define mate ;r; $example/Activation_cellml

#
#	Write out reference state for visualisation (View_Results.com)
#
FEM update solution from geometry
FEM UPDATE field 1 sub sol depvar 1
FEM UPDATE field 2 sub sol depvar 2
FEM UPDATE field 3 sub sol depvar 3
FEM export elem; $output/UnitCube as cube 
FEM export node; $output/UnitCube as cube

#
#	Define which nodal values are Dirichlet (essential) boundary conditions
#	and can be changed later on (see below 'FEM change disp ...' command)
#
FEM define init ;r; $example/Activation

#
#	Define grid for CellML and choose grid such that it coincides with
#	the gauss points of the elements
#
FEM define grid ;r; $example/Activation   gauss class 2
FEM update grid geometry
FEM group  grid element ALL_ELEMENTS as ALL_GRID_POINTS class 2

#
#	Define everything for CellML use
#	Note that the function fibrestress in MooneyRivlin+FibreStretch.xml is only
#	defined and not used in any calculations of the stress tensor. This is the
#	total (Cauchy) stress value which is added to the 2nd PKST and represents the
#	active contraction part. 

FEM define equa ;r; $example/Activation_cellml        class 2
FEM define mate ;d                                    class 2
FEM define cell ;r; $example/Activation_cellml        class 2
FEM define mate ;r; $example/Activation_cellml        cell class 2
FEM define init ;d                                    class 2
FEM define solv ;r; $example/Activation_cellml        class 2


#
#	By choosing option (3) in at the type of contraction model, we can subsequently
#	choose the type of stress (Cauchy or 2nd PKST-type), to which component of the
#	stress tensor you would like to add this stress (in body coordinates, 1,1 
#	represents the fibre direction, and the name of the CellML function which evaluates
#	the stress component. In this example, the function fibrestress.
#

FEM define acti ;r; $example/Activation_cellml


#
# 	The different load steps should demonstrate the different states, eg. only passive
#	or active and passive stress.
#

$LoadStepStart = 1
$MAX_LoadSteps = 3

for ( $LoadStep=$LoadStepStart; $LoadStep<$MAX_LoadSteps; $LoadStep++ )
{
   if ($LoadStep==1) {
      $disp  = 0.1	# displacement parameter
      $alpha = 0.0	# activation parameter
   }	     
   if ($LoadStep==2) {
      $disp  = 0.0	# displacement parameter
      $alpha = 0.2      # activation parameter
   }

   print "\n\ndisplacement: ",$disp,"   and alpha: ",$alpha,"\n\n"

#
#	Looks up in which cell array the variable 'alpha' is stored. 'Alpha' can be 
#	understood as a parameter which represents the percentage of active muscle fibres
#	within a muscle.
#
   FEM inquire cell_variable alpha return_variables v1,v2
   print "\n\n",$v1,"   ",$v2,"\n\n"

#
#	Update the CellML variable 'alpha' with a constant value
#
   FEM update grid $v1 $v2 const $alpha

#
#	Initialize all Dirichlet bc to 0.
#	This sets the boundary conditions which are equivalant to a translation of the 
#	nodes 2,4,6,8 (right side of the unit cube) by 0.1 (or 0.0 in the second load step) 
#
   FEM change disp init  
   FEM change disp translate nodes 2,4,6,8  by  $disp,0,0

#
#	Solve the problem
#
   FEM define solv;r; $example/Activation

   FEM solve iter 20

#
#	Compute the reaction forces necessary to obtain the deformed state, and
#	save the deformed state.
#
   fem define fit;r; $example/Activation   grid
   fem fit

   FEM UPDATE field 1 sub sol depvar 1
   FEM UPDATE field 2 sub sol depvar 2
   FEM UPDATE field 3 sub sol depvar 3

   fem eval reac nh 1 to 4 noview
   fem eval reac nh 2 to 5 noview
   fem eval reac nh 3 to 6 noview

   FEM export nodes; $output/DefUnitCube_${LoadStep} as cube

}

#
#	The results can be visualised by the command View_Result.com
#	The time slider allows one to look at the different load steps.
#

Files used by this example are:

Name                           Modified     Size


example_5d4.com                23-Aug-2006  5.8k
Activation.ipbase              23-Aug-2006   12k
Activation.ipequa              23-Aug-2006  2.0k
Activation.ipfit               13-Apr-2007  1.6k
Activation.ipgrid              23-Aug-2006  654 
Activation.ipinit              23-Aug-2006  4.5k
Activation.ipmate              23-Aug-2006  739 
Activation.ippara              23-Aug-2006  6.0k
Activation.ipsolv              16-Aug-2010  2.4k
Activation.ipsolv.old          13-Apr-2007  2.3k
Activation_cellml.ipacti       23-Aug-2006  605 
Activation_cellml.ipcell       23-Aug-2006  3.3k
Activation_cellml.ipequa       23-Aug-2006  1.6k
Activation_cellml.ipmatc       23-Aug-2006  3.0k
Activation_cellml.ipmate       23-Aug-2006  1.5k
Activation_cellml.ipsolv       13-Apr-2007  2.3k
MooneyRivlin+FibreStretch.xml  23-Aug-2006   19k
RunParameters.ippara           23-Aug-2006  5.9k
UnitCube.ipelem                23-Aug-2006  639 
UnitCube.ipelfb                23-Aug-2006  388 
UnitCube.ipelfd                23-Aug-2006  721 
UnitCube.ipfiel                23-Aug-2006   38k
UnitCube.ipnode                23-Aug-2006   15k
UnitCube30.ipfibr              23-Aug-2006   12k
View_Results.com               23-Aug-2006  3.8k

Download the entire example:

Name                      Modified     Size


examples_5_5d_5d4.tar.gz  17-Aug-2010   36k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmFailureSun Mar 6 00:16:42 20165
last breakTue Feb 24 03:22:00 20154
last successSun Apr 17 00:54:00 201110
cm-debugFailureSun Mar 6 00:40:30 201625
last breakTue Feb 24 03:40:00 201525
last successSat Apr 16 00:50:00 201139
mips-irix
cmSuccessSun Aug 19 02:16:49 200745
cm-debugSuccessWed Aug 15 02:54:58 2007228
cm-debug-clear-mallocSuccessSat Aug 18 03:23:32 2007259
cm-debug-clear-malloc7FailureTue Aug 21 05:58:55 200728
last breakMon Aug 28 01:39:00 2006102
cm64SuccessSun Aug 19 02:18:28 200747
cm64-debugSuccessTue Aug 21 03:35:16 2007256
rs6000-aix
cmSuccessWed Mar 4 01:12:59 20096
cm-debugSuccessMon Mar 2 01:21:12 200985
cm64SuccessWed Mar 4 01:12:25 20098
cm64-debugSuccessTue Mar 3 01:32:17 200980
x86_64-linux
cmFailureSun Mar 6 00:01:54 20163
last breakSun May 8 00:27:00 20114
last successSun Apr 17 00:27:00 20114
cm-debugFailureSun Mar 6 00:03:35 201612
last breakMon Apr 18 00:11:00 201116
last successSat Apr 16 00:20:00 201116

Testing status by file:

Html last generated: Sun Mar 6 05:50:27 2016

Input last modified: Mon Aug 16 11:19:15 2010

CMISS Help / Examples / 5 / 5d / 5d4