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 math
base_name = 'Cone_C40_28_Fixed'
model_name = base_name
index = 1
while model_name in mdb.models or 'Job_' + model_name in mdb.jobs:
    model_name = base_name + '_%02d' % index
    index += 1
model = mdb.Model(name=model_name, modelType=STANDARD_EXPLICIT)
length = 200.0
radius_bottom = 20.0
radius_top = 14.0
plate_width = 64.0
n_axial = 50
n_circ = 32
mesh_size = 4.0
imperfection = 0.1
duration = 0.04
material = model.Material(name='AA6061_T6')
material.Density(table=((2.7e-09,),))
material.Elastic(table=((70000.0, 0.33),))
eng_strain = 0.0598
true_stress = 279.98 * (1.0 + eng_strain)
plastic_strain = math.log(1.0 + eng_strain) - true_stress / 70000.0
material.Plastic(table=((200.0, 0.0), (true_stress, plastic_strain)))
tube = model.Part(name='Cone', dimensionality=THREE_D, type=DEFORMABLE_BODY)
for j in range(n_axial + 1):
    y = length * j / float(n_axial)
    radius = radius_bottom + (radius_top - radius_bottom) * y / length
    envelope = math.sin(math.pi * y / length)
    envelope *= math.exp(-((y - 24.0) / 24.0) ** 2)
    for i in range(n_circ):
        theta = 2.0 * math.pi * i / float(n_circ)
        r = radius + imperfection * envelope * math.cos(2.0 * theta)
        tube.Node(coordinates=(r * math.cos(theta), y, r * math.sin(theta)), label=j * n_circ + i + 1)
for j in range(n_axial):
    for i in range(n_circ):
        k = (i + 1) % n_circ
        indices = (j * n_circ + i, (j + 1) * n_circ + i, (j + 1) * n_circ + k, j * n_circ + k)
        tube.Element(nodes=tuple((tube.nodes[q] for q in indices)), elemShape=QUAD4, label=j * n_circ + i + 1)
tube.Set(name='AllElements', elements=tube.elements[:])
for (name, first) in (('BottomRing', 1), ('TopRing', n_axial * n_circ + 1)):
    tube.Set(name=name, nodes=tube.nodes.sequenceFromLabels(labels=tuple(range(first, first + n_circ))))
for j in range(n_axial):
    section_name = 'Wall_%02d' % j
    set_name = 'Band_%02d' % j
    thickness = 1.0 + (j + 0.5) / float(n_axial)
    model.HomogeneousShellSection(name=section_name, material='AA6061_T6', thickness=thickness, integrationRule=SIMPSON, numIntPts=5)
    tube.Set(name=set_name, elements=tube.elements.sequenceFromLabels(labels=tuple(range(j * n_circ + 1, (j + 1) * n_circ + 1))))
    tube.SectionAssignment(region=tube.sets[set_name], sectionName=section_name, offset=0.0, offsetType=MIDDLE_SURFACE, thicknessAssignment=FROM_SECTION)
tube.setElementType(regions=(tube.elements[:],), elemTypes=(mesh.ElemType(elemCode=S4R, elemLibrary=EXPLICIT),))
half_width = plate_width / 2.0
sketch = model.ConstrainedSketch(name='PlateSketch', sheetSize=200.0)
sketch.rectangle(point1=(-half_width, -half_width), point2=(half_width, half_width))
plate = model.Part(name='RigidPlate', dimensionality=THREE_D, type=DISCRETE_RIGID_SURFACE)
plate.BaseShell(sketch=sketch)
faces = plate.faces.getByBoundingBox(xMin=-half_width - 0.01, xMax=half_width + 0.01, yMin=-half_width - 0.01, yMax=half_width + 0.01, zMin=-0.01, zMax=0.01)
plate.Set(name='PlateFaces', faces=faces)
plate.Surface(name='PlateSurface', side1Faces=faces)
plate.setMeshControls(regions=faces, elemShape=QUAD, technique=STRUCTURED)
plate.seedPart(size=mesh_size, deviationFactor=0.1, minSizeFactor=0.1)
plate.setElementType(regions=(faces,), elemTypes=(mesh.ElemType(elemCode=R3D4, elemLibrary=EXPLICIT),))
plate.generateMesh()
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)
assembly.Instance(name='Cone-1', part=tube, dependent=ON)
for (side, height) in (('Bottom', 0.0), ('Top', length)):
    instance_name = side + 'Plate-1'
    assembly.Instance(name=instance_name, part=plate, dependent=ON)
    assembly.rotate(instanceList=(instance_name,), axisPoint=(0.0, 0.0, 0.0), axisDirection=(1.0, 0.0, 0.0), angle=90.0)
    assembly.translate(instanceList=(instance_name,), vector=(0.0, height, 0.0))
    feature = assembly.ReferencePoint(point=(0.0, height, 0.0))
    assembly.Set(name=side + 'RP', referencePoints=(assembly.referencePoints[feature.id],))
    model.RigidBody(name=side + 'RigidBody', refPointRegion=assembly.sets[side + 'RP'], bodyRegion=assembly.instances[instance_name].sets['PlateFaces'])
    arguments = dict(name=side + 'Tie', positionToleranceMethod=SPECIFIED, positionTolerance=0.05, adjust=OFF, tieRotations=ON, thickness=OFF)
    main_region = assembly.instances[instance_name].surfaces['PlateSurface']
    secondary_region = assembly.instances['Cone-1'].sets[side + 'Ring']
    try:
        model.Tie(main=main_region, secondary=secondary_region, **arguments)
    except TypeError:
        model.Tie(master=main_region, slave=secondary_region, **arguments)
