Example 5522: Dynamic active contraction of a high-order prolate LV using HMT

Created by: Martyn Nash and David Nickerson June 2000

WARNING: This example does not run correctly!

Dynamic HMT was not originally with compression in mind (model validation was based on tension experiments), but in any 3D incompressible model, some structures will experience compressive deformations. As a result, the history terms in the HMT fading memory description (Q), become unstable and hence the model solutions become unreliable.

Sets up an LV mechanics model coupled to a HMT cell model, from which the active fibre tension is derived. The LV, made up of two high order prolate elements, has been previously inflated (passively) to 1kPa. The solution phase for active contraction uses a [Ca]i pulse (defined using a time variable) to drive the HMT model, which computes the active tension at all computational points, causing the LV to contract.

The comfile run by this example is as follows:

# Example 5522 - Dynamic active contraction of a high-order prolate LV using HMT

# Created by:Martyn Nash and David Nickerson June 2000

#WARNING: This example does not run correctly!

# Dynamic HMT was not originally with compression in mind (model validation
# was based on tension experiments), but in any 3D incompressible model, some
# structures will experience compressive deformations.  As a result, the history
# terms in the HMT fading memory description (Q), become unstable and hence the
# model solutions become unreliable.

# Sets up an LV mechanics model coupled to a HMT cell model, from which
# the active fibre tension is derived.  The LV, made up of two high
# order prolate elements, has been previously inflated (passively) to
# 1kPa.  The solution phase for active contraction uses a [Ca]i pulse
# (defined using a time variable) to drive the HMT model, which computes
# the active tension at all computational points, causing the LV to
# contract.

$TRUE = 1
$FALSE = 0

$SOLVE_INFLATION = $FALSE
$EXPORT = $FALSE
$MECH_ERRTOL = 0.001

# Solve classes

$MECH = 1
$CELL = 2

# Model regions

$WALL = 1

fem define parameters;r;mech_hmt;example

# Set up the ventricular wall geometry
fem define coordinate;r;wall;example region $WALL #prolate spheriodal coords
fem define node;r;wall;example region $WALL    #6 nodes
fem define base;r;wall;example                 #lambda tricub; mu/theta trilin
fem define element;r;wall;example region $WALL #2 wall elements
fem define fibre;r;wall;example region $WALL   #fibre: endo  65 deg; epi -55 deg
fem define element;r;wall;example fibre region $WALL #sheet: endo -45 deg; epi  45 deg
fem group node 1..4,7,8 as ALL_NODES
fem group element 1..2 as ALL_ELEMENTS

# Set up the mechanics model
fem define equation;r;wall;example region $WALL class $MECH lock #Incomp finelas; FEM
fem define material;r;wall;example region $WALL class $MECH    #pole zero law
fem define active;r;mech_hmt;example region $WALL class $MECH  #passiv inflation
fem define initial;r;wall;example region $WALL class $MECH     #inflate to 1kPa
fem define solve;r;wall;example region $WALL class $MECH       #GMRES

# Set up the HMT cell model
fem define grid;r;mech_hmt;example region $WALL class $CELL  #5*5*5 grid scheme
fem group grid 1..225 as ALL_GRID_POINTS
fem update grid geometry region $WALL
fem update grid metric region $WALL
fem define equation;r;hmt;example region $WALL class $CELL lock  #HMT cell model
# Need to define these material parameters even though they are
# completely unused for this simulation
fem define material;r;hmt;example class $CELL
fem define cell;r;hmt;example class $CELL
fem define material;r;hmt;example cell class $CELL
fem update grid material class $CELL
fem define time;r;cai;example #time variable for the conc of intracellular Ca++
fem define initial;r;hmt;example class $CELL #use time-varying length/Cai
fem define solve;r;hmt;example class $CELL

# Solve the HMT cell model for initial setup
fem solve class $CELL to 0 

# Solve the mechanics model incrementally for passive inflation
if ($SOLVE_INFLATION) {
  fem solve class $MECH increment 0.1 iter 2 error $MECH_ERRTOL
  fem solve class $MECH increment 0.4 iter 2 error $MECH_ERRTOL
  fem solve class $MECH increment 0.5 iter 99 error $MECH_ERRTOL
  fem define initial;w;wall_infl region $WALL class $MECH #inflated
} else {
  fem define initial;r;wall_infl;example region $WALL class $MECH #inflated
  fem define solve;r;wall;example region $WALL class $MECH        #GMRES
}

if ($EXPORT) {
  fem export node;test as test
  fem export element;test as test
  fem export element;parameters as parameters field cell
  fem export element;numbers as numbers grid_numbers
}


########################################################################
# Solve the active contraction problem using HMT
########################################################################

# Output file names

