Example d42: Single rectus femoris knee flexion

This example illustrates flexion of a knee with a single rectus femoris muscle

Flexion solution Gauss point activation
Fig. 1 Flexion solution Fig. 2 Gauss point activation from nerve innervation


The comfile run by this example is as follows:

#################################################################################
#
# This example shows a single rectus muscle being loaded via flexion from motion
# capture. At angles of 10,30 and 60 degress an active tension is also added while 
# the muscle is lengthened (eccentric contraction). This is added via an ipacti
# file giving calcium values corresponding to a response curve triggered from a 
# single nerve activation. The Gauss points are then given calcium values according 
# to the response curve. 
#
# This example illustrates the modelling framework for quasi-static musculoskeletal 
# modelling using motion capture as the driving displacement boundary condition. The 
# perl script used to generate the kinematics (translations+rotations) is also included (jcs.pl).
# The parameters chosen are for stability and robustness and only a few load steps are
# simulated due to computation time. This example also uses the cellML framework for
# including material parameters and the 'pole-zero' law (soft tissue) and St-Venant-
# Kirchoff model (bone) are included as MathML (CellML files).
#
# At the chosen angles, the effect on passive stresses from active tension are analysed 
# at chosen Gauss points. 
#
# Also included are the perl scripts for generating the ipacti files according to the
# calcium activation curve, and for extracting the passive Gauss point stresses.
#
#################################################################################

# If the example path is not set, default to current directory
if(!defined $example){$example = "./";}

################################################ Motion capture kinematics

# Read in translations and rotations for current load step
$input_file1 = "$example/kinematics.txt";

# Open input file1
CORE::open (file1,$input_file1) or die "Error opening ";

################################################ Parameter setup

# Define parameters,regions,and coordinate system
fem def para;r;contact;example
fem def coord 3,1
fem def region;r;contact;example
fem def bases;r;contact;example 

$REG=1;

###################################################### Read in meshes

# list
@list = ('single_rectus','femur','tibia')
foreach $component (@list)
{ 
#Read in 
fem def nodes;r;$component;example reg $REG
fem def elements;r;$component;example reg $REG

if ($REG==1){fem change nodes translate by 0,0,7 reg $REG}

# Define fibre field (nj 4,5,6)
fem def fibre;r;$component;example reg $REG
fem def elements;r;$component;example fibre reg $REG

# Define field (nj 7,8,9)
if ($REG==1){fem def field;r;$component;example reg $REG}
if ($REG==1){fem def elements;r;$component;example field reg $REG}

# Node groups
if ($REG==1){fem group nodes 23113..23116 as TENDPRESCRIBE reg $REG}
if ($REG==1){fem group nodes 14001..14004 as MUSTOPFIX reg $REG}
if ($REG==1){fem group nodes 1..64,14001..14008,14013..14016,14021..14024,14029..14032,23109..23120  as MESHALL1 reg $REG}

# Element groups
if ($REG==1){fem group elements 24004,24008,14004 as TEN reg $REG}
if ($REG==1){fem group elements 24009..24010 as PATLIG reg $REG}
if ($REG==1){fem group elements 10..18 as CART reg $REG}
if ($REG==1){fem group elements 1..9,19..27 as BONE reg $REG}
if ($REG==1){fem group elements 14001..14003 as MUS reg $REG}

## Export initial mesh
#fem exp nodes;init_$component as $component reg $REG
#fem exp elements;init_$component as $component reg $REG
## Export gauss points for identifying activation points
#if ($REG==1){fem exp gauss;gauss_pts as gauss_pts reg $REG}

# Define equations (finite elasticity)
if ($REG==1){fem def equa;r;$component;example class 1 reg $REG lockregion}  
# Also need to define equations for other two regions
# to ensure they are defined as contact bodies.
if ($REG==2){fem def equa;r;$component;example class 1 reg $REG lockregion}  
if ($REG==3){fem def equa;r;$component;example class 1 reg $REG lockregion}  

# Define at grid points and read from cells
if ($REG==1){fem def mate;r;$component;example reg $REG} 

# Read initial conditions for single_rectus
if ($REG==1){fem def init;r;$component;example reg $REG}

# Define solve for mechanics problem
if ($REG==1){fem def solve;r;$component;example reg $REG}

if ($REG==1){fem define acti;r;$component;example reg $REG}

$REG=$REG+1;
}

# define contact parameters
fem def contact;r;contact;example

####################################################### Grid-based cell setup

