Holmes constitutive law constants optimisation, single prolate element

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:

	Status	Tested	Real time (s)
i686-linux
cm	Failure	Sun Mar 6 01:20:17 2016	73
	last break	Sun Aug 17 03:01:00 2014	168
	last success	Sun Nov 9 01:56:00 2003	187
cm-debug	Failure	Sun Mar 6 01:37:22 2016	140
	last break	Mon Aug 18 00:01:00 2014	247
mips-irix
cm	Failure	Thu Aug 16 06:11:39 2007	602
	last break	Fri Mar 17 01:39:00 2006	608
	last success	Sun Nov 9 03:44:00 2003	584
cm64	Failure	Fri Aug 17 05:07:56 2007	660
	last break	Fri Mar 17 03:47:00 2006	640
rs6000-aix
cm	Failure	Wed Mar 4 01:35:36 2009	65
	last break	Mon Aug 28 10:32:00 2006	66
	last success	Sun Nov 9 02:53:00 2003	105
cm64	Failure	Wed Mar 4 01:36:01 2009	70
	last break	Mon Aug 28 10:32:00 2006	69
x86_64-linux
cm	Failure	Sun Mar 6 00:05:21 2016	37
	last break	Sun Aug 17 00:01:00 2014	47
cm-debug	Failure	Sun Mar 6 00:06:48 2016	68
	last break	Mon Aug 18 00:01:00 2014	89

Testing status by file:

cmiss_test.log

i686-linux
Success	cm:	cmiss_test.log.retain.
Success	cm-debug:	cmiss_test.log.retain.
mips-irix
Success	cm:	cmiss_test.log.retain.
Success	cm64:	cmiss_test.log.retain.
rs6000-aix
Success	cm:	cmiss_test.log.retain.
Success	cm64:	cmiss_test.log.retain.
x86_64-linux
Success	cm:	cmiss_test.log.retain.
Success	cm-debug:	cmiss_test.log.retain.

cmiss_test.out

i686-linux
Failure	cm:	ndiff test: significant differences with generic answer.
Failure	cm-debug:	ndiff test: significant differences with generic answer.
mips-irix
Failure	cm:	ndiff test: significant differences with generic answer.
Failure	cm64:	ndiff test: significant differences with generic answer.
rs6000-aix
Failure	cm:	ndiff test: significant differences with generic answer.
Failure	cm64:	ndiff test: significant differences with generic answer.
x86_64-linux
Failure	cm:	ndiff test: significant differences with generic answer.
Failure	cm-debug:	ndiff test: significant differences with generic answer.

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