#### - Automatic generation of parametric hip models - ####
### Robert J. Cooper, University of Leeds, August 2017 ###

# basic definitions - leave alone
Hip = '01R'
centx = 0
centy = 0
centz = 0
modelx = 0
modely = 0
modelz = 0
HR = 25
numloftslices = 5
loftslices = [12.000000000000, 15.000000000000, 19.000000000000, 23.000000000000, 28.000000000000, ] 
ACT = 1*HR/22.95 
Frad = 1.1*HR
FR = 0.5*Frad	
bonethick = 1.5*HR/22.95
fcartoffset = 5*HR/22.95

## EDITABLE PARAMETERS ##

# (A) Femoral neck parameters - VARY these using matlab script and paste them in here
NRa = [23.816899208424, 22.364960237427, 20.550226483188, 19.133974596216, 23.950000000000, ]
NRb = [24.625000000000, 21.692000000000, 18.300000000000, 16.252500000000, 15.760000000000, ]
NRx = [1.072253460255, 1.034375000000, 0.714928993933, -0.150383733180, 0.000000000000, ]
NRy = [0.500000000000, 2.218224345902, 2.668151329163, 0.852868531952, 0.000000000000, ]
NA = [65.000000000000, 65.000000000000, 75.000000000000, 70.000000000000, 10.000000000000, ]
NL = [0.264400000000, -6.159900000000, -14.725600000000, -23.291200000000, -34.022400000000, ]
Cam_ang = [25.000000000000, 15.000000000000, ]
Cam_rad = [99.250000000000, 93.250000000000, ]
Cam_width = [45.000000000000, 31.000000000000, ]
CamRad = 96.250000000000 # *
Angle = 20.000000000000 # *
Width = 38.000000000000 # *


# (B) Acetabular parameters
CCD = 40 # * VARY THIS *
CE = 30  # * VARY THIS *
AV = 20  # * VARY THIS *

# (C) Seeds for meshing
paracartSeed = 0.0107*HR 
labrumSeed = 0.0214*HR 
parfemurseed = 0.085*HR 
fcartseed = 0.0425*HR 

## EDITABLE PARAMETERS ##


## MORE MODEL INPUTS ##

# labrum parameters - leave alone
eang = 19.8769 # gives Apc = 0.33
# eang = 23.9538 # gives Apc = 0.297, decrease bone 10%
# eang = 15.9024 # gives Apc = 0.363, increase bone 10%
# note Apc means percentage of sphere represented by acebone

# fix LH = value, in order to vary labrum / bone ratio with same coverage
# LH = -2.81297524651976 # fixed to value with Apc = 0.33

# keeps labrum length scaled to lablength when HR=25
lablength = 7 # baselength
# lablength = 6.3 # decrease lab 10%
# lablength = 7.7 # increase lab 10%
LH = (HR/25*lablength) - 0.5 * ( (Frad+ACT+bonethick)*sin(radians(eang)) + Frad*sin(radians(eang)) )

# further model parameters
angtheta = 29.463 # gives Hpc = 0.2541	
rotcent = 0
ManOp_OrthLab = 1
ManOp_OrthCart = 1
ManOp_FlexAng = 70
ManOp_irot = 0.0872665
ManOp_acartdelam = 0
ManOp_initialxyz = [-1.4, -0.5, 0.5]  

## END MODEL INPUTS ##


# import abaqus modules 
from part import *
from material import *
from section import *
from assembly import *
from step import *
from interaction import *
from load import *
from mesh import *
from optimization import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *

# stop jnl file from using getSequenceFromMask
session.journalOptions.setValues(replayGeometry=COORDINATE, recoverGeometry=COORDINATE)

# assign variable to shorten code
M = mdb.models['Model-1']

# define linearly elastic isotropic cartilage	
M.Material(name='cartilage-mat')
M.materials['cartilage-mat'].Elastic(table=((12.0, 0.4), ))

# articular layer 
M.Material(name='cartilage1-mat')
M.materials['cartilage1-mat'].Elastic(table=(
	(15.0, 15.0, 9.0, 
	0.4, 0.4, 0.4, 
	4.29, 4.29, 4.29), 
	), type=ENGINEERING_CONSTANTS)	

# middle layer - equivalent to isotropic 	
M.Material(name='cartilage2-mat')
M.materials['cartilage2-mat'].Elastic(table=(
	(12.0, 12.0, 12.0, 
	0.4, 0.4, 0.4, 
	4.29, 4.29, 4.29), 
	), type=ENGINEERING_CONSTANTS)	
	
# base layer adjacent to bone
M.Material(name='cartilage3-mat')
M.materials['cartilage3-mat'].Elastic(table=(
	(9.0, 9.0, 15.0, 
	0.4, 0.4, 0.4, 
	4.29, 4.29, 4.29), 
	), type=ENGINEERING_CONSTANTS)		
	
# make all layers the same if unselected orthotropic cartilage option
if ManOp_OrthCart == 0:
	M.materials['cartilage1-mat'].Elastic(table=(
		(12.0, 12.0, 12.0, 
		0.4, 0.4, 0.4, 
		4.29, 4.29, 4.29), 
		), type=ENGINEERING_CONSTANTS)	
		
	M.materials['cartilage3-mat'].Elastic(table=(
		(12.0, 12.0, 12.0, 
		0.4, 0.4, 0.4, 
		4.29, 4.29, 4.29), 
		), type=ENGINEERING_CONSTANTS)		

if ManOp_OrthLab == 1:
# define transversely isotropic labrum
	M.Material(name='labrum-mat')
	M.materials['labrum-mat'].Elastic(table=(
		(12.0, 12.0, 20.0, 
		0.4, 0.4, 0.4, 
		5.24, 5.24, 5.24), 
		), type=ENGINEERING_CONSTANTS)		
else:
	# define linearly elastic isotropic labrum	
	M.Material(name='labrum-mat')
	M.materials['labrum-mat'].Elastic(table=((20.0, 0.4), ))

