Example d43: Multiple Quadriceps knee flexion

This example illustrates flexion of a knee with the full quadriceps muscles

Initial setup Solution at 75 degrees
Fig. 1 Initial setup Fig. 2 Solution at 75 degrees


The comfile run by this example is as follows:

#################################################################################
#
# This example illustrates the modelling framework for quasi-static musculoskeletal 
# modelling using motion capture as the driving displacement boundary condition. 
# 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).
#
# Illustrated is the complete quadriceps knee simulation and the scripts used to get the
# patella kinematics and the surface contact are given.
#
# This example has been modified to run faster and be more stable for testing. Its main
# purpose is to serve as a template for other musculoskeletal models. In particular, more 
# stable material parameters and placing grid points at Gauss points are two main changes.
# Only 4*4*4 Gauss points are used for each element in this model rather than varying 
# amounts for bone and soft tissue.
#
#################################################################################

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

################################################ Read in information to prescribe tibia motion
# get input arguments
$input_file1 = "$example/kinematics.txt";

# Open input file1
CORE::open (file1,$input_file1) or die "Error opening ";
################################################paramater 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 = ('thesis_final_mesh','femur_cust','tibia_cust')
foreach $component (@list)
{ 
#Read in 
fem def nodes;r;$component;example reg $REG
fem def elements;r;$component;example 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}

if ($REG==1){fem group nodes 23113..23116 as TENDPRESCRIBE reg $REG}
if ($REG==1){fem group nodes 14001..14008,21067..21068,21071..21072,22005..22007,22132,22141..22142,23036..23038,23040,23042..23043 as MUSTOPFIX reg $REG}

#if ($REG==1){fem group nodes 21013,21015,21037..21038,21054..21055,22027,22029,22050 as MUSBACKFIXA reg $REG}
#if ($REG==1){fem group nodes 22052,22065,22068,22136,22138,22140,23013,23016..23017,23064..23066,23077..23079 as MUSBACKFIXB reg $REG}

# Extra fixed nodes on quadriceps (21013,21037,22138,22140,23064,23079)
if ($REG==1){fem group nodes 21013,21037,22138,22140,23064,23079 as MUSBACKFIX reg $REG}

if ($REG==1){fem group nodes 1..64,14001..14008,14013..14016,14021..14024,14029..14032,21009,21011,21013 as MESHALL1 reg $REG}
if ($REG==1){fem group nodes 21015,21033,21035,21037..21038,21050,21052,21054..21055,21067..21068 as MESHALL2 reg $REG}
if ($REG==1){fem group nodes 21071..21072,21074,21078..21080,22005..22007,22025,22027,22029..22030 as MESHALL3 reg $REG}
if ($REG==1){fem group nodes 22049..22050,22052..22053,22065,22068..22070,22089,22093,22095..22096,22132 as MESHALL4 reg $REG}
if ($REG==1){fem group nodes 22135..22142,23013..23014,23016..23017,23019..23020,23022,23026..23027 as MESHALL5 reg $REG}
if ($REG==1){fem group nodes 23036..23038,23040,23042..23043,23060..23061,23064..23067,23077..23079 as MESHALL6 reg $REG}
if ($REG==1){fem group nodes 23081..23083,23096..23120 as MESHALL7 reg $REG}

if ($REG==1){fem group elements 24001..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,21001..21004,22001..22008,23001..23008 as MUS reg $REG}

if ($REG==1){fem group elements 14001..14003 as RECGRP reg $REG}
if ($REG==1){fem group elements 21001..21004 as MEDGRP reg $REG}
if ($REG==1){fem group elements 22001..22008 as INTGRP reg $REG}
if ($REG==1){fem group elements 23001..23008 as LATGRP reg $REG}

# Export initial mesh
#fem exp nodes;init_$component as $component reg $REG
#fem exp elements;init_$component as $component reg $REG

####################################################### Define normal equations

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

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

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

# Define solve for mechanics problem
if ($REG==1){fem def solve;r;contact;example reg $REG}
if ($REG==1){fem define acti;r;multiple_muscle;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;multiple_muscle;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
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;w;grid_groups_multiple;example reg 1
#fem group grid;r;grid_groups_multiple;example reg 1

# Define equations (cell)
fem def equa;r;multiple_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;multiple_cell;example class 2   # Define cell variants and pole zero law parameters. Class 2 problem

# Define cell material
fem def mate;r;multiple_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;multiple_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

$cenx=0.0;
$ceny=0.0;
$cenz=0.0;
#for ($i=1;$i<71;$i++)
for ($i=1;$i<6;$i++)  # for testing
{
# Read in next set of angles
$li = <file1>;
@compnew = split(/\s+/,$li);  
 
#Make load step for tibia and prescribed 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

fem update solution from geometry substitute reg 1

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

#######################################Projection
# patella cartilage faces: 48,53,58,64,69,74,78,81,84
# quadriceps tendon faces: 229,234,239,248,252,256
fem group faces 48,53,58,64,69,74,78,81,84,229,234,239,248,252,256 as PATELLAFACE reg 1
fem group faces 355,359,363,367,371,375,379,383,419,423,427,431,435,439,443,447,483,487,491,495,499,503,551,555,559,563,619,623,627,683,687,691 as FEMURFACE reg 2

# Set contact Xi points on specified faces as 6x6 points
fem def xi;c contact_points faces PATELLAFACE points 10 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
# This problem has a contact penalty stiffness of 0.01 and will converge when the 
# residual is less than 0.1 or 6 iterations have been performed (whichever comes first).
# This has been chosen for speed and stability. 
fem solve error 0.1 iterate 6 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;example 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
}

# write out new rectus_patella geometry
#fem def nodes;w;def_simulate_rectus_patella_$i;example as simulate_rectus_patella reg 1
#fem exp nodes;def_simulate_rectus_patella_$i as simulate_rectus_patella reg 1

# export tibia rotated mesh
#fem exp nodes;def_tibia_cust$i as tibia_cust reg 3

# export femur rotated mesh
#fem exp nodes;def_femur_cust$i as femur_cust reg 2

}

