Example 722: Holmes constitutive law constants optimisation, single prolate element

This example optimises the constants C1 & C2 in the Holmes constitutive law for infarcted tissue. The constants are optimised to minimise the the sum of squares difference between the experimental strains and the model strains, at three transmural locations and four different pressures though the inflation phase. The experimental values used in this example are for normal myocardial tissue.


The comfile run by this example is as follows:

# Example 722

$TRUE = 1;  # Initialise logicals
$FALSE = 0;
$pressure = "time"; # A hack to use time variables to represent experimental pressure-strain data

#$RESIDUALS = "E11";          # Pick which strains to use for residuals 
$RESIDUALS = "E11 E22 E33";   # Be sure to set the correct number of residuals in the ipopti file

if($RESOLVE ne "TRUE") {      # Don't need to redo these commands when resolving during optimisation
fem;
fem define parameters;r;profib;example;
fem define coordinate;r;profib;example;
fem define base;r;profib;example;          
}

fem define node;r;profib;example;           
fem define element;r;profib;example;   
fem define fibre;r;profib;example;
fem define element;r;profib;example fibre;   
fem define equation;r;profib;example;  
fem define material;r;profib;example;
       
if($RESOLVE eq "TRUE") {         # Need to update the material parameters to those being 
  fem update material optimise;  # evaluated by the optimiser
}

set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n"; 
print "\033[0;30;42m                           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"; 
fem define initial;r;profib;example; 
fem define solve;r;profib;example;     

fem solve increment 0.5;   

fem solve increment 0.5;                         
if($RESOLVE eq "TRUE") {
		if($RESIDUALS eq "E11") {  # Residuals are calculated from the values loaded into YP by the following
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 1;  #load yp(1,7) & #load yp(1,8)
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 2;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 3;   
				fem evaluate $pressure variable circ_e11_normal_depth_5 $pressure 7.5 to YP 1;   # 1.0 kPa = 7.5 mmHg
				fem evaluate $pressure variable circ_e11_normal_depth_35 $pressure 7.5 to YP 2;   
				fem evaluate $pressure variable circ_e11_normal_depth_65 $pressure 7.5 to YP 3;   
		} elsif ($RESIDUALS eq "E11 E22 E33"){
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 1;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 2;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 3;   
				fem evaluate $pressure variable circ_e11_normal_depth_5 $pressure 7.5 to YP 1;   
				fem evaluate $pressure variable circ_e11_normal_depth_35 $pressure 7.5 to YP 2;   
				fem evaluate $pressure variable circ_e11_normal_depth_65 $pressure 7.5 to YP 3;   
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 4;   
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 5;   
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 6;   
				fem evaluate $pressure variable circ_e22_normal_depth_5 $pressure 7.5 to YP 4;   
				fem evaluate $pressure variable circ_e22_normal_depth_35 $pressure 7.5 to YP 5;   
				fem evaluate $pressure variable circ_e22_normal_depth_65 $pressure 7.5 to YP 6;   
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 7;   
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 8;   
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 9;   
				fem evaluate $pressure variable circ_e33_normal_depth_5 $pressure 7.5 to YP 7;   
				fem evaluate $pressure variable circ_e33_normal_depth_35 $pressure 7.5 to YP 8;   
				fem evaluate $pressure variable circ_e33_normal_depth_65 $pressure 7.5 to YP 9;  
		} 
}

fem solve increment 0.5;        
if($RESOLVE eq "TRUE") {
		if($RESIDUALS eq "E11") {
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 4;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 5;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 6;   
				fem evaluate $pressure variable circ_e11_normal_depth_5 $pressure 11.25 to YP 4; # 1.5 kPa = 11.25 mmHg
				fem evaluate $pressure variable circ_e11_normal_depth_35 $pressure 11.25 to YP 5; 
				fem evaluate $pressure variable circ_e11_normal_depth_65 $pressure 11.25 to YP 6; 
		} elsif ($RESIDUALS eq "E11 E22 E33"){
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 10;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 11;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 12;   
				fem evaluate $pressure variable circ_e11_normal_depth_5 $pressure 11.25 to YP 10; 
				fem evaluate $pressure variable circ_e11_normal_depth_35 $pressure 11.25 to YP 11; 
				fem evaluate $pressure variable circ_e11_normal_depth_65 $pressure 11.25 to YP 12; 
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 13;   
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 14;   
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 15;   
				fem evaluate $pressure variable circ_e22_normal_depth_5 $pressure 11.25 to YP 13; 
				fem evaluate $pressure variable circ_e22_normal_depth_35 $pressure 11.25 to YP 14; 
				fem evaluate $pressure variable circ_e22_normal_depth_65 $pressure 11.25 to YP 15; 
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 16;   
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 17;   
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 18;   
				fem evaluate $pressure variable circ_e33_normal_depth_5 $pressure 11.25 to YP 16; 
				fem evaluate $pressure variable circ_e33_normal_depth_35 $pressure 11.25 to YP 17; 
				fem evaluate $pressure variable circ_e33_normal_depth_65 $pressure 11.25 to YP 18; 
		} 
}

