Example 811d: Deformation of the porcine ventricular model

Deformation of the porcine ventricular model


The comfile run by this example is as follows:

$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



Files used by this example are:

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

Download the entire example:

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


CMISS Help / Examples / 8 / 81 / 811 / 811d