A simple cylinder with default fibre axes (0deg) initially inflated to 1kPa
and then a cavity volume mesh is coupled in and there is first isovolumic
contraction followed by ejection.
# If the example path is not set, default to current directory
if (! defined $example) {
$example = "./";
}
# Drop off the trailing / in the example path
$chopped = chop $example;
if ($chopped ne "/") {
$example .= $chopped;
}
# Set up some pretty colour highlighting for messages
read com;"$example/formats.com";
# Define the function used to define activity levels
# usage: DefineCalciumLevel(activation_level)
read com;"$example/define_active.com";
# Define the function used to define cavity stiffness
# usage: DefineCavityStiffness(stiffness,region)
# usage: ($newStiffness) = CalculateNewCavityStiffness(currentStiffness)
read com;"$example/define_cavity_stiffness.com";
set echo;
set num_threads 1;
# Region
$WALL = "1";
$CAVITY = "2";
# Classes
$TISSUE = "1";
# Output file names
$name = "cylinder";
$undeformed = $name;
$deformed = $name."_def";
# Solution parameters
$MAX_ITERS = 20;
$TOL = 0.001;
$TOL = 1.0e-6;
$CAVITY_STIFFNESS = 1000.0;
# Do we solve inflation or read it in ?
if (! defined $SOLVE_INFLATION) {
$SOLVE_INFLATION = 1;
}
# Do we isovolumic contraction or read it in ?
if (! defined $SOLVE_ISOVOL) {
$SOLVE_ISOVOL = 1;
}
# Do we solve ejection or read it in ?
if (! defined $SOLVE_EJECTION) {
$SOLVE_EJECTION = 1;
}
# cut down run time for nightly testing
if ($TESTING) {
$SOLVE_INFLATION=0;
$SOLVE_ISOVOL=1;
$SOLVE_EJECTION=0;
$maxCai1 = 0.1;
$dCai = 0.1;
}
$CAVITY_PRESSURE_NY = 2681;
$cylinderInputFilename = "$example/cylinder_large_refined";
$cavityInputFile = "$example/cavity_large_refined";
$cylinderSaveName = "cylinder_large_refined";
$cavitySaveName = "cavity_large_refined";
fem define parameters;r;$example/cylinder;
fem define coor;r;$example/3drc;
# Define all the bases (unit scale factors)
# 1 Geometry
fem define bases;r;$example/3CubicHermite_4x4x4Gauss;
# 2 Faces
fem define ;add bases;r;$example/2CubicHermite_4x4Gauss;
# 3 Fibre angle
fem define ;add bases;r;$example/2Linear-CubicHermite_4x4x4Gauss;
# 4 Imbrication angle
fem define ;add bases;r;$example/3Linear_4x4x4Gauss;
# 5 Sheet angle
fem define ;add bases;r;$example/Linear-2CubicHermite_4x4x4Gauss;
# 6 Faces
fem define ;add bases;r;$example/Linear-CubicHermite_4x4x4Gauss;
# 7 Faces
fem define ;add bases;r;$example/2Linear_4x4Gauss;
# 8 Faces
fem define ;add bases;r;$example/CubicHermite-Linear_4x4x4Gauss;
# 9 Pressure
fem define ;add bases;r;$example/3Linear-2PressAuxXi3_4x4x4Gauss;
# 10 Cavity geometry
fem define ;add bases;r;$example/2CubicHermite-Linear_4x4x4Gauss;
# 11 Cavity pressure
fem define ;add bases;r;$example/PressAuxXi3_4x4x4Gauss;
# Define the wall mesh
fem define nodes;r;$cylinderInputFilename reg $WALL;
fem define elements;r;$cylinderInputFilename reg $WALL;
###fem define fibre;r;$cylinderInputFilename reg $WALL;
fem define fibre;d reg $WALL;
fem define element;r;$cylinderInputFilename fibre reg $WALL;
# Groups used to define boundary conditions
fem group node 2,4 as fixed_x region $WALL;
fem group node 1,3 as fixed_y region $WALL;
fem group node 1..4,21,25,29,33 as fixed_z region $WALL;
fem group element 1..8 as endo region $WALL;
fem group element 9..16 as epi region $WALL;
fem group element 1..16 as wall region $WALL;
# Export initial undefomed geometry
fem export node;$undeformed as cylinder region $WALL;
fem export elem;$undeformed as cylinder region $WALL;
# Define static analysis of 3D finite elasticity
fem define equation;r;$example/finelas_tch_incomp class $TISSUE reg $WALL;
# Want to use the pole-zero law without specifying shear components
fem define material;r;$example/orth_incomp_active class $TISSUE reg $WALL;
# Initialise calcium activation to 0
&DefineCalciumLevel(0);
# Fix the base in-plane and apply a pressure BC to the inner surface
fem define initial;r;$example/cylinder class $TISSUE reg $WALL;
# Solve with Newton-Raphson and LU Decomposition
fem define solve;r;$example/newton class $TISSUE reg $WALL;
# Inflation solution
$dP = 0.5; # kPa
$maxPressure = 1.0; # kPa
$pressure = 0.0;
print "\n";
print "$Bgreen$Fblue $FBnone\n";
print "$Bgreen$Fblue Solving inflation $FBnone\n";
print "$Bgreen$Fblue $FBnone\n";
print "\n";
$file = $cylinderSaveName."_endinflation";
if ($SOLVE_INFLATION) {
read com;$example/solve_inflation;
# Save deformed configurations for end inflation
fem define initial;w;$file class $TISSUE region $WALL;
} else {
fem define initial;r;"$example/$file" class $TISSUE region $WALL;
fem define solve;r;$example/newton class $TISSUE reg $WALL;
#fem solve increment 0.0 iter $MAX_ITERS error $TOL class $TISSUE region $WALL;
fem export node;$file field as cylinder class $TISSUE region $WALL;
fem export elem;$file field as cylinder class $TISSUE region $WALL;
}
print "\n";
print "$Bgreen$Fblue $FBnone\n";
print "$Bgreen$Fblue Finished inflation $FBnone\n";
print "$Bgreen$Fblue $FBnone\n";
print "\n";
#
# Set up cavity region
#
# Check that the current coupled solution is converged (as expected!)
# i.e. check that the sum of squared residuals at the bottom of the file
# is a very small number!!
fem evaluate residuals;bob wrt geom_params;
fem list map;bob solution;
#fem list node;bob0 solution region $WALL;
# Cavity region and geometry
fem define region;r;$example/coupled;
fem define coord;r;$example/cavity region $CAVITY;
fem define node;r;$cavityInputFile region $CAVITY;
fem define element;r;$cavityInputFile region $CAVITY;
# The BC's for the cavity/wall interface nodes must match
fem group node 2,4 as cavity_fixed_x region $CAVITY;
fem group node 1,3 as cavity_fixed_y region $CAVITY;
fem group node 1..4,21,25,29,33 as cavity_fixed_z region $WALL;
fem group node 54 as cavity_fixed_nodes region $CAVITY;
fem group node 53 as cavity_free_node region $CAVITY;
$VALVE_NODE = "53";
$VALVE_NODES = "1..4,21,25,29,33";
# Set up equations, material properties and initial conditions
# for cavity region
fem define equation;r;$example/coupled region $WALL,$CAVITY lock;
&DefineCavityStiffness(1000,$CAVITY);
fem define initial;r;$example/cavity region $CAVITY;
# Define cavity-to-wall coupling and solution procedure information
fem define coupling;r;$example/coupled;
fem define solve;r;$example/coupled_lu_nosparse coupled region $WALL,$CAVITY;
#fem define solve;r;$example/coupled_gmres_nosparse coupled region $WALL,$CAVITY;
# Transfer wall deformation from the inflation problem to
# geometric dependent variables for the cavity region
fem update solution coupled source_region $WALL class $TISSUE;
# Set up the reference state for the cavity region to be the deformed
# configuration from the inflation problem. To set up the
# 'inflated cavity', the x-coord ('in 1' represents the 1st RC coord)
# for the central cavity node ($VALVE_NODE) is set to the average of the
# x-coords for the adjacent interface nodes ($VALVE_NODES - this can
# be a list of nodes)
fem update solution cavity_reference average $VALVE_NODE in 3 node $VALVE_NODES region $CAVITY;
#
# List undeformed (pre-inflation) and inflated cavity volumes
#fem list element total region $CAVITY;
#fem list element deformed total region $CAVITY;
#fem list element deformed total basis 10 variable cavityVolume region $CAVITY;
#print "$Bcyan$Fred Cavity volume = $cavityVolume $FBnone\n";
# Export initial cavity geometry
fem export node;cavity_initial as cavity region all;
fem export elem;cavity_initial as cavity region $CAVITY;
fem export node;cavity_initial_def field as cavity region all;
fem export elem;cavity_initial_def field as cavity region $CAVITY;
# Check that the current coupled solution is converged (as expected!)
# i.e. check that the sum of squared residuals at the bottom of the file
# is a very small number!!
fem evaluate residuals;coupled_residuals wrt geom_params coupled;
fem list map;bob2 solution coupled;
# Save deformed configurations for end diastole
$file = $cylinderSaveName."_enddiast";
fem define initial;w;$file region $WALL,$CAVITY;
#
# ISOVOLUMIC CONTRACTION
#
if (! defined $dCai) {
$dCai = 0.05;
}
if (! defined $maxCai1) {
$maxCai1 = 0.2;
}
if (! defined $maxCai2) {
$maxCai2 = 0.9;
}
$maxCai = $maxCai1;
$Cai = $dCai;
print "\n";
print "$Bgreen$Fblue $FBnone\n";
print "$Bgreen$Fblue Solving isovolumic contraction $FBnone\n";
print "$Bgreen$Fblue $FBnone\n";
print "\n";
$file = $cylinderSaveName."_endisovolcontr";
if ($SOLVE_ISOVOL) {
read com;$example/solve_contraction;
# Save deformed configurations for end diastole
fem define initial;w;$file region $WALL,$CAVITY;
} else {
fem define initial;r;"$example/$file" region $WALL,$CAVITY;
&DefineCalciumLevel($maxCai);
fem solve coupled increment 0.0 iter $MAX_ITERS error $TOL;
fem export nodes;bob field as cylinder class $TISSUE region $WALL;
fem export nodes;fred field as cavity class $TISSUE region all;
}
print "\n";
print "$Bgreen$Fblue $FBnone\n";
print "$Bgreen$Fblue Finished isovolumic contraction $FBnone\n";
print "$Bgreen$Fblue $FBnone\n";
print "\n";
# Check for iso-volumic contraction
#fem list element deformed total basis 10 variable cavityVolume region $CAVITY;
#print "$Bcyan$Fred Cavity volume = $cavityVolume $FBnone\n";
#
# EJECTION
#
# Decrease cavity 'stiffness parameter' and solve for ejection
&DefineCavityStiffness(10,$CAVITY);
# Continue from previous calcium activation level
$Cai = $maxCai + $dCai;
$maxCai = $maxCai2;
print "\n";
print "$Bgreen$Fblue $FBnone\n";
print "$Bgreen$Fblue Solving ejection $FBnone\n";
print "$Bgreen$Fblue $FBnone\n";
print "\n";
if ($SOLVE_EJECTION) {
read com;$example/solve_contraction;
}
print "\n";
print "$Bgreen$Fblue $FBnone\n";
print "$Bgreen$Fblue Finished ejection $FBnone\n";
print "$Bgreen$Fblue $FBnone\n";
print "\n";
# Save deformed configurations for end diastole
$file = $cylinderSaveName."_endeject";
fem define initial;w;$file region $WALL,$CAVITY;
# Determine cavity volume at end-ejection.
# NOTE: to calc cavity volume, the deformed position of the central
# cavity node must be recalc'ed from the averaged adjacent wall nodes.
# After doing this the current solution is meaningless, ie. best to save
# deformed configurations for both regions before this (as above).
# ?????????????????????????????????????????????????????????????
#fem update solution cavity_reference average 6 in 1 node 2 region $CAVITY
#fem list element deformed total region $CAVITY
# ?? Check ejection fraction
#fem list element deformed total basis 10 variable ejectedVolume region $CAVITY;
#print "$Bcyan$Fred Ejected volume = $ejectedVolume (from $cavityVolume) $FBnone\n";
