Example 5d5: Steady state tension-length-Ca relation in CellML

This example duplicates example 5d2 but uses CellML to define the active tension model. Currently this example is for demonstrations and development purposes in order to share the concepts being implemented in order to solve these types of simulations.

There are a few restrictions currently:


The comfile run by this example is as follows:

# Example 5d5

set echo;

# 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 directory to store results
if (! defined $output) {
  $output = ".";
}
if (! -d $output) {
	mkdir $output,0700;
}

$MECH_MAX_ITERS = 20;
$MECH_ERRTOL = 0.001;

sub updateGridPoints {
	my $cellClass = shift;
	my $tissueClass = shift;
	my $lMechIdx = shift;
	my $lGridArray = shift;
	my $lGridIdx = shift;

	fem update grid geometry deformed class $cellClass from_class $tissueClass;
	fem update grid metric deformed class $cellClass from_class $tissueClass;
	fem update grid material class $cellClass;
	# The need for this has been removed from UPGRID() - might need to go back...
	#fem update gauss strain fibre components class $tissueClass;
	#fem update grid extension_ratio no_ze_calc component $lMechIdx $lGridArray $lGridIdx class $tissueClass;
}

# Set up the geometry 1*1*1 element 4*4*4mm tri-cubic Hermite basis, 
# 0 fibre angle
fem define parameter;r;$example/minimal;
fem define coordinates;r;$example/emech;
fem define node;r;$example/emech-1x1x1-cubic;
fem define base;r;$example/emech-cubic-linear;
fem define element;r;$example/emech-1x1x1-cubic;
fem define fibre;r;$example/emech-0deg;
fem define element;r;$example/emech-1x1x1-linear fibre;

# Need to calculate the undeformed nodal derivatives
fem update node derivative 1 linear;
fem update node derivative 2 linear;
fem update node derivative 3 linear;

# Node/element groups to set BC's
fem group node external as bdy_nodes;
fem group node 1,10,19,28 as left_end;
fem group node 2,11,20,29 as right_end;
fem group node left_end as corners;
fem group node 29 as wall_x;
fem group element 1 as ALL_ELEMENTS;

# Set up the mechanics
fem define equation;r;$example/mech-cubic;
fem define material;r;$example/mech-cubic;
# fix the corners group and allow wall_x group to slide in x-axis
fem define initial;r;$example/mech-cubic-fixed;
fem define solve;r;$example/newton;

# Export the undeformed geometry
fem export node;$output/slab as slab;
fem export element;$output/slab as slab;
# ... and define the deformed field ...
fem export node;$output/deformed_0 field as slab;
fem export element;$output/deformed field as slab;

# We want to use a CellML model to define the active component of the stress tensor (in the fibre direction).
fem define active;r;$example/cellml;

#
# Now set up a second solve class which we will use to manage the CellML model.
#
# Define the grid to be aligned with the gauss points of the elements.
fem define grid;r;$example/cellml gauss class 2;
fem update grid geometry class 2;
fem update grid metric class 2;
fem group grid element ALL_ELEMENTS as ALL_GRID_POINTS class 2;
fem define equation;r;$example/cellml class 2;
fem define material;d class 2;
# We'll use the Nash steady state tension-length-Ca relation as it is defined internally in cm.
fem define cell;r;$example/nash class 2;
fem define material;r;$example/nash cell class 2;
# define the time variable used to mimic the Ca_actn parameter value from example 5d2
fem define time;r;$example/5d2;
# The variable number comes from the output of the CellML model parsing which shows that Cai is PARAM(4). There 3 derived entries and 1 (reserved) state entries so we want the variable number to be 1+3+4 = 8
fem define initial;r;$example/nash class 2;
fem define solve;r;$example/nash class 2;
fem update grid material class 2;

# export the grid point numbers
fem export element;$output/grid as slab grid_numbers;

# Set the various indices for the cellular (YQS) and gauss (YG) arrays
# Tension index for grid (YQS) and Gauss (YG) variables
fem inquire cell_variable T return_variables T_GRIDARRAY,T_GRIDIDX;
$T_GAUSSARRAY = "YG";
$T_GAUSSIDX = 1;
# !!!! These need to be greater than 6 so that when a update gauss
# !!!! stress/strain is done we don't loose anything
# use this to avoid repeating the evaluation of tension at gauss points
$T_FROM_CELLS_GAUSSIDX = 7;
# Fibre extension ratio index for grid (RCQS) and mechanics variables
# !!!! Note that lambda must be spatially varying?
fem inquire cell_variable lambda return_variables L_GRIDARRAY,L_GRIDIDX;
$L_MECHIDX = 1;
# [Ca]i index for grid (YQS)
fem inquire cell_variable CaiC return_variables CAI_GRIDARRAY,CAI_GRIDIDX;
# Extension ratio index for grid (YQS)
fem inquire cell_variable lambdaC return_variables LC_GRIDARRAY,LC_GRIDIDX;

# initialise the cellular extension ratios
&updateGridPoints(2,1,$L_MECHIDX,$L_GRIDARRAY,$L_GRIDIDX);

# initialise the solve class for the CellML model
fem solve to 0 class 2;

# Open the history files for tension, extension ratio, and [Ca]i
$ITER_UNIT = 21;
$CONV_UNIT = 22;
fem open history;$output/cellml unit $CONV_UNIT write variables yqs niqslist $T_GRIDIDX,$LC_GRIDIDX,$CAI_GRIDIDX binary class 2;
fem open history;$output/cellml-iter unit $ITER_UNIT write variables yqs niqslist $T_GRIDIDX,$LC_GRIDIDX,$CAI_GRIDIDX binary class 2;

