Example 5d3: Active contraction of LV wedge (Gauss point based)

This example uses the old approach to active tension modelling to simlate the active contraction and relaxation of two elements taken from the LV free wall of the RC heart model. In this simulation, the initial geometry is deformed purely by the generation of active tension within the cardiac tissue.

The level of activation is controlled via the ipacti file, which sets the level of activation to use for the current solution. This animation shows the results from the model when the activation level rises from zero to 1.5 and then returns to zero.

Cool animation

draw.com can be used to visualise the results of this simulation.

The comfile run by this example is as follows:

# An example which solves the active contraction of 2 elements taken from
# the LV free wall of the RC heart. It uses the old technique for active 
# contractions (Gauss point based).

# Set up some cool text 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";

sub DefineCalciumLevel {
    # A function to define the current calcium activation factor - done
    # by writing out a temporary IPACTI file with the given level of
    # activation and then reading it in.

    # Parameters
    ($actn) = @_;

    # Create the IPACTI file
    CORE::open(IPACTI,">tmp_ipacti_file.ipacti") || die "Cannot open temporary ipacti file: $!\n";
    print IPACTI " CMISS Version 1.21 ipacti File Version 2\n";
    print IPACTI " Heading: Temp. ipacti file for activation level of $actn\n";
    print IPACTI " \n";
    print IPACTI " Specify type of contraction model [2]:\n";
    print IPACTI "   (1) SS tension-length-Ca relation (OLD)\n";
    print IPACTI "   (2) Hunter/McCulloch/ter Keurs (coupled to cell/grid model)\n";
    print IPACTI "   (3) Unused\n";
    print IPACTI "   (4) Outer hair cell stiffness\n";
    print IPACTI "    1\n";
    print IPACTI " Enter max isometric tension at ext.ratio=1 (Tref) [100kPa]: 0.10000D+03\n";
    print IPACTI " Enter non-dimensional slope parameter (beta) [1.45]: 0.14500D+01\n";
    print IPACTI " Enter c50 for [Ca]i saturation curve (0<c<1) [0.5]: 0.50000D+00\n";
    print IPACTI " Enter Hill coeff. for [Ca]i saturation curve (h) [3.0]: 0.30000D+01\n";
    print IPACTI " Enter max [Ca]i (Ca_max) [1]: 0.10000D+01\n";
    print IPACTI " Specify whether the initial calcium level [Ca]i is [1]:\n";
    print IPACTI "  (1) Constant spatially\n";
    print IPACTI "  (2) Piecewise constant (defined by elements)\n";
    print IPACTI "  (3) Defined by Gauss points\n";
    print IPACTI "   1\n";
    print IPACTI " Enter initial calcium level [Ca]i [0]: $actn\n";
    print "Calcium = $actn\n";
    CORE::close IPACTI;
    # Define the current level of activation
    fem define active;r;tmp_ipacti_file;
    # And delete the temporary file
    unlink "tmp_ipacti_file.ipacti";

} # sub DefineCalciumLevel

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

# Set the number of threads to use
if (! defined $num_threads) {
    $num_threads = 1;
}
set num_threads $num_threads;

# Set the output directory and make sure it exists
if (! defined $output) {
    #$output = "output-$ENV{'HOSTTYPE'}";
    $output = ".";
}
if (! -d $output) {
    (mkdir $output,0700) || die "Unable to create output directory ($output): $!\n";
}

# Decide what information to output
$output_iterations = 0;
$export_gauss = 0;
$list_strains = 0;
@gauss_points = (1,38,39,42,43,22,23,26,27,64);

# Set up the geometry 
fem define parameter;r;$example/minimal;
fem define coordinates;r;$example/emech;
fem define node;r;$example/wedge;
# Need a tri-cubic Hermite for geometry and fibres (and bi-cubic Hermite for 
# faces) and a tri-linear (and bi-linear for faces) and a auxilary basis for 
# pressure
fem define base;r;$example/emech-cubic-linear;
fem define element;r;$example/wedge;
# Use the fitted fibres and assume a 90 degree sheet angle everywhere
fem define fibre;r;$example/wedge;
fem define element;r;$example/wedge fibre;

# Node/element groups to set BC's
fem group node external as bdy_nodes;
# Fixed in x,y,z and all derivatives
fem group node 17,21,82,85 as fixed_all;
# Fixed in z and all derivatives
fem group node 19,22,84,86,31,33,35,36 as fixed_deriv;
fem group element 10,29 as ALL_ELEMENTS;

# Set up the mechanics
fem define equation;r;$example/mech-cubic;
fem define material;r;$example/mech-cubic-carey;
# fix the corners group
fem define initial;r;$example/mech-cubic-fixed;
fem define solve;r;$example/newton;

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

