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