Example 5d1: Active contraction of a 3D slab using cellular steady-state HMT

This example demonstrates the use of a cellular mechanics model used to generate the active contraction for use in a finite elasticity model. It models a 3D slab of tissue, intially at rest, which is subjected to a uniform [Ca]i transient which drives the contraction. The slab is fixed at one end (in the plane of the element face) and the opposite end is allowed to slide as the muscle contracts.

draw.com can be used to visualise the results of this simulation, as shown in the animation below - which shows the results of this model for a 200ms simulation.


The comfile run by this example is as follows:

#
# 3D slab with HMT used for active tension calculation
#
# Uses the steady state version of HMT at the cellular (grid point)
# level to update the active tension at the gauss points for 
# integration in a finite elasticity solution of a 3-dimensional slab
# of cardiac (pig?) tissue
#

# Set up some colour colour highlighting...
# To turn off all formating
$FBnone = "\033[0m";
# Foreground colurs
$Fblack   = "\033[30m";
$Fred     = "\033[31m";
$Fgreen   = "\033[32m";
$Fyellow  = "\033[33m";
$Fblue    = "\033[34m";
$Fmagneta = "\033[35m";
$Fcyan    = "\033[36m";
$Fwhite   = "\033[37m";
# Background colours
$Bblack   = "\033[40m";
$Bred     = "\033[41m";
$Bgreen   = "\033[42m";
$Byellow  = "\033[43m";
$Bblue    = "\033[44m";
$Bmagneta = "\033[45m";
$Bcyan    = "\033[46m";
$Bwhite   = "\033[47m";

# 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;
}

#
# Main section ...
#

set echo;

$FORMAT = "binary";

# for this example, we do not want to export data for every iteration
$output_iterations = 0;
# we do want to export cellular membrane potentials
$export_potentials = 0;
# we do not want to export any Gauss point information
$export_gauss = 0;
# and we do not want to list out any strains
$list_strains = 0;
# although, if we did we would want them at these Gauss points
@gauss_points = (1,38,39,42,43,22,23,26,27,64);

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

# Enable to solution of the finite elasticty model
$SOLVE_DISTRIBUTED_MECHANICS = 1;
# Set to false to not solve the cell model
$SOLVE_CELL_MODEL = 1;
# We want to use the steady state version of HMT
$STEADY_STATE = 1;
# Set this to 1 to fix both ends of the slab
$ISOMETRIC = 0;

# The different [Ca]i transients - cai.iptime is from Martyn Nash,
# cai_resting.iptime is a constant resting level of [Ca]i and 
# cai_new.iptime is an approximation of the transient presented in
# Stuyvers et. al. Prog. Biophys. Mol. Biol. 69:1998 (Fig. 1), and 
# cai_test.iptime is a scaled down version of this transient - i.e.
# one that works...

#$cai = "cai";
#$cai = "cai_resting";
#$cai = "cai_new";
$cai = "cai_test";

# The two class definitions
$MECHANICS = 1;
$CELL = 2;

# Parameters for the distributed mechanics solution
$MECH_ERRTOL = 0.01;
$MAX_MECH_ITER = 20;

# Set up the geometry 1*8*1 element 4*32*4mm bicubic-Hermite basis, 
# 0 fibre angle
fem define parameter;r;$example/emech;
fem define coordinates;r;$example/emech;
fem define node;r;$example/emech-1x8x1-cubic;
# Need a tri-cubic Hermite for geometry (and bi-cubic Hermite for faces)
# and a tri-linear for fibres (and bi-linear for faces)
# and a auxilary basis for pressure
fem define base;r;$example/emech-cubic-linear;
fem define element;r;$example/emech-1x8x1-cubic;
# 0deg fibre angles everywhere, using a tri-linear interpolation, and
# need to define sheets etc...
fem define fibre;r;$example/emech-0deg;
fem define element;r;$example/emech-1x8x1-linear fibre;

