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 part, material, section, assembly, step, interaction, load, mesh, job
import builtins as bi
import math
import os
import re
(H, YC, T, NS) = (47.23, 18.0, 0.5, 20)
GROWTH = 1.18
angles = []
for octant in range(8):
    for degree in (0.0, 3.0, 6.0, 10.0, 15.0, 21.0, 28.0, 36.0, 45.0):
        local = degree if octant % 2 == 0 else 45.0 - degree
        angles.append((octant * 45.0 + local) % 360.0)
angles = bi.sorted(set(angles))
N = len(angles)
(pts, quads, kinds, extlin, outlines, outer_tags) = ([], [], [], [], [], [])
(lookup, weights) = ({}, {})

def node(point):
    key = tuple([round(value, 9) for value in point])
    if key not in lookup:
        lookup[key] = len(pts)
        pts.append(key)
    return lookup[key]

def add_quad(points, kind):
    quads.append(tuple([node(point) for point in points]))
    kinds.append(kind)
    return len(quads) - 1

def ring(inner, outer, divisions, kind):
    (rows, cells) = ([], [])
    for j in range(divisions + 1):
        factor = j / float(divisions)
        if kind == 'rock':
            factor = (GROWTH ** j - 1.0) / (GROWTH ** divisions - 1.0)
        rows.append([(p[0] + (q[0] - p[0]) * factor, p[1] + (q[1] - p[1]) * factor) for (p, q) in zip(inner, outer)])
    for j in range(divisions):
        row = []
        for i in range(len(inner)):
            following = (i + 1) % len(inner)
            row.append(add_quad((rows[j][i], rows[j + 1][i], rows[j + 1][following], rows[j][following]), kind))
        cells.append(row)
    return (rows, cells)
for (side, rx, rise, drop) in ((1, 5.88, 6.7, 1.3), (-1, 4.5, 5.8, 1.2)):
    cx = -side * (rx + 1.2)
    (inner, outer, core) = ([], [], [])
    for degree in angles:
        theta = math.radians(degree)
        (cosine, sine) = (math.cos(theta), math.sin(theta))
        local_angle = math.atan2(sine, side * cosine)
        ry = rise if sine >= 0.0 else drop
        inner.append((cx + rx * cosine, YC + ry * sine))
        scale = bi.max(bi.abs(cosine), bi.abs(sine))
        core.append((cx + 0.4 * rx * cosine / scale, YC + 0.4 * bi.min(rise, drop) * sine / scale))
        (x, z) = (cx + (rx + T) * cosine, YC + (ry + T) * sine)
        (top_join, bottom_join) = (22.4, 16.2)
        target = top_join if local_angle >= 0.0 else bottom_join
        if bi.abs(local_angle) <= math.pi / 4.0 + 1e-09:
            x = 0.0
            z = (top_join + bottom_join) / 2.0
            z += (top_join - bottom_join) / 2.0 * local_angle / (math.pi / 4.0)
        elif bi.abs(local_angle) < math.pi / 2.0:
            factor = (math.pi / 2.0 - bi.abs(local_angle)) / (math.pi / 4.0)
            x += side * factor * (bi.abs(cx) - (rx + T) / math.sqrt(2.0))
            sign = 1.0 if local_angle > 0.0 else -1.0
            z += factor * (target - (YC + (ry + T) * sign / math.sqrt(2.0)))
        outer.append((x, z))
    (rows, cells) = ring(inner, outer, 3, 'liner')
    if side == 1:
        exposed = [i for i in range(N) if 45.0 <= angles[i] < 315.0]
    else:
        exposed = [i for i in range(N) if angles[i] >= 225.0]
        exposed += [i for i in range(N) if angles[i] < 135.0]
    for i in exposed:
        extlin.append(cells[-1][i])
        outlines.append(outer[i])
        outer_tags.append((side, angles[i]))
    for j in range(4):
        for i in range(N):
            weights[node(rows[j][i])] = (side, angles[i], j / 3.0)
    ring(core, inner, 4, 'core')
    xx = bi.sorted(set([round(point[0], 9) for point in core]))
    zz = bi.sorted(set([round(point[1], 9) for point in core]))
    for i in range(len(xx) - 1):
        for j in range(len(zz) - 1):
            add_quad(((xx[i], zz[j]), (xx[i + 1], zz[j]), (xx[i + 1], zz[j + 1]), (xx[i], zz[j + 1])), 'core')
