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.
# 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;
}
#Testing example 722 FEM update material optimise; FEM li material constitutive;
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 13-Sep-2001 80 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
Name Modified Size
examples_7_72_722.tar.gz 17-Aug-2010 11k
| Status | Tested | Real time (s) | |
| i686-linux | |||
| cm | Failure | Thu Jan 10 02:00:46 2013 | 155 |
| last break | Sun Apr 16 02:10:00 2006 | 210 | |
| last success | Sun Nov 9 01:56:00 2003 | 187 | |
| 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 |
| i686-linux | |||
| Success | cm: | 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. |
| i686-linux | |||
| Failure | cm: | 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. |
Html last generated: Thu Jan 10 05:53:51 2013
Input last modified: Mon Aug 16 11:21:48 2010