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 io
modelName = 'BallPlateImpactForce'
jobName = 'BallPlateImpactForceJob'
plateWidth = 24.0
plateThickness = 0.6
ballRadius = 3.0
initialGap = 0.2
impactSpeed = 300000.0
stepTime = 8e-05
density = 7.85e-09
elasticModulus = 210000.0
poissonRatio = 0.3
plateMeshSize = 0.6
ballMeshSize = 0.8
thicknessDivisions = 3
historyIntervals = 1600

def desktopDirectory():
    if os.name == 'nt':
        import ctypes
        buffer = ctypes.create_unicode_buffer(32768)
        result = ctypes.windll.shell32.SHGetFolderPathW(None, 16, None, 0, buffer)
        if result == 0 and buffer.value:
            return buffer.value
    return os.path.join(os.path.expanduser('~'), 'Desktop')
desktop = desktopDirectory()
if not os.path.isdir(desktop):
    os.makedirs(desktop)
if jobName in mdb.jobs:
    del mdb.jobs[jobName]
if modelName in mdb.models:
    del mdb.models[modelName]
model = mdb.Model(name=modelName)
halfWidth = plateWidth / 2.0
sketch = model.ConstrainedSketch(name='PlateProfile', sheetSize=plateWidth * 2.0)
sketch.rectangle(point1=(-halfWidth, -halfWidth), point2=(halfWidth, halfWidth))
plate = model.Part(name='Plate', dimensionality=THREE_D, type=DEFORMABLE_BODY)
plate.BaseSolidExtrude(sketch=sketch, depth=plateThickness)
del model.sketches['PlateProfile']
steel = model.Material(name='PlateSteel')
steel.Density(table=((density,),))
steel.Elastic(table=((elasticModulus, poissonRatio),))
steel.Plastic(table=((300.0, 0.0), (340.0, 0.02), (380.0, 0.05), (420.0, 0.1), (460.0, 0.2)))
steel.DuctileDamageInitiation(table=((0.08, 0.333333, 0.0),))
steel.ductileDamageInitiation.DamageEvolution(type=DISPLACEMENT, softening=LINEAR, table=((0.02,),))
model.HomogeneousSolidSection(name='PlateSection', material='PlateSteel')
plate.SectionAssignment(region=regionToolset.Region(cells=plate.cells[:]), sectionName='PlateSection')
plate.seedPart(size=plateMeshSize, deviationFactor=0.1, minSizeFactor=0.1)
thicknessEdges = []
for edge in plate.edges:
    ids = edge.getVertices()
    if len(ids) == 2:
        p0 = plate.vertices[ids[0]].pointOn[0]
        p1 = plate.vertices[ids[1]].pointOn[0]
        if abs(p0[0] - p1[0]) < 1e-06 and abs(p0[1] - p1[1]) < 1e-06 and (abs(p0[2] - p1[2]) > 1e-06):
            thicknessEdges.append(edge)
plate.seedEdgeByNumber(edges=tuple(thicknessEdges), number=thicknessDivisions, constraint=FIXED)
plate.setMeshControls(regions=plate.cells[:], elemShape=HEX, technique=STRUCTURED)
plate.setElementType(regions=(plate.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D8R, elemLibrary=EXPLICIT, hourglassControl=ENHANCED, elemDeletion=ON, maxDegradation=0.99),))
plate.generateMesh()
plate.Set(name='PlateElements', elements=plate.elements[:])
sketch = model.ConstrainedSketch(name='BallProfile', sheetSize=ballRadius * 6.0)
sketch.ConstructionLine(point1=(0.0, -2.0 * ballRadius), point2=(0.0, 2.0 * ballRadius))
sketch.ArcByCenterEnds(center=(0.0, 0.0), point1=(0.0, ballRadius), point2=(0.0, -ballRadius), direction=CLOCKWISE)
sketch.Line(point1=(0.0, -ballRadius), point2=(0.0, ballRadius))
ball = model.Part(name='Ball', dimensionality=THREE_D, type=DEFORMABLE_BODY)
ball.BaseSolidRevolve(sketch=sketch, angle=360.0, flipRevolveDirection=OFF)
del model.sketches['BallProfile']
material = model.Material(name='BallSteel')
material.Density(table=((density,),))
material.Elastic(table=((elasticModulus, poissonRatio),))
model.HomogeneousSolidSection(name='BallSection', material='BallSteel')
ball.SectionAssignment(region=regionToolset.Region(cells=ball.cells[:]), sectionName='BallSection')
ball.seedPart(size=ballMeshSize, deviationFactor=0.1, minSizeFactor=0.1)
ball.setMeshControls(regions=ball.cells[:], elemShape=TET, technique=FREE)
ball.setElementType(regions=(ball.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D4, elemLibrary=EXPLICIT),))
ball.generateMesh()
volume = 0.0
for element in ball.elements:
    points = [node.coordinates for node in element.getNodes()]
    (a, b, c) = [[points[j][i] - points[0][i] for i in range(3)] for j in range(1, 4)]
    determinant = a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0])
    volume += abs(determinant) / 6.0
ballMass = density * volume
if ballMass <= 0.0:
    raise RuntimeError('Invalid projectile mass.')