$SOLNSDIR = "solutions/"
$OUTFILE = "${SOLNSDIR}mech_hmt"
$MATRIXFILE = "${SOLNSDIR}matrices_yq_yqs"
$MECHSOLFILE = "${SOLNSDIR}mechsol"
$HISTFILE = "${SOLNSDIR}hmt"
$SIGFILEBASE = "${SOLNSDIR}hmt_"

# Set the reference tension (in kPa) and tolerance for the 
# maximum and minimum absolute differences (percentages)
# (note this is specified at a percentage of T_REF)
$T_REF = 100.0

$TDIFF_TOLMIN = 1
$TDIFF_TOLMAX = 10
$PCT_TDIFF = 0

# Other variables

$T_MAXABSDIFF = 0.0
$SOLVEMECH = $FALSE
$MAX_ITER_MECH = 10
$ITER_MECH2 = 0
$MAX_ITER_MECH2 = 3
$ITER_CELL = 0
$MAX_ITER_CELL = 10

# Signals to export
# Vm = 1, Cab = 2, z = 3, L = 4, Tension = 5, IsoTension = 6, Cai = 7
$SIGNALS = "1..7"

# File format
$FORMAT = "binary"

# 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

# Fibre extension ratio index for grid (YQS) and mechanics variables
fem inquire cell_variable ExtensionRatio return_variables L_GRIDARRAY,L_GRIDIDX

$L_MECHIDX = 1


# Time and stepping variables

$TIME = 0.0
$TEND = 7.0
$TSTEP = 1.0

# Open the HMT cell history file
#fem open history;$HISTFILE write unit 21 variables yqs niqslist $SIGNALS $FORMAT class $CELL

# Compute and store fibre extension ratio for HMT grid problem
fem update grid extension_ratio component $L_MECHIDX $L_GRIDARRAY $L_GRIDIDX region $WALL class $MECH

# Interpolate initial active tension from grid point variables to 
# Gauss pts
fem update gauss gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_GAUSSIDX;

# Set up checkpoint matrix
# Write the YQ (grid contiuum vars) and YQS (cell vars) matrices
#fem write matrix;$MATRIXFILE//_//format($LASTMECHTIME,F8.4) matrices YQ,YQS $FORMAT

$LASTMECHTIME = $TIME

# Echo commands for command file debugging
set echo

########################################################################
# MAIN SOLUTION LOOP
########################################################################
while(($TIME < $TEND) && ($ITER_MECH2 < $MAX_ITER_MECH2)) {
  $TIME = $TIME + $TSTEP

  print " \n"
  print " *******************************************\n"

  print " TIME = ".sprintf("%16.4f",${TIME})." msec\n"


  fem evaluate time variable cai time $TIME

# Solve the HMT cell model for the current time
  fem solve restart to $TIME class $CELL

# Get maximum absolute tension difference over all Gauss pts
  fem compare gaussvars gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_GAUSSIDX region $WALL return maxabsdiff to_variable T_MAXABSDIFF 


########################################################################
# Determine if the mechanics problem needs resolving/updating

  $PCT_TDIFF = 100.0 * $T_MAXABSDIFF / $T_REF

  print " \n"
  print " Maximum tension difference: ".sprintf("%8.2f",${PCT_TDIFF})." %\n"


  if($PCT_TDIFF < $TDIFF_TOLMIN) {
#   Tension has not changed appreciably so don't solve mechanics

#   Write the grid variables to the history file
#    fem write history unit 21 time $TIME variables yqs $FORMAT class $CELL

#    TSTEP = $TSTEP * 2.0

    $SOLVEMECH=$FALSE


########################################################################

  } elsif($PCT_TDIFF > $TDIFF_TOLMAX) {

#   Tension has changed too much so restart from previous cell checkpoint
#   time step, but using a smaller time step (don't solve mechanics)

    print " >>>INFO: Tension has changed by too much. Re-starting from previous checkpoint.\n"

    print " \n"

#   Restart cell problem from last checkpoint
#   Read the YQ (grid contiuum vars) and YQS (cell vars) matrices
		fem read matrix;${MATRIXFILE}."_".sprintf("%8.4f",$LASTMECHTIME) matrices YQ,YQS $FORMAT

    $TIME = $TIME - $TSTEP

    fem solve from $LASTMECHTIME to $TIME class $CELL

    $TSTEP = $TSTEP / 2.0

    print " \n"

		print " Time step halved: TSTEP = ".sprintf("%16.8f",${TSTEP})." msec\n"


    print " \n"
    print " *******************************************\n"

		print " TIME = ".sprintf("%16.4f",${TIME})." msec\n"

    print " \n"

    fem evaluate time variable cai time $TIME

    $ITER_CELL = $ITER_CELL + 1

    if($ITER_CELL == $MAX_ITER_CELL) {
      $SOLVEMECH=$TRUE
    } else {
      $SOLVEMECH=$FALSE
    }

########################################################################

  } else {

#   Tension has changed by a suitable amount to now resolve mechanics

    $SOLVEMECH=$TRUE
  }



########################################################################
  if($SOLVEMECH) {
    if($LASTMECHTIME == $TIME) {
      $ITER_MECH2 = $ITER_MECH2 + 1
    } else {
      $ITER_MECH2 = 0
    }

#   Checkpoint the cell problem
    $LASTMECHTIME=$TIME


#   Write the YQ (grid contiuum vars) and YQS (cell vars) matrices
#    fem write matrix;$MATRIXFILE//_//format($LASTMECHTIME,F8.4) matrices YQ,YQS $FORMAT 

#   Interpolate active tension from grid point variables to Gauss pts
    fem update gauss gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_GAUSSIDX region $WALL

#   Solve the mechanics problem using the updated active tension
    fem solve class $MECH increment 0.0 iterate $MAX_ITER_MECH error $MECH_ERRTOL

#   Save mechanics solution
#    fem define initial;w;$MECHSOLFILE//_//format($TIME,F8.4) region $WALL class $MECH

#   Compute and store fibre extension ratio for HMT grid problem
    fem update grid extension_ratio component $L_MECHIDX $L_GRIDARRAY $L_GRIDIDX region $WALL class $MECH

    $ITER_CELL = 0
  }
}