# Need to calculate the undeformed nodal derivatives (using the linear 
# flag since we are modelling a square)
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 9,18,27,36 as right_end;
if ($ISOMETRIC) {
    fem group node left_end,right_end as corners;
} else {
    fem group node left_end as corners;
}
fem group node 18 as wall_x;
fem group element 1..8 as ALL_ELEMENTS;

# Set up the mechanics
fem define equation;r;$example/mech-cubic class $MECHANICS lock;
# using Carey's new material parameters
fem define material;r;$example/mech-cubic class $MECHANICS;
# use the cellular HMT active tension
fem define active;r;$example/mech-cubic class $MECHANICS;
# fix the corners group and allow wall_x group to slide in x-axis
fem define initial;r;$example/mech-cubic-fixed class $MECHANICS;
# use the Newton solver with GMRES
fem define solve;r;$example/newton class $MECHANICS;

# Set up the HMT cellular model - used to calculate the active tension
# Define a 9x9x9 grid scheme for use in all elements and calculate 
# their spatial location and metrics
fem define grid;r;$example/cell-9x9x9 class $CELL;
fem update grid geometry;
fem update grid metric;
# Could use this for coupled models...
#fem group grid xi1=0 oneoff as stimulus_site;
# Always a useful group to have
fem group grid element ALL_ELEMENTS as ALL_GRID_POINTS;
# Define the HMT cardiac mechanics model
fem define equation;r;$example/cell-hmt class $CELL;
# Need to define these material parameters even though they are
# completely unused for this simulation
fem define material;r;$example/cell-hmt class $CELL;
# Define the model parameters - if we wnat to use the steady-state 
# version we need to set the STEADY_STATE flag
if ($STEADY_STATE) {
    fem define cell;r;$example/cell-hmt-ss class $CELL;
} else {
    fem define cell;r;$example/cell-hmt class $CELL;
}
# Define the spatially varying cellular parameters - there are none
fem define material;r;$example/cell-hmt cell class $CELL;
# and update the grid point materials
fem update grid material class $CELL;
# Define the [Ca]i transient (defined above)
fem define time;r;$example."/".$cai;
# Define the initial condition for the cellular model - setting the
# [Ca]i variable in the cell model to be defined via the above time
# variable
fem define initial;r;$example/cell-hmt class $CELL;
# Define a simple Euler interagtion for the cellular ODE's
fem define solve;r;$example/cell-hmt class $CELL;

# Export the undeformed geometry
fem export node;$output/slab as tissue;
fem export element;$output/slab as tissue;
# ... and define the deformed field ...
fem export element;$output/deformed field as tissue;
# and the grid point numbers
fem export element;$output/numbers as tissue grid_numbers;

# Initialise the cell integration
fem solve class $CELL to 0;

# 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 Tension 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
$T_FROM_CELLS_GAUSSIDX = 7; ## use this to avoid repeating the evaluation of tension at gauss points
# Fibre extension ratio index for grid (YQS) and mechanics variables
fem inquire cell_variable ExtensionRatio return_variables L_GRIDARRAY,L_GRIDIDX;
$L_MECHIDX = 1;
# and previous fibre extension ratio index for grid (YQS) and mechanics variables
fem inquire cell_variable ExtensionRatioPrevious return_variables L_GRIDARRAY_PREV,L_GRIDIDX_PREV;
if (($L_GRIDARRAY ne $L_GRIDARRAY_PREV) || ($L_GRIDIDX == $L_GRIDIDX_PREV)) {
    die "Extension ratio current ($L_GRIDARRAY,$L_GRIDIDX) and previous ($L_GRIDARRAY_PREV,$L_GRIDIDX_PREV) arrays do not match!!";
}
# [Ca]i index for grid (YQS)
fem inquire cell_variable Cai return_variables CAI_GRIDARRAY,CAI_GRIDIDX;

