from abaqus import mdb
from abaqusConstants import OFF
import os as _sw_os
_sw_before_jobs = {str(n): mdb.jobs[n] for n in mdb.jobs.keys()}
_sw_written_jobs = []
def _sw_record_write(_job, *args, **kwargs):
    _answer = _job.writeInput(*args, **kwargs)
    _name = str(getattr(_job, 'name', ''))
    if not _name or _name not in mdb.jobs.keys():
        _names = [str(n) for n in mdb.jobs.keys() if mdb.jobs[n] is _job]
        if len(_names) != 1:
            raise RuntimeError('Cannot identify the job which wrote the INP')
        _name = _names[0]
    _sw_written_jobs.append((_name, _sw_os.path.abspath(_name + '.inp')))
    return _answer

from abaqus import *
from abaqusConstants import *
import mesh
import regionToolset
import os
import time
name = 'CutCompare_' + time.strftime('%Y%m%d_%H%M%S')
while name in mdb.models or name in mdb.jobs:
    name += 'A'
(E, NU, SY) = (30.0, 0.3, 0.6)
(EP0, UF, RHO) = (0.25, 0.4, 1e-09)
(LX, LY, H) = (24.0, 6.0, 12.0)
(TRAVEL, PERIOD) = (13.0, 0.06)
vumat_path = os.path.abspath(name + '.for')
source = "\n      subroutine vumat(\n     1 nblock,ndir,nshr,nstatev,nfieldv,nprops,jInfoArray,\n     2 stepTime,totalTime,dtArray,cmname,coordMp,charLength,\n     3 props,density,strainInc,relSpinInc,\n     4 tempOld,stretchOld,defgradOld,fieldOld,\n     5 stressOld,stateOld,enerInternOld,enerInelasOld,\n     6 tempNew,stretchNew,defgradNew,fieldNew,\n     7 stressNew,stateNew,enerInternNew,enerInelasNew)\n      include 'vaba_param.inc'\n      dimension jInfoArray(*),dtArray(*),props(nprops),\n     1 density(nblock),coordMp(nblock,*),charLength(nblock),\n     2 strainInc(nblock,ndir+nshr),relSpinInc(nblock,nshr),\n     3 tempOld(nblock),stretchOld(nblock,ndir+nshr),\n     4 defgradOld(nblock,ndir+2*nshr),fieldOld(nblock,nfieldv),\n     5 stressOld(nblock,ndir+nshr),stateOld(nblock,nstatev),\n     6 enerInternOld(nblock),enerInelasOld(nblock),\n     7 tempNew(nblock),stretchNew(nblock,ndir+nshr),\n     8 defgradNew(nblock,ndir+2*nshr),fieldNew(nblock,nfieldv),\n     9 stressNew(nblock,ndir+nshr),stateNew(nblock,nstatev),\n     1 enerInternNew(nblock),enerInelasNew(nblock)\n      character*80 cmname\n      integer k,i\n      double precision g,bulk,sy,ep0,uf,tr,pm,q,dp,fac\n      double precision ep,up,d,du,work,eff(6),dev(6),oldD\n      g=props(1)/(2.d0*(1.d0+props(2)))\n      bulk=props(1)/(3.d0*(1.d0-2.d0*props(2)))\n      sy=props(3)\n      ep0=props(4)\n      uf=props(5)\n      do k=1,nblock\n        do i=1,nstatev\n          stateNew(k,i)=stateOld(k,i)\n        enddo\n        do i=1,6\n          stressNew(k,i)=0.d0\n        enddo\n        enerInternNew(k)=enerInternOld(k)\n        enerInelasNew(k)=enerInelasOld(k)\n        if (totalTime.eq.0.d0) stateNew(k,3)=1.d0\n        if (stateNew(k,3).gt.0.5d0) then\n          ep=stateOld(k,1)\n          oldD=stateOld(k,2)\n          up=stateOld(k,4)\n          tr=sum(strainInc(k,1:3))\n          do i=1,6\n            eff(i)=stateOld(k,4+i)+2.d0*g*strainInc(k,i)\n          enddo\n          do i=1,3\n            eff(i)=eff(i)+(bulk-2.d0*g/3.d0)*tr\n          enddo\n          pm=sum(eff(1:3))/3.d0\n          dev=eff\n          dev(1:3)=eff(1:3)-pm\n          q=sqrt(1.5d0*(sum(dev(1:3)**2)\n     1      +2.d0*sum(dev(4:6)**2)))\n          dp=max(0.d0,(q-sy)/(3.d0*g))\n          if (totalTime.eq.0.d0) dp=0.d0\n          fac=1.d0\n          if (q.gt.1.d-14) fac=max(0.d0,1.d0-3.d0*g*dp/q)\n          eff=fac*dev\n          eff(1:3)=eff(1:3)+pm\n          du=charLength(k)*max(0.d0,\n     1       ep+dp-max(ep,ep0))\n          ep=ep+dp\n          up=up+du\n          d=max(oldD,min(1.d0,up/uf))\n          stateNew(k,1)=ep\n          stateNew(k,2)=d\n          stateNew(k,4)=up\n          if (d.ge.1.d0) stateNew(k,3)=0.d0\n          do i=1,6\n            stateNew(k,4+i)=eff(i)\n            stressNew(k,i)=(1.d0-d)*eff(i)\n          enddo\n          work=0.d0\n          do i=1,6\n            fac=0.5d0\n            if (i.gt.3) fac=1.d0\n            work=work+fac*(stressOld(k,i)+stressNew(k,i))\n     1           *strainInc(k,i)\n          enddo\n          enerInternNew(k)=enerInternOld(k)+work/density(k)\n          enerInelasNew(k)=enerInelasOld(k)\n     1      +(1.d0-0.5d0*(oldD+d))*sy*dp/density(k)\n        endif\n      enddo\n      return\n      end\n"
with open(vumat_path, 'w') as handle:
    handle.write(source)
