Example 811e: Cardiac cycle simulation using the porcine ventricular model

Simulation of the whole cardiac cyle. The various phases of the cardiac cycle are simulated with appropriate boundary conditions.

The model is first passively inflated during the diastasis and atrial contraction phase.

Isovolumic contraction is simulated by incrementing active fibre tension (incrementing the CA parameter), while at each increment the cavity pressures are adjusted to keep the cavity volumes constant.

The rapid ejection phase is simulated by incrementing the pressure. At each pressure increment a specified volume is given in a file. The CA parameter (active fiber tension) is adjusted to obtain the desired volume. The pressure-volume relationship may be changed to e.g. a Windkessel-model by adjusting the values in the file.

The reduced ejection phase is simulated by keeping CA constant and reducing the pressure slightly.

Isovolumic relaxation is simulated similarly to IVC.

The rapid filling phase is simulated by reducing the CA and pressures down to 0.

Note that the MRP01.ipline and MRP01.ipnode files in this example are slightly different to the same files in examples 811a-d. A couple of scalefactors and derivatives are changed as I believe they were wrong and caused non-overlapping element surfaces. In addition the LV endocardial apex is made slightly rounder in an attempt to improve numerical stability.

Created by Espen Remme May 2004.
 


The comfile run by this example is as follows:

# Example 811e

# If the example path is not set, default to current directory
if(!defined $example)
{
  $example = "./";
}

$TRUE = 1;  # Initialise logicals
$FALSE = 0;
#
use MySubs; # file manipulation & plotting routines 
#
do "formats.com" # Cool colours

set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n"; 
print "\033[0;30;42m                       Execution Flow                         \033[0m\n"; 
print "\033[0;30;42m ============================================================ \033[0m\n"; 
#==========================================================================================================
set echo;
#
set num_threads 1; 
#
# Cell model 
#
$CELL_MODEL="FM_GAUSS";
#$CELL_MODEL="FM";
#$CELL_MODEL="HMT";
#
#
# Output 
#
$OUTPUT_DIR = "output_".$CELL_MODEL."/";
if( ! -d ${OUTPUT_DIR})
{
    mkdir ${OUTPUT_DIR};
}
$output = "output/";
if( ! -d ${output})
{
    mkdir ${output};
}
#
# Solve phases
#
$SOLVE_INFLATION = $TRUE;
#$SOLVE_INFLATION = $FALSE;
#
$SOLVE_ACTIVE_CONTRACTION = $TRUE;
#$SOLVE_ACTIVE_CONTRACTION = $FALSE;

$CAVITY_REGIONS = $TRUE;

$SOLVE_EJECTION = $TRUE;
#$SOLVE_EJECTION = $FALSE;
#
set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n"; 
print "\033[0;30;42m                    Regions & Classes                         \033[0m\n"; 
print "\033[0;30;42m ============================================================ \033[0m\n"; 
#==========================================================================================================
set echo;
#
# Regions
#
$WALL = 1;
$BASE = 2;
$LV_CAVITY = 3;
$RV_CAVITY = 4;
#
# Classes
#
$MECHANICS = 1;
$CELLULAR = 2;
$STRESS_FIT = 3;