assembly.engineeringFeatures.PointMassInertia(name='ImpactorMass', region=assembly.sets['TopRP'], mass=0.07, alpha=0.0, composite=0.0)
model.ExplicitDynamicsStep(name='Impact', previous='Initial', timePeriod=duration, improvedDtMethod=ON)
model.EncastreBC(name='FixedBottom', createStepName='Initial', region=assembly.sets['BottomRP'])
model.DisplacementBC(name='TopGuide', createStepName='Initial', region=assembly.sets['TopRP'], u1=0.0, u2=UNSET, u3=0.0, ur1=0.0, ur2=0.0, ur3=0.0)
model.Velocity(name='InitialImpactVelocity', region=assembly.sets['TopRP'], velocity1=0.0, velocity2=-7000.0, velocity3=0.0, omega=0.0)
contact_property = model.ContactProperty('CrushContactProperty')
contact_property.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON)
contact_property.TangentialBehavior(formulation=PENALTY, directionality=ISOTROPIC, table=((0.15,),), maximumElasticSlip=FRACTION, fraction=0.005)
contact = model.ContactExp(name='GeneralContact', createStepName='Initial')
contact.includedPairs.setValuesInStep(stepName='Initial', useAllstar=ON)
contact.contactPropertyAssignments.appendInStep(stepName='Initial', assignments=((GLOBAL, SELF, 'CrushContactProperty'),))
for name in list(model.fieldOutputRequests.keys()):
    del model.fieldOutputRequests[name]
for name in list(model.historyOutputRequests.keys()):
    del model.historyOutputRequests[name]
model.FieldOutputRequest(name='DeformationFrames', createStepName='Impact', variables=('S', 'PEEQ', 'U', 'V', 'STH'), numIntervals=160)
for (name, set_name, variables) in (('TopMotion', 'TopRP', ('U2', 'V2', 'A2', 'RF2')), ('BottomReaction', 'BottomRP', ('RF2',))):
    model.HistoryOutputRequest(name=name, createStepName='Impact', variables=variables, region=assembly.sets[set_name], timeInterval=duration / 2000.0)
model.HistoryOutputRequest(name='EnergyBalance', createStepName='Impact', variables=('ALLKE', 'ALLIE', 'ALLAE', 'ALLVD', 'ALLFD', 'ETOTAL'), timeInterval=duration / 1000.0)
model.setValues(description='Paper-based C40-28, L200, t1-2. Units mm N s tonne. S4R cone, R3D4 plates, nominal mesh 4. Assumptions: diameters on midsurface, 64 square plates, friction 0.15, nominal 0.10 geometric imperfection, 50 thickness bands, converted bilinear tensile endpoint, duration 0.040 s. No eigenmode import or mass scaling.')
assembly.regenerate()
job_name = 'Job_' + model_name
job = mdb.Job(name=job_name, model=model_name, description='C40-28 explicit axial impact', type=ANALYSIS, numCpus=4, numDomains=4, parallelizationMethodExplicit=DOMAIN, multiprocessingMode=DEFAULT, explicitPrecision=DOUBLE, nodalOutputPrecision=SINGLE, memory=80, memoryUnits=PERCENTAGE)
_sw_record_write(job, consistencyChecking=ON)
if 'Viewport: 1' in session.viewports:
    viewport = session.viewports['Viewport: 1']
    viewport.setValues(displayedObject=assembly)
    viewport.view.fitView()
print('MODEL_CREATED: ' + model_name)
print('INPUT_WRITTEN: ' + job_name + '.inp')
print('TUBE_ELEMENTS: %d' % len(tube.elements))
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')