##########################################################################
# This problem shows a cylinder tied (or glued) to an elastic substrate
##########################################################################
# If the example path is not set, default to current directory
if (!defined $example) {$example = "./";}
# Drop off the trailing / in the example path
$chopped = chop $example;
if ($chopped ne "/") {$example .= $chopped;}
# Define parameters, regions, coordinate system and bases
fem def para;r;hertz;example
fem def region;r;two;example
fem def coord 3,1
fem def bases;r;collapsed;example
# Read in half-space described by trilinear elements
fem def nodes;r;half;example reg 1
fem def elem;r;half;example reg 1
# Define fibre directions
fem def fibre;r;half;example reg 1
fem def elem;r;half;example fibre reg 1
# Define field to store pressure solution
fem def field;r;half;example reg 1
fem def elem;r;half;example field reg 1
# Read in rigid cylinder described by bicubic-linear elements
fem def nodes;r;cylinder;example reg 2
fem def elem;r;cylinder;example reg 2
# Define fibre directions
fem def fibre;r;cylinder;example reg 2
fem def elem;r;cylinder;example fibre reg 2
## Export slave (half-space)
#fem export nodes;half as half reg 1
#fem export elem;half as half reg 1
## Export Master (cylinder)
#fem export nodes;cylinder as cylinder reg 2
#fem export elem;cylinder as cylinder reg 2
# Node groups for fixed boundary conditions
fem group nodes 101..108,233,235,297..299,301,303,305,427..431,433,435,437,439,441,443,445,687..694 as fixedx reg 1
fem group nodes 105..110,113..114,117..118,121..122,125..126,129..130,133..134,137..138,141..142,145..146,149..150,153..154,157..158,161..162 as fixedy1 reg 1
fem group nodes 165..166,169..170,173..174,177..178,181..182,185..186,189..190,193..194,197..198,201..202,205..206,209..210,213..214,217..218,221..222,225..226,229..230 as fixedy2 reg 1
fem group nodes 101..694 as fixedz reg 1
fem group nodes 1..20 as cylinderfixedz reg 2
fem group nodes 1..20 as cylinderfixedy reg 2
fem group nodes 1..20 as cylinderfixedx reg 2
# Define equations
fem def equa;r;half;example reg 1
fem def equa;r;cylinder;example reg 2
# Define material properties
fem def mat;r;half;example reg 1
fem def mat;r;cylinder;example reg 2
# Define initial fixed boundary conditions
fem def init;r;half;example reg 1
fem def init;r;cylinder;example reg 2
## Export deformed Master
#fem export nodes;cylinder_def as cylinder_def reg 2 offset 2000
#fem export elem;cylinder_def as cylinder_def reg 2 offset_elem 2000
# define contact parameters
fem def contact;r;contact;example
# Define solve for reg 1
fem def solve;r;half;example reg 1
#######################################Define cell problem
# Set up a grid at gauss scheme
# Trilinear basis to interpolate between grid points
fem def grid;r;half;example gauss class 2
# Sets up the grid geometry
fem update grid geometry
# Group grid points in the elements 101 to 356 as ALL_GRID
fem group grid element 101..356 as ALL_GRID
# Define a cellml equation
# Options chosen allow use of structures for cardiac electrical activation
fem def equa;r;cellml;example class 2
# Needed to setup
fem def mate;d class 2
# Define cellml file
fem define cell;r;half;example class 2
# Define cell material
fem def mate;r;half;example cell class 2
# Default init conditions for cell class 2
fem def init;d class 2
# Setup solve for cellml class 2
fem def solv;r;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
for($k=1;$k<18;$k++) # load steps
{
if (($k>1)&&($k<5))
{
fem change nodes rotate by -3 about 100,-10,10 axis 0,0,1 nodes cylinderfixedy reg 2
}
if (($k>4)&&($k<11))
{
fem change nodes rotate by 3 about 100,-10,10 axis 0,0,1 nodes cylinderfixedy reg 2
}
if (($k>10)&&($k<14))
{
fem change nodes rotate by -3 about 100,-10,10 axis 0,0,1 nodes cylinderfixedy reg 2
}
if (($k>13)&&($k<16))
{
fem change nodes translate by 0,5,0 nodes cylinderfixedy reg 2
}
if (($k>15)&&($k<18))
{
fem change nodes translate by 0,-5,0 nodes cylinderfixedy reg 2
}
$j=0;
$CONVERGED=0;
while ($CONVERGED==0)
{
$j++;
######################################Projection
# Set contact Xi points on specified faces as 8x8 points
fem def xi;c contact_tied faces 678,746,848,865,979,988,997,1006 points 4
# Define data at xi positions
fem def data;c from_xi reg 1
# Update slave info
fem update data field to slave
# Project onto target face
fem def xi;c closest_face faces 1077,1081,1085,1089 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
if ($j>1)
{
FEM up geom from sol YP_index 3 to 1..3 reg 1
}
######################################Mechanics problem
# Solve finite elasticity/contact problem
fem solve error 0.00001 iterate 25 reg 1
FEM up geom from sol to 1..3 reg 1
}
# Place initial geometry YP(ny,3) into XP
FEM up geom from sol YP_index 3 to 1..3 reg 1
## Stresses at gauss points
#fem up gauss stress fibre reg 1
#fem def gauss;w;stress_$k as half number 3 reg 1
# Update XP from YP
FEM up geom from sol to 1..3 reg 1
################################################
# Place initial geometry YP(ny,6) into XP
FEM up geom from sol YP_index 6 to 7..9 reg 1
## Write out new geometry
#fem def nodes;w;half_def_$k as half reg 1
#fem def elem;w;half_def_$k as half reg 1
#fem exp nodes;half_def_$k as half reg 1
## Export deformed Master
#fem export nodes;cylinder_def_$k as cylinder_def reg 2 offset 2000
#fem export elem;cylinder_def_$k as cylinder_def reg 2 offset_elem 2000
}
#load steps
