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