# -*- coding: utf-8 -*-
from abaqus import mdb
from abaqusConstants import *
import mesh, regionToolset, math, os, time
# mm, N, s, tonne, MPa. X forward, Y up, Z lateral.
# Reference-image interpretation: straight travel through a crank angle.
# This is a reduced engineering model, not a validated biological model.
pelican_body_length, pelican_body_width, pelican_body_height = 500.,220.,240.
head_radius, beak_length, beak_width, neck_length = 75.,300.,65.,300.
leg_length, wheel_radius, wheel_width, wheelbase = 350.,300.,40.,900.
frame_radius, crank_radius, pedal_length = 12.5,150.,100.
seat_height, body_x, crank_y, lateral = 780.,390.,440.,90.
duration, drive_start, drive_ramp, crank_speed = 0.6,0.15,0.25,0.8
front_teeth, rear_teeth = 40.,20.
gear_ratio, mu, slip_fraction = front_teeth/rear_teeth,0.7,0.005
body_n, small_n, wheel_n, tire_depth = 12,6,96,40.
name = 'PelicanDrive_' + time.strftime('%Y%m%d_%H%M%S')
folder = os.path.abspath(name)
if not os.path.isdir(folder):
    os.makedirs(folder)
os.chdir(folder)
m = mdb.Model(name=name,modelType=STANDARD_EXPLICIT)
a = m.rootAssembly
a.DatumCsysByDefault(CARTESIAN)
for mat,rho in (('Steel',7.85e-9),('Tissue',1.e-9),
                ('Rubber',1.1e-9),('Keratin',1.2e-9)):
    m.Material(name=mat)
    m.materials[mat].Density(table=((rho,),))
    m.HomogeneousSolidSection(name=mat,material=mat)
m.materials['Steel'].Elastic(table=((210000.,0.3),))
m.materials['Steel'].Plastic(table=((355.,0.),(450.,0.1)))
m.materials['Keratin'].Elastic(table=((1000.,0.35),))
# Tissue initial E=10 MPa, nu=0.45; compressibility avoids locking.
for mat,c,d in (('Tissue',1.724137931,0.06),('Rubber',1.,0.04)):
    m.materials[mat].Hyperelastic(materialType=ISOTROPIC,testData=OFF,
        type=NEO_HOOKE,volumetricResponse=VOLUMETRIC_DATA,table=((c,d),))
J,P,occupied = {},{},{}
sequence = [0]
def uid(prefix):
    sequence[0] += 1
    return prefix+str(sequence[0])
def add(x,y):
    return tuple(x[i]+y[i] for i in range(3))
def sub(x,y):
    return tuple(x[i]-y[i] for i in range(3))
def mul(x,s):
    return tuple(v*s for v in x)
def dot(x,y):
    return sum(x[i]*y[i] for i in range(3))
def unit(x):
    return mul(x,1./math.sqrt(dot(x,x)))
def cross(x,y):
    return (x[1]*y[2]-x[2]*y[1],x[2]*y[0]-x[0]*y[2],
            x[0]*y[1]-x[1]*y[0])
def joint(label,xyz):
    f = a.ReferencePoint(point=xyz)
    J[label] = a.Set(name=label,referencePoints=(a.referencePoints[f.id],))
    P[label] = xyz
    return label
def hinge(label,parent):
    joint(label,P[parent])
    # Independent planar rotation; lateral DOF handled by guide BCs.
    for d in (1,2):
        m.Equation(name=uid('Pin'),terms=((1.,label,d),(-1.,parent,d)))
    return label
def couple(label,inst,labels):
    reg = regionToolset.Region(nodes=inst.nodes.sequenceFromLabels(
        labels=tuple(labels)))
    m.Coupling(name=uid('C'),controlPoint=J[label],surface=reg,
        influenceRadius=WHOLE_SURFACE,couplingType=KINEMATIC,
        u1=ON,u2=ON,u3=ON,ur1=ON,ur2=ON,ur3=ON)
