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