# Close the HMT cell history file
#fem close history unit 21 $FORMAT class $CELL

Additional testing commands:

fem list initial region $WALL class $MECH

Files used by this example are:

Name                Modified     Size


example_5522.com    09-Apr-2003  9.2k
cai.iptime          30-Jul-2001  1.0k
hmt.ipcell          31-Jul-2000  3.8k
hmt.ipequa          26-May-2003  1.2k
hmt.ipexpo          19-Jul-2000  659 
hmt.ipinit          19-Jul-2000  548 
hmt.ipmatc          07-Apr-2003  356 
hmt.ipmate          03-May-2001  1.9k
hmt.ipsolv          03-Oct-2001  900 
hmt.ipsolv.old      19-Jul-2000  854 
hmt_all.ipexpo      19-Jul-2000  657 
mech_hmt.ipacti     19-Jul-2000  331 
mech_hmt.ipgrid     06-Mar-2003  610 
mech_hmt.ippara     12-Nov-2002  5.9k
mech_hmt_setup.com  24-Aug-2000  4.2k
mech_hmt_solve.com  24-Aug-2000   13k
solutions/          27-Feb-2004    - 
test_output.com     19-Jul-2000   43 
wall.ipbase         19-Jul-2000   12k
wall.ipcoor         19-Jul-2000  710 
wall.ipelem         19-Jul-2000  1.2k
wall.ipelfb         19-Jul-2000  618 
wall.ipequa         02-May-2004  2.0k
wall.ipfibr         30-Jan-2001  1.7k
wall.ipinit         12-Dec-2002  1.6k
wall.ipmate         12-Dec-2002  6.3k
wall.ipnode         19-Jul-2000  4.8k
wall.ippara         12-Nov-2002  5.9k
wall.ipsolv         16-Aug-2010  2.3k
wall.ipsolv.old     13-Apr-2007  2.1k
wall_infl.ipinit    12-Dec-2002   15k

Download the entire example:

Name                           Modified     Size


examples_5_55_552_5522.tar.gz  17-Aug-2010   71k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmFailureSun Apr 6 01:06:01 20032
last breakSun Feb 2 01:05:00 20032
cm-debugFailureWed Apr 9 01:07:11 200322
last breakTue Apr 8 01:10:00 200322
last successSun Sep 9 22:11:00 200184
mips-irix
cmFailureSun Apr 6 02:12:44 20037
last breakSun Feb 2 02:11:00 200317
cm-debugFailureWed Apr 9 01:31:57 200313
last breakTue Apr 8 01:44:00 200318
last successMon Sep 10 01:24:00 2001134
cm-debug-clear-mallocFailureSun Apr 6 22:32:19 200314
last breakSun Mar 2 22:28:00 200315
last successSun Sep 9 22:31:00 2001142
rs6000-aix
cmFailureSun Apr 6 03:19:36 20031
last breakSun Feb 2 03:25:00 20032
cm-debugFailureWed Apr 9 02:21:56 20032
last breakTue Apr 8 02:25:00 20033
cm64FailureTue Feb 18 20:21:19 20032
last breakFri Jan 31 09:48:00 20036

Testing status by file:

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

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

CMISS Help / Examples / 5 / 55 / 552 / 5522