Example 5535: Passive inflation/isovolumic contraction of a trilinear prolate spheroid with varying fibre angles

A single prolate. Geometry is interpolated isochorically (Focus^3.Cosh[lambda].Sinh^2[lambda]). Hydrostatic pressure is interpolated quadratically in Xi3 (element based) and the pole-zero law defines the incompressible material behaviour. The prolate is inflated then contracted isovolumically with differing fibre angles. The relationship between fibre angles and stress is plotted.


The comfile run by this example is as follows:

# Example 5535

# A single axisymetric prolate LV model is used to investigate
# the local feedback mechanisms affecting the fibre 
# orientation.
# (Tissue Remodelling with Micro Structurally Based
#  Material Laws, Peter Hunter & Theo Arts 1997)
#
# Muscle Fibre Orientation
# PROBLEM_1: Sheet angles are all zero, fibre angles are zero at
#            the apex nodes, -90 at the endocardial base, and the 
#            epicardial base node fibre angle is varied from 20 to 90 degs.
# PROBLEM_2: Endocardial sheet angles are -45 epicardial are 45, 
#            Endocardial fibre angles are 65 epicardial
#            are varied from -20 to -90 degs.

$TRUE = 1
$FALSE = 0  

# Choose which problem to run
$PROBLEM_1 = $TRUE
$PROBLEM_2 = $FALSE

$MIN_ANGLE = 20
$MAX_ANGLE = 90
$INC_ANGLE = 10


system "rm plot.dat"
set echo