set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n"; 
print "\033[0;30;42m                            Load Model                        \033[0m\n"; 
print "\033[0;30;42m ============================================================ \033[0m\n"; 
#==========================================================================================================
set echo;
#
set echo;
#==========================================================================================================
print "\033[0;30;42m ==================================================================== \033[0m\n"; 
print "\033[0;30;42m         Nodes, Elements, Fibres & Sheets for Pig heart Model         \033[0m\n"; 
print "\033[0;30;42m ==================================================================== \033[0m\n"; 
#==========================================================================================================
set echo;
#
fem define parameters;r;${example}MRP01_refined;
fem define coordinates;r;${example}MRP01;
fem define region;r;${example}MRP01;
fem define base;r;${example}MRP01;                                 # 1..8 Geometry and fibres
fem define;add bases;r;${example}3Linear-2PressAuxXi3_4x4x4Gauss;  # 9 Pressure
#
fem define node;r;${example}MRP01;
fem define nodes;r;${example}MRP01_base region 2;
fem define elem;r;${example}MRP01;
fem define elements;r;${example}MRP01_base region 2;
#
fem define base;r;${example}MRP01_readse;                          #redefine bases to read in scale factors 
fem define;add bases;r;${example}3Linear-2PressAuxXi3_4x4x4Gauss; 
fem define line;r;${example}MRP01;
fem define line;r;${example}MRP01_base region 2;
#
#fem update nodes hanging nodes 144,146 element 47 direction 1;
# 
fem define fibre;r;${example}MRP01;
fem define elem;r;${example}MRP01 fibre;
#
# Node Groups
#
fem group nodes 3,4,7,8,18,20,32,34,89..94,105..106,115..120,121,122,123,134..136,109,110 as FIXED;
fem group nodes 9,13,21,23,35,37,41..42,45..46,49..50,53..54,81..82,95..96 as EQUATOR_LV_ENDO;
fem group nodes 22,24,36,38,83..84,11..12,15..16,97..98,101..102 as EQUATOR_LV_EPI;
fem group nodes 1..8,10,14,17..20,25..28,30..34,39..40,43..44,47..48,51..52,55..64,66,69,72,75,85..94,99,100,103..110,120..133,135..142,65,67..68,70..71,73..74,76..80,111..119,134 as REMAINDER_INITIAL;
# Nodes for material properties
fem group nodes 29,30 as apexNodes;
fem group nodes 25..28,39,40,87,88,103,104,124..133,137..142 as subApexNodes;
fem group nodes 13..16,23,24,37,38,41,43,45,47,49,51,53,55,66,69,72,75,81,83,95,97,99,101 as ring2;
fem group nodes 9..12,21,22,35,36,42,44,46,48,50,52,54,56,65,67,68,70,71,73,74,76,82,84,96,98,100,102 as ring3;
fem group nodes 1,2,5,6,17,19,31,33,57..64,77..80,85,86,107..114 as ring4;
fem group nodes ring2,ring3,ring4 as normalNodes;
fem group nodes 3,4,7,8,18,20,32,34,89..94,105,106,115..123,134..136 as baseNodes;
#
# Element Groups
#
fem group elements 14..21,32..34,37,43,51..52 as SEPTUM;
fem group elements 22,26,23..25,27,36,45..50,55 as RV_FREE_WALL;
#
#fem group elements 8,13,31,39,53..54,59..63,76..80,22,23,71,72,1,36,75 as ISO;
#fem group elements 2..7,9..12,14..21,24..30,32..35,37,38,40..52,55,57..58,64,68..70,73,74,65,66,67 as ANISO;
fem group elements 2,3,5..7,10..12,14..30,35..47,55,64..75 as ISO;
fem group elements 8,13,31,39,53,54,59..63,76..80,1,4,9,41,32..34,48..52,57,58 as ANISO;


fem group elements 4..13,28,29,30,31,35,39,40,41,42,53,54,57,59,61..62,65..69  as LV_FREE_WALL_INITIAL;
#
set echo;

fem group nodes REMAINDER_INITIAL as REMAINDER;
fem group elements LV_FREE_WALL_INITIAL as LV_FREE_WALL;



set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n"; 
print "\033[0;30;42m                     Equations                                \033[0m\n"; 
print "\033[0;30;42m ============================================================ \033[0m\n"; 
#==========================================================================================================
set echo;
#
fem define equation;r;${example}finelas_tch_incomp lock region $WALL;

# Removed mapping as cmiss didn't seem to handle it correctly - EWR
#fem define mapping;r;${example}MRP01;

