Example 5f2: Isovolumic contraction of cylinder mesh in RC coordinates

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.

Initial cylinder geometry (blue) and inflated geometry (brown). Addition of the undeformed cavity volume mesh (green lines, red surface).
Update the cavity mesh to match the inflated cyclinder deformation (gold lines). Active contraction of the cylinder tissue while keeping the cavity volume constant.
Simulating ejection by allowing one node from the cavity mesh to move.


The comfile run by this example is as follows:


# 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";


Files used by this example are:

Name                                        Modified     Size

5f2.com 27-Apr-2006 10k 2CubicHermite-Linear_4x4x4Gauss.ipbase 20-Apr-2006 1.9k 2CubicHermite_4x4Gauss.ipbase 20-Apr-2006 1.6k 2Linear-CubicHermite_4x4x4Gauss.ipbase 20-Apr-2006 1.8k 2Linear_4x4Gauss.ipbase 20-Apr-2006 1.2k 3CubicHermite_4x4x4Gauss.ipbase 20-Apr-2006 1.9k 3Linear-2PressAuxXi3_4x4x4Gauss.ipbase 20-Apr-2006 1.9k 3Linear_4x4x4Gauss.ipbase 20-Apr-2006 1.4k 3drc.ipcoor 20-Apr-2006 570 CubicHermite-Linear_4x4x4Gauss.ipbase 20-Apr-2006 1.5k Linear-2CubicHermite_4x4x4Gauss.ipbase 20-Apr-2006 1.9k Linear-CubicHermite_4x4x4Gauss.ipbase 20-Apr-2006 1.5k PressAuxXi3_4x4x4Gauss.ipbase 20-Apr-2006 1.1k cavity.ipcoor 20-Apr-2006 565 cavity.ipinit 20-Apr-2006 4.2k cavity_large_refined.ipelem 20-Apr-2006 5.2k cavity_large_refined.ipnode 20-Apr-2006 32k coupled.ipcoup 20-Apr-2006 477 coupled.ipregi 20-Apr-2006 111 coupled.irequa 20-Apr-2006 3.9k coupled_gmres_nosparse.irsolv 16-Aug-2010 3.0k coupled_gmres_nosparse.irsolv.old 13-Apr-2007 2.9k coupled_lu_nosparse.irsolv 16-Aug-2010 2.6k coupled_lu_nosparse.irsolv.old 13-Apr-2007 2.4k cylinder.ipinit 20-Apr-2006 2.8k cylinder.ippara 20-Apr-2006 5.9k cylinder_large_refined.ipelem 20-Apr-2006 12k cylinder_large_refined.ipelfb 20-Apr-2006 7.5k cylinder_large_refined.ipfibr 20-Apr-2006 21k cylinder_large_refined.ipnode 20-Apr-2006 86k cylinder_large_refined_endinflation.ipinit 20-Apr-2006 214k define_active.com 20-Apr-2006 1.8k define_cavity_stiffness.com 20-Apr-2006 1.6k finelas_tch_incomp.ipequa 20-Apr-2006 2.0k formats.com 20-Apr-2006 727 gmres_nosparse.ipsolv 16-Aug-2010 2.2k gmres_nosparse.ipsolv.old 13-Apr-2007 2.1k newton.ipsolv 16-Aug-2010 2.2k newton.ipsolv.old 13-Apr-2007 2.1k orth_incomp.ipmate 20-Apr-2006 5.7k orth_incomp_active.ipmate 20-Apr-2006 6.2k solve_contraction.com 20-Apr-2006 1.7k solve_inflation.com 20-Apr-2006 1.3k view_cavity.com 20-Apr-2006 1.6k view_deformation.com 20-Apr-2006 2.9k view_inflation.com 20-Apr-2006 2.1k view_mesh.com 20-Apr-2006 1.5k

Download the entire example:

Name                      Modified     Size

examples_5_5f_5f2.tar.gz 17-Aug-2010 1.4M

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmFailureSun Mar 6 01:09:24 201646
last breakTue Sep 2 01:05:00 201448
last successSun Apr 17 01:16:00 2011123
cm-debugFailureSun Mar 6 02:01:59 2016282
last breakTue Sep 2 03:30:00 2014281
last successSat Apr 16 03:10:00 2011440
mips-irix
cmSuccessThu Aug 16 05:21:19 2007449
cm-debugSuccessWed Aug 15 04:41:34 20072400
cm-debug-clear-mallocFailureFri Apr 28 14:04:38 200623128
last breakFri Apr 28 07:34:00 2006211
cm-debug-clear-malloc7FailureFri Apr 28 14:05:05 200623289
last breakFri Apr 28 07:36:00 200623289
cm64SuccessThu Aug 16 05:29:25 2007476
cm64-debugSuccessFri Jul 27 06:42:36 20072661
cm64-debug-clear-mallocSuccessThu Apr 27 15:05:02 20062657
rs6000-aix
cmSuccessWed Mar 4 01:26:01 200953
cm-debugSuccessMon Mar 2 02:36:01 2009888
cm64SuccessWed Mar 4 01:27:07 200954
cm64-debugSuccessTue Mar 3 02:39:56 2009824
x86_64-linux
cmFailureSun Mar 6 00:04:55 201628
last breakTue Aug 12 00:03:00 201434
last successSun Apr 17 00:33:00 201137
cm-debugFailureSun Mar 6 00:09:35 2016128
last breakTue Aug 12 00:03:00 2014168
last successSat Apr 16 00:43:00 2011175

Testing status by file:


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

Input last modified: Mon Aug 16 11:19:16 2010


CMISS Help / Examples / 5 / 5f / 5f2