# Initialise the distributed mechanics
fem update gauss gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_FROM_CELLS_GAUSSIDX;
fem update gauss gauss_field destination_field $T_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
fem solve increment 0.0 iterate $MECH_MAX_ITERS error $MECH_ERRTOL;

if (!$CONVERGED) {
	die "\n\nSomething went wrong with the initial mechanics solve\n\n";
}

# Set the initial fibre extension ratio for cell problem
#fem update grid extension_ratio component $L_MECHIDX $L_GRIDARRAY $L_GRIDIDX;

# the "time" loop, matching the Ca_actn steps from 5d2

$iterations = 0;
for ($time=0;$time<=4;$time++) {
	print "time = $time\n";
	$CONVERGED = 0;
	$iter = 0;
	# Update the cellular extension ratios at the grid points
	&updateGridPoints(2,1,$L_MECHIDX,$L_GRIDARRAY,$L_GRIDIDX);
	# ...and solve the cell model for the current step
	fem solve class 2 restart to $time;
	fem write history unit $ITER_UNIT time $iterations variables yqs binary class 2;
	$iterations++;
		
	# Update the active tension at the Gauss points
	#fem update gauss gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_FROM_CELLS_GAUSSIDX;
	#fem update gauss gauss_field destination_field $T_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
	# ...and solve the finite elasticity with the new active tensions
	fem solve increment 0.0 iterate $MECH_MAX_ITERS error $MECH_ERRTOL;
	$iter++;

	if ($CONVERGED) {
		print "\nConverged solution at time: $time; after $iter iterations\n";
		# Export the deformed nodes - need field on the export node to get the deformed geometry
		fem export node;$output/deformed_$time field as slab;
		# Export the converged Gauss point strains
		#fem update gauss strain fibre components;
		#fem export gauss;$output/gauss_strain_$time yg as slab;
		# Export the converged Gauss point stresses first need to 
		# put the active tension back into YG (the strains from the
		# previous update have overwritten it)
		#fem update gauss gauss_field destination_field $T_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
		#fem update gauss stress fibre components;
		#fem export gauss;$output/gauss_stress_$time yg as slab;
		#fem update gauss gauss_field destination_field $T_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
		
		# ...and save cellular solution
		fem write history unit $CONV_UNIT time $time variables yqs binary class 2;
	} else {
		print "\n\nMechanics solve unsuccessful at time: $time; $iter iterations\n\n";
		last;
	}
}

# Close the cell history files
fem close history unit $CONV_UNIT binary class 2;
fem close history unit $ITER_UNIT binary class 2;

# Compute the CMISS signal files from the history file
# [Ca]i
fem evaluate electrode;$output/cai history $output/cellml from grid yqs iy $CAI_GRIDIDX binary class 2;
fem evaluate electrode;$output/cai-iter history $output/cellml-iter from grid yqs iy $CAI_GRIDIDX binary class 2;
# tension
fem evaluate electrode;$output/tension history $output/cellml from grid yqs iy $T_GRIDIDX binary class 2;
fem evaluate electrode;$output/tension-iter history $output/cellml-iter from grid yqs iy $T_GRIDIDX binary class 2;
# and extension ratio
fem evaluate electrode;$output/extension_ratio history $output/cellml from grid yqs iy $LC_GRIDIDX binary class 2;
fem evaluate electrode;$output/extension_ratio-iter history $output/cellml-iter from grid yqs iy $LC_GRIDIDX binary class 2;

# Export to UnEMAP signal files
fem define export;r;$example/nash;
fem export signal;$output/cai-iter electrode signal $output/cai-iter;
fem export signal;$output/tension-iter electrode signal $output/tension-iter;
fem export signal;$output/extension_ratio-iter electrode signal $output/extension_ratio-iter;
fem export signal;$output/cai electrode signal $output/cai;
fem export signal;$output/tension electrode signal $output/tension;
fem export signal;$output/extension_ratio electrode signal $output/extension_ratio;

# Check for isovolumic contraction
fem list element total;
fem list element deformed total;

Files used by this example are:

Name                       Modified     Size

example_5d5.com 13-Jul-2009 8.8k 5d2.iptime 13-Jul-2009 1.1k cellml.ipacti 13-Jul-2009 461 cellml.ipequa 13-Jul-2009 1.5k cellml.ipgrid 13-Jul-2009 628 emech-0deg.ipfibr 13-Jul-2009 6.9k emech-1x1x1-cubic.ipelem 13-Jul-2009 573 emech-1x1x1-cubic.ipnode 13-Jul-2009 15k emech-1x1x1-linear.ipelfb 13-Jul-2009 377 emech-cubic-linear.ipbase 13-Jul-2009 7.4k emech.ipcoor 13-Jul-2009 570 mech-cubic-fixed.ipinit 13-Jul-2009 4.0k mech-cubic.ipequa 13-Jul-2009 2.0k mech-cubic.ipmate 13-Jul-2009 6.2k minimal.ippara 13-Jul-2009 5.9k nash.cml 13-Jul-2009 9.0k nash.ipcell 13-Jul-2009 2.1k nash.ipexpo 13-Jul-2009 658 nash.ipinit 13-Jul-2009 660 nash.ipmatc 13-Jul-2009 747 nash.ipsolv 13-Jul-2009 1.8k newton.ipsolv 16-Aug-2010 2.7k newton.ipsolv.old 13-Jul-2009 2.5k

Download the entire example:

Name                      Modified     Size

examples_5_5d_5d5.tar.gz 17-Aug-2010 15k

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

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


CMISS Help / Examples / 5 / 5d / 5d5