m = mdb.Model(name=name)
umat = m.Material(name='CustomJ2Damage')
umat.Density(table=((RHO,),))
umat.Depvar(n=10, deleteVar=3)
umat.UserMaterial(mechanicalConstants=(E, NU, SY, EP0, UF))
builtin = m.Material(name='BuiltinJ2Damage')
builtin.Density(table=((RHO,),))
builtin.Elastic(table=((E, NU),))
builtin.Plastic(table=((SY, 0.0),))
builtin.DuctileDamageInitiation(table=((EP0, 0.0, 0.0),))
builtin.ductileDamageInitiation.DamageEvolution(type=DISPLACEMENT, softening=LINEAR, table=((UF,),))
toolmat = m.Material(name='RigidToolMaterial')
toolmat.Density(table=((7.8e-09,),))
toolmat.Elastic(table=((210000.0, 0.3),))
for (sec, material) in (('CustomSection', 'CustomJ2Damage'), ('BuiltinSection', 'BuiltinJ2Damage'), ('ToolSection', 'RigidToolMaterial')):
    m.HomogeneousSolidSection(name=sec, material=material)
a = m.rootAssembly
a.DatumCsysByDefault(CARTESIAN)
m.ExplicitDynamicsStep(name='Cut', previous='Initial', timePeriod=PERIOD)
m.SmoothStepAmplitude(name='Travel', timeSpan=STEP, data=((0.0, 0.0), (0.035, 1.0), (0.04, 1.0), (0.055, 0.0), (0.06, 0.0)))
for (label, offset, section) in (('VUMAT', -20.0, 'CustomSection'), ('BUILTIN', 20.0, 'BuiltinSection')):
    sk = m.ConstrainedSketch(name='Profile', sheetSize=100.0)
    sk.rectangle(point1=(-LX / 2.0, 0.0), point2=(LX / 2.0, LY))
    p = m.Part(name='Block_' + label, dimensionality=THREE_D, type=DEFORMABLE_BODY)
    p.BaseSolidExtrude(sketch=sk, depth=H)
    del m.sketches['Profile']
    p.Set(name='ALL', cells=p.cells[:])
    p.SectionAssignment(region=p.sets['ALL'], sectionName=section)
    p.setMeshControls(regions=p.cells[:], elemShape=HEX, technique=STRUCTURED)
    p.setElementType(regions=(p.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D8R, elemLibrary=EXPLICIT, hourglassControl=ENHANCED, elemDeletion=ON, maxDegradation=1.0),))
    p.seedPart(size=1.0, deviationFactor=0.1, minSizeFactor=0.1)
    p.generateMesh()
    bi = a.Instance(name='Block_' + label, part=p, dependent=ON)
    a.translate(instanceList=(bi.name,), vector=(offset, 0.0, 0.0))
    a.regenerate()
    a.Set(name='MAT_' + label, elements=bi.elements[:])
    for (side, x) in (('L', offset - LX / 2.0), ('R', offset + LX / 2.0)):
        fs = bi.faces.getByBoundingBox(xMin=x - 0.001, xMax=x + 0.001, yMin=-0.001, yMax=LY + 0.001, zMin=-0.001, zMax=H + 0.001)
        if len(fs) != 1:
            raise RuntimeError('Clamp face selection failed.')
        a.Set(name=label + side, faces=fs)
        m.EncastreBC(name=label + side, createStepName='Initial', region=a.sets[label + side])
    sk = m.ConstrainedSketch(name='Profile', sheetSize=60.0)
    pts = ((-0.3, 0.0), (0.3, 0.0), (2.0, 10.0), (-2.0, 10.0))
    for i in range(4):
        sk.Line(point1=pts[i], point2=pts[(i + 1) % 4])
    p = m.Part(name='Knife_' + label, dimensionality=THREE_D, type=DEFORMABLE_BODY)
    p.BaseSolidExtrude(sketch=sk, depth=LY + 4.0)
    del m.sketches['Profile']
    p.Set(name='ALL', cells=p.cells[:])
    p.SectionAssignment(region=p.sets['ALL'], sectionName='ToolSection')
    p.setMeshControls(regions=p.cells[:], elemShape=HEX, technique=SWEEP)
    p.setElementType(regions=(p.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D8R, elemLibrary=EXPLICIT),))
    p.seedPart(size=2.0, deviationFactor=0.1, minSizeFactor=0.1)
    p.generateMesh()
    ki = a.Instance(name='Knife_' + label, part=p, dependent=ON)
    a.rotate(instanceList=(ki.name,), axisPoint=(0.0, 0.0, 0.0), axisDirection=(1.0, 0.0, 0.0), angle=90.0)
    a.translate(instanceList=(ki.name,), vector=(offset, LY + 2.0, H + 0.1))
    a.regenerate()
    rp = a.ReferencePoint(point=(offset, LY / 2.0, H + 10.1))
    a.Set(name='RP_' + label, referencePoints=(a.referencePoints[rp.id],))
    m.RigidBody(name='Rigid_' + label, refPointRegion=a.sets['RP_' + label], bodyRegion=regionToolset.Region(cells=ki.cells[:]))
    bc = m.DisplacementBC(name='Motion_' + label, createStepName='Initial', region=a.sets['RP_' + label], u1=0.0, u2=0.0, u3=0.0, ur1=0.0, ur2=0.0, ur3=0.0)
    bc.setValuesInStep(stepName='Cut', u3=-TRAVEL, amplitude='Travel')
    m.HistoryOutputRequest(name='Force_' + label, createStepName='Cut', variables=('U3', 'RF3'), region=a.sets['RP_' + label], numIntervals=300)