# Grid based problem setup
# Ensure that grid points in each xi direction match 
# gauss points in xi directions from ipbase file
fem def grid;r;rectus_cell;example gauss class 2  # place grid points at Gauss points
fem update grid geometry  # defines geometric connectivity of grid points. Class 2 problem
#fem exp node;check as check reg 1 # export grid for checking
#fem exp elem;check as check grid_numbers reg 1 # export grid for checking

# Only have to do once (grid groups)
fem group grid element TEN as TENGRP reg 1 # group grid points within elements
fem group grid element PATLIG as PATLIGGRP reg 1 # group grid points within elements
fem group grid element CART as CARTGRP reg 1 # group grid points within elements
fem group grid element BONE as BONEGRP reg 1 # group grid points within elements
fem group grid element MUS as MUSGRP reg 1 # group grid points within elements
fem group grid element TEN,PATLIG,MUS as ALL_SOFT reg 1 # group grid points within elements
fem group grid element CART,BONE as ALL_HARD reg 1 # group grid points within elements
#fem group grid element CART,BONE,TEN,PATLIG,MUS as ALL_GRID reg 1 # group grid points within elements
#fem group grid;w;grid_groups reg 1
#fem group grid;r;grid_groups_rectus;example reg 1

# Define equations (cell)
fem def equa;r;rectus_cell;example class 2  # Define cell equation. Not actually used. Just defined to setup structures

# Needed to setup
fem def mate;d class 2  # set up conductivities- default

# Define cell variants and material properties
fem def cell;r;rectus_cell;example class 2   # Define cell variants and pole zero law parameters. Class 2 problem

# Define cell material
fem def mate;r;rectus_cell;example cell class 2  # Define material properties for cell

# Default init conditions for cell class 2
fem def init;d class 2 

# Setup solve for cellml class 2
fem def solv;r;rectus_cell;example class 2 

#runs the cellml time integration solver stuff once to initialise cellml solver things. 
#The "0" class 2 to 0 says you have a 0 time step for the integration
#fem solve class 2 to 0 


####################################################### Prescribe boundary kinematics

# Set centre of rotation
$cenx=0.0;
$ceny=0.0;
$cenz=0.0;