def rigid(label,inst):
    m.RigidBody(name=uid('Rigid'),refPointRegion=J[label],
        bodyRegion=regionToolset.Region(elements=inst.elements[:]))
def build(label,points,cells,material,radius=None,orient=None):
    p = m.Part(name=label,dimensionality=THREE_D,type=DEFORMABLE_BODY)
    for i,xyz in enumerate(points):
        p.Node(coordinates=xyz,label=i+1)
    for i,conn in enumerate(cells):
        p.Element(nodes=tuple(p.nodes[q-1] for q in conn),
            elemShape=HEX8 if radius is None else LINE2,label=i+1)
    reg = p.Set(name='All',elements=p.elements[:])
    p.setElementType(regions=(p.elements[:],),elemTypes=(mesh.ElemType(
        elemCode=C3D8R if radius is None else B31,elemLibrary=EXPLICIT),))
    sec = material
    if radius is not None:
        prof,sec = uid('Profile'),uid('Section')
        m.CircularProfile(name=prof,r=radius)
        m.BeamSection(name=sec,integration=DURING_ANALYSIS,
            poissonRatio=0.3,profile=prof,material=material,
            temperatureVar=LINEAR,consistentMassMatrix=False)
        p.assignBeamSectionOrientation(region=reg,method=N1_COSINES,n1=orient)
    p.SectionAssignment(region=reg,sectionName=sec)
    return a.Instance(name=label,part=p,dependent=ON)
def beam(label,j1,j2,r,n=3):
    v = sub(P[j2],P[j1])
    xyz = [add(P[j1],mul(v,float(i)/n)) for i in range(n+1)]
    orient = (0.,0.,1.) if abs(unit(v)[2])<0.9 else (0.,1.,0.)
    inst = build(label,xyz,[(i+1,i+2) for i in range(n)],'Steel',r,orient)
    couple(j1,inst,[1])
    couple(j2,inst,[n+1])
def ellipsoid(label,center,radii,n=small_n,axis=(1.,0.,0.),mat='Tissue'):
    ex = unit(axis)
    guide = (0.,0.,1.) if abs(ex[2])<0.9 else (0.,1.,0.)
    ey = unit(cross(guide,ex))
    ez = cross(ex,ey)
    points,grid,skin,cells = [],{},[],[]
    for k in range(n+1):
        for j in range(n+1):
            for i in range(n+1):
                x,y,z = [2.*t/n-1. for t in (i,j,k)]
                u = x*math.sqrt(1-y*y/2-z*z/2+y*y*z*z/3)
                v = y*math.sqrt(1-z*z/2-x*x/2+z*z*x*x/3)
                w = z*math.sqrt(1-x*x/2-y*y/2+x*x*y*y/3)
                points.append(add(center,add(mul(ex,radii[0]*u),
                    add(mul(ey,radii[1]*v),mul(ez,radii[2]*w)))))
                grid[i,j,k] = len(points)
                if i in (0,n) or j in (0,n) or k in (0,n):
                    skin.append(len(points))
    offsets = ((0,0,0),(1,0,0),(1,1,0),(0,1,0),
               (0,0,1),(1,0,1),(1,1,1),(0,1,1))
    for k in range(n):
        for j in range(n):
            for i in range(n):
                cells.append(tuple(grid[i+u,j+v,k+w] for u,v,w in offsets))
    return build(label,points,cells,mat),skin
def patch(obj,label,count=9):
    inst,skin = obj
    used = occupied.setdefault(inst.name,set())
    labels = [q for q in skin if q not in used]
    labels.sort(key=lambda q: dot(sub(inst.nodes[q-1].coordinates,P[label]),
                                  sub(inst.nodes[q-1].coordinates,P[label])))
    if len(labels)<count:
        raise RuntimeError('Insufficient coupling patch: '+inst.name)
    used.update(labels[:count])
    couple(label,inst,labels[:count])
