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
base_name = 'Cantilever3D'
model_name = base_name
index = 1
while model_name in mdb.models or model_name + '_Job' in mdb.jobs:
    model_name = '{}_{}'.format(base_name, index)
    index += 1
job_name = model_name + '_Job'
length = 200.0
width = 20.0
height = 20.0
young_modulus = 210000.0
poisson_ratio = 0.3
tip_force = 1000.0
traction = tip_force / (width * height)
model = mdb.Model(name=model_name)
sketch = model.ConstrainedSketch(name='BeamProfile', sheetSize=500.0)
sketch.rectangle(point1=(0.0, -width / 2.0), point2=(length, width / 2.0))
beam = model.Part(name='Beam', dimensionality=THREE_D, type=DEFORMABLE_BODY)
beam.BaseSolidExtrude(sketch=sketch, depth=height)
del model.sketches['BeamProfile']
material = model.Material(name='Steel')
material.Elastic(table=((young_modulus, poisson_ratio),))
material.Density(table=((7.85e-09,),))
model.HomogeneousSolidSection(name='BeamSection', material='Steel', thickness=None)
beam.Set(name='AllCells', cells=beam.cells[:])
beam.SectionAssignment(region=beam.sets['AllCells'], sectionName='BeamSection')
beam.setMeshControls(regions=beam.cells[:], elemShape=HEX, technique=STRUCTURED)
beam.setElementType(regions=(beam.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D20R, elemLibrary=STANDARD),))
beam.seedPart(size=5.0, deviationFactor=0.1, minSizeFactor=0.1)
beam.generateMesh()
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)
instance = assembly.Instance(name='Beam-1', part=beam, dependent=ON)
assembly.regenerate()
fixed_faces = instance.faces.findAt(((0.0, 0.0, height / 2.0),))
tip_faces = instance.faces.findAt(((length, 0.0, height / 2.0),))
assembly.Set(name='FixedEnd', faces=fixed_faces)
assembly.Set(name='TipFace', faces=tip_faces)
assembly.Surface(name='TipSurface', side1Faces=tip_faces)
tolerance = 1e-05
tip_nodes = instance.nodes.getByBoundingBox(xMin=length - tolerance, xMax=length + tolerance, yMin=-tolerance, yMax=tolerance, zMin=height / 2.0 - tolerance, zMax=height / 2.0 + tolerance)
if len(tip_nodes) != 1:
    raise RuntimeError('Expected one node at the free-end face center.')
assembly.Set(name='TipCenter', nodes=tip_nodes)
model.StaticStep(name='Bending', previous='Initial', nlgeom=OFF, timePeriod=1.0, initialInc=0.1, maxInc=0.1, minInc=1e-06, maxNumInc=100)
model.EncastreBC(name='Clamp', createStepName='Initial', region=assembly.sets['FixedEnd'])
model.SurfaceTraction(name='TipLoad', createStepName='Bending', region=assembly.surfaces['TipSurface'], magnitude=traction, directionVector=((0.0, 0.0, 0.0), (0.0, 0.0, -1.0)), distributionType=UNIFORM, traction=GENERAL, follower=OFF, resultant=OFF)
model.fieldOutputRequests['F-Output-1'].setValues(variables=('S', 'U', 'RF', 'E'))
model.HistoryOutputRequest(name='TipDisplacement', createStepName='Bending', variables=('U3',), region=assembly.sets['TipCenter'])
job = mdb.Job(name=job_name, model=model_name, type=ANALYSIS, numCpus=4, numDomains=4)
_sw_record_write(job, consistencyChecking=ON)
if len(session.viewports):
    viewport = session.viewports[list(session.viewports.keys())[0]]
    viewport.setValues(displayedObject=assembly)
    viewport.view.setValues(session.views['Iso'])
    viewport.view.fitView()
second_moment = width * height ** 3 / 12.0
tip_deflection = tip_force * length ** 3 / (3.0 * young_modulus * second_moment)
nominal_stress = tip_force * length * height / (2.0 * second_moment)
print('Model: {}'.format(model_name))
print('Job: {}'.format(job_name))
print('Input file: {}.inp'.format(job_name))
print('Euler-Bernoulli tip U3: {:.6f}'.format(-tip_deflection))
print('Nominal root bending stress: {:.6f}'.format(nominal_stress))
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')