M = len(outlines)
corners = [(-60.0, H), (-60.0, 0.0), (60.0, 0.0), (60.0, H)]

def polar(point):
    return math.atan2(point[1] - YC, point[0])
anchors = []
for corner in corners:
    angle = polar(corner)
    distances = [bi.abs(math.atan2(math.sin(polar(point) - angle), math.cos(polar(point) - angle))) for point in outlines]
    anchors.append(distances.index(bi.min(distances)))
far = [None] * M
for j in range(4):
    start = anchors[j]
    divisions = (anchors[(j + 1) % 4] - start) % M
    assert divisions > 0
    for k in range(divisions):
        factor = k / float(divisions)
        (p, q) = (corners[j], corners[(j + 1) % 4])
        far[(start + k) % M] = (p[0] * (1.0 - factor) + q[0] * factor, p[1] * (1.0 - factor) + q[1] * factor)
(rows, rock_cells) = ring(outlines, far, NS, 'rock')
for j in range(NS + 1):
    for i in range(M):
        (side, angle) = outer_tags[i]
        weights[node(rows[j][i])] = (side, angle, 1.0 - j / float(NS))
stations = [0.0, 8.0, 16.0, 22.0, 26.0, 28.0, 28.5, 29.0, 30.0, 32.0, 34.0, 35.0, 36.0, 38.0, 40.0, 41.0, 41.5, 42.0, 44.0, 48.0, 54.0, 62.0, 70.0]
(NP, NQ, NZ) = (len(pts), len(quads), len(stations) - 1)
(base, serial) = ('TwinArch_GradedMesh', 1)
name = base
while name in mdb.models or name + '_Job' in mdb.jobs:
    name = base + '_%02d' % serial
    serial += 1
model = mdb.Model(name=name)
model.setValues(description='Approximate section; graded mesh; independent NV/V domains; immediate support; SI.')
rock = model.Material(name='Rock')
rock.Density(table=((18000.0 / 9.81,),))
rock.Elastic(table=((850000000.0, 0.32),))
rock.MohrCoulombPlasticity(table=((22.0, 0.1),))
rock.mohrCoulombPlasticity.MohrCoulombHardening(table=((180000.0, 0.0),))
concrete = model.Material(name='C30_Elastic')
concrete.Density(table=((25000.0 / 9.81,),))
concrete.Elastic(table=((30000000000.0, 0.2),))
for (sec, mat) in (('RockSection', 'Rock'), ('LiningSection', 'C30_Elastic')):
    model.HomogeneousSolidSection(name=sec, material=mat, thickness=None)
a = model.rootAssembly
a.DatumCsysByDefault(CARTESIAN)
model.GeostaticStep(name='Geostatic', previous='Initial', nlgeom=ON, timePeriod=1.0, initialInc=0.1, minInc=1e-06, maxInc=0.5, maxNumInc=100)
previous = 'Geostatic'
for step_name in ('ExcavateSupport', 'SupportBalance', 'Load_0p4MPa', 'Load_2p4MPa'):
    model.StaticStep(name=step_name, previous=previous, nlgeom=ON, timePeriod=1.0, initialInc=0.01, minInc=1e-08, maxInc=0.05, maxNumInc=1500, matrixStorage=UNSYMMETRIC, stabilizationMethod=DISSIPATED_ENERGY_FRACTION, stabilizationMagnitude=1e-05)
    previous = step_name
prop = model.ContactProperty('RockLining')
prop.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON, constraintEnforcementMethod=DEFAULT)
prop.TangentialBehavior(formulation=PENALTY, directionality=ISOTROPIC, table=((0.3,),), maximumElasticSlip=FRACTION, fraction=0.005)

def element_labels(indices):
    return tuple([k * NQ + i + 1 for k in range(NZ) for i in indices])

def void_gap(index, length):
    (side, angle, weight) = weights.get(index, (0, 0.0, 0.0))
    deviation = bi.abs(angle - 90.0)
    cap = bi.max(0.0, bi.abs(length - 35.0) - 6.0)
    if side != 1 or deviation > 15.0 or cap > 1.0:
        return 0.0
    return weight * math.sqrt(bi.max(0.0, 1.0 - (deviation / 15.0) ** 2)) * math.sqrt(bi.max(0.0, 1.0 - cap ** 2))