def limb(label,j1,j2,r,n=small_n):
    v = sub(P[j2],P[j1])
    obj = ellipsoid(label,mul(add(P[j1],P[j2]),0.5),
        (math.sqrt(dot(v,v))/2,r,r),n,v)
    patch(obj,j1)
    patch(obj,j2)
rear = joint('RearAxle',(0.,wheel_radius,0.))
front = joint('FrontAxle',(wheelbase,wheel_radius,0.))
bb = joint('BottomBracket',(400.,crank_y,0.))
seat = joint('Seat',(body_x,seat_height,0.))
steer = joint('Steering',(780.,740.,0.))
for i,(u,v) in enumerate(((bb,seat),(seat,steer),(steer,bb))):
    beam('Frame%d'%i,u,v,frame_radius)
for side in (-1,1):
    r = joint('RearStay%s'%side,(0.,wheel_radius,side*35.))
    f = joint('Fork%s'%side,(wheelbase,wheel_radius,side*35.))
    for i,(u,v) in enumerate(((r,bb),(r,seat),(f,steer),(rear,r),(front,f))):
        beam('Support%s_%d'%(side,i),u,v,9.,2)
handle = joint('HandleStem',(820.,850.,0.))
beam('Stem',steer,handle,10.,2)
grips = {}
for side in (-1,1):
    grips[side] = joint('Grip%s'%side,(820.,850.,side*120.))
    beam('Handle%s'%side,handle,grips[side],9.,2)
wheels,treads = [],[]
for wi,parent in enumerate((rear,front)):
    hub = hinge('WheelRP%d'%wi,parent)
    wheels.append(hub)
    cx,cy,cz = P[hub]
    points,grid,cells,outer = [],{},[],[]
    for k in range(3):
        for j in range(wheel_n):
            theta = 2*math.pi*j/wheel_n
            for i in range(3):
                r = wheel_radius-tire_depth+i*tire_depth/2
                points.append((cx+r*math.cos(theta),cy+r*math.sin(theta),
                               -wheel_width/2+k*wheel_width/2))
                grid[i,j,k] = len(points)
    for k in range(2):
        for j in range(wheel_n):
            jj = (j+1)%wheel_n
            for i in range(2):
                cells.append((grid[i,j,k],grid[i+1,j,k],grid[i+1,jj,k],
                    grid[i,jj,k],grid[i,j,k+1],grid[i+1,j,k+1],
                    grid[i+1,jj,k+1],grid[i,jj,k+1]))
                if i==1:
                    outer.append(len(cells))
    tire = build('Tire%d'%wi,points,cells,'Rubber')
    treads.append(a.Surface(name='Tread%d'%wi,face4Elements=
        tire.elements.sequenceFromLabels(labels=tuple(outer))))
    rigid(hub,ellipsoid('Hub%d'%wi,P[hub],(22.,22.,25.),4,mat='Steel')[0])
    for j in range(0,wheel_n,8):
        q = joint('Rim%d_%d'%(wi,j),points[grid[0,j,1]-1])
        couple(q,tire,[grid[0,t,k] for t in (j,(j+1)%wheel_n) for k in range(3)])
        beam('Spoke%d_%d'%(wi,j),hub,q,2.5,2)
crank = hinge('CrankRP',bb)
feet = []
for side in (-1,1):
    root = joint('CrankRoot%s'%side,(400.,crank_y,side*lateral))
    theta = math.radians(65.)+(math.pi if side<0 else 0.)
    tip = joint('CrankTip%s'%side,(400.+crank_radius*math.cos(theta),
                crank_y+crank_radius*math.sin(theta),side*lateral))
    beam('Shaft%s'%side,crank,root,12.,1)
    beam('CrankArm%s'%side,root,tip,12.,2)
    pedal = hinge('PedalRP%s'%side,tip)
    rigid(pedal,ellipsoid('Pedal%s'%side,P[pedal],(pedal_length/2,8.,30.),
                          4,mat='Steel')[0])
    foot = ellipsoid('Foot%s'%side,add(P[pedal],(15.,26.,0.)),(65.,18.,32.))
    patch(foot,pedal)
    ankle = joint('Ankle%s'%side,add(P[pedal],(0.,44.,0.)))
    patch(foot,ankle)
    feet.append((side,hinge('AnkleHinge%s'%side,ankle)))
