# -*- coding: mbcs -*-
# Abaqus/CAE Python script
# Two steel plates impact analysis - fast explicit model

from abaqus import *
from abaqusConstants import *
import regionToolset
import mesh

# ------------------------------------------------------------
# 0. Clean model database
# ------------------------------------------------------------
MODEL_NAME = 'TwoPlateImpact'
JOB_NAME = 'TwoPlateImpact_Job'

if MODEL_NAME in mdb.models.keys():
    del mdb.models[MODEL_NAME]

model = mdb.Model(name=MODEL_NAME)

# ------------------------------------------------------------
# 1. Parameters (mm, N, s, tonne)
# ------------------------------------------------------------
L = 100.0
W = 80.0
T = 8.0

gap = 1.0

# Moving plate begins above fixed plate
z_fixed = 0.0
z_move = T + gap

# Impact velocity: 3000 mm/s = 3 m/s
impact_velocity = -3000.0

# Mesh size: relatively coarse for fast calculation
mesh_size = 10.0

# Short explicit event
step_time = 0.003

# ------------------------------------------------------------
# 2. Geometry
# ------------------------------------------------------------
def create_plate(part_name):
    sketch = model.ConstrainedSketch(
        name='__profile__',
        sheetSize=300.0
    )
    sketch.rectangle(
        point1=(-L/2.0, -W/2.0),
        point2=(L/2.0, W/2.0)
    )

    part = model.Part(
        name=part_name,
        dimensionality=THREE_D,
        type=DEFORMABLE_BODY
    )
    part.BaseSolidExtrude(
        sketch=sketch,
        depth=T
    )
    del model.sketches['__profile__']
    return part

fixedPart = create_plate('FixedPlate')
movingPart = create_plate('MovingPlate')

# ------------------------------------------------------------
# 3. Material - structural steel with simple plasticity
# ------------------------------------------------------------
steel = model.Material(name='Steel')
steel.Density(table=((7.85e-9,),))
steel.Elastic(table=((210000.0, 0.30),))

# Mild steel simplified elastic-plastic behavior
steel.Plastic(
    table=(
        (250.0, 0.0),
        (300.0, 0.05),
        (350.0, 0.15)
    )
)

model.HomogeneousSolidSection(
    name='SteelSection',
    material='Steel',
    thickness=None
)

for part in (fixedPart, movingPart):
    region = regionToolset.Region(cells=part.cells)
    part.SectionAssignment(
        region=region,
        sectionName='SteelSection'
    )

# ------------------------------------------------------------
# 4. Assembly
# ------------------------------------------------------------
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)

fixedInst = assembly.Instance(
    name='FixedPlate-1',
    part=fixedPart,
    dependent=ON
)

movingInst = assembly.Instance(
    name='MovingPlate-1',
    part=movingPart,
    dependent=ON
)

# Fixed plate spans Z = 0...T
# Moving plate spans Z = T+gap...2T+gap
assembly.translate(
    instanceList=('MovingPlate-1',),
    vector=(0.0, 0.0, z_move)
)

# ------------------------------------------------------------
# 5. Explicit dynamic step
# ------------------------------------------------------------
model.ExplicitDynamicsStep(
    name='Impact',
    previous='Initial',
    timePeriod=step_time,
    improvedDtMethod=ON
)

# ------------------------------------------------------------
# 6. Contact
# ------------------------------------------------------------
contactProp = model.ContactProperty('SteelContact')

contactProp.NormalBehavior(
    pressureOverclosure=HARD,
    allowSeparation=ON,
    constraintEnforcementMethod=DEFAULT
)

contactProp.TangentialBehavior(
    formulation=PENALTY,
    directionality=ISOTROPIC,
    table=((0.15,),),
    maximumElasticSlip=FRACTION,
    fraction=0.005
)

model.ContactExp(
    name='GeneralContact',
    createStepName='Impact'
)

model.interactions['GeneralContact'].includedPairs.setValuesInStep(
    stepName='Impact',
    useAllstar=ON
)

model.interactions['GeneralContact'].contactPropertyAssignments.appendInStep(
    stepName='Impact',
    assignments=((GLOBAL, SELF, 'SteelContact'),)
)

# ------------------------------------------------------------
# 7. Sets
# ------------------------------------------------------------
assembly.Set(
    name='FIXED_ALL',
    cells=fixedInst.cells
)

assembly.Set(
    name='MOVING_ALL',
    cells=movingInst.cells
)

# ------------------------------------------------------------
# 8. Boundary condition
# ------------------------------------------------------------
model.EncastreBC(
    name='BC_FixedPlate',
    createStepName='Initial',
    region=assembly.sets['FIXED_ALL']
)

# ------------------------------------------------------------
# 9. Initial velocity of moving plate
# ------------------------------------------------------------
model.Velocity(
    name='InitialImpactVelocity',
    region=assembly.sets['MOVING_ALL'],
    velocity1=0.0,
    velocity2=0.0,
    velocity3=impact_velocity,
    omega=0.0
)

# ------------------------------------------------------------
# 10. Mesh
# ------------------------------------------------------------
elemType1 = mesh.ElemType(
    elemCode=C3D8R,
    elemLibrary=EXPLICIT,
    kinematicSplit=AVERAGE_STRAIN,
    hourglassControl=DEFAULT,
    distortionControl=DEFAULT
)

for part in (fixedPart, movingPart):
    part.setElementType(
        regions=(part.cells,),
        elemTypes=(elemType1,)
    )

    part.seedPart(
        size=mesh_size,
        deviationFactor=0.1,
        minSizeFactor=0.1
    )

    # Structured/sweep hex mesh for simple plate geometry
    part.setMeshControls(
        regions=part.cells,
        elemShape=HEX,
        technique=STRUCTURED
    )

    part.generateMesh()

# Regenerate after meshing
assembly.regenerate()

# ------------------------------------------------------------
# 11. Output requests
# ------------------------------------------------------------
# Reduce output frequency to keep ODB small and calculation fast
if 'F-Output-1' in model.fieldOutputRequests.keys():
    model.fieldOutputRequests['F-Output-1'].setValues(
        variables=('S', 'PEEQ', 'LE', 'U', 'V'),
        numIntervals=30
    )

if 'H-Output-1' in model.historyOutputRequests.keys():
    del model.historyOutputRequests['H-Output-1']

model.HistoryOutputRequest(
    name='EnergyOutput',
    createStepName='Impact',
    variables=(
        'ALLKE',
        'ALLIE',
        'ALLAE',
        'ETOTAL'
    ),
    numIntervals=100
)

# ------------------------------------------------------------
# 12. Job
# ------------------------------------------------------------
job = mdb.Job(
    name=JOB_NAME,
    model=MODEL_NAME,
    description='Fast explicit impact analysis of two steel plates',
    type=ANALYSIS,
    explicitPrecision=SINGLE,
    nodalOutputPrecision=SINGLE,
    memory=90,
    memoryUnits=PERCENTAGE,
    multiprocessingMode=DEFAULT,
    numCpus=4,
    numDomains=4
)

# Write exactly one input file; do not submit here
job.writeInput(consistencyChecking=OFF)

print('Model created successfully.')
print('Job name: %s' % JOB_NAME)
print('Input file written: %s.inp' % JOB_NAME)
print('Submit this job later in Abaqus/CAE Job Manager.')