# Load steps
#for ($i=1;$i<128;$i++)
for ($i=1;$i<11;$i++) # for testing only
{
# Read in next set of angles and translations
$li = <file1>;
@compnew = split(/\s+/,$li);  

# Make load step for tibia and prescribed (attached) tendon nodes

# Rot about x-axis
fem change nodes rotate by $compnew[0] about $cenx,$ceny,$cenz axis 1,0,0 reg 3
fem change nodes rotate by $compnew[0] node TENDPRESCRIBE about $cenx,$ceny,$cenz axis 1,0,0 reg 1

# Rot about y-axis
fem change nodes rotate by $compnew[1] about $cenx,$ceny,$cenz axis 0,1,0 reg 3
fem change nodes rotate by $compnew[1] node TENDPRESCRIBE about $cenx,$ceny,$cenz axis 0,1,0 reg 1

# Rot about z-axis
fem change nodes rotate by $compnew[2] about $cenx,$ceny,$cenz axis 0,0,1 reg 3
fem change nodes rotate by $compnew[2] node TENDPRESCRIBE about $cenx,$ceny,$cenz axis 0,0,1 reg 1

# Translation
fem change nodes translate by $compnew[3],$compnew[4],$compnew[5] reg 3
fem change nodes translate by $compnew[3],$compnew[4],$compnew[5] node TENDPRESCRIBE reg 1

# Offset tibia motion so femur is stationary (reverse order)
fem change nodes translate by -$compnew[9],-$compnew[10],-$compnew[11] reg 3
fem change nodes rotate by -$compnew[8] about $cenx,$ceny,$cenz axis 0,0,1 reg 3
fem change nodes rotate by -$compnew[7] about $cenx,$ceny,$cenz axis 0,1,0 reg 3
fem change nodes rotate by -$compnew[6] about $cenx,$ceny,$cenz axis 1,0,0 reg 3
fem change nodes translate by -$compnew[9],-$compnew[10],-$compnew[11] node TENDPRESCRIBE reg 1
fem change nodes rotate by -$compnew[8] node TENDPRESCRIBE about $cenx,$ceny,$cenz axis 0,0,1 reg 1
fem change nodes rotate by -$compnew[7] node TENDPRESCRIBE about $cenx,$ceny,$cenz axis 0,1,0 reg 1
fem change nodes rotate by -$compnew[6] node TENDPRESCRIBE about $cenx,$ceny,$cenz axis 1,0,0 reg 1

# Update boundary conditions to solution file
fem update solution from geometry substitute reg 1

# activation files
# 10 degress activation
if(($i>5)&&($i<32))
{
$counter=$i-5;
fem define acti;r;$counter;example reg 1
}

# 30 degress activation
if(($i>54)&&($i<81))
{
$counter=$i-54;
fem define acti;r;$counter;example reg 1
}

# 60 degress activation
if(($i>101)&&($i<128))
{
$counter=$i-101;
fem define acti;r;$counter;example reg 1
}

$CONVERGED=0;
while ($CONVERGED==0) 
##for ($j=1;$j<1;$j++)
{

####################################### Contact routines

# Master and slave face groups
#fem group faces 48,53,58,64,69,74,78,81,84,134 as PATELLAFACE reg 1
fem group faces 64,69,74,78,81,84 as PATELLAFACE reg 1
#fem group faces 233,237,241,245,249,253,257,297,301,305,309,313,317,321,361,365,369,373,377,381,385,433,437,441,497,501,505,561,565,569,625,629,633 as FEMURFACE reg 2
fem group faces 233,237,241,245,249,253,257,297,301,305,309,313,317,321,361,365,369,373,377,381,433,437,441,497,501,505,561,565 as FEMURFACE reg 2

# Set contact Xi points on specified faces as 4x4 points
fem def xi;c contact_points faces PATELLAFACE points 4 reg 1

# Define data at xi positions 
fem def data;c from_xi 

# Update slave info
fem update data field to slave

# Project onto target face
fem def xi;c closest_face faces FEMURFACE reg 2

# Store projection gap in data fields
fem update data field from gap 

# Update master info
fem update data field to master 

# Place initial geometry YP(ny,3) into XP
fem up geom from sol YP_index 3 to 1..3 reg 1

##################################### Mechanics problem

# Solve finite elasticity/contact problem
fem solve error 0.1 iterate 15 reg 1

##Place nodal reaction forces YP(ny,4) into XP nj 7..9
##fem up geom from sol YP_index 4 to 7..9 reg 1
##write out nodal reaction field
##fem def field;w;forces_$i reg 1

## Place contact forces YP(ny,6) into XP nj 7..9
#fem up geom from sol YP_index 6 to 7..9 reg 1

# update XP from YP
fem up geom from sol to 1..3 reg 1

} # Newton iterations

################################################ Stresses (Gauss point)

# Place initial geometry YP(ny,3) into XP
fem up geom from sol YP_index 3 to 1..3 reg 1

# Export stress solutions at Gauss points

# Fibre passive stresses at Gauss points  
fem up gauss stress fibre passive principal reg 1
#fem def gauss;w;passive_stress_rectus_patella_$i as single_rectus number 3 reg 1

# Update XP from YP
fem up geom from sol to 1..3 reg 1

################################################ Export solution

# Write out new single_rectus geometry
#fem def nodes;w;def_single_rectus$i as single_rectus reg 1
#fem exp nodes;def_single_rectus$i as single_rectus reg 1

# Export femur rotated mesh
#fem exp nodes;def_femur$i as femur reg 2

# Export tibia rotated mesh
#fem exp nodes;def_tibia$i as tibia reg 3

} # Load steps

# Close the files.
CORE::close file1;

Additional testing commands:



Files used by this example are:

Name                                Modified     Size