# Open the history files for tension, extension ratio, and [Ca]i
# One for the solutions for each iteration and one for the 
# coverged solutions
$ITERATIONS_UNIT = 21;
$CONVERGED_UNIT = 22;
fem open history;$output/cell-hmt_iter write unit $ITERATIONS_UNIT variables yqs niqslist $T_GRIDIDX,$L_GRIDIDX,$CAI_GRIDIDX binary class $CELL;
fem open history;$output/cell-hmt write unit $CONVERGED_UNIT variables yqs niqslist $T_GRIDIDX,$L_GRIDIDX,$CAI_GRIDIDX binary class $CELL;

# Initialise the distributed mechanics
# ?? Need to initialise YG ??
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;
if ($SOLVE_DISTRIBUTED_MECHANICS) {
    fem solve class $MECHANICS increment 0.0 iterate $MAX_MECH_ITER error $MECH_ERRTOL;
} else {
    $CONVERGED = 1;
}

#####################################################################
# NOTE:                                                             #
#   $CONVERGED is set in the CMISS routine NONLIN, either to 1 or 0 #
#     depending on whether convergence was achieved or not in the   #
#     non-linear solution.                                          #
#####################################################################


# save the intial converged solution
if ($CONVERGED) {
    if ($output_iterations && $export_gauss) {
        # Export the intial strains ??
        fem update gauss strain fibre components;
        fem export gauss;$output/gauss_strain_0 yg as tissue;
        # and 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_0 yg as tissue;
        # Need to replace the stress values that have been put into
        # YG(1..) with the initial active tension values
        fem update gauss gauss_field destination_field $T_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
    }
} else {
    die "\n\nSomething is very, very wrong - see ya...\n\n";
}

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

#
# The main time loop....
#
$tend = 5.0;
$tstep = 1.0;
$finish = 0;
$current_time = 0;
$previous_time = 0;
$history_time = 0;
##set echo;
for ($time=0;($time<=$tend)&&(!$finish);$time+=$tstep) {
    $previous_time = $current_time;
    $current_time = $time;

    # Solve the current step
    $iter = 0;
    $CONVERGED = 0;
    while ($iter < $MAX_MECH_ITER) {
	$iter++;
	print "\nIteration: $iter\n\n";
        # Update the cellular extension ratios at the grid points
	fem update grid extension_ratio component $L_MECHIDX $L_GRIDARRAY $L_GRIDIDX class $MECHANICS;
        # Solve the cell model for the current step
	fem solve class $CELL from $previous_time to $current_time;
        # Write to the iterations history file
        fem write history unit $ITERATIONS_UNIT time $history_time variables yqs $FORMAT class $CELL;
        $history_time++;
        # 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;
        # Solve the finite elasticity with the new active tensions
	fem solve class $MECHANICS increment 0.0 iterate 1 error $MECH_ERRTOL;
        # Stop iterating if we have reached a converged solution
	if ($CONVERGED) {
	    last;
	}
    }
    if ($CONVERGED) {
        print $Fcyan."Solve was successful, convergence was found in $iter iteration(s)$FBnone\n";
        # Write to the converged history file
        fem write history unit $CONVERGED_UNIT time $time variables yqs $FORMAT class $CELL;
        # Export the deformed nodes - need field on the export node to 
        # get the deformed geometry
        fem export node;$output/deformed_$time field as tissue;
        if ($export_gauss) {
            # Export the converged Gauss point strains
            fem update gauss strain fibre components;
            fem export gauss;$output/gauss_strain_$time yg as tissue;
            # 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 tissue;
            fem update gauss gauss_field destination_field $T_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
        }
        if ($list_strains) {
            # List out the strain at the specified Gauss points
            $gauss_string = "";
            foreach $gauss_point (@gauss_points) {
                fem list strain;$output/gauss_point_${gauss_point}_strain_$time gauss at $gauss_point fibre;
                $gauss_string .= "$gauss_point,";
            }
            chop $gauss_string;
            # and the grid point closest to them
            fem evaluate strain;$output/grid_point_strains_$time from cell_extension_ratios index $L_GRIDIDX at gauss $gauss_string nearest class $CELL;
        }
    } else {
        print $Bred.$Fcyan."Solve was unsuccessful, convergence was not found with $iter iterations$FBnone\n";
        $finish = 1;
    }
} # end of time loop