center = (body_x,seat_height+pelican_body_height/2,0.)
body = ellipsoid('Body',center,(pelican_body_length/2,
    pelican_body_height/2,pelican_body_width/2),body_n)
patch(body,seat,16)
rigid(seat,ellipsoid('Saddle',(body_x,seat_height-10.,0.),(85.,10.,65.),
                     4,mat='Rubber')[0])
for side,ankle in feet:
    root = joint('WingRoot%s'%side,add(center,(-20.,30.,side*95.)))
    patch(body,root)
    patch(ellipsoid('Wing%s'%side,add(center,(-55.,30.,side*110.)),
                    (205.,75.,24.),8),root)
    hip = joint('Hip%s'%side,(body_x,seat_height+35.,side*65.))
    patch(body,hip)
    knee = joint('Knee%s'%side,(body_x+100.,seat_height+35.-leg_length*0.40,
                              side*lateral))
    limb('Thigh%s'%side,hinge('HipHinge%s'%side,hip),knee,24.)
    limb('Shank%s'%side,hinge('KneeHinge%s'%side,knee),ankle,19.)
    shoulder = joint('Shoulder%s'%side,add(center,(130.,35.,side*75.)))
    patch(body,shoulder)
    elbow = joint('Elbow%s'%side,(650.,875.,side*105.))
    limb('Arm%s'%side,shoulder,elbow,24.)
    limb('Forearm%s'%side,elbow,grips[side],18.)
base = joint('NeckBase',add(center,(175.,65.,0.)))
patch(body,base)
mid = joint('NeckMid',add(P[base],(30.,neck_length*0.48,0.)))
top = joint('NeckTop',add(P[base],(-20.,neck_length,0.)))
for i,(u,v,r) in enumerate(((base,mid,40.),(mid,top,34.))):
    limb('Neck%d'%i,u,v,r,8)
hc = add(P[top],(25.,45.,0.))
head = ellipsoid('Head',hc,(head_radius,head_radius,head_radius*0.82),8)
patch(head,top)
billroot = joint('BillRoot',add(hc,(head_radius*0.75,-12.,0.)))
patch(head,billroot)
for label,offset,radii,mat in (
    ('Beak',(beak_length/2,-25.,0.),(beak_length/2,18.,beak_width/2),'Keratin'),
    ('Pouch',(100.,-45.,0.),(125.,43.,beak_width*0.48),'Tissue')):
    patch(ellipsoid(label,add(P[billroot],offset),radii,6,mat=mat),billroot)
ground = build('Ground',[(-600.,-50.,-500.),(2200.,-50.,-500.),
    (2200.,0.,-500.),(-600.,0.,-500.),(-600.,-50.,500.),
    (2200.,-50.,500.),(2200.,0.,500.),(-600.,0.,500.)],
    [(1,2,3,4,5,6,7,8)],'Steel')
road = a.Surface(name='Road',face5Elements=ground.elements[:])
m.EncastreBC(name='GroundFixed',createStepName='Initial',
    region=regionToolset.Region(nodes=ground.nodes[:]))
m.ExplicitDynamicsStep(name='Ride',previous='Initial',timePeriod=duration)
m.SmoothStepAmplitude(name='GravityRamp',data=((0.,0.),(0.10,1.),(duration,1.)))
m.SmoothStepAmplitude(name='DriveRamp',data=((0.,0.),(drive_start,0.),
    (drive_start+drive_ramp,1.),(duration,1.)))
