Example 5d6: Active contraction of a cube using HMT 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 5d6

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.01;


sub hmtCaiTransient {
    # Routine to define a time variable describing the [Ca]i transient
    # from the HMT paper
    $Cai_init = 0.0001; # mM
    $Cai_max = 0.001;   # mM
    $Cai_Tinit = 1.0;   # ms
    $Tau = 60;          # ms
    $numPoints=100;
    $variableName = "cai";

    CORE::open(IPTIME,">tmp-cai-HMT.iptime") || die "Cannot open iptime file: $!\n";
    print IPTIME " CMISS Version 1.21 iptime File Version 3\n";
    print IPTIME " Heading:\n";
    print IPTIME "\n";
    print IPTIME " Enter the name of the time variable [EXIT]: $variableName\n";
    print IPTIME "\n";
    print IPTIME " Do you want this time variable to be periodic [N]? N\n";
    print IPTIME "\n";
    print IPTIME " Enter the number of time points to be set [2]: $numPoints\n";
    print IPTIME " Enter the time variable value before time point 1 [0.0]: $Cai_init\n";
    for ($point=1;$point<=$numPoints;$point++){
        $time=1000/$numPoints*($point-1);
        if ($time < $Cai_Tinit) {
            $Cai = $Cai_init;
        } else {
            $t = $time - $Cai_Tinit;
            $Cai = $Cai_init+($Cai_max-$Cai_init)*($t/$Tau)*exp(1-$t/$Tau);
        }
        print IPTIME " Enter the time for time point $point [0.0]: $time\n";
        print IPTIME " Enter the time variable value at time point $point [0.0]: $Cai\n";
    }
    print IPTIME " Enter the time variable value after time point $numPoints [0.0]: $Cai_init\n";
    print IPTIME "\n";
    print IPTIME " The type of interpolation is [2]:\n";
    print IPTIME "   (1) Constant\n";
    print IPTIME "   (2) Linear Lagrange\n";
    print IPTIME "  *(3) Quadratic Lagrange\n";
    print IPTIME "  *(4) Cubic Lagrange\n";
    print IPTIME "  *(5) Cubic Hermite\n";
    print IPTIME "    2\n";
    print IPTIME " Enter the name of the time variable [EXIT]: EXIT\n";
    CORE::close IPTIME;
    fem define time;r;tmp-cai-HMT;
    unlink "tmp-cai-HMT.iptime";
} # sub hmtCaiTransient

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/hmt class 2;
fem define material;r;$example/hmt cell class 2;
# define the time variable used to define the [Ca]i transient
&hmtCaiTransient();
# 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/hmt 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<=500;$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_5d6.com 18-Aug-2009 10k 5d2.iptime 18-Aug-2009 1.1k cellml.ipacti 18-Aug-2009 462 cellml.ipequa 18-Aug-2009 1.5k cellml.ipgrid 18-Aug-2009 628 emech-0deg.ipfibr 18-Aug-2009 6.9k emech-1x1x1-cubic.ipelem 18-Aug-2009 573 emech-1x1x1-cubic.ipnode 18-Aug-2009 15k emech-1x1x1-linear.ipelfb 18-Aug-2009 377 emech-cubic-linear.ipbase 18-Aug-2009 7.4k emech.ipcoor 18-Aug-2009 570 hmt.cml 18-Aug-2009 23k hmt.ipcell 18-Aug-2009 3.9k hmt.ipinit 18-Aug-2009 662 hmt.ipmatc 18-Aug-2009 945 mech-cubic-fixed.ipinit 18-Aug-2009 4.0k mech-cubic.ipequa 18-Aug-2009 2.0k mech-cubic.ipmate 18-Aug-2009 6.2k minimal.ippara 18-Aug-2009 5.9k nash.ipexpo 18-Aug-2009 658 nash.ipsolv 18-Aug-2009 1.8k newton.ipsolv 16-Aug-2010 2.7k newton.ipsolv.old 18-Aug-2009 2.5k

Download the entire example:

Name                      Modified     Size

examples_5_5d_5d6.tar.gz 17-Aug-2010 17k

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

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


CMISS Help / Examples / 5 / 5d / 5d6