for ($ANGLE=$MIN_ANGLE;$ANGLE<=$MAX_ANGLE;$ANGLE+=$INC_ANGLE)
{

	print "\n\n\n\$ANGLE = $ANGLE\n";

  fem
  fem define parameters;r;profib;example
  fem define coordinate;r;profib;example #Defines 3 GLOBAL PROLATE SPHEROIDAL
  #                                      !  coordinates.
  #                                      !  Transmural interpolation is in
  #                                      !  Focus^3.Cosh(Lamda).Sinh^2(Lamda).
  fem define node;r;;example              #Focus = 37.5. Reads 4 NODES, 3 COORDS,
  #                                      !  1 VERSION per node per nj, 0
  #                                      !  DERIVATIVES for all coords. Nodal
  #                                      !  coordinates, Xj, are (lamda,mu,theta):
  #                                      !  1=(0.43,0,0);2=(0.43,120,0);
  #                                      !  3=(0.72,0,0); 4=(0.72,120,0).
  fem define base;r;;example              #Defines 3 basis functions:
  #                                      !  1) Lagrange/Hermite with 3 Xi coords,
  #                                      !     Linear Lagrange interpolant and 1
  #                                      !     Gauss point in Xi1 dirctn (theta-
  #                                      !     circumferential) but 3 in Xi2 (mu-
  #                                      !     longitudinal) and Xi3 (lamda-
  #                                      !     transmural);
  #                                      !  2) An auxiliary basis with 5 auxiliary
  #                                      !     element parameters, pressure basis
  #                                      !     same number of Gauss points as for
  #                                      !     basis function 1, and polynomial
  #                                      !     degrees of: 0 for Xi1 and Xi2 for
  #                                      !     all 5 parameters; 0 in Xi3 for
  #                                      !     parameter 1; 1 for parameter 2; 2
  #                                      !     for parameter 3; -1 for parameter
  #                                      !     4 (Xi3=0 face pressure bc); & -2
  #                                      !     for parameter 5 (Xi3=1 face
  #                                      !     pressure bc).
  #                                      !  3) 2 Xi coords with linear Lagrange
  #                                      !     interpolant, and 3 Gauss points
  #                                      !     in Xi1 and Xi2 directions.
  fem define element;r;;example           #Defines 1 prolate spheroidal element.
  set echo
  if($PROBLEM_1)
  {
    fem define fibre;r;problem_1/profib$ANGLE;example
  } elsif($PROBLEM_2)
  {
    fem define fibre;r;problem_2/profib$ANGLE;example
  }
  set echo
  #                                      #Defines fibre angle orientation wrt Xi1
  #                                      !  direction (theta) in degrees. 
  #                                      !  degrees at outer nodes 3,4.
  fem define element;r;profib;example fibre    #Defines fibre elements.
  fem define window                      #Defines window 
  fem draw lines                         #Makes line segments visible on window.
  fem define equation;r;;example          #Static 3D finite elasticity problem with
  #                                      !  geometric coords as the dependent
  #                                      !  variables. Solution is by Galerkin
  #                                      !  finite elements in a nonlinear
  #                                      !  analysis of a nonlinear equation.
  #                                      !  Incompressible material with basis 1
  #                                      !  in all elements for dep vars 1,2,3
  #                                      !  and basis 2 in all elements for dep
  #                                      !  var 4 (hydrostatic pressure).
  fem define material;r;;example          #Stresses in constitutive law are
  #                                      !  referred to body coords.  Defines a
  #                                      !  hyperelastic constitutive law for
  #                                      !  which the strain energy function, W,
  #                                      !  is a function of fibre/transverse
  #                                      !  strains that is pole-zero in strains
  #                                      !  with 16 constant constitutive law
  #                                      !  parameters. The Lagrangian strain
  #                                      !  components are referred to the
  #                                      !  undeformed fibre coordinate system
  #                                      !  where the 1-direction is aligned with
  #                                      !  fibre axis, 2-direction is the sheet
  #                                      !  axis (in the plane of the sheet,
  #                                      !  perpendicular to the fibres), and the
  #                                      !  3-direction is the sheet-normal coordinate.

  # 
  #
  # INFLATION
  #
  #

  fem define initial;r;;example           #Boundary pressure increments are entered
  #                                      !  and hydrostatic pressure is matched
  #                                      !  across boundaries with adjacent Xi3
  #                                      !  faces. Initial displacements are all
  #                                      !  zero. For dependent variable 2
  #                                      !  (mu-dirn) fix all nodes (1..4);
  #                                      !  Dependent variable 3 (theta-dirn),
  #                                      !  fix node 2. No force bcs.
  #                                      !  Dependent variable 4 (hyd press) fix
  #                                      !  element 1, parameter 4 (Xi3=0 face
  #                                      !  inside prolate) with increment
  #                                      !  magnitude 1, and parameter 5 (Xi3=1
  #                                      !  face outside prolate) with increment
  #                                      !  magnitude 0. No force bcs.
  fem define solve;r;profib;example      #  Read in solution information.
  fem solve increment 0.2 iter 8         #Solve the problem.
  #                                      !Incrementing of the ventricular
  #                                      !  pressure is neccessary to avoid
  #                                      !  unstable solutions. This depends on
  #                                      !  the constitutive law and pressure
  #                                      !  the constitutive law and pressure
  #                                      !  magnitude, as the wall is gradually
  #                                      !  stiffened.
  fem solve increment 0.3 iter 8                         
  fem solve increment 0.5        
  #                              
  #                              
  #                              
  fem draw lines def dotted      #Draw deformed mesh on window.
  fem def initial;w;profib_def
  fem eval reac
  # ############################################ 
  # Inflated solution is converged OK! 
  # ############################################ 

  # 
  #
  # COUPLING SET UP 
  #
  #
  $WALL = 1
  $CAVITY = 2
  # 
  # Set up cavity region geometry 
  fem define region;r;coupled;example 
  fem define coord;r;cavity;example region $CAVITY
  #
  fem define;add base;r;cavity;example 
  fem define node;r;cavity;example region $CAVITY
  fem define element;r;cavity;example region $CAVITY
  # 
  # Set up equations, material properties and initial conditions 
  # for cavity region 
  fem define equation;r;coupled;example region $WALL,$CAVITY lock
  fem define material;r;cavity_isovol;example region $CAVITY 
  fem define initial;r;cavity;example region $CAVITY
  # 
  # Define cavity-to-wall coupling and solution procedure information 
  fem define coupling;r;coupled;example 
  fem define solve;r;coupled;example 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
  # 
  # 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 (node 6) is set to the average of the 
  # x-coords for the adjacent interface nodes (node 2 - this could 
  # be a list of nodes) 
  fem update solution cavity_reference average 6 in 1 node 2 region $CAVITY
  # 
  # List undeformed (pre-inflation) and inflated cavity volumes 
  fem list element total region $CAVITY
  fem list element deformed total region $CAVITY
  # ############################################ 
  # cavity volume inflated from 40ml to 60ml. 
  # ############################################ 
  # 
  fem draw lines region $WALL,$CAVITY
  fem draw lines dot def region $WALL,$CAVITY
  # 
  # Check that the current coupled solution is converged (as expected!) 
  fem evaluate residuals wrt geom_params coupled
  # ############################################ 
  # inflated state is a converged solution for the coupled problem! 
  # ############################################ 
 
  #
  # 
  # ISOVOLUMIC CONTRACTION 
  # 
  # 

  # Increment intracell calcium conc and solve for isovolumic contraction 
  fem define material;r;profib_active;example 
  if($PROBLEM_1)
  {
    fem define active;r;problem_1/profib_contract_Ca20;example
  } elsif($PROBLEM_2)
  {
    fem define active;r;problem_2/profib_contract_Ca10;example
  }

  fem solve coupled increment 0.0
  # ############################################ 
  # Converged solution for the coupled wall/cavity problem.
  # ############################################ 
  #fem define initial;r;profib_endisovolcontr region WALL,CAVITY 
  # 
  # Save deformed configurations for end ejection 
  fem define initial;w;profib_endisovolcontr region $WALL,$CAVITY 
  # 
  fem draw lines deformed dotted rgb=blue region $WALL,$CAVITY
  # 
  # Check that cavity volume has not changed 
  fem list element deformed total region $CAVITY
  # ############################################ 
  # This contraction phase was indeed isovolumic! 
  # ############################################ 

  # 
  #
  # POSTPROCESSING, VISUALISATION & PLOTTING
  #  
  #
  fem export nodes;heart$ANGLE as heart
  fem export elements;heart$ANGLE as heart
  fem export gauss;heart$ANGLE as gauss
  fem evaluate field;fibre_angles$ANGLE geometric gauss # get fibre & sheet angles
  fem list stress;stresses$ANGLE fibre # get stresses
  system "echo \"\# Fibre angle ${ANGLE} degrees\" >> plot.dat"
  system "awk -v angle=${ANGLE} -f ${example}extract_fibre_angles_stress.awk ${example}elements.txt fibre_angles${ANGLE}.opfiel stresses${ANGLE}.opstre >> plot.dat"
}
system "gnuplot ${example}angle_stress.gnu"