set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n"; 
print "\033[0;30;42m                     Material Properties                      \033[0m\n"; 
print "\033[0;30;42m ============================================================ \033[0m\n"; 
#==========================================================================================================
set echo;
#fem define material;r;${example}polezero region $WALL;
fem define material;r;${example}polezero_gauss region $WALL;
MySubs::changeCalcium(0);
fem define active;r;${example}."calcium0xx" region $WALL;

set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n"; 
print "\033[0;30;42m                      Initial Conditions                      \033[0m\n"; 
print "\033[0;30;42m ============================================================ \033[0m\n"; 
#==========================================================================================================
set echo;
fem define initial;r;${example}MRP01 region $WALL;


set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n"; 
print "\033[0;30;42m                       Solution Method                        \033[0m\n"; 
print "\033[0;30;42m ============================================================ \033[0m\n"; 
#==========================================================================================================
set echo;
#
#fem define solve;r;newton region $WALL; 
fem define solve;r;${example}lu region $WALL; 

$name = "heartOriginal"
$file = $output.$name
fem export nodes;$file as $name offset 1000
fem export elements;$file as $name offset_elem 1000


set echo; 
$ITERS = 15;
$TOL = 0.001;

if ( $SOLVE_INFLATION ) 
{
  read commands;${example}solveInflation; 
} else 
{
    #==========================================================================================================
    print "\033[0;30;42m ============================================================ \033[0m\n"; 
    print "\033[0;30;42m                    End-Diastolic State                       \033[0m\n"; 
    print "\033[0;30;42m ============================================================ \033[0m\n"; 
    #==========================================================================================================
    #
    print "\033[0;30;43m                Reading in end-diastolic state                \033[0m\n";
    $FILE = "full_heart_press_1";
    $PRESS = 1;
    $FILE = "full_heart_press_0.9";
    $PRESS = 0.9;
    fem define initial;r;${OUTPUT_DIR}.${FILE} region $WALL;
#    fem define solve;r;newton region $WALL class $MECHANICS;
    fem define solve;r;${example}lu region $WALL;
    #fem solve increment 0.0 iter 20 error $TOL ;
    fem eval resid wrt geom reg $WALL
}
$name = "heartInflated"
$file = $output.$name
fem export nodes;$file field as $name offset 2000
fem export elements;$file field as $name offset_elem 2000
set echo; 

#die("")

set echo; 
if ( $CAVITY_REGIONS ) 
{           
    #==========================================================================================================
    print "\033[0;30;42m ============================================================ \033[0m\n"; 
    print "\033[0;30;42m                    Set up Cavity Regions                     \033[0m\n"; 
    print "\033[0;30;42m ============================================================ \033[0m\n"; 
    #==========================================================================================================
    #
    fem define base;r;${example}MRP01;                                 # Redefine bases to calculate scale factors 
    fem define;add bases;r;${example}3Linear-2PressAuxXi3_4x4x4Gauss;  
    #
    # Define geometry for LV and RV cavities
    #
    fem define region;r;${example}cavities;
    #
    fem define coord;r;${example}lv_cavity region $LV_CAVITY;
    fem define node;r;${example}lv_cavity region $LV_CAVITY;
    fem define element;r;${example}lv_cavity region $LV_CAVITY;
    #
    fem define coord;r;${example}rv_cavity region $RV_CAVITY;
    fem define node;r;${example}rv_cavity region $RV_CAVITY;
    fem define element;r;${example}rv_cavity region $RV_CAVITY;
    #
    # Calculate the cavity volumes
    #
    fem list elements deformed total cavity_volume region $LV_CAVITY; # Defines and sets the var $LV_CAVITY_VOLUME
    fem list elements deformed total cavity_volume region $RV_CAVITY; # Defines and sets the var $RV_CAVITY_VOLUME
    print "\nLV cavity volume = ${LV_CAVITY_VOLUME}\n";
    print "RV cavity volume = ${RV_CAVITY_VOLUME}\n\n";
    $DIASTOLIC_LV_CAVITY_VOLUME = $LV_CAVITY_VOLUME;
    $DIASTOLIC_RV_CAVITY_VOLUME = $RV_CAVITY_VOLUME;
    #
    # Read back in scale factors for the ventricular wall elements
    #
    fem define base;r;${example}MRP01_readse;                          # Redefine bases to read in scale factors 
    fem define;add bases;r;${example}3Linear-2PressAuxXi3_4x4x4Gauss; 
    fem define line;r;${example}MRP01;
    $LV_PRESS = $PRESS;
    $RV_PRESS = $PRESS*0.2;
}
set echo;