# create sections
M.HomogeneousSolidSection(material='cartilage-mat', name=
    'cartilage-sec', thickness=None)
M.HomogeneousSolidSection(material='cartilage1-mat', name=
    'cartilage1-sec', thickness=None)
M.HomogeneousSolidSection(material='cartilage2-mat', name=
    'cartilage2-sec', thickness=None)
M.HomogeneousSolidSection(material='cartilage3-mat', name=
    'cartilage3-sec', thickness=None)	
M.HomogeneousSolidSection(material='labrum-mat', name=
    'labrum-sec', thickness=None)


def makesimpleace():
	## - Acetabular bone part - ##	
	
	# make parace from sketches  
	M.ConstrainedSketch(name='__profile__', sheetSize=200.0)
	M.sketches['__profile__'].ConstructionLine(point1=(0.0, 
		-100.0), point2=(0.0, 100.0))
	M.sketches['__profile__'].Arc3Points(
		point1=((Frad+ACT+bonethick)*cos(radians(90)),(Frad+ACT+bonethick)*sin(radians(90))), 
		point2=((Frad+ACT+bonethick)*cos(radians(eang)),(Frad+ACT+bonethick)*sin(radians(eang))), 
		point3=((Frad+ACT+bonethick)*cos(radians(45)),(Frad+ACT+bonethick)*sin(radians(45))))
	M.sketches['__profile__'].Arc3Points(
		point1=((Frad+ACT)*cos(radians(90)),(Frad+ACT)*sin(radians(90))), 
		point2=((Frad+ACT)*cos(radians(eang)),(Frad+ACT)*sin(radians(eang))), 
		point3=((Frad+ACT)*cos(radians(45)),(Frad+ACT)*sin(radians(45))))
	M.sketches['__profile__'].Line(
		point1=((Frad+ACT)*cos(radians(eang)),(Frad+ACT)*sin(radians(eang))), 
		point2=((Frad+ACT+bonethick)*cos(radians(eang)),(Frad+ACT+bonethick)*sin(radians(eang))))		
	M.Part(dimensionality=THREE_D, name='parace', type=
		ANALYTIC_RIGID_SURFACE)
	M.parts['parace'].AnalyticRigidSurfRevolve(sketch=M.sketches['__profile__'])
	del M.sketches['__profile__']

	# instance acebone
	M.rootAssembly.Instance(dependent=ON, name='parace-1'
		, part=M.parts['parace'])

	# create acebone surface	
	M.rootAssembly.Surface(name='parace-sur', side1Faces=
		M.rootAssembly.instances['parace-1'].faces.findAt(
		((0, Frad+ACT, 0), )))	

	# no mesh required because it is analytical rigid

	## - Labrum part - ##	
		
	# labrum 
	M.ConstrainedSketch(name='__profile__', sheetSize=200.0)
	M.sketches['__profile__'].ConstructionLine(point1=(0.0, 
		-100.0), point2=(0.0, 100.0))
	M.sketches['__profile__'].Line(
		point1=(Frad*cos(radians(eang)), Frad*sin(radians(eang))), 
		point2=(Frad, -LH))
	M.sketches['__profile__'].Line(
		point1=(Frad, -LH), 
		point2=((Frad+ACT+bonethick)*cos(radians(eang)), (Frad+ACT+bonethick)*sin(radians(eang))))
	M.sketches['__profile__'].Line(
		point1=((Frad+ACT+bonethick)*cos(radians(eang)), (Frad+ACT+bonethick)*sin(radians(eang))), 
		point2=(Frad*cos(radians(eang)), Frad*sin(radians(eang))))
	M.Part(dimensionality=THREE_D, name='labrum', type=DEFORMABLE_BODY)
	M.parts['labrum'].BaseSolidRevolve(angle=360.0, flipRevolveDirection=OFF, sketch=M.sketches['__profile__'])
	del M.sketches['__profile__']

	# instance labrum
	M.rootAssembly.Instance(dependent=ON, name='labrum-1', 
		part=M.parts['labrum'])	
		
	# assign material orientation to labrum
	M.parts['labrum'].MaterialOrientation(
		additionalRotationField='', additionalRotationType=ROTATION_NONE, angle=0.0
		, axis=AXIS_3, fieldName='', localCsys=
		M.parts['labrum'].datums[M.parts['labrum'].DatumCsysByTwoLines(CYLINDRICAL, 
		line1=M.parts['labrum'].datums[M.parts['labrum'].DatumAxisByPrincipalAxis(principalAxis=XAXIS).id], 
		line2=M.parts['labrum'].datums[M.parts['labrum'].DatumAxisByPrincipalAxis(principalAxis=ZAXIS).id], 
		name='Datum csys-1').id], orientationType=SYSTEM, 
		region=Region(cells=M.parts['labrum'].cells.getByBoundingSphere( (0,0,0), 500)), stackDirection=STACK_ORIENTATION)		
	
	# assign section to labrum	
	M.parts['labrum'].SectionAssignment(offset=0.0, offsetField=
		'', offsetType=MIDDLE_SURFACE, region=Region(
		cells=M.parts['labrum'].cells.getByBoundingSphere( (0,0,0), 500)),
		sectionName='labrum-sec', thicknessAssignment=FROM_SECTION)	
		
	# labrum all surface
	M.rootAssembly.Surface(name='labrum-sur', side1Faces=
		M.rootAssembly.instances['labrum-1'].faces.getByBoundingSphere( (0,0,0), 500))

	# mesh labrum
	M.parts['labrum'].seedPart(deviationFactor=0.1, 
		minSizeFactor=0.1, size=labrumSeed)
	M.parts['labrum'].generateMesh()

	# give labrum elements enhanced hourglass control
	M.parts['labrum'].setElementType(elemTypes=(ElemType(
		elemCode=C3D8R, elemLibrary=STANDARD, secondOrderAccuracy=OFF, 
		kinematicSplit=AVERAGE_STRAIN, hourglassControl=ENHANCED, 
		distortionControl=DEFAULT), ElemType(elemCode=C3D6, elemLibrary=STANDARD), 
		ElemType(elemCode=C3D4, elemLibrary=STANDARD)), 
		regions=(M.parts['labrum'].cells.findAt(((Frad, -LH,0), )), ))			
	
	## - Acetabular cartilage part - ##	

	# make acart layers from from sketches
	acfact1 = 0.333
	acfact2 = 0.666
	acfact3 = 1
	
	M.ConstrainedSketch(name='__profile__', sheetSize=200.0)
	M.sketches['__profile__'].ConstructionLine(point1=(0.0, 
		-100.0), point2=(0.0, 100.0))
	M.sketches['__profile__'].Arc3Points(
		point1=((Frad+acfact1*ACT)*cos(radians(90)),(Frad+acfact1*ACT)*sin(radians(90))), 
		point2=((Frad+acfact1*ACT)*cos(radians(eang)),(Frad+acfact1*ACT)*sin(radians(eang))), 
		point3=((Frad+acfact1*ACT)*cos(radians(45)),(Frad+acfact1*ACT)*sin(radians(45))))
	M.sketches['__profile__'].Arc3Points(
		point1=((Frad)*cos(radians(90)),(Frad)*sin(radians(90))), 
		point2=((Frad)*cos(radians(eang)),(Frad)*sin(radians(eang))), 
		point3=((Frad)*cos(radians(45)),(Frad)*sin(radians(45))))
	M.sketches['__profile__'].Line(
		point1=((Frad+acfact1*ACT)*cos(radians(90)),(Frad+acfact1*ACT)*sin(radians(90))), 
		point2=((Frad)*cos(radians(90)),(Frad)*sin(radians(90))))
	M.sketches['__profile__'].Line(
		point1=((Frad)*cos(radians(eang)),(Frad)*sin(radians(eang))), 
		point2=((Frad+acfact1*ACT)*cos(radians(eang)),(Frad+acfact1*ACT)*sin(radians(eang))))
	M.Part(dimensionality=THREE_D, name='acart1', type=
		DEFORMABLE_BODY)
	M.parts['acart1'].BaseSolidRevolve(angle=360.0, 
		flipRevolveDirection=OFF, sketch=
		M.sketches['__profile__'])
	del M.sketches['__profile__']
	
	M.ConstrainedSketch(name='__profile__', sheetSize=200.0)
	M.sketches['__profile__'].ConstructionLine(point1=(0.0, 
		-100.0), point2=(0.0, 100.0))
	M.sketches['__profile__'].Arc3Points(
		point1=((Frad+acfact2*ACT)*cos(radians(90)),(Frad+acfact2*ACT)*sin(radians(90))), 
		point2=((Frad+acfact2*ACT)*cos(radians(eang)),(Frad+acfact2*ACT)*sin(radians(eang))), 
		point3=((Frad+acfact2*ACT)*cos(radians(45)),(Frad+acfact2*ACT)*sin(radians(45))))
	M.sketches['__profile__'].Arc3Points(
		point1=((Frad+acfact1*ACT)*cos(radians(90)),(Frad+acfact1*ACT)*sin(radians(90))), 
		point2=((Frad+acfact1*ACT)*cos(radians(eang)),(Frad+acfact1*ACT)*sin(radians(eang))), 
		point3=((Frad+acfact1*ACT)*cos(radians(45)),(Frad+acfact1*ACT)*sin(radians(45))))
	M.sketches['__profile__'].Line(
		point1=((Frad+acfact2*ACT)*cos(radians(90)),(Frad+acfact2*ACT)*sin(radians(90))), 
		point2=((Frad+acfact1*ACT)*cos(radians(90)),(Frad+acfact1*ACT)*sin(radians(90))))
	M.sketches['__profile__'].Line(
		point1=((Frad+acfact1*ACT)*cos(radians(eang)),(Frad+acfact1*ACT)*sin(radians(eang))), 
		point2=((Frad+acfact2*ACT)*cos(radians(eang)),(Frad+acfact2*ACT)*sin(radians(eang))))
	M.Part(dimensionality=THREE_D, name='acart2', type=
		DEFORMABLE_BODY)
	M.parts['acart2'].BaseSolidRevolve(angle=360.0, 
		flipRevolveDirection=OFF, sketch=
		M.sketches['__profile__'])
	del M.sketches['__profile__']	
	
	M.ConstrainedSketch(name='__profile__', sheetSize=200.0)
	M.sketches['__profile__'].ConstructionLine(point1=(0.0, 
		-100.0), point2=(0.0, 100.0))
	M.sketches['__profile__'].Arc3Points(
		point1=((Frad+acfact3*ACT)*cos(radians(90)),(Frad+acfact3*ACT)*sin(radians(90))), 
		point2=((Frad+acfact3*ACT)*cos(radians(eang)),(Frad+acfact3*ACT)*sin(radians(eang))), 
		point3=((Frad+acfact3*ACT)*cos(radians(45)),(Frad+acfact3*ACT)*sin(radians(45))))
	M.sketches['__profile__'].Arc3Points(
		point1=((Frad+acfact2*ACT)*cos(radians(90)),(Frad+acfact2*ACT)*sin(radians(90))), 
		point2=((Frad+acfact2*ACT)*cos(radians(eang)),(Frad+acfact2*ACT)*sin(radians(eang))), 
		point3=((Frad+acfact2*ACT)*cos(radians(45)),(Frad+acfact2*ACT)*sin(radians(45))))
	M.sketches['__profile__'].Line(
		point1=((Frad+acfact3*ACT)*cos(radians(90)),(Frad+acfact3*ACT)*sin(radians(90))), 
		point2=((Frad+acfact2*ACT)*cos(radians(90)),(Frad+acfact2*ACT)*sin(radians(90))))
	M.sketches['__profile__'].Line(
		point1=((Frad+acfact2*ACT)*cos(radians(eang)),(Frad+acfact2*ACT)*sin(radians(eang))), 
		point2=((Frad+acfact3*ACT)*cos(radians(eang)),(Frad+acfact3*ACT)*sin(radians(eang))))
	M.Part(dimensionality=THREE_D, name='acart3', type=
		DEFORMABLE_BODY)
	M.parts['acart3'].BaseSolidRevolve(angle=360.0, 
		flipRevolveDirection=OFF, sketch=
		M.sketches['__profile__'])
	del M.sketches['__profile__']	
	
	# make horse part from sketch
	fangle = acos(Frad/FR*2/5)
	M.ConstrainedSketch(name='__profile__', sheetSize=200.0)
	M.sketches['__profile__'].Line(
		point1=(-FR*cos(fangle), -Frad*4/5), 
		point2=(-FR*cos(fangle), -FR*sin(fangle)))
	M.sketches['__profile__'].Line(
		point1=(FR*cos(fangle), -Frad*4/5), 
		point2=(FR*cos(fangle), -FR*sin(fangle)))
	M.sketches['__profile__'].Line(
		point1=(-FR*cos(fangle), -Frad*4/5), 
		point2=(FR*cos(fangle), -Frad*4/5))
	M.sketches['__profile__'].Arc3Points(
		point1=(-FR*cos(fangle), -FR*sin(fangle)), 
		point2=(FR*cos(fangle), -FR*sin(fangle)), 
		point3=(0.0, FR))
	M.Part(dimensionality=THREE_D, name='horse', type=
		DEFORMABLE_BODY)
	M.parts['horse'].BaseSolidExtrude(depth=3*Frad, sketch=
		M.sketches['__profile__'])
	del M.sketches['__profile__']
		
	# make and instance paracart1 by cutting acart1 with horse
	M.rootAssembly.Instance(dependent=ON, name='acart1-1', 
		part=M.parts['acart1'])
	M.rootAssembly.Instance(dependent=ON, name='horse-1', 
		part=M.parts['horse'])
	M.rootAssembly.rotate(angle=90.0, axisDirection=(1.0, 0.0, 
		0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('horse-1', ))
	M.rootAssembly.translate(instanceList=('horse-1', ), 
		vector=(0.0, 1.5*Frad, -0.3*Frad))
	M.rootAssembly.InstanceFromBooleanCut(cuttingInstances=(
		M.rootAssembly.instances['horse-1'], ), 
		instanceToBeCut=M.rootAssembly.instances['acart1-1'], 
		name='paracart1', originalInstances=DELETE)

	# assign section to paracart1	
	M.parts['paracart1'].SectionAssignment(offset=0.0, offsetField=
		'', offsetType=MIDDLE_SURFACE, region=Region(
		cells=M.parts['paracart1'].cells.getByBoundingSphere( (0,0,0), 500)),
		sectionName='cartilage1-sec', thicknessAssignment=FROM_SECTION)		
		
	# partition paracart1 in XY plane
	M.parts['paracart1'].PartitionCellByDatumPlane(cells=
		M.parts['paracart1'].cells.findAt(
		(((Frad+acfact1*ACT)*cos(radians(45)), (Frad+acfact1*ACT)*sin(radians(45)), 0), )), 
		datumPlane=M.parts['paracart1'].datums
		[M.parts['paracart1'].DatumPlaneByPrincipalPlane(offset=0.0, principalPlane=XYPLANE).id])

	# make and instance paracart1 by cutting acart1 with horse
	M.rootAssembly.Instance(dependent=ON, name='acart2-1', 
		part=M.parts['acart2'])
	M.rootAssembly.Instance(dependent=ON, name='horse-1', 
		part=M.parts['horse'])
	M.rootAssembly.rotate(angle=90.0, axisDirection=(1.0, 0.0, 
		0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('horse-1', ))
	M.rootAssembly.translate(instanceList=('horse-1', ), 
		vector=(0.0, 1.5*Frad, -0.3*Frad))
	M.rootAssembly.InstanceFromBooleanCut(cuttingInstances=(
		M.rootAssembly.instances['horse-1'], ), 
		instanceToBeCut=M.rootAssembly.instances['acart2-1'], 
		name='paracart2', originalInstances=DELETE)

	# assign section to paracart2	
	M.parts['paracart2'].SectionAssignment(offset=0.0, offsetField=
		'', offsetType=MIDDLE_SURFACE, region=Region(
		cells=M.parts['paracart2'].cells.getByBoundingSphere( (0,0,0), 500)),
		sectionName='cartilage2-sec', thicknessAssignment=FROM_SECTION)		
		
	# partition paracart2 in XY plane
	M.parts['paracart2'].PartitionCellByDatumPlane(cells=
		M.parts['paracart2'].cells.findAt(
		(((Frad+acfact2*ACT)*cos(radians(45)), (Frad+acfact2*ACT)*sin(radians(45)), 0), )), 
		datumPlane=M.parts['paracart2'].datums
		[M.parts['paracart2'].DatumPlaneByPrincipalPlane(offset=0.0, principalPlane=XYPLANE).id])

	# make and instance paracart1 by cutting acart1 with horse
	M.rootAssembly.Instance(dependent=ON, name='acart3-1', 
		part=M.parts['acart3'])
	M.rootAssembly.Instance(dependent=ON, name='horse-1', 
		part=M.parts['horse'])
	M.rootAssembly.rotate(angle=90.0, axisDirection=(1.0, 0.0, 
		0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('horse-1', ))
	M.rootAssembly.translate(instanceList=('horse-1', ), 
		vector=(0.0, 1.5*Frad, -0.3*Frad))
	M.rootAssembly.InstanceFromBooleanCut(cuttingInstances=(
		M.rootAssembly.instances['horse-1'], ), 
		instanceToBeCut=M.rootAssembly.instances['acart3-1'], 
		name='paracart3', originalInstances=DELETE)

	# assign section to paracart3	
	M.parts['paracart3'].SectionAssignment(offset=0.0, offsetField=
		'', offsetType=MIDDLE_SURFACE, region=Region(
		cells=M.parts['paracart3'].cells.getByBoundingSphere( (0,0,0), 500)),
		sectionName='cartilage3-sec', thicknessAssignment=FROM_SECTION)		
		
	# partition paracart3 in XY plane
	M.parts['paracart3'].PartitionCellByDatumPlane(cells=
		M.parts['paracart3'].cells.findAt(
		(((Frad+acfact3*ACT)*cos(radians(45)), (Frad+acfact3*ACT)*sin(radians(45)), 0), )), 
		datumPlane=M.parts['paracart3'].datums
		[M.parts['paracart3'].DatumPlaneByPrincipalPlane(offset=0.0, principalPlane=XYPLANE).id])

	# merge layers
	M.rootAssembly.InstanceFromBooleanMerge(domain=GEOMETRY, 
		instances=(M.rootAssembly.instances['paracart1-1'], 
		M.rootAssembly.instances['paracart2-1'], 
		M.rootAssembly.instances['paracart3-1']), 
		keepIntersections=ON, name='paracart', originalInstances=DELETE)	
		
	# give acart elements enhanced hourglass control
	M.parts['paracart'].setElementType(elemTypes=(ElemType(
		elemCode=C3D8R, elemLibrary=STANDARD, secondOrderAccuracy=OFF, 
		kinematicSplit=AVERAGE_STRAIN, hourglassControl=ENHANCED, 
		distortionControl=DEFAULT), ElemType(elemCode=C3D6, elemLibrary=STANDARD), 
		ElemType(elemCode=C3D4, elemLibrary=STANDARD)), 
		regions=(M.parts['paracart'].cells.getByBoundingSphere( (0,0,0), 500), ))	
		
	# assign material orientation to acart using equation for points lying on sphere
	M.parts['paracart'].MaterialOrientation(
		additionalRotationField='', additionalRotationType=ROTATION_NONE, angle=0.0
		, axis=AXIS_1, flipNormalDirection=False, flipPrimaryDirection=False, 
		normalAxisDefinition=SURFACE, normalAxisDirection=AXIS_3, normalAxisRegion=
		Region(side1Faces=M.parts['paracart'].faces.findAt(
		((Frad*cos(radians(60)),Frad*sin(radians(60))*sin(radians(30)),Frad*sin(radians(60))*cos(radians(30))), ), 
		)), orientationType=DISCRETE, 
		primaryAxisDefinition=VECTOR, primaryAxisDirection=AXIS_1, 
		primaryAxisVector=(0.0, 0.0, 1.0), 
		region=Region(cells=M.parts['paracart'].cells.getByBoundingSphere( (0,0,0), 500)), stackDirection=STACK_ORIENTATION)			
		
	# create paracart all surface 
	M.rootAssembly.Surface(name='paracart-sur', side1Faces=
		M.rootAssembly.instances['paracart-1'].faces.getByBoundingSphere( (0,0,0), 500))
		
	# create paracart mesh
	M.parts['paracart'].seedPart(deviationFactor=0.1, 
		minSizeFactor=0.1, size=paracartSeed)
	M.parts['paracart'].generateMesh()			

	# Make acetabular axes match	
	M.rootAssembly.rotate(angle=90.0, axisDirection=(0.0, 0.0, 
		1.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('paracart-1', 'labrum-1', 'parace-1'))
	M.rootAssembly.rotate(angle=90.0, axisDirection=(0.0, 1.0, 
		0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('paracart-1', 'labrum-1', 'parace-1'))		
	# Set CE angle
	M.rootAssembly.rotate(angle=-90+CE, axisDirection=(0.0, 1.0, 
		0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('paracart-1', 'labrum-1', 'parace-1'))

	if Hip[-1] == 'R':
		# RIGHT AV angle	
		M.rootAssembly.rotate(angle=-360+AV, axisDirection=(1.0, 0.0, 
			0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('paracart-1', 'labrum-1', 'parace-1'))			

	if Hip[-1] == 'L':			
		# LEFT AV angle	 	
		M.rootAssembly.rotate(angle=180, axisDirection=(0.0, 0.0, 
			1.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('paracart-1', 'labrum-1', 'parace-1'))		
		M.rootAssembly.rotate(angle=AV, axisDirection=(1.0, 0.0, 
			0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('paracart-1', 'labrum-1', 'parace-1'))		

		
# function to create and instance femur 
# includes cartilage layer  	
def makeparfemur():	
	
	## - Head part - ##	
	
	# make rigidhead top from sketches 
	M.ConstrainedSketch(name='__profile__', sheetSize=200.0)

	M.sketches['__profile__'].ConstructionLine(point1=(0.0, 
		-100.0), point2=(0.0, 100.0))
		
	M.sketches['__profile__'].Arc3Points(
		point1=(0.0, HR), 
		point2=(HR*cos(radians(angtheta)), HR*sin(radians(angtheta))), 
		point3=(HR*cos(radians(-80)), -HR*sin(radians(-80))))	
		
	M.Part(dimensionality=THREE_D, name='headtop', type=
		DISCRETE_RIGID_SURFACE)
		
	M.parts['headtop'].BaseShellRevolve(angle=360.0, 
		flipRevolveDirection=OFF, sketch=M.sketches['__profile__'])
	del M.sketches['__profile__']

	# instance head top
	M.rootAssembly.Instance(dependent=ON, name='headtop-1', 
		part=M.parts['headtop'])
	
	# Make axes match Bergmann's axes - this will now match with ipfemur	
	M.rootAssembly.rotate(angle=90.0, axisDirection=(0.0, 0.0, 
		1.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('headtop-1', ))
	M.rootAssembly.rotate(angle=90.0, axisDirection=(0.0, 1.0, 
		0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('headtop-1', ))		
		
	# make rigidhead neck by lofting through wires
	
	# create wire list
	wires = []
	for ii in range(0,numloftslices):
		wires.append('neckw%d' % ii)
	
	for ii in range(0,numloftslices):	
		# sketch wires
		M.ConstrainedSketch(name='__profile__', sheetSize=200.0)
		M.sketches['__profile__'].EllipseByCenterPerimeter(
			axisPoint1=(NRa[ii], 0.0), axisPoint2=(0.0, NRb[ii]), center=(0.0, 0.0))
		M.Part(dimensionality=THREE_D, name=wires[ii], type=
			DISCRETE_RIGID_SURFACE)
		M.parts[wires[ii]].BaseWire(sketch=M.sketches['__profile__'])
		del M.sketches['__profile__']
	
		# instance neck wires 
		M.rootAssembly.Instance(dependent=ON, name=wires[ii]+'-1', 
		part=M.parts[wires[ii]])
		
		# shift to correct neck length positions 
		M.rootAssembly.translate(instanceList=(wires[ii]+'-1', ), 
			vector=(0.0, 0.0, NL[ii]))	
			
		# rotate to correct orientation 	
		M.rootAssembly.rotate(angle=NA[ii], axisDirection=(0.0, 0.0, 
			1.0), axisPoint=(0.0, 0.0, 0.0), instanceList=(wires[ii]+'-1', ))	

		# shift to correct centre positions 
		M.rootAssembly.translate(instanceList=(wires[ii]+'-1', ), 
		vector=(NRx[ii], NRy[ii], 0.0))	

	# merge to create parfemur
	M.rootAssembly.InstanceFromBooleanMerge(domain=GEOMETRY, 
		instances=(M.rootAssembly.instances['headtop-1'], 
		M.rootAssembly.instances[wires[0]+'-1'],
		M.rootAssembly.instances[wires[1]+'-1'], 
		M.rootAssembly.instances[wires[2]+'-1'], 
		M.rootAssembly.instances[wires[3]+'-1'], 
		M.rootAssembly.instances[wires[-1]+'-1']), 
		name='parfemur', originalInstances=DELETE)
		
	# loft to complete head
	M.parts['parfemur'].ShellLoft(loftsections=(
			(M.parts['parfemur'].edges.findAt(coordinates=(0.0, HR*cos(radians(angtheta)), HR*sin(radians(angtheta)))), ),
			(M.parts['parfemur'].edges.findAt(coordinates=(NRa[0]*cos(radians(NA[0]))+NRx[0], NRa[0]*sin(radians(NA[0]))+NRy[0], NL[0])), ), 
			(M.parts['parfemur'].edges.findAt(coordinates=(NRa[1]*cos(radians(NA[1]))+NRx[1], NRa[1]*sin(radians(NA[1]))+NRy[1], NL[1])), ),
			(M.parts['parfemur'].edges.findAt(coordinates=(NRa[2]*cos(radians(NA[2]))+NRx[2], NRa[2]*sin(radians(NA[2]))+NRy[2], NL[2])), ),
			(M.parts['parfemur'].edges.findAt(coordinates=(NRa[3]*cos(radians(NA[3]))+NRx[3], NRa[3]*sin(radians(NA[3]))+NRy[3], NL[3])), ),
			(M.parts['parfemur'].edges.findAt(coordinates=(NRa[-1]*cos(radians(NA[-1]))+NRx[-1], NRa[-1]*sin(radians(NA[-1]))+NRy[-1], NL[-1])), )),
			startCondition=NONE, endCondition=NONE)	
			# head top, ellipses adjusted by rotation and centre - order of ellipses obviously matters! 

	# Set CCD angle
	if Hip[-1] == 'L':
		M.rootAssembly.rotate(angle=CCD, axisDirection=
			(0.0, 1.0, 0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('parfemur-1', ))
	if Hip[-1] == 'R':
		M.rootAssembly.rotate(angle=-CCD, axisDirection=
			(0.0, 1.0, 0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('parfemur-1', ))			
		
	# create parameterised fcart as deformable copy of femur
	M.Part(name='parfcart', objectToCopy=M.parts['parfemur'])
	M.parts['parfcart'].setValues(space=THREE_D, type=DEFORMABLE_BODY)
		
	DatPFC = M.parts['parfcart'].DatumPlaneByPrincipalPlane(offset=-fcartoffset, 
		principalPlane=XYPLANE).id 	

	M.parts['parfcart'].PartitionFaceByDatumPlane(datumPlane=
		M.parts['parfcart'].datums[DatPFC], 
		faces=M.parts['parfcart'].faces)	
	
	# instancing parfcart with CCD rotation	
	M.rootAssembly.Instance(dependent=ON, name='parfcart-1'
		, part=M.parts['parfcart'])	

	if Hip[-1] == 'L':	
		M.rootAssembly.rotate(angle=CCD, axisDirection=(0.0, 1.0, 
			0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('parfcart-1', ))	
	if Hip[-1] == 'R':	
		M.rootAssembly.rotate(angle=-CCD, axisDirection=(0.0, 1.0, 
			0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('parfcart-1', ))				
			
	# set mesh controls to quad
	M.parts['parfcart'].setMeshControls(elemShape=QUAD, regions=M.parts['parfcart'].faces)
	M.parts['parfemur'].setMeshControls(elemShape=QUAD, regions=M.parts['parfemur'].faces)	
	
	# mesh parfemur
	M.parts['parfemur'].seedPart(deviationFactor=0.1, 
		minSizeFactor=0.1, size=parfemurseed)
	M.parts['parfemur'].generateMesh()	
	
	# find some coordinates on parfcart faces we want to mesh 
	# - i.e. above our partition which is definitley below y=0 
	rpf = M.parts['parfcart'].faces.getClosest(((0,HR,0),(0,0,HR),))	

	# mesh required region for fcart - requires seeding to work
	M.parts['parfcart'].seedPart(deviationFactor=0.1, 
		minSizeFactor=0.1, size=fcartseed)
	M.parts['parfcart'].generateMesh(regions=
		M.parts['parfcart'].faces.findAt((rpf[0][1], ), (rpf[1][1], ), ))	

	# mesh offset for cartilage	
	M.parts['parfcart'].generateMeshByOffset(meshType=SOLID, 
		numLayers=2, offsetDirection=OUTWARD, region=Region(
		side1Elements=M.parts['parfcart'].elements), totalThickness=1.0)
		
	# delete geometry features
	M.parts['parfcart'].deleteFeatures(('Part From Merge-1', ))

	# give elements a set		
	M.parts['parfcart'].Set(elements=
		M.parts['parfcart'].elements, name='parfcartOffset-1')
		
	# assign cartilage section	
	M.parts['parfcart'].SectionAssignment(offset=0.0, 
		offsetField='', offsetType=MIDDLE_SURFACE, region=Region(
		elements=M.parts['parfcart'].elements), 
		sectionName='cartilage-sec', thicknessAssignment=FROM_SECTION)		
		

# function to set up initial interactions
def setinteractions():
	M.rootAssembly.regenerate()	

	# surface for parameterised fcart	
	M.rootAssembly.Surface(
		face1Elements=M.rootAssembly.instances['parfcart-1'].elements, 
		face2Elements=M.rootAssembly.instances['parfcart-1'].elements,
		face3Elements=M.rootAssembly.instances['parfcart-1'].elements,
		face4Elements=M.rootAssembly.instances['parfcart-1'].elements,
		face5Elements=M.rootAssembly.instances['parfcart-1'].elements,
		face6Elements=M.rootAssembly.instances['parfcart-1'].elements,
		name='parfcart-sur')	
	
	# create surface for parameterised femur by selecting all faces 
	M.rootAssembly.Surface(name='parfemur-sur', side1Faces=
		M.rootAssembly.instances['parfemur-1'].faces) 	

	# define frictionless surface interaction property
	M.ContactProperty('IntProp-1')
	M.interactionProperties['IntProp-1'].TangentialBehavior(
		formulation=FRICTIONLESS)
	M.interactionProperties['IntProp-1'].NormalBehavior(
		allowSeparation=ON, constraintEnforcementMethod=DEFAULT, 
		pressureOverclosure=HARD)
		
	# define friction propert for acart / ace interaction	
	M.ContactProperty('IntProp-2')
	M.interactionProperties['IntProp-2'].NormalBehavior(
		allowSeparation=ON, constraintEnforcementMethod=DEFAULT, 
		pressureOverclosure=HARD)
	M.interactionProperties['IntProp-2'].TangentialBehavior(
		dependencies=0, directionality=ISOTROPIC, elasticSlipStiffness=None, 
		formulation=PENALTY, fraction=0.005, maximumElasticSlip=FRACTION, 
		pressureDependency=OFF, shearStressLimit=None, slipRateDependency=OFF, 
		table=((0.1, ), ), temperatureDependency=OFF)

	# surface interactions for parameterised model	
	M.SurfaceToSurfaceContactStd(adjustMethod=OVERCLOSED, 
		clearanceRegion=None, createStepName='Initial', datumAxis=None, 
		initialClearance=OMIT, interactionProperty='IntProp-1', master=
		M.rootAssembly.surfaces['parfcart-sur'], name='Int-1', 
		slave=M.rootAssembly.surfaces['paracart-sur'], 
		sliding=FINITE, thickness=ON)
		
	M.SurfaceToSurfaceContactStd(adjustMethod=OVERCLOSED, 
		clearanceRegion=None, createStepName='Initial', datumAxis=None, 
		initialClearance=OMIT, interactionProperty='IntProp-1', master=
		M.rootAssembly.surfaces['parfcart-sur'], name='Int-L1', 
		slave=M.rootAssembly.surfaces['labrum-sur'], 
		sliding=FINITE, thickness=ON)		

	M.SurfaceToSurfaceContactStd(adjustMethod=OVERCLOSED, 
		clearanceRegion=None, createStepName='Initial', datumAxis=None, 
		initialClearance=OMIT, interactionProperty='IntProp-1', master=
		M.rootAssembly.surfaces['parfemur-sur'], name='Int-2', 
		slave=M.rootAssembly.surfaces['paracart-sur'], 
		sliding=FINITE, thickness=ON)	
		
	M.SurfaceToSurfaceContactStd(adjustMethod=OVERCLOSED, 
		clearanceRegion=None, createStepName='Initial', datumAxis=None, 
		initialClearance=OMIT, interactionProperty='IntProp-1', master=
		M.rootAssembly.surfaces['parfemur-sur'], name='Int-L2', 
		slave=M.rootAssembly.surfaces['labrum-sur'], 
		sliding=FINITE, thickness=ON)			

	#  surface tie constraints for parameterised model
	M.Tie(adjust=OFF, master=
		M.rootAssembly.surfaces['parfemur-sur'], name=
		'Constraint-1', positionToleranceMethod=COMPUTED, slave=
		M.rootAssembly.surfaces['parfcart-sur'], thickness=ON, 
		tieRotations=ON)	

	M.Tie(adjust=OFF, master=
		M.rootAssembly.surfaces['parace-sur'], name=
		'Constraint-2', positionToleranceMethod=COMPUTED, slave=
		M.rootAssembly.surfaces['paracart-sur'], thickness=ON, 
		tieRotations=ON)	

	# swap tie for interaction acart / ace
	if ManOp_acartdelam == 1:
		del M.constraints['Constraint-2']
		M.SurfaceToSurfaceContactStd(adjustMethod=OVERCLOSED, 
			clearanceRegion=None, createStepName='Initial', datumAxis=None, 
			initialClearance=OMIT, interactionProperty='IntProp-2', master=
			M.rootAssembly.surfaces['parace-sur'], name='Int-A', 
			slave=M.rootAssembly.surfaces['paracart-sur'], sliding=
			FINITE, thickness=ON, tied=OFF)				

	M.Tie(adjust=ON, master=
		M.rootAssembly.surfaces['parace-sur'], name=
		'Constraint-3', positionToleranceMethod=COMPUTED, slave=
		M.rootAssembly.surfaces['labrum-sur'], thickness=ON
		, tieRotations=ON)
		
	M.Tie(adjust=ON, master=
		M.rootAssembly.surfaces['labrum-sur'], name=
		'Constraint-4', positionToleranceMethod=COMPUTED, slave=
		M.rootAssembly.surfaces['paracart-sur'], thickness=
		ON, tieRotations=ON)

# create some BCs 		
def setBCsPin():

	# initial rotation / translation to start simulation
	M.rootAssembly.rotate(angle=ManOp_FlexAng, axisDirection=(1.0, 0.0, 
		0.0), axisPoint=(0.0, 0.0, 0.0), instanceList=('parfemur-1', 'parfcart-1', ))		

	# set reference points with keys
	RPheadP = M.parts['parfemur'].ReferencePoint(point=(0.0, 0.0, 0.0*HR)).id
	RPaceboneP = M.parts['parace'].ReferencePoint(point=(0.0,Frad+bonethick+ACT, 0.0)).id

	# create steps after initial
	M.StaticStep(maxNumInc=15, name='Move', previous='Initial', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='FlexRot1', previous='Move', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='FlexRot2', previous='FlexRot1', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='FlexRot3', previous='FlexRot2', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='FlexRot4', previous='FlexRot3', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='Rot1', previous='FlexRot4', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='Rot2', previous='Rot1', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='Rot3', previous='Rot2', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='Rot4', previous='Rot3', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='Rot5', previous='Rot4', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='Rot6', previous='Rot5', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.StaticStep(maxNumInc=15, name='Rot7', previous='Rot6', matrixSolver=DIRECT, matrixStorage=UNSYMMETRIC)
	M.steps['Move'].setValues(nlgeom=ON)
	
	# BC1 fixes acetabulum
	M.DisplacementBC(amplitude=UNSET, createStepName='Initial', 
		distributionType=UNIFORM, fieldName='', localCsys=None, name='BC-1', 
		region=Region(referencePoints=(
		M.rootAssembly.instances['parace-1'].referencePoints[RPaceboneP], 
		)), u1=SET, u2=SET, u3=SET, ur1=SET, ur2=SET, ur3=SET)	
		
	# BC2 - translation to establish contact
	# propagated as this acts as pin BC
	M.DisplacementBC(amplitude=UNSET, createStepName='Move', 
		distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None, name=
		'BC-2', region=Region(referencePoints=(
		M.rootAssembly.instances['parfemur-1'].referencePoints[RPheadP], 
		)), u1=ManOp_initialxyz[0] , u2=ManOp_initialxyz[1], u3=ManOp_initialxyz[2], 
		ur1=UNSET, ur2=UNSET, ur3=UNSET)

	# BC3 - constrain rotations
	M.DisplacementBC(amplitude=UNSET, createStepName='Move', 
		distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None, name=
		'BC-3', region=Region(referencePoints=(M.rootAssembly.instances['parfemur-1'].referencePoints[RPheadP], 
		)), u1=UNSET, u2=UNSET, u3=UNSET, ur1=0.0, ur2=0.0, ur3=0.0)
	M.boundaryConditions['BC-3'].deactivate('FlexRot1')
	
	# BC4 - rotate in x axis for flexion
	M.VelocityBC(amplitude=UNSET, createStepName='FlexRot1', 
		distributionType=UNIFORM, fieldName='', localCsys=None, name='BC-4', 
		region=Region(referencePoints=(M.rootAssembly.instances['parfemur-1'].referencePoints[RPheadP], 
		)), v1=UNSET, v2=UNSET, v3=UNSET, vr1=ManOp_irot, vr2=0.0, vr3=0.0)
	M.boundaryConditions['BC-4'].deactivate('Rot1')	

	# BC5 - rotate in y axis for impingement
	if Hip[-1] == 'R':
		M.VelocityBC(amplitude=UNSET, createStepName='Rot1', 
			distributionType=UNIFORM, fieldName='', localCsys=None, name='BC-5', 
			region=Region(referencePoints=(M.rootAssembly.instances['parfemur-1'].referencePoints[RPheadP], 
			)), v1=UNSET, v2=UNSET, v3=UNSET, vr1=0.0, vr2=-ManOp_irot, vr3=0.0)
	if Hip[-1] == 'L':
		M.VelocityBC(amplitude=UNSET, createStepName='Rot1', 
			distributionType=UNIFORM, fieldName='', localCsys=None, name='BC-5', 
			region=Region(referencePoints=(M.rootAssembly.instances['parfemur-1'].referencePoints[RPheadP], 
			)), v1=UNSET, v2=UNSET, v3=UNSET, vr1=0.0, vr2=ManOp_irot, vr3=0.0)	
		
	# solver control options 
	timeAlgoParam=(0, 0, 0, 0, 10.0, 0, 12.0, 20.0, 0, 0, 50.0)    
	M.steps['Rot1'].control.setValues(allowPropagation=OFF, resetDefaultValues=OFF, discontinuous=OFF, timeIncrementation= timeAlgoParam)
	M.steps['Rot1'].solverControl.setValues(allowPropagation=OFF, resetDefaultValues=OFF, maxIterations=10)	
	
	M.rootAssembly.regenerate()	

	
def reqrestart():
	# request restart files for start of steps 
	# step - output - restart requests
	M.steps['Rot5'].Restart(frequency=1, numberIntervals=0, 
		overlay=ON, timeMarks=OFF)	
		
def runjob():

	# import python stuff
	import datetime
	from collections import deque

	# record start time
	with open("record.txt.tmp", mode='a') as file:
		file.write('\n Let\'s go: %s \n' % (datetime.datetime.now()))
		file.close 

	# define and run job	
	myJob = mdb.Job(name=Hip+'job', model='Model-1')
	myJob.submit()
	myJob.waitForCompletion()

	# record finish time
	with open("record.txt.tmp", mode='a') as file:
		file.write('\n Done: %s \n' % (datetime.datetime.now()))
		file.close 			
		
#### Manual steps are defined in protocol ####

# call functions
makesimpleace()
makeparfemur()
setinteractions()
setBCsPin()

#reqrestart()	
#runjob()