Additional testing commands:

fem list initial region $WALL,$CAVITY

Files used by this example are:

Name                             Modified     Size

example_5535.com 07-Jul-2003 11k angle_stress.gnu 10-Apr-2000 1.3k cavity.ipbase 10-Apr-2000 1.1k cavity.ipcoor 10-Apr-2000 707 cavity.ipelem 10-Apr-2000 404 cavity.ipinit 12-Dec-2002 1.4k cavity.ipnode 10-Apr-2000 1.3k cavity_eject.ipmate 10-Apr-2000 386 cavity_isovol.ipmate 10-Apr-2000 390 coupled.ipcoup 21-Aug-2002 437 coupled.ipcoup.old 10-Apr-2000 386 coupled.ipregi 10-Apr-2000 94 coupled.irequa 02-May-2004 3.8k coupled.irsolv 16-Aug-2010 3.1k coupled.irsolv.old 13-Apr-2007 2.9k elements.txt 10-Apr-2000 2 extract_fibre_angles_stress.awk 10-Apr-2000 2.4k extract_fibre_wall_strains.awk 10-Apr-2000 2.7k problem_1/ 03-Mar-2004 - problem_2/ 03-Mar-2004 - profib.ipacti 03-Mar-2004 735 profib.ipbase 10-Apr-2000 4.2k profib.ipcoor 10-Apr-2000 707 profib.ipelem 10-Apr-2000 409 profib.ipelfb 10-Apr-2000 349 profib.ipelfd 10-Apr-2000 276 profib.ipequa 02-May-2004 2.0k profib.ipfibr 30-Jan-2001 1.3k profib.ipfiel 10-Apr-2000 601 profib.ipfit 13-Apr-2007 1.5k profib.ipinit 12-Dec-2002 1.8k profib.ipmate 12-Dec-2002 5.7k profib.ipnode 10-Apr-2000 1.2k profib.ippara 12-Nov-2002 5.9k profib.ipsolv 16-Aug-2010 2.3k profib.ipsolv.old 13-Apr-2007 2.1k profib.ipwind 10-Apr-2000 262 profib_active.ipmate 12-Dec-2002 6.2k profib_contract.ipacti 03-Mar-2004 717 profib_def.ipinit 12-Dec-2002 6.2k profib_endeject.irinit 12-Dec-2002 11k profib_endisovolcontr.irinit 12-Dec-2002 11k prol.ippara 12-Nov-2002 5.9k test_output.com 20-Nov-2001 38 view.com 20-Nov-2001 3.2k

Download the entire example:

Name                           Modified     Size

examples_5_55_553_5535.tar.gz 17-Aug-2010 24k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmFailureSun Mar 6 00:31:22 201615
last breakTue Aug 12 00:37:00 201428
last successSun Apr 17 01:02:00 201131
cm-debugFailureSun Mar 6 00:46:00 201630
last breakTue Aug 12 00:40:00 201447
last successThu Apr 7 01:10:00 201154
mips-irix
cmSuccessSun Aug 19 03:01:05 200799
cm-debugSuccessWed Aug 15 03:11:35 2007260
cm-debug-clear-mallocSuccessSat Aug 18 03:44:28 2007298
cm-debug-clear-malloc7SuccessMon Aug 20 03:45:44 2007301
cm64SuccessSun Aug 19 03:03:38 2007116
cm64-debugSuccessTue Aug 21 04:13:40 2007267
cm64-debug-clear-mallocSuccessFri May 20 10:36:11 2005146
rs6000-aix
cmSuccessWed Mar 4 01:15:21 200913
cm-debugSuccessMon Mar 2 01:28:40 200981
cm64SuccessWed Mar 4 01:16:01 200913
cm64-debugSuccessTue Mar 3 01:32:30 200975
x86_64-linux
cmFailureSun Mar 6 00:02:17 20167
last breakTue Aug 12 00:19:00 201410
last successSun Apr 17 00:30:00 201111
cm-debugFailureSun Mar 6 00:03:52 201614
last breakTue Aug 12 00:40:00 201418
last successFri Apr 8 00:13:00 201119

Testing status by file:


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

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


CMISS Help / Examples / 5 / 55 / 553 / 5535