def create_part(part_name, indices, section_name, with_void=False):
    used = bi.sorted(set([value for i in indices for value in quads[i]]))
    (node_labels, coordinates, connectivity) = ([], [], [])
    for (k, length) in enumerate(stations):
        layer = [(p[0], p[1] + (void_gap(i, length) if with_void else 0.0)) for (i, p) in enumerate(pts)]
        for index in indices:
            q = quads[index]
            for j in range(4):
                (p, u, v) = (layer[q[j]], layer[q[(j + 1) % 4]], layer[q[(j - 1) % 4]])
                determinant = (u[0] - p[0]) * (v[1] - p[1]) - (u[1] - p[1]) * (v[0] - p[0])
                assert determinant > 1e-10, (part_name, k, index)
        node_labels.extend([k * NP + i + 1 for i in used])
        coordinates.extend([(layer[i][0], -length, layer[i][1]) for i in used])
    for k in range(NZ):
        for index in indices:
            connectivity.append(tuple([k * NP + v + 1 for v in quads[index]] + [(k + 1) * NP + v + 1 for v in quads[index]]))
    p = model.PartFromNodesAndElements(name=part_name, dimensionality=THREE_D, type=DEFORMABLE_BODY, nodes=(tuple(node_labels), tuple(coordinates)), elements=(('C3D8R', element_labels(indices), tuple(connectivity)),))
    p.Set(name='ALL', elements=p.elements[:])
    p.SectionAssignment(region=p.sets['ALL'], sectionName=section_name)
    return p
lining_ids = [i for i in range(NQ) if kinds[i] == 'liner']
fill_ids = [i for i in range(NQ) if kinds[i] != 'rock']
rock_ids = [i for i in range(NQ) if kinds[i] == 'rock']
top_ids = [rock_cells[-1][i] for i in range(M) if bi.abs(far[i][1] - H) < 1e-08 and bi.abs(far[(i + 1) % M][1] - H) < 1e-08]
assert top_ids
model.Gravity(name='Gravity', createStepName='Geostatic', comp3=-9.81)
for (case, offset, cavity) in (('NV', 0.0, False), ('V', 160.0, True)):
    ground = create_part('GROUND_' + case, list(range(NQ)), 'RockSection', cavity)
    lining = create_part('LINING_' + case, lining_ids, 'LiningSection')
    for (tag, indices) in (('EXCAVATION', fill_ids), ('RETAINED', rock_ids)):
        ground.Set(name=tag, elements=ground.elements.sequenceFromLabels(labels=element_labels(indices)))
    ground.Surface(name='ROCK_CONTACT', face6Elements=ground.elements.sequenceFromLabels(labels=element_labels(rock_cells[0])))
    ground.Surface(name='TOP', face4Elements=ground.elements.sequenceFromLabels(labels=element_labels(top_ids)))
    lining.Surface(name='OUTER', face4Elements=lining.elements.sequenceFromLabels(labels=element_labels(extlin)))
    (gn, ln) = ('GROUND_' + case + '-1', 'LINING_' + case + '-1')
    a.Instance(name=gn, part=ground, dependent=ON)
    a.Instance(name=ln, part=lining, dependent=ON)
    if offset:
        a.translate(instanceList=(gn, ln), vector=(offset, 0.0, 0.0))
    (g, l) = (a.instances[gn], a.instances[ln])
    model.ModelChange(name='Install_' + case, createStepName='Geostatic', region=l.sets['ALL'], regionType=ELEMENTS, activeInStep=OFF)
    model.interactions['Install_' + case].setValuesInStep(stepName='ExcavateSupport', activeInStep=ON, includeStrain=OFF)
    model.ModelChange(name='Excavate_' + case, createStepName='ExcavateSupport', region=g.sets['EXCAVATION'], regionType=ELEMENTS, activeInStep=OFF)
    model.GeostaticStress(name='InitialStress_' + case, region=g.sets['ALL'], stressMag1=0.0, vCoord1=H, stressMag2=-18000.0 * H, vCoord2=0.0, lateralCoeff1=0.32 / 0.68, lateralCoeff2=0.32 / 0.68)
    for (prefix, instance) in (('G', g), ('L', l)):
        for (value, tag) in ((0.0, 'FRONT'), (-70.0, 'BACK')):
            nds = instance.nodes.getByBoundingBox(xMin=offset - 61, xMax=offset + 61, yMin=value - 1e-07, yMax=value + 1e-07, zMin=-1, zMax=H + 1)
            assert len(nds) > 0
            bn = case + '_' + prefix + '_' + tag
            a.Set(name=bn, nodes=nds)
            model.DisplacementBC(name=bn, createStepName='Initial', region=a.sets[bn], u2=0.0)
    for (tag, x0, x1, z0, z1, dof) in (('LEFT', offset - 60 - 1e-07, offset - 60 + 1e-07, -1, H + 1, 'u1'), ('RIGHT', offset + 60 - 1e-07, offset + 60 + 1e-07, -1, H + 1, 'u1'), ('BASE', offset - 61, offset + 61, -1e-07, 1e-07, 'u3')):
        nds = g.nodes.getByBoundingBox(xMin=x0, xMax=x1, yMin=-71, yMax=1, zMin=z0, zMax=z1)
        assert len(nds) > 0
        bn = case + '_' + tag
        a.Set(name=bn, nodes=nds)
        model.DisplacementBC(name=bn, createStepName='Initial', region=a.sets[bn], **{dof: 0.0})
    args = dict(name='Contact_' + case, createStepName='ExcavateSupport', sliding=FINITE, interactionProperty='RockLining', adjustMethod=NONE, initialClearance=OMIT)
    try:
        model.SurfaceToSurfaceContactStd(main=l.surfaces['OUTER'], secondary=g.surfaces['ROCK_CONTACT'], **args)
    except TypeError:
        model.SurfaceToSurfaceContactStd(master=l.surfaces['OUTER'], slave=g.surfaces['ROCK_CONTACT'], **args)
    model.Pressure(name='Pressure_' + case, createStepName='Load_0p4MPa', region=g.surfaces['TOP'], magnitude=400000.0)
    model.loads['Pressure_' + case].setValuesInStep(stepName='Load_2p4MPa', magnitude=2400000.0)
    for (side, cx, z) in (('L', -7.08, 24.7), ('R', 5.7, 23.8)):
        distances = [(bi.sum([(nd.coordinates[0] - offset - cx) ** 2, (nd.coordinates[1] + 35.0) ** 2, (nd.coordinates[2] - z) ** 2]), int(nd.label)) for nd in l.nodes]
        label = bi.min(distances)[1]
        tag = case + '_CROWN_' + side
        a.Set(name=tag, nodes=l.nodes.sequenceFromLabels(labels=(label,)))
        model.HistoryOutputRequest(name=tag, createStepName='ExcavateSupport', variables=('U1', 'U2', 'U3'), region=a.sets[tag], frequency=1)