set echo;
if ( $SOLVE_ACTIVE_CONTRACTION ) 
{
  read commands;${example}solveIVC; 

  $name = "heartIVC"
  $file = $output.$name
  fem export nodes;$file field as $name offset 3000
  fem export elements;$file field as $name offset_elem 3000
} else 
{
}

set echo;


set echo;
if ( $SOLVE_EJECTION ) 
{
  read commands;${example}solveEjection; 
}

if(1)
{
  read commands;${example}solveReduced; 

  $name = "heartEndEjection"
  $file = $output.$name
  fem export nodes;$file field as $name offset 5000
  fem export elements;$file field as $name offset_elem 5000
}

set echo;


read com;solveIVR;
#read com;solveIVR2;

read com;solveRecoil;




Files used by this example are:

Name                                    Modified     Size

example_811e.com 04-May-2004 12k 3Linear-2PressAuxXi3_4x4x4Gauss.ipbase 03-May-2004 1.9k MRP01.ipbase 03-May-2004 11k MRP01.ipcoor 03-May-2004 688 MRP01.ipelem 03-May-2004 100k MRP01.ipelfb 03-May-2004 67k MRP01.ipfibr 03-May-2004 415k MRP01.ipinit 03-May-2004 70k MRP01.ipline 03-May-2004 455k MRP01.ipmap 03-May-2004 146k MRP01.ipnode 03-May-2004 461k MRP01.ipregi 03-May-2004 142 MRP01_base.ipelem 03-May-2004 12k MRP01_base.ipline 03-May-2004 36k MRP01_base.ipnode 03-May-2004 121k MRP01_readse.ipbase 03-May-2004 11k MRP01_refined.ippara 03-May-2004 5.9k MySubs.pm 03-May-2004 4.4k calcium0xx.ipacti 03-May-2004 1.5k cavities.ipregi 03-May-2004 130 cavityAuxiliaryBasis.ipbase 03-May-2004 1.1k coupled.ipcoup 03-May-2004 437 export.com 03-May-2004 592 finelas_tch_incomp.ipequa 04-May-2004 2.1k formats.com 03-May-2004 726 initLVCavity.ipinit 03-May-2004 31k initRVCavity.ipinit 03-May-2004 31k lu.ipsolv 16-Aug-2010 2.4k lu.ipsolv.old 13-Apr-2007 2.2k lv_cavity.ipcoor 03-May-2004 688 lv_cavity.ipelem 03-May-2004 42k lv_cavity.ipnode 03-May-2004 157k polezero_gauss.ipmate 03-May-2004 19k rv_cavity.ipcoor 03-May-2004 688 rv_cavity.ipelem 03-May-2004 12k rv_cavity.ipnode 03-May-2004 125k solveEjection.com 03-May-2004 5.0k solveIVC.com 03-May-2004 3.6k solveIVR.com 03-May-2004 4.2k solveInflation.com 03-May-2004 3.0k solveRecoil.com 03-May-2004 3.5k solveReduced.com 03-May-2004 2.1k

Download the entire example:

Name                           Modified     Size

examples_8_81_811_811e.tar.gz 17-Aug-2010 204k

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

Input last modified: Mon Aug 16 11:22:36 2010


CMISS Help / Examples / 8 / 81 / 811 / 811e