prop = m.ContactProperty('Grip')
prop.NormalBehavior(pressureOverclosure=HARD,allowSeparation=ON)
prop.TangentialBehavior(formulation=PENALTY,table=((mu,),),
    maximumElasticSlip=FRACTION,fraction=slip_fraction)
contact = m.ContactExp(name='RollingContact',createStepName='Ride')
contact.includedPairs.setValuesInStep(stepName='Ride',useAllstar=OFF,
    addPairs=tuple((s,road) for s in treads))
contact.contactPropertyAssignments.appendInStep(stepName='Ride',
    assignments=((GLOBAL,SELF,'Grip'),))
for label in J:
    if not label.startswith('Rim'):
        m.DisplacementBC(name='Guide_'+label,createStepName='Initial',
            region=J[label],u3=0.,ur1=0.,ur2=0.)
# Ideal engaged chain: rear relative angle = gear_ratio * crank relative angle.
# No tooth contact, chain slack, derailment or chain stress is represented.
m.Equation(name='ChainTransmission',terms=((1.,wheels[0],6),
    (-1.,rear,6),(-gear_ratio,crank,6),(gear_ratio,bb,6)))
m.VelocityBC(name='CrankMotor',createStepName='Ride',region=J[crank],
    vr3=-crank_speed,amplitude='DriveRamp')
# Wheel translations and front-wheel rotation are NOT prescribed.
m.Gravity(name='Weight',createStepName='Ride',comp2=-9810.,
          amplitude='GravityRamp')
m.fieldOutputRequests['F-Output-1'].setValues(
    variables=('S','LE','U','V','A','RF'),numIntervals=40)
for key in list(m.historyOutputRequests.keys()):
    del m.historyOutputRequests[key]
m.HistoryOutputRequest(name='Energy',createStepName='Ride',
    variables=('ALLIE','ALLKE','ALLAE','ALLVD','ALLWK','ETOTAL'),
    timeInterval=duration/200.)
for label in [crank]+wheels:
    m.HistoryOutputRequest(name='Motion_'+label,createStepName='Ride',
        region=J[label],variables=('U1','U2','UR3','V1','VR3','RM3'),
        timeInterval=duration/200.)
a.regenerate()
job = mdb.Job(name=name,model=name,numCpus=4,numDomains=4,
    explicitPrecision=DOUBLE,memory=80,memoryUnits=PERCENTAGE)
job.writeInput(consistencyChecking=ON)
mdb.saveAs(pathName=os.path.join(folder,name+'.cae'))
print('CAE and INP: '+folder)
print('SOUTHWEITION_JOB_PREPARED')
import os as _sw_os
try:
    print('SOUTHWEITION_CWD ' + _sw_os.getcwd())
    def _sw_job_model(_job):
        _model = getattr(_job, 'model', '')
        try:
            return str(_model.name)
        except Exception:
            return str(_model or '')
    _sw_models = list(mdb.models.keys())
    _sw_latest_model = _sw_models[-1] if _sw_models else ''
    _sw_name = ''
    for _cand in reversed(list(mdb.jobs.keys())):
        if _sw_latest_model and _sw_job_model(mdb.jobs[_cand]) == _sw_latest_model:
            _sw_name = str(_cand)
            break
    if not _sw_name:
        _sw_names = list(mdb.jobs.keys())
        if _sw_names:
            _sw_name = str(_sw_names[-1])
    if _sw_name:
        try:
            mdb.jobs[_sw_name].setValues(numCpus=4, numDomains=4)
        except Exception:
            pass
        mdb.jobs[_sw_name].writeInput(consistencyChecking=OFF)
        print('SOUTHWEITION_JOB ' + str(_sw_name))
        print('SOUTHWEITION_MODEL ' + str(_sw_latest_model))
        print('SOUTHWEITION_INP ' + _sw_os.path.abspath(str(_sw_name) + '.inp'))
    print('SOUTHWEITION_JOB_PREPARED')
except Exception:
    print('SOUTHWEITION_JOB_PREPARED')