# Loop through the activation levels desired
$counter = 0;
$iteration = 0;
$max = 30.0;
$MAX_ITERS = 20;
for ($activation_level=0;$activation_level<=$max;$activation_level+=0.5) {
#    if ($activation_level > $max/2) {
#	$counter++;
#	$actn = $max/2 - $counter;
#    } else {
#	$actn = $activation_level;
#    }
    if (($activation_level <= 5) || ($activation_level > 25)) {
        $actn = 0.0;
    } elsif ($activation_level <= $max/2) {
        $actn = $activation_level - 5;
    } else {
        $counter+=0.5;
        $actn = $max/2 - $counter - 5;
    }
    $actn *= 0.1;
    print $Fred.$Byellow."Activation level: $actn$FBnone\n";
    # Define the current activation level
    &DefineCalciumLevel($actn);
    # Solve the current step
    $CONVERGED = 0;
    $iter = 0;
    while (! $CONVERGED) {
        $iter++;
        $iteration++;
        fem solve increment 0.0 iterate 1 error 0.001;
        if ($output_iterations) {
            if ($export_gauss) {
                # Export the strains
                fem update gauss strain fibre components;
                fem export gauss;$output/iter_gauss_strain_$iteration yg as wedge;
                # and stresses
                fem update gauss stress fibre components;
                fem export gauss;$output/iter_gauss_stress_$iteration yg as wedge;
                fem update gauss stress fibre active components nominal;
                fem export gauss;$output/iter_gauss_active_stress_$iteration yg as wedge;
            }
            if ($list_strains) {
                # List out the Gauss point strains
                foreach $gauss_point (@gauss_points) {
                    fem list strain;$output/gauss_point_."$gauss_point"._strain_iterations_$iteration gauss at $gauss_point fibre;
                }
            }
            # and nodes
            fem export node;$output/iter_deformed_$iteration field as wedge;
        }
        if ($iter > $MAX_ITERS) {
            die "Too many iterations";
        }
    }
    
    if ($CONVERGED) {
        print $Fyellow.$Bred." Convergence achieved in $iter iteration(s) $FBnone\n";
        if ($export_gauss) {
            # Export the strains
            fem update gauss strain fibre components;
            fem export gauss;$output/gauss_strain_$activation_level yg as wedge;
            # and stresses
            fem update gauss stress fibre components;
            fem export gauss;$output/gauss_stress_$activation_level yg as wedge;
        }
        # and nodes
        fem export node;$output/deformed_$activation_level field as wedge;
        # List out the Gauss point strains
        if ($list_strains) {
            foreach $gauss_point (@gauss_points) {
                fem list strain;$output/gauss_point_."$gauss_point"._strain_$activation_level gauss at $gauss_point fibre;
            }
        }
    } else {
        die $Fred.$Byellow." Failed to converge... $FBnone\n";
    }
}

Files used by this example are:

Name                       Modified     Size


example_5d3.com            07-Mar-2004  7.4k
draw.com                   01-Oct-2001  9.5k
emech-cubic-linear.ipbase  01-Oct-2001  7.4k
emech.ipcoor               01-Oct-2001  570 
mech-cubic-carey.ipmate    05-Dec-2002  6.2k
mech-cubic-fixed.ipinit    05-Dec-2002  5.4k
mech-cubic.ipequa          02-May-2004  2.0k
minimal.ippara             12-Nov-2002  5.9k
newton.ipsolv              16-Aug-2010  2.7k
newton.ipsolv.old          13-Apr-2007  2.5k
wedge.ipelem               01-Oct-2001  1.2k
wedge.ipelfb               01-Oct-2001  624 
wedge.ipfibr               01-Oct-2001   19k
wedge.ipnode               01-Oct-2001   23k

Download the entire example:

Name                      Modified     Size


examples_5_5d_5d3.tar.gz  17-Aug-2010  179k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmFailureSun Mar 6 00:33:41 201618
last breakTue Feb 24 03:34:00 201519
last successSun Jul 3 02:31:00 200576
mips-irix
cmFailureTue Aug 21 06:03:35 2007186
last breakThu Apr 13 03:51:00 2006229
last successSun Jul 3 01:43:00 200573
cm64FailureMon Aug 20 05:36:41 2007213
last breakThu Apr 13 04:04:00 2006197
last successFri May 20 10:49:00 200579
rs6000-aix
cmFailureWed Mar 4 01:29:40 200924
last breakSat Apr 14 02:47:00 200725
last successSun Jul 3 03:10:00 200527
cm64FailureWed Mar 4 01:30:03 200927
last breakFri Nov 3 02:32:00 200647
last successSun Jul 3 03:11:00 200528

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