# Close the cell history files
fem close history unit $ITERATIONS_UNIT $FORMAT class $CELL;
fem close history unit $CONVERGED_UNIT $FORMAT class $CELL;

# Compute the CMISS signal files from the history file
##set echo;
# [Ca]i
fem evaluate electrode;$output/cai history $output/cell-hmt from grid yqs iy $CAI_GRIDIDX $FORMAT class $CELL;
fem evaluate electrode;$output/cai_iter history $output/cell-hmt_iter from grid yqs iy $CAI_GRIDIDX $FORMAT class $CELL;
# tension
fem evaluate electrode;$output/tension history $output/cell-hmt from grid yqs iy $T_GRIDIDX $FORMAT class $CELL;
fem evaluate electrode;$output/tension_iter history $output/cell-hmt_iter from grid yqs iy $T_GRIDIDX $FORMAT class $CELL;
# and extension ratio
fem evaluate electrode;$output/extension_ratio history $output/cell-hmt from grid yqs iy $L_GRIDIDX $FORMAT class $CELL;
fem evaluate electrode;$output/extension_ratio_iter history $output/cell-hmt_iter from grid yqs iy $L_GRIDIDX $FORMAT class $CELL;

# Export to UnEMAP signal files
fem define export;r;$example/cell;
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;
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;

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



Files used by this example are:

Name                       Modified     Size

example_5d1.com 10-May-2003 14k cai.iptime 30-Jul-2001 1.0k cai_new.iptime 30-Jul-2001 1.3k cai_resting.iptime 12-Apr-2001 728 cai_test.iptime 27-Sep-2001 1.3k cell-9x9x9.ipgrid 06-Mar-2003 610 cell-hmt-ss.ipcell 10-May-2003 3.5k cell-hmt.ipcell 10-May-2003 3.5k cell-hmt.ipequa 26-May-2003 1.2k cell-hmt.ipinit 12-Apr-2001 548 cell-hmt.ipmatc 10-May-2003 356 cell-hmt.ipmate 12-Apr-2001 1.8k cell-hmt.ipsolv 03-Oct-2001 900 cell-hmt.ipsolv.old 12-Apr-2001 854 cell.ipexpo 12-Apr-2001 659 draw.com 10-May-2003 10k emech-0deg.ipfibr 12-Apr-2001 6.9k emech-1x8x1-cubic.ipelem 12-Apr-2001 3.8k emech-1x8x1-cubic.ipnode 12-Apr-2001 67k emech-1x8x1-linear.ipelfb 12-Apr-2001 2.3k emech-cubic-linear.ipbase 12-Apr-2001 7.4k emech.ipcoor 12-Apr-2001 570 emech.ippara 12-Nov-2002 5.9k mech-cubic-fixed.ipinit 10-May-2003 4.0k mech-cubic.ipacti 12-Apr-2001 253 mech-cubic.ipequa 02-May-2004 2.0k mech-cubic.ipmate 10-May-2003 6.2k newton.ipsolv 16-Aug-2010 2.7k newton.ipsolv.old 13-Apr-2007 2.5k

Download the entire example:

Name                      Modified     Size

examples_5_5d_5d1.tar.gz 17-Aug-2010 197k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmSuccessSun Mar 6 00:59:36 2016271
mips-irix
cmSuccessSun Aug 19 04:57:53 20072403
cm64SuccessSun Jul 29 05:23:28 20072543
rs6000-aix
cmSuccessWed Mar 4 01:47:17 2009332
cm64SuccessWed Mar 4 01:48:38 2009326

Testing status by file:


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 / 5d1