# Close the files.
CORE::close file1;


Additional testing commands:



Files used by this example are:

Name                                   Modified     Size

example_d43.com 23-Apr-2008 11k backup.ipacti 28-Feb-2005 1.0k contact.ipbase 28-Feb-2005 10k contact.ipcont 25-Sep-2008 628 contact.ippara 21-Feb-2008 6.0k contact.ipregi 28-Feb-2005 93 contact.ipsolv 16-Aug-2010 2.7k contact.ipsolv.old 13-Apr-2007 2.5k contact.ircoor 28-Feb-2005 1.8k contact_force_ratio.dat 28-Feb-2005 4.7k contact_new.ippara 21-Feb-2008 6.1k def_simulate_rectus_patella_10.exnode 28-Feb-2005 334k def_simulate_rectus_patella_30.exnode 28-Feb-2005 334k def_simulate_rectus_patella_70.exnode 28-Feb-2005 334k def_tibia_cust10.exnode 28-Feb-2005 118k def_tibia_cust30.exnode 28-Feb-2005 118k def_tibia_cust70.exnode 28-Feb-2005 118k femur_cust.ipelem 28-Feb-2005 453k femur_cust.ipelfb 28-Feb-2005 230k femur_cust.ipequa 28-Feb-2005 2.2k femur_cust.ipfibr 28-Feb-2005 616k femur_cust.ipnode 28-Feb-2005 895k fit_gauss.com 28-Feb-2005 1.2k grid_groups.ipgrgp 28-Feb-2005 125k grid_groups_multiple.ipgrgp 28-Feb-2005 46k kinematics.txt 28-Feb-2005 16k material.ipcell 28-Feb-2005 6.2k multiple_cell.ipcell 28-Feb-2005 9.4k multiple_cell.ipequa 28-Feb-2005 1.6k multiple_cell.ipmatc 28-Feb-2005 13k multiple_cell.ipsolv 13-Apr-2007 2.5k multiple_muscle.ipacti 28-Feb-2005 1.0k multiple_muscle.ipgrid 28-Feb-2005 622 multiple_muscle.ipmate 28-Feb-2005 1.7k new.ipfibr 28-Feb-2005 287k plot_results.gp 28-Feb-2005 2.0k polezero_nearly_incompressible.cml 28-Feb-2005 52k st_venant_kirchoff.cml 28-Feb-2005 10k tendon_force_ratio.dat 28-Feb-2005 4.7k test.ipfit 13-Apr-2007 2.2k test_output.com 28-Feb-2005 0 thesis.ipcell 28-Feb-2005 6.2k thesis.ipmatc 28-Feb-2005 2.1k thesis.ipmate.bak 28-Feb-2005 6.5k thesis_final_mesh.ipelem 28-Feb-2005 31k thesis_final_mesh.ipelfb 28-Feb-2005 16k thesis_final_mesh.ipelfd 28-Feb-2005 16k thesis_final_mesh.ipequa 28-Feb-2005 2.2k thesis_final_mesh.ipfibr 28-Feb-2005 287k thesis_final_mesh.ipfiel 28-Feb-2005 287k thesis_final_mesh.ipinit 28-Feb-2005 8.8k thesis_final_mesh.ipnode 28-Feb-2005 332k thesiscell.ipequa 28-Feb-2005 1.0k tibia_cust.ipelem 28-Feb-2005 96k tibia_cust.ipelfb 28-Feb-2005 47k tibia_cust.ipequa 28-Feb-2005 2.2k tibia_cust.ipfibr 28-Feb-2005 132k tibia_cust.ipnode 28-Feb-2005 197k view.com 28-Feb-2005 5.9k

Download the entire example:

Name                      Modified     Size

examples_d_d4_d43.tar.gz 17-Aug-2010 802k

Testing status by version:

StatusTestedReal time (s)
rs6000-aix
cm64SuccessSun Mar 13 05:45:20 20051081

Testing status by file:


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

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


CMISS Help / Examples / d / d4 / d43