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