rp = ball.ReferencePoint(point=(0.0, 0.0, 0.0))
ball.Set(name='BallRP', referencePoints=(ball.referencePoints[rp.id],))
ball.Set(name='BallElements', elements=ball.elements[:])
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)
plateInstance = assembly.Instance(name='PlateInstance', part=plate, dependent=ON)
ballInstance = assembly.Instance(name='BallInstance', part=ball, dependent=ON)
assembly.translate(instanceList=('BallInstance',), vector=(0.0, 0.0, plateThickness + initialGap + ballRadius))
assembly.regenerate()
model.RigidBody(name='RigidBall', refPointRegion=ballInstance.sets['BallRP'], bodyRegion=ballInstance.sets['BallElements'])
labels = [node.label for node in plateInstance.nodes if abs(abs(node.coordinates[0]) - halfWidth) < 1e-06 or abs(abs(node.coordinates[1]) - halfWidth) < 1e-06]
assembly.Set(name='PlateBoundary', nodes=plateInstance.nodes.sequenceFromLabels(labels=tuple(labels)))
model.EncastreBC(name='ClampPlate', createStepName='Initial', region=assembly.sets['PlateBoundary'])
model.Velocity(name='InitialVelocity', region=ballInstance.sets['BallRP'], velocity1=0.0, velocity2=0.0, velocity3=-impactSpeed, omega=0.0)
model.ExplicitDynamicsStep(name='Impact', previous='Initial', timePeriod=stepTime, nlgeom=ON, linearBulkViscosity=0.06, quadBulkViscosity=1.2)
prop = model.ContactProperty('ImpactProperty')
prop.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON)
prop.TangentialBehavior(formulation=FRICTIONLESS)
contact = model.ContactExp(name='Contact', createStepName='Initial')
contact.includedPairs.setValuesInStep(stepName='Initial', useAllstar=ON)
contact.contactPropertyAssignments.appendInStep(stepName='Initial', assignments=((GLOBAL, SELF, 'ImpactProperty'),))
model.fieldOutputRequests['F-Output-1'].setValues(variables=('S', 'U', 'V', 'STATUS'), numIntervals=60)
model.FieldOutputRequest(name='Damage', createStepName='Impact', region=plateInstance.sets['PlateElements'], variables=('PEEQ', 'SDEG', 'STATUS'), numIntervals=60)
model.historyOutputRequests['H-Output-1'].setValues(variables=('ALLKE', 'ALLIE', 'ALLAE', 'ETOTAL'), numIntervals=200)
model.HistoryOutputRequest(name='BallAcceleration', createStepName='Impact', region=ballInstance.sets['BallRP'], variables=('A3', 'V3', 'U3'), timeInterval=stepTime / float(historyIntervals))
job = mdb.Job(name=jobName, model=modelName, type=ANALYSIS, description='BallMass=%.16g; Fz=BallMass*A3' % ballMass, numCpus=4, numDomains=4, explicitPrecision=DOUBLE, nodalOutputPrecision=FULL, parallelizationMethodExplicit=DOMAIN)
_sw_record_write(job, consistencyChecking=ON)
mdb.saveAs(pathName=os.path.join(desktop, 'BallPlateImpactForce.cae'))
exportCode = "from abaqus import *\nimport os\nimport io\n\nballMass = __MASS__\njobName = __JOB__\noutputPath = __OUTPUT__\nmatches = [\n    odb for key, odb in session.odbs.items()\n    if os.path.splitext(os.path.basename(odb.name))[0].lower()\n    == jobName.lower()\n]\nif len(matches) != 1:\n    raise RuntimeError('Open the completed matching result in CAE first.')\nodb = matches[0]\nstep = odb.steps['Impact']\nif not step.frames or step.frames[-1].frameValue < __END__ * 0.999:\n    raise RuntimeError('The requested analysis has not reached its end time.')\nregions = [\n    region for region in step.historyRegions.values()\n    if 'A3' in region.historyOutputs\n    and 'V3' in region.historyOutputs\n    and 'U3' in region.historyOutputs\n]\nif len(regions) != 1:\n    raise RuntimeError('Expected exactly one projectile acceleration history.')\ndata = [\n    (float(t), ballMass * float(a))\n    for t, a in regions[0].historyOutputs['A3'].data\n]\nif not data:\n    raise RuntimeError('Empty acceleration history.')\nwith io.open(outputPath, 'w', encoding='utf-8-sig', newline='') as stream:\n    stream.write(u'time_s,force_z_N,force_magnitude_N\\n')\n    for t, force in data:\n        stream.write(u'%.16g,%.16g,%.16g\\n' % (t, force, abs(force)))\nxyName = 'BallContactForceZ'\nif xyName in session.xyDataObjects:\n    del session.xyDataObjects[xyName]\nxy = session.XYData(name=xyName, data=tuple(data))\nplotName = 'BallContactForcePlot'\nif plotName in session.xyPlots:\n    del session.xyPlots[plotName]\nplot = session.XYPlot(name=plotName)\nchart = plot.charts[list(plot.charts.keys())[0]]\nchart.setValues(curvesToPlot=(session.Curve(xyData=xy),))\nsession.viewports[session.currentViewportName].setValues(displayedObject=plot)\nprint(outputPath)\n"
exportCode = exportCode.replace('__MASS__', repr(ballMass))
exportCode = exportCode.replace('__JOB__', ascii(jobName) if 'ascii' in dir(__builtins__) else repr(jobName))
exportCode = exportCode.replace('__OUTPUT__', repr(os.path.join(desktop, 'BallForceCurve.csv')))
exportCode = exportCode.replace('__END__', repr(stepTime))
exportPath = os.path.join(desktop, 'ExportBallForce.py')
with io.open(exportPath, 'w', encoding='utf-8') as stream:
    stream.write(exportCode)
print(exportPath)
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')