prop = m.ContactProperty('Friction')
prop.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON)
prop.TangentialBehavior(formulation=PENALTY, table=((0.1,),), maximumElasticSlip=FRACTION, fraction=0.005)
ct = m.ContactExp(name='Contact', createStepName='Initial')
ct.includedPairs.setValuesInStep(stepName='Initial', useAllstar=ON)
ct.contactPropertyAssignments.appendInStep(stepName='Initial', assignments=((GLOBAL, SELF, 'Friction'),))
for key in list(m.fieldOutputRequests.keys()):
    del m.fieldOutputRequests[key]
m.FieldOutputRequest(name='Common', createStepName='Cut', variables=('S', 'U', 'V', 'STATUS'), numIntervals=80)
for (label, variables) in (('VUMAT', ('SDV',)), ('BUILTIN', ('PEEQ', 'SDEG'))):
    m.FieldOutputRequest(name='State_' + label, createStepName='Cut', variables=variables, region=a.sets['MAT_' + label], numIntervals=80)
m.HistoryOutputRequest(name='Energy', createStepName='Cut', variables=('ALLIE', 'ALLKE', 'ALLAE', 'ETOTAL'), numIntervals=200)
for part in m.parts.values():
    if len(part.elements) == 0:
        raise RuntimeError('Empty mesh: ' + part.name)
job = mdb.Job(name=name, model=name, type=ANALYSIS, userSubroutine=vumat_path, numCpus=4, numDomains=4, explicitPrecision=DOUBLE, nodalOutputPrecision=FULL)
_sw_record_write(job, consistencyChecking=ON)
if len(session.viewports):
    vp = session.viewports[session.currentViewportName]
    vp.setValues(displayedObject=a)
    vp.view.setValues(session.views['Iso'])
    vp.view.fitView()
print('VUMAT specimen: x = -20; built-in specimen: x = +20')
print('Job: ' + name)
print('INP: ' + os.path.abspath(name + '.inp'))
print('VUMAT: ' + vumat_path)
print('Demo parameters only; calibration and validation required.')
print('SOUTHWEITION_JOB_PREPARED')

if not _sw_written_jobs:
    _sw_candidates = [n for n in mdb.jobs.keys() if n not in _sw_before_jobs or mdb.jobs[n] is not _sw_before_jobs[n]]
    if len(_sw_candidates) != 1:
        raise RuntimeError('Cannot identify this script job: call writeInput() on exactly the intended job')
    _sw_record_write(mdb.jobs[_sw_candidates[0]], consistencyChecking=OFF)
_sw_name, _sw_inp = _sw_written_jobs[-1]
if not _sw_os.path.isfile(_sw_inp):
    raise RuntimeError('This script did not produce its INP: ' + _sw_inp)
print('SOUTHWEITION_CWD ' + _sw_os.path.dirname(_sw_inp))
print('SOUTHWEITION_JOB ' + _sw_name)
print('SOUTHWEITION_MODEL ' + str(getattr(mdb.jobs[_sw_name], 'model', '')))
print('SOUTHWEITION_INP ' + _sw_inp)
print('SOUTHWEITION_JOB_PREPARED')