example_d42.com 23-Apr-2008 11k 0.1passive_10d_g1.txt 18-Feb-2005 390 0.1passive_10d_g2.txt 18-Feb-2005 390 0.1passive_10d_g3.txt 18-Feb-2005 390 0.1passive_30d_g1.txt 18-Feb-2005 390 0.1passive_30d_g2.txt 18-Feb-2005 390 0.1passive_30d_g3.txt 18-Feb-2005 390 0.1passive_60d_g1.txt 18-Feb-2005 390 0.1passive_60d_g2.txt 18-Feb-2005 390 0.1passive_60d_g3.txt 18-Feb-2005 390 1.ipacti 18-Feb-2005 2.0k 10.ipacti 18-Feb-2005 2.1k 11.ipacti 18-Feb-2005 2.1k 12.ipacti 18-Feb-2005 2.1k 13.ipacti 18-Feb-2005 2.1k 14.ipacti 18-Feb-2005 2.1k 15.ipacti 18-Feb-2005 2.1k 16.ipacti 18-Feb-2005 2.1k 17.ipacti 18-Feb-2005 2.1k 18.ipacti 18-Feb-2005 2.1k 19.ipacti 18-Feb-2005 2.1k 2.ipacti 18-Feb-2005 2.0k 20.ipacti 18-Feb-2005 2.1k 21.ipacti 18-Feb-2005 2.1k 22.ipacti 18-Feb-2005 2.1k 23.ipacti 18-Feb-2005 2.1k 24.ipacti 18-Feb-2005 2.1k 25.ipacti 18-Feb-2005 2.1k 26.ipacti 18-Feb-2005 2.1k 3.ipacti 18-Feb-2005 2.0k 4.ipacti 18-Feb-2005 2.0k 5.ipacti 18-Feb-2005 2.1k 6.ipacti 18-Feb-2005 2.1k 7.ipacti 18-Feb-2005 2.1k 8.ipacti 18-Feb-2005 2.1k 9.ipacti 18-Feb-2005 2.1k abd.dat 18-Feb-2005 2.3k angles_dynamic.txt 18-Feb-2005 12k backup.ipcell 18-Feb-2005 9.3k backup.ipmatc 18-Feb-2005 13k both.ipcell 18-Feb-2005 9.3k both.ipmatc 18-Feb-2005 10k cmiss_angles.txt 18-Feb-2005 5.3k cmiss_rot_tra.bak 18-Feb-2005 16k cmiss_rot_tra.txt 18-Feb-2005 16k contact.ipbase 18-Feb-2005 10k contact.ipcont 25-Sep-2008 628 contact.ippara 21-Feb-2008 6.1k contact.ipregi 18-Feb-2005 93 create_ipacti_muscle_paper.pl 18-Feb-2005 8.8k extract_gauss_stress_passive.pl 18-Feb-2005 6.2k femur.ipelem 18-Feb-2005 453k femur.ipelfb 18-Feb-2005 230k femur.ipequa 18-Feb-2005 2.2k femur.ipfibr 18-Feb-2005 616k femur.ipnode 18-Feb-2005 895k flex.dat 18-Feb-2005 2.1k flexion9.trc.txt 18-Feb-2005 358k gauss_pts.exdata 18-Feb-2005 155k grid_groups.ipgrgp 18-Feb-2005 53k grid_groups_rectus.ipgrgp 18-Feb-2005 272k hertz.ippara 21-Feb-2008 5.9k jcs.pl 18-Feb-2005 63k kinematics.txt 18-Feb-2005 9.8k pole.ipcell 18-Feb-2005 7.0k pole.ipmatc 18-Feb-2005 7.6k polezero_nearly_incompressible.cml 18-Feb-2005 52k rectus_cell.ipcell 18-Feb-2005 9.4k rectus_cell.ipcell.bak 18-Feb-2005 9.3k rectus_cell.ipcellbak 18-Feb-2005 9.4k rectus_cell.ipequa 18-Feb-2005 1.6k rectus_cell.ipgrid 18-Feb-2005 622 rectus_cell.ipmatc 18-Feb-2005 13k rectus_cell.ipmatc.bak 18-Feb-2005 13k rectus_cell.ipmatcbak 18-Feb-2005 18k rectus_cell.ipsolv 13-Apr-2007 2.5k rot.dat 18-Feb-2005 2.1k saint.ipcell 18-Feb-2005 2.5k saint.ipmatc 18-Feb-2005 2.7k single_rectus.ipacti 18-Feb-2005 921 single_rectus.ipelem 18-Feb-2005 18k single_rectus.ipelfb 18-Feb-2005 9.3k single_rectus.ipelfd 18-Feb-2005 9.2k single_rectus.ipequa 18-Feb-2005 2.2k single_rectus.ipfibr 18-Feb-2005 147k single_rectus.ipfiel 18-Feb-2005 150k single_rectus.ipinit 18-Feb-2005 5.1k single_rectus.ipmate 18-Feb-2005 1.7k single_rectus.ipnode 18-Feb-2005 173k single_rectus.ipsolv 16-Aug-2010 2.7k single_rectus.ipsolv.old 13-Apr-2007 2.5k st_venant_kirchoff.cml 18-Feb-2005 10k static.trc.txt 18-Feb-2005 171k test_output.com 18-Feb-2005 0 testboth.ipcell 18-Feb-2005 9.3k testboth.ipmatc 18-Feb-2005 9.8k tibia.ipelem 18-Feb-2005 96k tibia.ipelfb 18-Feb-2005 47k tibia.ipequa 18-Feb-2005 2.2k tibia.ipfibr 18-Feb-2005 132k tibia.ipnode 18-Feb-2005 197k view.com 18-Feb-2005 5.8k

Download the entire example:

Name                      Modified     Size

examples_d_d4_d42.tar.gz 17-Aug-2010 816k

Testing status by version:

StatusTestedReal time (s)
rs6000-aix
cm64SuccessSun Mar 13 05:27:19 2005686

Testing status by file:


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

Input last modified: Mon Aug 16 11:35:06 2010


CMISS Help / Examples / d / d4 / d42