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 = 'BuriedTunnel3D'
name = base_name
index = 1
while name in mdb.models or name + '_Job' in mdb.jobs:
    name = '{}_{}'.format(base_name, index)
    index += 1
model = mdb.Model(name=name)
length = 10.0
cover = 10.0
ri = 4.0
ro = 4.3
spring = 0.8
floor = 0.3
radius = 0.3
outer_height = spring + ro
offset = -(cover + outer_height)
soil_width = 40.0
soil_depth = 30.0

def outer_profile(sketch, dy):
    sketch.ArcByCenterEnds(center=(0.0, spring + dy), point1=(ro, spring + dy), point2=(-ro, spring + dy), direction=COUNTERCLOCKWISE)
    points = ((-ro, spring + dy), (-ro, dy), (ro, dy), (ro, spring + dy))
    for i in range(3):
        sketch.Line(point1=points[i], point2=points[i + 1])
lining_sketch = model.ConstrainedSketch(name='LiningProfile', sheetSize=100.0)
outer_profile(lining_sketch, 0.0)
lining_sketch.ArcByCenterEnds(center=(0.0, spring), point1=(ri, spring), point2=(-ri, spring), direction=COUNTERCLOCKWISE)
for side in (-1.0, 1.0):
    lining_sketch.Line(point1=(side * ri, spring), point2=(side * ri, floor + radius))
lining_sketch.ArcByCenterEnds(center=(-ri + radius, floor + radius), point1=(-ri, floor + radius), point2=(-ri + radius, floor), direction=COUNTERCLOCKWISE)
lining_sketch.Line(point1=(-ri + radius, floor), point2=(ri - radius, floor))
lining_sketch.ArcByCenterEnds(center=(ri - radius, floor + radius), point1=(ri - radius, floor), point2=(ri, floor + radius), direction=COUNTERCLOCKWISE)
lining = model.Part(name='Lining', dimensionality=THREE_D, type=DEFORMABLE_BODY)
lining.BaseSolidExtrude(sketch=lining_sketch, depth=length)
soil_sketch = model.ConstrainedSketch(name='SoilProfile', sheetSize=100.0)
soil_sketch.rectangle(point1=(-soil_width / 2.0, -soil_depth), point2=(soil_width / 2.0, 0.0))
outer_profile(soil_sketch, offset)
soil = model.Part(name='Soil', dimensionality=THREE_D, type=DEFORMABLE_BODY)
soil.BaseSolidExtrude(sketch=soil_sketch, depth=length)
properties = (('LiningPlaceholder', 1000000000.0, 0.2, 1.0), ('SoilPlaceholder', 10000000.0, 0.3, 1.0))
for (material_name, modulus, poisson, density) in properties:
    material = model.Material(name=material_name)
    material.Elastic(table=((modulus, poisson),))
    material.Density(table=((density,),))
    model.HomogeneousSolidSection(name=material_name + 'Section', material=material_name)
for (part, material_name, size) in ((lining, 'LiningPlaceholder', 0.15), (soil, 'SoilPlaceholder', 1.5)):
    part.Set(name='AllCells', cells=part.cells[:])
    part.SectionAssignment(region=part.sets['AllCells'], sectionName=material_name + 'Section')
    part.setMeshControls(regions=part.cells[:], elemShape=TET, technique=FREE)
    part.setElementType(regions=(part.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D10, elemLibrary=STANDARD),))
    part.seedPart(size=size, deviationFactor=0.1, minSizeFactor=0.1)
    part.generateMesh()
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)
lining_instance = assembly.Instance(name='Lining-1', part=lining, dependent=ON)
soil_instance = assembly.Instance(name='Soil-1', part=soil, dependent=ON)
assembly.translate(instanceList=('Lining-1',), vector=(0.0, offset, 0.0))
assembly.regenerate()
interface_points = ((0.0, offset + outer_height, length / 2.0), (-ro, offset + spring / 2.0, length / 2.0), (ro, offset + spring / 2.0, length / 2.0), (0.0, offset, length / 2.0))
queries = tuple(((point,) for point in interface_points))
assembly.Surface(name='LiningOuter', side1Faces=lining_instance.faces.findAt(*queries))
assembly.Surface(name='SoilCavity', side1Faces=soil_instance.faces.findAt(*queries))
model.Tie(name='PlaceholderBond', main=assembly.surfaces['SoilCavity'], secondary=assembly.surfaces['LiningOuter'], positionToleranceMethod=COMPUTED, adjust=OFF, tieRotations=ON, thickness=ON)
tolerance = 1e-05
bottom = soil_instance.faces.getByBoundingBox(yMin=-soil_depth - tolerance, yMax=-soil_depth + tolerance)
assembly.Set(name='SoilBottom', faces=bottom)
model.DisplacementBC(name='BottomFixed', createStepName='Initial', region=assembly.sets['SoilBottom'], u1=0.0, u2=0.0, u3=0.0)
for (side, x) in (('Left', -soil_width / 2.0), ('Right', soil_width / 2.0)):
    faces = soil_instance.faces.getByBoundingBox(xMin=x - tolerance, xMax=x + tolerance)
    region = assembly.Set(name=side, faces=faces)
    model.DisplacementBC(name=side + 'Roller', createStepName='Initial', region=region, u1=0.0)
for instance in (lining_instance, soil_instance):
    for (label, z) in (('Front', 0.0), ('Back', length)):
        faces = instance.faces.getByBoundingBox(zMin=z - tolerance, zMax=z + tolerance)
        set_name = instance.name + '_' + label
        region = assembly.Set(name=set_name, faces=faces)
        model.DisplacementBC(name=set_name, createStepName='Initial', region=region, u3=0.0)
model.StaticStep(name='PlaceholderStep', previous='Initial')
model.Gravity(name='PlaceholderGravity', createStepName='PlaceholderStep', comp2=-1.0)
job = mdb.Job(name=name + '_Job', model=name, numCpus=4, numDomains=4)
_sw_record_write(job)
viewport_name = 'Viewport: 1'
if viewport_name not in session.viewports:
    session.Viewport(name=viewport_name)
viewport = session.viewports[viewport_name]
viewport.makeCurrent()
viewport.setValues(displayedObject=assembly)
viewport.assemblyDisplay.setValues(mesh=OFF)
viewport.view.setValues(session.views['Iso'])
viewport.view.fitView()
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')