for key in list(model.fieldOutputRequests.keys()):
    del model.fieldOutputRequests[key]
model.FieldOutputRequest(name='Fields', createStepName='Geostatic', variables=('S', 'LE', 'U', 'RF', 'PE', 'PEEQ', 'CSTRESS', 'CDISP'), frequency=1)
model.HistoryOutputRequest(name='Energy', createStepName='Geostatic', variables=('ALLIE', 'ALLSE', 'ALLAE', 'ALLSD', 'ALLWK'), frequency=1)
a.regenerate()
jobname = name + '_Job'
analysis_job = mdb.Job(name=jobname, model=name, type=ANALYSIS, numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE)
_sw_record_write(analysis_job, consistencyChecking=ON)
inppath = os.path.join(os.getcwd(), jobname + '.inp')
with open(inppath, 'r') as stream:
    inptext = stream.read()
types = re.findall('\\*Element\\s*,[^\\n]*type\\s*=\\s*([A-Za-z0-9]+)', inptext, re.I)
assert types and bi.all([value.upper() == 'C3D8R' for value in types])
mdb.saveAs(pathName=os.path.join(os.getcwd(), name + '.cae'))
if len(session.viewports):
    vp = session.viewports[list(session.viewports.keys())[0]]
    vp.setValues(displayedObject=model.parts['GROUND_NV'])
    vp.partDisplay.setValues(mesh=ON)
    vp.view.setValues(cameraPosition=(90.0, 220.0, 100.0), cameraTarget=(0.0, -35.0, 23.0), cameraUpVector=(0.0, 0.0, 1.0))
    vp.view.fitView()
print('MODEL:', name)
print('TOTAL_ELEMENTS:', 2 * (NQ + len(lining_ids)) * NZ)
print('INP:', inppath)
print('NV: no cavity; V: left crown cavity, length 14, angle 30, depth 1')
print('U3_BASELINE: end of SupportBalance; JOB_NOT_SUBMITTED')
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')