fem solve increment 0.5;        
if($RESOLVE eq "TRUE") {
		if($RESIDUALS eq "E11") {
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 7;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 8;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 9;   
				fem evaluate $pressure variable circ_e11_normal_depth_5 $pressure 15 to YP 7;    # 2.0 kPa = 15 mmHg
				fem evaluate $pressure variable circ_e11_normal_depth_35 $pressure 15 to YP 8;    
				fem evaluate $pressure variable circ_e11_normal_depth_65 $pressure 15 to YP 9;    
		} elsif ($RESIDUALS eq "E11 E22 E33"){
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 19;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 20;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 21;   
				fem evaluate $pressure variable circ_e11_normal_depth_5 $pressure 15 to YP 19;    
				fem evaluate $pressure variable circ_e11_normal_depth_35 $pressure 15 to YP 20;    
				fem evaluate $pressure variable circ_e11_normal_depth_65 $pressure 15 to YP 21;    
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 22;   
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 23;   
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 24;   
				fem evaluate $pressure variable circ_e22_normal_depth_5 $pressure 15 to YP 22;    
				fem evaluate $pressure variable circ_e22_normal_depth_35 $pressure 15 to YP 23;    
				fem evaluate $pressure variable circ_e22_normal_depth_65 $pressure 15 to YP 24;    
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 25;   
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 26;   
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 27;   
				fem evaluate $pressure variable circ_e33_normal_depth_5 $pressure 15 to YP 25;    
				fem evaluate $pressure variable circ_e33_normal_depth_35 $pressure 15 to YP 26;    
				fem evaluate $pressure variable circ_e33_normal_depth_65 $pressure 15 to YP 27;    
		} 
}

fem solve increment 0.5;        
if($RESOLVE eq "TRUE") {
		if($RESIDUALS eq "E11") {
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 10;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 11;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 12;   
				fem evaluate $pressure variable circ_e11_normal_depth_5 $pressure 18.75 to YP 10; # 2.5 kPa = 18.75 mmHg
				fem evaluate $pressure variable circ_e11_normal_depth_35 $pressure 18.75 to YP 11; 
				fem evaluate $pressure variable circ_e11_normal_depth_65 $pressure 18.75 to YP 12; 
		} elsif ($RESIDUALS eq "E11 E22 E33"){
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 28;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 29;   
				fem li strain e11 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 30;   
				fem evaluate $pressure variable circ_e11_normal_depth_5 $pressure 18.75 to YP 28; 
				fem evaluate $pressure variable circ_e11_normal_depth_35 $pressure 18.75 to YP 29; 
				fem evaluate $pressure variable circ_e11_normal_depth_65 $pressure 18.75 to YP 30; 
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 31;   
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 32;   
				fem li strain e22 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 33;   
				fem evaluate $pressure variable circ_e22_normal_depth_5 $pressure 18.75 to YP 31; 
				fem evaluate $pressure variable circ_e22_normal_depth_35 $pressure 18.75 to YP 32; 
				fem evaluate $pressure variable circ_e22_normal_depth_65 $pressure 18.75 to YP 33; 
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.5 element 1 wall to YP 34;   
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.65 element 1 wall to YP 35;   
				fem li strain e33 xi_point xi_1 0.5 xi_2 0.5 xi_3 0.35 element 1 wall to YP 36;   
				fem evaluate $pressure variable circ_e33_normal_depth_5 $pressure 18.75 to YP 34; 
				fem evaluate $pressure variable circ_e33_normal_depth_35 $pressure 18.75 to YP 35; 
				fem evaluate $pressure variable circ_e33_normal_depth_65 $pressure 18.75 to YP 36; 
		} 
}

