Deformation of the porcine ventricular model
$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";
#
# LV free wall local refinement
# !!! only for inflation
$REFINE=$FALSE;
$REFINE_INFARCT=$FALSE;
#
# Output
#
if (!$REFINE || !$REFINE_INFARCT)
{
$OUTPUT_DIR = "output_".$CELL_MODEL."/";
} else {
$OUTPUT_DIR = "output_refined/";
}
if( ! -d ${OUTPUT_DIR})
{
mkdir ${OUTPUT_DIR};
}
#
# Solve phases
#
$SOLVE_INFLATION = $TRUE;
#
$SOLVE_ACTIVE_CONTRACTION = $TRUE;
#
$SOLVE_EJECTION = $TRUE;
#
if ( $SOLVE_ACTIVE_CONTRACTION )
{
$CAVITY_REGIONS = $TRUE;
} else {
$CAVITY_REGIONS = $FALSE;
}
$GRID_BASED_MECHANICS_PROPERTIES = $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;
#
read commands;example_811b;
#
# 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;
#
# 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 4..13,28,29,30,31,35,39,40,41,42,53,54,57,59,61..62,65..69 as LV_FREE_WALL_INITIAL;
#
fem group element all_elements as ALL_ELEMENTS;
#
set echo;
if($REFINE){
fem refine Xi 1 element 4 at 0.7;
fem refine Xi 1 element 5 at 0.7;
fem refine Xi 1 element 6 at 0.7;
fem refine Xi 1 element 7 at 0.7;
#
fem refine Xi 3 element 10 at 0.5;
fem refine Xi 3 element 11 at 0.5;
fem refine Xi 3 element 29 at 0.5;
fem refine Xi 3 element 28 at 0.5;
fem refine Xi 3 element 30 at 0.5;
fem refine Xi 3 element 12 at 0.5;
#
fem refine Xi 3 element 5 at 0.5;
fem refine Xi 3 element 6 at 0.5;
fem refine Xi 3 element 7 at 0.5;
#
fem group nodes 158..161,166,167,168..183 as REMAINDER_REFINED;
fem group nodes REMAINDER_INITIAL,REMAINDER_REFINED as REMAINDER;
#
fem group elements 89..92 as LV_FREE_WALL_REFINED;
fem group elements LV_FREE_WALL_INITIAL,LV_FREE_WALL_REFINED as LV_FREE_WALL;
#
fem export nodes;${OUTPUT_DIR}."MRP01_refined" as MRP01 region all;
fem export elements;${OUTPUT_DIR}."MRP01_refined" as MRP01 region $WALL;
fem export elements;${OUTPUT_DIR}."MRP01_base_refined" as base_skel region $BASE;
} elsif ($REFINE_INFARCT) {
fem refine Xi 1 element 4,5,6,7 at 0.7;
#
fem refine Xi 3 element 10,11,29,28,30,12 at 0.5;
#
fem refine Xi 3 element 5,6,7 at 0.5;
#
fem refine Xi 1 element 10,11,12,93,94,98 at 0.5;
#fem group nodes 1..8,10,14,17..20,25..28,30..34,39..40,43..44,47..48,51..52,55..80,85..94 as REMAINDER_REFINED1;
#fem group nodes 99..100,103..142,158,160,166..183,185..187,189..191,193,195 as REMAINDER_REFINED2;
#fem group elements 102..104 as LV_FREE_WALL_REFINED;
#
fem refine Xi 1 element 10,11,12,93,94,98 at 0.5;
#fem group nodes 1..8,10,14,17..20,25..28,30..34,39..40,43..44,47..48,51..52,55..80,85..94 as REMAINDER_REFINED1;
#fem group nodes 99..100,103..142,158,160,166..183,185..187,189..191,193,195,198..199,201..203,207 as REMAINDER_REFINED2;
#fem group elements 102..104,108..110 as LV_FREE_WALL_REFINED;
#
fem refine Xi 1 element 102,103,104,105,106,107 at 0.5;
fem group nodes 1..8,10,14,17..20,25..28,30..34,39..40,43..44,47..48,51..52,55..80,85..94,99..100,103..142 as REMAINDER_REFINED1;
fem group nodes 158,160,166..183,185..187,189..191,193,195,197..199,201..203,205,207,209..211,213..215,217,219 as REMAINDER_REFINED2;
fem group elements 102..104,108..110,114..116 as LV_FREE_WALL_REFINED;
fem group nodes REMAINDER_REFINED1,REMAINDER_REFINED2 as REMAINDER;
# fem group nodes REMAINDER_INITIAL,REMAINDER_REFINED1,REMAINDER_REFINED2 as REMAINDER;
#
# fem group elements 102..104,108..110,114..116,120..122,126..128,132..134,138..140 as LV_FREE_WALL_REFINED;
fem group elements LV_FREE_WALL_INITIAL,LV_FREE_WALL_REFINED as LV_FREE_WALL;
#
fem export nodes;${OUTPUT_DIR}."MRP01_refined" as MRP01 region all;
fem export elements;${OUTPUT_DIR}."MRP01_refined" as MRP01 region $WALL;
fem export elements;${OUTPUT_DIR}."MRP01_base_refined" as base_skel region $BASE;
fem group elements 103,106 as INFARCT_ELEMENTS;
#fem define fibre;w;MRP01_refined;
fem define fibre;r;MRP01_refined;
} else {
fem group nodes REMAINDER_INITIAL as REMAINDER;
fem group elements LV_FREE_WALL_INITIAL as LV_FREE_WALL;
}
if ($GRID_BASED_MECHANICS_PROPERTIES) {
set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n";
print "\033[0;30;42m Grid \033[0m\n";
print "\033[0;30;42m ============================================================ \033[0m\n";
#==========================================================================================================
set echo;
#
fem define grid;r;MRP01_refined class $MECHANICS region $WALL;
fem update grid geometry class $CELLULAR region $WALL;
#
fem export element;${OUTPUT_DIR}."MRP01_grid" grid_numbers as grid class $CELLULAR;
#
#fem group grid grid 1 as NORMAL_CELLS;
#fem group grid within elements INFARCT_ELEMENTS as INFARCT_CELLS region $WALL;
#fem group grid within elements ISO as ISO_CELLS region $WALL;
#fem group grid;w;MRP01_refined region $WALL;
fem group grid;r;MRP01_refined region $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;finelas_tch_incomp lockregion $WALL;
fem define mapping;r;MRP01;
if ($GRID_BASED_MECHANICS_PROPERTIES) {
fem define equation;r;cellular class $CELLULAR region $WALL; # Cellular equations (set, not used)
}
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;
if (!$GRID_BASED_MECHANICS_PROPERTIES) {
fem define material;r;polezero region $WALL;
} else {
fem define cell;r;polezero class $CELLULAR region $WALL;# Cell variant definitions
# fem define material;r;polezero_infarct cell class $CELLULAR region $WALL;
fem define material;r;polezero_infarct_edited cell class $CELLULAR region $WALL;
fem define material;r;polezero_grid region $WALL; # Mechanics materials set from cells
}
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;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;lu region $WALL;
set echo;
$ITERS = 15;
$TOL = 0.001;
if ( $SOLVE_INFLATION ) {
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n";
print "\033[0;30;42m Solve Inflation \033[0m\n";
print "\033[0;30;42m ============================================================ \033[0m\n";
#==========================================================================================================
print "\033[0;30;43m Increase cavity pressures incrementally to simulate diastole \033[0m\n";
#
set out;${OUTPUT_DIR}steps_inflate on;
$NAME1 = "full_heart_press_";
$MAXIMUM_INCREM = 9;
$PRESS = 0.0;
for $i ( 0..$MAXIMUM_INCREM )
{
if ($i < 1)
{
$INCREM = 0.0;
} else {
$INCREM = 0.1;
# $INCREM = 0.3333;
# $INCREM = 0.16666;
}
#
$PRESS = $PRESS + $INCREM;
if ($PRESS < 1.3)
{
$PRESS_mmHg = sprintf "%1.1f",${PRESS}*7.500637554192106;
} else {
$PRESS_mmHg = sprintf "%2.1f",${PRESS}*7.500637554192106;
}
$FILE = ${NAME1}.${PRESS};
# $FILE = ${NAME1}.${PRESS_mmHg};
fem solve increment $INCREM iter $ITERS error $TOL class $MECHANICS;
fem define initial;w;${OUTPUT_DIR}.${FILE};
fem export nodes;${OUTPUT_DIR}.${FILE} as heart;
fem export nodes;${OUTPUT_DIR}${FILE}_def field as heart;
fem export elements;${OUTPUT_DIR}.${FILE} as heart;
fem export elements;${OUTPUT_DIR}${FILE}_def field as wall;
#
fem update gauss strain fibre components region $WALL;
fem export gauss;${OUTPUT_DIR}.${FILE}."_gauss_strain" yg as gauss_strain;
fem update gauss stress fibre components region $WALL;
fem export gauss;${OUTPUT_DIR}.${FILE}."_gauss_stress" yg as gauss_stress;
if ($GRID_BASED_MECHANICS_PROPERTIES) {
fem export element;${OUTPUT_DIR}.${FILE}."_grid_fields" field cell grid_numbers as grid class $CELLULAR;
fem export gauss;${OUTPUT_DIR}.${FILE}."_gauss_material" parameters class $MECHANICS region $WALL;
}
#
for $ne ( 1..80 )
{
if ($ne != 56 ){
$xi1 = 0.5;
$xi2 = 0.5;
MySubs::Create_File_strain_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_gauss",$ne);
MySubs::Create_File_strain_gauss("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_gauss",$ne);
MySubs::Create_File_stress_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_gauss",$ne);
#
MySubs::Create_File_strain_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_strain_xi("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_stress_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_xi",$ne,$xi1,$xi2);
}
}
#
if ($REFINE)
{
for $ne ( 89..101 )
{
$xi1 = 0.5;
$xi2 = 0.5;
MySubs::Create_File_strain_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_gauss",$ne);
MySubs::Create_File_strain_gauss("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_gauss",$ne);
MySubs::Create_File_stress_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_gauss",$ne);
#
MySubs::Create_File_strain_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_strain_xi("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_stress_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_xi",$ne,$xi1,$xi2);
}
MySubs::Combine_Element_Files(${OUTPUT_DIR}.${FILE}."_ne_11_fibre_strain_gauss",${OUTPUT_DIR}.${FILE}."_ne_94_fibre_strain_gauss",${OUTPUT_DIR}.${FILE}."_ne_11.94_fibre_strain_gauss")
#
MySubs::Combine_Element_Files(${OUTPUT_DIR}.${FILE}."_ne_11_fibre_strain_xi",${OUTPUT_DIR}.${FILE}."_ne_94_fibre_strain_xi",${OUTPUT_DIR}.${FILE}."_ne_11.94_fibre_strain_xi")
}
if ($REFINE_INFARCT)
{
for $ne ( 89..119 )
{
$xi1 = 0.5;
$xi2 = 0.5;
MySubs::Create_File_strain_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_gauss",$ne);
MySubs::Create_File_strain_gauss("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_gauss",$ne);
MySubs::Create_File_stress_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_gauss",$ne);
#
MySubs::Create_File_strain_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_strain_xi("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_stress_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_xi",$ne,$xi1,$xi2);
}
MySubs::Combine_Element_Files(${OUTPUT_DIR}.${FILE}."_ne_11_fibre_strain_gauss",${OUTPUT_DIR}.${FILE}."_ne_94_fibre_strain_gauss",${OUTPUT_DIR}.${FILE}."_ne_11.94_fibre_strain_gauss")
#
MySubs::Combine_Element_Files(${OUTPUT_DIR}.${FILE}."_ne_11_fibre_strain_xi",${OUTPUT_DIR}.${FILE}."_ne_94_fibre_strain_xi",${OUTPUT_DIR}.${FILE}."_ne_11.94_fibre_strain_xi")
}
}
#
if (!$REFINE)
{
# MySubs::Create_Slice_EXR_FILE();
# MySubs::Create_Section_EXR_FILE();
}
if (!$REFINE_INFARCT)
{
# MySubs::Create_Slice_EXR_FILE();
# MySubs::Create_Section_EXR_FILE();
}
} 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;
fem define initial;r;${OUTPUT_DIR}.${FILE} region $WALL;
# fem define solve;r;newton region $WALL class $MECHANICS;
fem define solve;r;lu region $WALL;
fem solve increment 0.0 iter 20 error $TOL ;
}
set echo;
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;MRP01; # Redefine bases to calculate scale factors
fem define;add bases;r;3Linear-2PressAuxXi3_4x4x4Gauss;
#
# Define geometry for LV and RV cavities
#
fem define region;r;cavities;
#
fem define coord;r;lv_cavity region $LV_CAVITY;
fem define node;r;lv_cavity region $LV_CAVITY;
fem define element;r;lv_cavity region $LV_CAVITY;
#
fem define coord;r;rv_cavity region $RV_CAVITY;
fem define node;r;rv_cavity region $RV_CAVITY;
fem define element;r;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;MRP01_readse; # Redefine bases to read in scale factors
fem define;add bases;r;3Linear-2PressAuxXi3_4x4x4Gauss;
fem define line;r;MRP01;
$LV_PRESS = $PRESS;
$RV_PRESS = $PRESS*0.2;
}
set echo;
set echo;
if ( $SOLVE_ACTIVE_CONTRACTION ) {
if($CELL_MODEL eq "FM_GAUSS") {
fem define material;r;polezero_active region $WALL;
fem define active;r;active0_0 region $WALL;
# Isovolumically contract the ventricles by incrementally
# increasing intracellular calcium concentration and solving.
# If the cavity volume decreases, increment the cavity pressure
# and resolve for the same activation level.
#
set output;${OUTPUT_DIR}steps_isovolumic_contract on;
#
$NAME1 = "active0_";
$NAME2 = "fullheart_def_active_";
$ITER = 20;
$ERROR_TOLERANCE = 0.001;
$CA = 0;
$MAX_CA = 24;
#
while ($CA <= $MAX_CA )
{
print "\033[0;30;43m Calcium Level = ${CA} \033[0m\n";
$FILE=${NAME2}.${CA};
#
fem define active;r;${NAME1}.${CA} region $WALL;
fem solve increment 0.0 iter $ITER error $ERROR_TOLERANCE;
fem define initial;w;${OUTPUT_DIR}.${FILE} region $WALL
#
fem list elements deformed total cavity_volume region $LV_CAVITY;
fem list elements deformed total cavity_volume region $RV_CAVITY;
print "\nLV cavity volume = ${LV_CAVITY_VOLUME}\n";
print "RV cavity volume = ${RV_CAVITY_VOLUME}\n\n";
#
if ($LV_CAVITY_VOLUME < ($DIASTOLIC_LV_CAVITY_VOLUME-$DIASTOLIC_LV_CAVITY_VOLUME*0.01) )
{
$INCREM = 0.1;
fem update pressure boundary elements LV_FREE_WALL auxillary 1 increment $INCREM;
fem update pressure boundary elements SEPTUM auxillary 1 increment $INCREM;
$LV_PRESS = $LV_PRESS+$INCREM;
print "\nIncreasing LV cavity pressure by $INCREM to $LV_PRESS kPa\n\n";
} elsif ($RV_CAVITY_VOLUME < ($DIASTOLIC_RV_CAVITY_VOLUME-$DIASTOLIC_RV_CAVITY_VOLUME*0.01) )
{
$INCREM = 0.1;
fem update pressure boundary elements SEPTUM auxillary 2 increment $INCREM;
fem update pressure boundary elements RV_FREE_WALL auxillary 1 increment $INCREM;
$RV_PRESS = $RV_PRESS+$INCREM*0.2;
print "\nIncreasing RV cavity pressure by $INCREM to $RV_PRESS kPa\n\n";
} else {
$CA++;
}
#
fem export nodes;${OUTPUT_DIR}.${FILE} field as heart region $WALL;
fem export elements;${OUTPUT_DIR}.${FILE} field as wall region $WALL;
fem export nodes;${OUTPUT_DIR}."lv_cavity_ref_".${CA} as lv_cavity region all;
fem export elements;${OUTPUT_DIR}."lv_cavity_ref_".${CA} as lv_cavity region $LV_CAVITY;
for $ne ( 1..80 )
{
if ($ne != 56 ){
$xi1 = 0.5;
$xi2 = 0.5;
MySubs::Create_File_strain_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_gauss",$ne);
MySubs::Create_File_strain_gauss("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_gauss",$ne);
MySubs::Create_File_stress_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_gauss",$ne);
#
MySubs::Create_File_strain_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_strain_xi("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_stress_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_xi",$ne,$xi1,$xi2);
}
}
}
#
print "Total LV cavity pressure = $LV_PRESS\n";
print "Total RV cavity pressure = $RV_PRESS\n";
print "End-diastolic LV cavity volume = ${DIASTOLIC_LV_CAVITY_VOLUME}\n";
print "End-isovolumic contraction LV cavity volume = ${LV_CAVITY_VOLUME}\n\n";
print "End-diastolic RV cavity volume = ${DIASTOLIC_RV_CAVITY_VOLUME}\n";
print "End-isovolumic contraction RV cavity volume = ${RV_CAVITY_VOLUME}\n\n";
#
} else {
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n";
print "\033[0;30;42m Solving Active Contraction \033[0m\n";
print "\033[0;30;42m ============================================================ \033[0m\n";
#==========================================================================================================
fem define material;r;polezero_active region $WALL;
fem define active;r;$CELL_MODEL
#
# Set up the HMT cellular model - used to calculate the active tension
# Define a 3x3x3 grid scheme for use in all elements and calculate
# their spatial location and metrics
#
fem define grid;r;cell-3x3x3 gauss class $CELLULAR;
fem update grid geometry;
fem update grid metric;
#
fem group grid element ALL_ELEMENTS as ALL_GRID_POINTS;
#
# Define the cardiac cell mechanics model
#
fem define equation;r;cell-$CELL_MODEL class $CELLULAR;
fem define material;r;cell-FM class $CELLULAR; # unused
#
# Define the model parameters, there are no spatially varying cellular parameters
#
fem define cell;r;cell-$CELL_MODEL class $CELLULAR;
fem define material;r;cell-$CELL_MODEL cell class $CELLULAR;
fem update grid material class $CELLULAR;
#
# Define the [Ca]i transient
#
$cai = "cai";
fem define time;r;cai-$CELL_MODEL;
#
# 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;cell-$CELL_MODEL class $CELLULAR;
#
# Define a simple Euler interagtion for the cellular ODE's
#
fem define solve;r;cell-$CELL_MODEL class $CELLULAR;
#
# Cai index for grid (YQS) variables
#
fem inquire cell_variable Cai return_variables CAI_GRIDARRAY,CAI_GRIDIDX;
#
# Fibre extension ratio index for grid (YQS) and mechanics variables
#
fem inquire cell_variable ExtensionRatio return_variables L_GRIDARRAY,L_GRIDIDX;
$L_MECHIDX = 1;
#
# Tension index for grid (YQS) and Gauss (YG) variables
#
fem inquire cell_variable IsometricTension return_variables T_GRIDARRAY,T_GRIDIDX;
$T_GAUSSARRAY = "YG";
$T_GAUSSIDX = 1;
#
# Open a history file to store the cellular extension ratios
#
fem open history;${OUTPUT_DIR}."cell-ext" write variables yqs niqslist $L_GRIDIDX,$T_GRIDIDX,$CAI_GRIDIDX binary class $CELLULAR unit 21;
#
fem update gauss strain fibre components;
fem update grid extension_ratio component $L_MECHIDX $L_GRIDARRAY $L_GRIDIDX class $MECHANICS;
fem solve class $CELLULAR to 0;
fem update gauss gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_GAUSSIDX;
fem solve class $MECHANICS increment 0.0 iterate 10 error 0.001;
#
#
# Loop through time, solving the coupled active dynamic HTM with coupled mecahnics
#
$RESOLVE=$FALSE;
$PRESSURE_INCREMENT=$FALSE;
$TIME = 0.0;
if($CELL_MODEL eq "HMT") {
$TIME_MAX=120;
} elsif ($CELL_MODEL eq "FM") {
$TIME_MAX=30; # Time is actually activation level
}
$dt=1;
while ($TIME <= $TIME_MAX) {
$iter = 0;
$CONVERGED = 0;
while ($iter < 1) {
$iter++;
print "\nIteration: $iter\n\n";
if(!$PRESSURE_INCREMENT){
fem update gauss strain fibre components;
fem update grid extension_ratio component $L_MECHIDX $L_GRIDARRAY $L_GRIDIDX class $MECHANICS;
}
$local_time = $TIME;
$local_tend = $TIME+$dt;
$FILE = "fullheart_def_active_".${TIME};
#
if($RESOLVE){
$TIME_last=$TIME-$dt;
$FILE_LAST="fullheart_def_active_".${TIME_last};
fem read matrix;${OUTPUT_DIR}.${FILE_LAST}."_YQ_YQS" matrices YQ,YQS binary;
}
#
fem solve class $CELLULAR from $local_time to $local_tend;
fem update gauss gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_GAUSSIDX ;
fem solve increment 0.0 iterate 10 error 0.001;
print "\n\n\n convergence iterations $CONVERGENCE_ITERATIONS \n\n\n";
#
fem list elements deformed total cavity_volume region $LV_CAVITY;
fem list elements deformed total cavity_volume region $RV_CAVITY;
print "\nLV cavity volume = ${LV_CAVITY_VOLUME}\n";
print "RV cavity volume = ${RV_CAVITY_VOLUME}\n\n";
#
# Increase cavity pressures to maintain constant volume
#
if ($CONVERGED) {
$PRESSURE_INCREMENT=$FALSE;
if ($LV_CAVITY_VOLUME < ($DIASTOLIC_LV_CAVITY_VOLUME-$DIASTOLIC_LV_CAVITY_VOLUME*0.01) )
{
$INCREM = 0.1;
fem update pressure boundary elements LV_FREE_WALL auxillary 1 increment $INCREM;
fem update pressure boundary elements SEPTUM auxillary 1 increment $INCREM;
$LV_PRESS = $LV_PRESS+$INCREM;
print "\nIncreasing LV cavity pressure by $INCREM to $LV_PRESS kPa\n\n";
$RESOLVE=$TRUE;
$PRESSURE_INCREMENT=$TRUE;
last;
}
if ($RV_CAVITY_VOLUME < ($DIASTOLIC_RV_CAVITY_VOLUME-$DIASTOLIC_RV_CAVITY_VOLUME*0.01) )
{
$INCREM = 0.1;
fem update pressure boundary elements SEPTUM auxillary 2 increment $INCREM;
fem update pressure boundary elements RV_FREE_WALL auxillary 1 increment $INCREM;
$RV_PRESS = $RV_PRESS+$INCREM*0.2;
print "\nIncreasing RV cavity pressure by $INCREM to $RV_PRESS kPa\n\n";
$RESOLVE=$TRUE;
$PRESSURE_INCREMENT=$TRUE;
last;
}
}
#
# Increase time if all has gone well, otherwise try again
#
if ($CONVERGED) {
if (($CONVERGENCE_ITERATIONS == 1) || ($CELL_MODEL eq "FM" )) {
$LAST_CONVERGED_TIME=$TIME;
$TIME += $dt;
$RESOLVE=$FALSE
last;
} else {
$RESOLVE=$TRUE
last;
}
}
}
if ($CONVERGED == 0) {
print "\n\nCouldn't Converge, Exiting.\n\n\n";
fem export nodes;${OUTPUT_DIR}.${FILE} field as heart region $WALL;
fem export elements;${OUTPUT_DIR}.${FILE} field as wall region $WALL;
fem export nodes;${OUTPUT_DIR}."lv_cavity_ref_".${TIME} as lv_cavity region all;
fem export elements;${OUTPUT_DIR}."lv_cavity_ref_".${TIME} as lv_cavity region $LV_CAVITY;
last;
}
#
fem export nodes;${OUTPUT_DIR}.${FILE} field as heart region $WALL;
fem export elements;${OUTPUT_DIR}.${FILE} field as wall region $WALL;
fem export nodes;${OUTPUT_DIR}."lv_cavity_ref_".${TIME} as lv_cavity region all;
fem export elements;${OUTPUT_DIR}."lv_cavity_ref_".${TIME} as lv_cavity region $LV_CAVITY;
#
# Update grid and get new extension ratio values
#
fem update grid geometry deformed class $CELLULAR from_class $MECHANICS;
fem update grid metric deformed class $CELLULAR from_class $MECHANICS;
fem update grid material class $CELLULAR;
#
if (!$RESOLVE) {
# A converged solution was obtained, so we are done for this time step..
# Store the current cell state for the next set of iterations
fem write matrix;${OUTPUT_DIR}.${FILE}."_YQ_YQS" matrices YQ,YQS binary;
fem write history time $TIME variables yqs binary class $CELLULAR unit 21;
}
}
#
fem close history binary class $CELLULAR unit 21;
#
# Export the history file to signal files
#
fem evaluate electrode;${OUTPUT_DIR}."/cell-ext" history ${OUTPUT_DIR}."/cell-ext" from grid yqs iy $L_GRIDIDX binary class $CELLULAR;
fem evaluate electrode;${OUTPUT_DIR}."/cell-tension" history ${OUTPUT_DIR}."/cell-ext" from grid yqs iy $T_GRIDIDX binary class $CELLULAR;
fem evaluate electrode;${OUTPUT_DIR}."/cell-cai" history ${OUTPUT_DIR}."/cell-ext" from grid yqs iy $CAI_GRIDIDX binary class $CELLULAR;
#
fem define export;r;cell-3x3x3;
#
fem export signal;${OUTPUT_DIR}."/cell-ext" electrode signal ${OUTPUT_DIR}."/cell-ext";
fem export signal;${OUTPUT_DIR}."/cell-tension" electrode signal ${OUTPUT_DIR}."/cell-tension";
fem export signal;${OUTPUT_DIR}."/cell-cai" electrode signal ${OUTPUT_DIR}."/cell-cai";
#
print "Total LV cavity pressure = $LV_PRESS\n";
print "Total RV cavity pressure = $RV_PRESS\n";
print "End-diastolic LV cavity volume = ${DIASTOLIC_LV_CAVITY_VOLUME}\n";
print "End-isovolumic contraction LV cavity volume = ${LV_CAVITY_VOLUME}\n\n";
print "End-diastolic RV cavity volume = ${DIASTOLIC_RV_CAVITY_VOLUME}\n";
print "End-isovolumic contraction RV cavity volume = ${RV_CAVITY_VOLUME}\n\n";
}
}
set echo;
set echo;
if ( $SOLVE_EJECTION ) {
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n";
print "\033[0;30;42m Solving Ejection \033[0m\n";
print "\033[0;30;42m ============================================================ \033[0m\n";
#==========================================================================================================
#
# Decrement the cavity pressures to simulate ejection
#
set output;${OUTPUT_DIR}steps_ejection on;
#
$FILE = "fullheart_def_active_30_press_";
$INCREM = -0.1;
$ITER = 20;
$ERROR_TOLERANCE = 0.001;
$PRESS=$LV_PRESS*10;
while ($PRESS > 0.0 )
{
fem solve increment $INCREM iter $ITER error $ERROR_TOLERANCE;
fem define initial;w;${OUTPUT_DIR}.${FILE}.${PRESS} region $WALL
#
fem list elements deformed total cavity_volume region $LV_CAVITY;
print "\nLV cavity volume = ${LV_CAVITY_VOLUME}\n\n";
#
#
fem export nodes;${OUTPUT_DIR}.${FILE}.${PRESS} field as heart region $WALL;
fem export elements;${OUTPUT_DIR}.${FILE}.${PRESS} field as wall region $WALL;
for $ne ( 1..80 )
{
if ($ne != 56 ){
$xi1 = 0.5;
$xi2 = 0.5;
MySubs::Create_File_strain_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_gauss",$ne);
MySubs::Create_File_strain_gauss("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_gauss",$ne);
MySubs::Create_File_stress_gauss("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_gauss",$ne);
#
MySubs::Create_File_strain_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_strain_xi("wall",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_wall_strain_xi",$ne,$xi1,$xi2);
MySubs::Create_File_stress_xi("fibre",${OUTPUT_DIR}.${FILE}."_ne_".$ne."_fibre_stress_xi",$ne,$xi1,$xi2);
}
}
$PRESS--;
}
set out off;
}
set echo;
quit
Name Modified Size
example_811d.com 11-Sep-2002 33k 3Linear-2PressAuxXi3_4x4x4Gauss.ipbase 11-Sep-2002 1.9k FM.ipacti 11-Sep-2002 258 HMT.ipacti 11-Sep-2002 258 MRP01.exelem 11-Sep-2002 528k MRP01.exnode 11-Sep-2002 340k MRP01.ipbase 11-Sep-2002 11k MRP01.ipcoor 11-Sep-2002 688 MRP01.ipelem 11-Sep-2002 100k MRP01.ipelfb 11-Sep-2002 67k MRP01.ipfibr 11-Sep-2002 408k MRP01.ipinit 05-Dec-2002 31k MRP01.ipline 11-Sep-2002 455k MRP01.ipmap 11-Sep-2002 146k MRP01.ipnode 11-Sep-2002 461k MRP01.ippara 12-Nov-2002 5.9k MRP01.ipregi 11-Sep-2002 142 MRP01_base.exelem 11-Sep-2002 42k MRP01_base.ipelem 11-Sep-2002 12k MRP01_base.ipline 11-Sep-2002 36k MRP01_base.ipnode 11-Sep-2002 121k MRP01_readse.ipbase 11-Sep-2002 11k MRP01_refined.ipgrid 06-Mar-2003 1.5k MRP01_refined.ippara 12-Nov-2002 5.9k MySubs.pm 11-Sep-2002 40k active0_0.ipacti 03-Mar-2004 777 active0_0.opstra 11-Sep-2002 37k active0_1.ipacti 03-Mar-2004 777 active0_1.opstra 11-Sep-2002 37k active0_10.ipacti 03-Mar-2004 777 active0_10.opstra 11-Sep-2002 37k active0_11.ipacti 03-Mar-2004 777 active0_11.opstra 11-Sep-2002 37k active0_12.ipacti 03-Mar-2004 777 active0_12.opstra 11-Sep-2002 37k active0_13.ipacti 03-Mar-2004 777 active0_13.opstra 11-Sep-2002 37k active0_14.ipacti 03-Mar-2004 777 active0_14.opstra 11-Sep-2002 37k active0_15.ipacti 03-Mar-2004 777 active0_15.opstra 11-Sep-2002 37k active0_16.ipacti 03-Mar-2004 777 active0_16.opstra 11-Sep-2002 37k active0_17.ipacti 03-Mar-2004 777 active0_17.opstra 11-Sep-2002 37k active0_18.ipacti 03-Mar-2004 777 active0_18.opstra 11-Sep-2002 37k active0_19.ipacti 03-Mar-2004 777 active0_19.opstra 11-Sep-2002 37k active0_2.ipacti 03-Mar-2004 777 active0_2.opstra 11-Sep-2002 37k active0_20.ipacti 03-Mar-2004 777 active0_20.opstra 11-Sep-2002 37k active0_21.ipacti 03-Mar-2004 777 active0_21.opstra 11-Sep-2002 37k active0_22.ipacti 03-Mar-2004 777 active0_22.opstra 11-Sep-2002 37k active0_23.ipacti 03-Mar-2004 777 active0_23.opstra 11-Sep-2002 37k active0_24.ipacti 03-Mar-2004 777 active0_24.opstra 11-Sep-2002 37k active0_25.ipacti 03-Mar-2004 777 active0_26.ipacti 03-Mar-2004 777 active0_27.ipacti 03-Mar-2004 777 active0_28.ipacti 03-Mar-2004 777 active0_29.ipacti 03-Mar-2004 777 active0_3.ipacti 03-Mar-2004 777 active0_3.opstra 11-Sep-2002 37k active0_30.ipacti 03-Mar-2004 777 active0_4.ipacti 03-Mar-2004 777 active0_4.opstra 11-Sep-2002 37k active0_5.ipacti 03-Mar-2004 777 active0_5.opstra 11-Sep-2002 37k active0_6.ipacti 03-Mar-2004 777 active0_6.opstra 11-Sep-2002 37k active0_7.ipacti 03-Mar-2004 777 active0_7.opstra 11-Sep-2002 37k active0_8.ipacti 03-Mar-2004 777 active0_8.opstra 11-Sep-2002 37k active0_9.ipacti 03-Mar-2004 777 active0_9.opstra 11-Sep-2002 37k cai-FM.iptime 11-Sep-2002 956 cai-HMT.iptime 11-Sep-2002 120k cavities.ipregi 11-Sep-2002 130 cell-3x3x3.ipexpo 11-Sep-2002 661 cell-3x3x3.ipgrid 06-Mar-2003 610 cell-9x9x9.ipexpo 11-Sep-2002 658 cell-9x9x9.ipgrid 06-Mar-2003 610 cell-FM.ipcell 11-Sep-2002 1.4k cell-FM.ipequa 26-May-2003 1.4k cell-FM.ipinit 11-Sep-2002 545 cell-FM.ipmatc 11-Sep-2002 360 cell-FM.ipmate 11-Sep-2002 1.8k cell-FM.ipsolv 13-Apr-2007 1.8k cell-HMT.ipcell 11-Sep-2002 3.8k cell-HMT.ipequa 26-May-2003 1.2k cell-HMT.ipinit 11-Sep-2002 548 cell-HMT.ipmatc 11-Sep-2002 360 cell-HMT.ipmate 11-Sep-2002 1.8k cell-HMT.ipsolv 13-Apr-2007 1.8k cell.ipcell 11-Sep-2002 3.3k cellular.ipequa 26-May-2003 1.0k example_811b.com 11-Sep-2002 1.4k finelas_tch_incomp.ipequa 02-May-2004 2.1k formats.com 11-Sep-2002 726 lu.ipsolv 16-Aug-2010 2.8k lu.ipsolv.old 13-Apr-2007 2.6k lv_cavity.ipcoor 11-Sep-2002 688 lv_cavity.ipelem 11-Sep-2002 42k lv_cavity.ipnode 11-Sep-2002 158k newton.ipsolv 16-Aug-2010 2.7k newton.ipsolv.old 13-Apr-2007 2.5k polezero.cell.xml 11-Sep-2002 2.5k polezero.ipcell 11-Sep-2002 3.1k polezero.ipmate 11-Sep-2002 6.7k polezero.xml 11-Sep-2002 4.9k polezero_active.ipmate 11-Sep-2002 6.7k polezero_grid.ipmate 11-Sep-2002 6.4k polezero_infarct.ipmatc 11-Sep-2002 2.3k polezero_infarct_edited.ipmatc 11-Sep-2002 8.1M rv_cavity.ipcoor 11-Sep-2002 688 rv_cavity.ipelem 11-Sep-2002 12k rv_cavity.ipnode 11-Sep-2002 125k view_contract.com 11-Sep-2002 5.3k
Name Modified Size
examples_8_81_811_811d.tar.gz 17-Aug-2010 807k
Html last generated: Sun Mar 6 05:50:33 2016
Input last modified: Mon Aug 16 11:22:36 2010