if($RESOLVE eq "TRUE") {
  $ITERATION=$ITERATION + 1;
  print "\033[0;30;43m                    Iteration = ${ITERATION}                  \033[0m\n"; 
} else {
  print "\033[0;30;43m               Inflated solution is converged OK!             \033[0m\n"; 
}

set echo;

if($RESOLVE ne "TRUE") {

set echo;
#==========================================================================================================
print "\033[0;30;42m ============================================================ \033[0m\n"; 
print "\033[0;30;42m               Material Paramter Optimisation                 \033[0m\n"; 
print "\033[0;30;42m ============================================================ \033[0m\n"; 
#==========================================================================================================
set echo;
fem define $pressure;r;wall_e11_normal;example;      # Curves defining the experimental pressure-strain data
fem define $pressure;r;wall_e22_normal;example;      # at various depths for normal (healthy) tissue
fem define $pressure;r;wall_e33_normal;example;
fem define $pressure;r;wall_e12_normal;example;
fem define $pressure;r;wall_e13_normal;example;
fem define $pressure;r;wall_e23_normal;example;

$ITERATION = 1;
$RESOLVE = "TRUE";
#fem define optimise;r;C1_C2_e11;example;           # Pick which set of residuals to use
fem define optimise;r;C1_C2_e11e22e33;example;       
optimise; 
}


Additional testing commands:

#Testing example 722
FEM update material optimise;
FEM li material constitutive;

Files used by this example are:

Name                    Modified     Size

example_722.com 13-Sep-2001 11k C1_C2_e11.ipopti 16-Sep-2001 1.0k C1_C2_e11e22e33.ipopti 05-Jul-2005 1.0k profib.ipbase 13-Sep-2001 4.2k profib.ipcoor 13-Sep-2001 707 profib.ipelem 13-Sep-2001 409 profib.ipelfb 13-Sep-2001 349 profib.ipequa 02-May-2004 2.0k profib.ipfibr 13-Sep-2001 1.3k profib.ipinit 17-Dec-2002 1.7k profib.ipmate 13-Oct-2003 4.2k profib.ipnode 13-Sep-2001 1.2k profib.ippara 12-Nov-2002 5.9k profib.ipsolv 16-Aug-2010 2.3k profib.ipsolv.old 13-Apr-2007 2.1k test_output.com 16-Aug-2014 81 wall_e11_normal.iptime 13-Sep-2001 3.0k wall_e12_normal.iptime 13-Sep-2001 3.0k wall_e13_normal.iptime 13-Sep-2001 3.0k wall_e22_normal.iptime 13-Sep-2001 3.0k wall_e23_normal.iptime 13-Sep-2001 3.0k wall_e33_normal.iptime 13-Sep-2001 3.0k

Download the entire example:

Name                      Modified     Size

examples_7_72_722.tar.gz 17-Aug-2014 12k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmFailureSun Mar 6 01:20:17 201673
last breakSun Aug 17 03:01:00 2014168
last successSun Nov 9 01:56:00 2003187
cm-debugFailureSun Mar 6 01:37:22 2016140
last breakMon Aug 18 00:01:00 2014247
mips-irix
cmFailureThu Aug 16 06:11:39 2007602
last breakFri Mar 17 01:39:00 2006608
last successSun Nov 9 03:44:00 2003584
cm64FailureFri Aug 17 05:07:56 2007660
last breakFri Mar 17 03:47:00 2006640
rs6000-aix
cmFailureWed Mar 4 01:35:36 200965
last breakMon Aug 28 10:32:00 200666
last successSun Nov 9 02:53:00 2003105
cm64FailureWed Mar 4 01:36:01 200970
last breakMon Aug 28 10:32:00 200669
x86_64-linux
cmFailureSun Mar 6 00:05:21 201637
last breakSun Aug 17 00:01:00 201447
cm-debugFailureSun Mar 6 00:06:48 201668
last breakMon Aug 18 00:01:00 201489

Testing status by file:


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

Input last modified: Sat Aug 16 10:27:04 2014


CMISS Help / Examples / 7 / 72 / 722