# -*- coding: utf-8 -*-
from abaqus import mdb
from abaqusConstants import *
import mesh
import math
import os
import tempfile
import traceback
try:
    import builtins as py
except ImportError:
    import __builtin__ as py

# ============================================================
# 仅建立橡胶：主簧、两个止挡和连续外环共节点连接。
# mm、N、MPa；外环厚度1.2、轴向厚度30为待核实参数。
# 空腔R3/R10过渡为样条近似，不宣称严格复刻全部论文尺寸。
# 只写入CAE/INP，不执行求解，不读取ODB，不使用网络请求。
# WinError 10054可能来自外部调度连接；本脚本不能修复该连接。
# OUT_ROOT可改成已有的本地磁盘目录，例如 r'E:\AbaqusRuns'。
# ============================================================

OUT_ROOT = os.path.join(tempfile.gettempdir(), 'AbaqusRubberRuns')
if not os.path.isdir(OUT_ROOT):
    os.makedirs(OUT_ROOT)
OUTDIR = tempfile.mkdtemp(prefix='Connected_', dir=OUT_ROOT)
BASE = 'Rubber_' + os.path.basename(OUTDIR).split('_')[-1]
os.chdir(OUTDIR)

RO, RI = 34.495, 16.5
RING_T = 1.2
RC = RO-RING_T
DEPTH, NZ = 30.0, 6
H, HF = 1.8, 0.6
FMAX = 6000.0
PERIOD = 0.1
LOG = os.path.join(OUTDIR, 'build_progress.txt')

def report(message):
    with py.open(LOG, 'a') as stream:
        stream.write(message+'\n')
    print(message)

def polar(angle):
    t = math.radians(angle)
    return (RC*math.cos(t), RC*math.sin(t))

def build():
    report('1/7 Starting fresh model: '+BASE)
    report('Output directory: '+OUTDIR)
    if not 0.0 < RING_T < 2.0:
        raise ValueError('This traced profile requires 0 < RING_T < 2 mm.')
    m = mdb.Model(name=BASE)
    sk = m.ConstrainedSketch(name='RubberProfile', sheetSize=160.0)
    for radius in (RO, RI):
        sk.CircleByCenterPerimeter(center=(0.0, 0.0),
                                   point1=(radius, 0.0))
    a, b, c, d = [polar(t) for t in (125.0, 110.0, 70.0, 55.0)]
    cx = 8.0
    cy = 35.5-math.sqrt(9.5**2-cx**2)
    u, v = (-cx, cy), (cx, cy)
    if cy >= b[1]:
        raise ValueError('Cavity transition incompatible with ring thickness.')

    # 两个完全位于外圆内的闭合空腔，保留连续橡胶外环。
    for sign in (1.0, -1.0):
        def tr(q):
            return (q[0], sign*q[1])
        cw = CLOCKWISE if sign > 0.0 else COUNTERCLOCKWISE
        ccw = COUNTERCLOCKWISE if sign > 0.0 else CLOCKWISE
        for start, end in ((a, b), (c, d)):
            sk.ArcByCenterEnds(center=(0.0, 0.0),
                point1=tr(start), point2=tr(end), direction=cw)
        lm = ((b[0]+u[0])/2.0, (b[1]+u[1])/2.0)
        rm = (-lm[0], lm[1])
        lower = [d, (17.5, 24.6), (15.4, 22.0), (10.5, 20.65),
                 (5.0, 20.5), (0.0, 20.5), (-5.0, 20.5),
                 (-10.5, 20.65), (-15.4, 22.0), (-17.5, 24.6), a]
        for points in ([b, lm, u], [v, rm, c], lower):
            sk.Spline(points=py.tuple([tr(q) for q in points]))
        sk.ArcByCenterEnds(center=tr((0.0, 35.5)),
            point1=tr(u), point2=tr(v), direction=ccw)

    src = m.Part(name='SectionTemplate', dimensionality=THREE_D,
                 type=DEFORMABLE_BODY)
    src.BaseShell(sketch=sk)
    if py.len(src.faces) != 1:
        raise RuntimeError('Expected one connected rubber cross-section.')
    report('2/7 Connected profile created; generating quadrilateral mesh.')
    src.seedPart(size=H, deviationFactor=0.04, minSizeFactor=0.1)
    for edge in src.edges:
        x, y, z = edge.pointOn[0]
        if math.hypot(x, y) > RO-0.05 or py.abs(y) > 19.0:
            src.seedEdgeBySize(edges=(edge,), size=HF,
                deviationFactor=0.03, minSizeFactor=0.1, constraint=FINER)
    src.setMeshControls(regions=src.faces[:], elemShape=QUAD,
                         technique=FREE, algorithm=MEDIAL_AXIS)
    src.setElementType(regions=(src.faces[:],),
        elemTypes=(mesh.ElemType(elemCode=S4R, elemLibrary=STANDARD),))
    src.generateMesh()
    if not py.len(src.elements):
        raise RuntimeError('Section mesh is empty.')

    # 检查四边形与网格连通性，不允许独立橡胶岛。
    parent = {}
    def root(n):
        while parent[n] != n:
            parent[n] = parent[parent[n]]
            n = parent[n]
        return n
    for element in src.elements:
        nodes = py.list(element.getNodes())
        if py.len(nodes) != 4:
            raise RuntimeError('Non-quadrilateral element detected.')
        for node in nodes:
            parent.setdefault(node.label, node.label)
        first = root(nodes[0].label)
        for node in nodes[1:]:
            parent[root(node.label)] = first
    if py.len(py.set([root(n) for n in py.list(parent.keys())])) != 1:
        raise RuntimeError('Rubber mesh has disconnected components.')

    report('3/7 Extruding section mesh into HEX8 elements.')
    p = m.Part(name='Rubber3D', dimensionality=THREE_D,
               type=DEFORMABLE_BODY)
    nodes = py.list(src.nodes)
    nn = py.len(nodes)
    lookup = {}
    for j, node in py.enumerate(nodes):
        lookup[node.label] = j
    for k in py.range(NZ+1):
        z = DEPTH*py.float(k)/NZ-DEPTH/2.0
        for j, node in py.enumerate(nodes):
            p.Node(coordinates=(node.coordinates[0], node.coordinates[1], z),
                   label=k*nn+j+1)
    label = 0
    for k in py.range(NZ):
        for element in src.elements:
            q = py.list(element.getNodes())
            area = 0.0
            for j in py.range(4):
                x1, y1 = q[j].coordinates[:2]
                x2, y2 = q[(j+1) % 4].coordinates[:2]
                area += x1*y2-x2*y1
            if py.abs(area) < 1.e-10:
                raise RuntimeError('Degenerate section element.')
            if area < 0.0:
                q.reverse()
            ids = [lookup[node.label] for node in q]
            conn = [p.nodes[k*nn+j] for j in ids]
            conn.extend([p.nodes[(k+1)*nn+j] for j in ids])
            label += 1
            p.Element(nodes=py.tuple(conn), elemShape=HEX8, label=label)
    del m.parts['SectionTemplate']

    report('4/7 Assigning NR45 and explicit material orientation.')
    rubber = p.Set(name='RUBBER', elements=p.elements[:])
    m.Material(name='NR45').Hyperelastic(
        materialType=ISOTROPIC, testData=OFF, type=POLYNOMIAL, n=2,
        volumetricResponse=VOLUMETRIC_DATA,
        table=((-0.438, 1.537, 0.0, 0.132, 0.0, 0.0, 0.0),))
    m.HomogeneousSolidSection(name='RubberSection', material='NR45')
    p.SectionAssignment(region=rubber, sectionName='RubberSection')
    p.setElementType(regions=(p.elements[:],),
        elemTypes=(mesh.ElemType(elemCode=C3D8H, elemLibrary=STANDARD),))
    datum = p.DatumCsysByThreePoints(name='MaterialAxes',
        coordSysType=CARTESIAN, origin=(0.0, 0.0, 0.0),
        point1=(1.0, 0.0, 0.0), point2=(0.0, 1.0, 0.0))
    p.MaterialOrientation(region=rubber, orientationType=SYSTEM,
        localCsys=p.datums[datum.id], axis=AXIS_3,
        additionalRotationType=ROTATION_NONE, stackDirection=STACK_3)
    for name, radius in (('INNER_BOND', RI), ('OUTER_BOND', RO)):
        labels = []
        for node in p.nodes:
            r = math.hypot(node.coordinates[0], node.coordinates[1])
            if py.abs(r-radius) < 1.e-4:
                labels.append(node.label)
        if not labels:
            raise RuntimeError('Missing interface: '+name)
        p.Set(name=name, nodes=p.nodes.sequenceFromLabels(
              labels=py.tuple(labels)))

    report('5/7 Defining interface constraints and cavity contact.')
    assembly = m.rootAssembly
    assembly.DatumCsysByDefault(CARTESIAN)
    inst = assembly.Instance(name='RUBBER', part=p, dependent=ON)
    ref = assembly.ReferencePoint(point=(0.0, 0.0, 0.0))
    assembly.Set(name='DRIVER',
        referencePoints=(assembly.referencePoints[ref.id],))
    m.Coupling(name='InnerDriver', controlPoint=assembly.sets['DRIVER'],
        surface=inst.sets['INNER_BOND'], influenceRadius=WHOLE_SURFACE,
        couplingType=KINEMATIC, localCsys=None,
        u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)
    m.DisplacementBC(name='OuterFixed', createStepName='Initial',
        region=inst.sets['OUTER_BOND'], u1=0.0, u2=0.0, u3=0.0)
    m.DisplacementBC(name='DriverGuide', createStepName='Initial',
        region=assembly.sets['DRIVER'], u3=0.0,
        ur1=0.0, ur2=0.0, ur3=0.0)
    m.ContactProperty('RubberContact')
    m.interactionProperties['RubberContact'].NormalBehavior(
        pressureOverclosure=HARD, allowSeparation=ON)
    contact = m.ContactStd(name='CavityContact', createStepName='Initial')
    contact.includedPairs.setValuesInStep(stepName='Initial', useAllstar=ON)
    contact.contactPropertyAssignments.appendInStep(stepName='Initial',
        assignments=((GLOBAL, SELF, 'RubberContact'),))

    report('6/7 Preparing one complete quasi-static fatigue cycle.')
    m.StaticStep(name='Preload', previous='Initial', nlgeom=ON,
        timePeriod=1.0, initialInc=0.01, minInc=1.e-8,
        maxInc=0.05, maxNumInc=2000)
    m.StaticStep(name='FatigueCycle', previous='Preload', nlgeom=ON,
        timePeriod=PERIOD, initialInc=PERIOD/200.0, minInc=1.e-10,
        maxInc=PERIOD/80.0, maxNumInc=5000)
    data = []
    for j in py.range(161):
        t = PERIOD*j/160.0
        data.append((t, math.cos(2.0*math.pi*t/PERIOD)))
    m.TabularAmplitude(name='CycleAmplitude', timeSpan=STEP,
                       data=py.tuple(data))
    m.ConcentratedForce(name='CycleForce', createStepName='Preload',
        region=assembly.sets['DRIVER'], cf1=FMAX)
    m.loads['CycleForce'].setValuesInStep(
        stepName='FatigueCycle', amplitude='CycleAmplitude')
    m.fieldOutputRequests['F-Output-1'].setValues(
        variables=('S', 'NE', 'LE', 'U', 'RF'), numIntervals=4)
    m.fieldOutputRequests['F-Output-1'].setValuesInStep(
        stepName='FatigueCycle', numIntervals=40)
    m.HistoryOutputRequest(name='DriverHistory', createStepName='Preload',
        variables=('U1', 'U2', 'RF1', 'RF2'),
        region=assembly.sets['DRIVER'], frequency=1)
    m.HistoryOutputRequest(name='EnergyHistory', createStepName='Preload',
        variables=('ALLSE', 'ALLWK'), region=MODEL, frequency=1)
    m.keywordBlock.synchVersions(storeNodesAndElements=False)
    for i, block in py.enumerate(py.list(m.keywordBlock.sieBlocks)):
        lines = block.splitlines()
        if lines and lines[0].strip().upper().startswith('*ELEMENT OUTPUT'):
            if 'DIRECTIONS=' not in lines[0].upper().replace(' ', ''):
                lines[0] += ', DIRECTIONS=YES'
                m.keywordBlock.replace(i, '\n'.join(lines))

    report('7/7 Creating Job, saving CAE and writing INP.')
    job = mdb.Job(name=BASE, model=BASE, type=ANALYSIS,
        numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE)
    cae = os.path.join(OUTDIR, BASE+'.cae')
    mdb.saveAs(pathName=cae)
    job.writeInput(consistencyChecking=ON)
    with py.open(BASE+'.inp', 'r') as stream:
        text = stream.read().upper()
    if '*ORIENTATION' not in text or 'ORIENTATION=' not in text:
        raise RuntimeError('Material orientation missing from INP.')
    report('COMPLETE: '+cae)
    report('JOB: '+BASE)
    report('HEX elements: %d; connected components: 1' % py.len(p.elements))
    report('Submit this Job later; import FatigueCycle S/NE into fe-safe.')

try:
    build()
except Exception:
    with py.open(LOG, 'a') as stream:
        stream.write(traceback.format_exc())
    raise
print('SOUTHWEITION_JOB_PREPARED')
import os as _sw_os
try:
    print('SOUTHWEITION_CWD ' + _sw_os.getcwd())
    def _sw_job_model(_job):
        _model = getattr(_job, 'model', '')
        try:
            return str(_model.name)
        except Exception:
            return str(_model or '')
    _sw_models = list(mdb.models.keys())
    _sw_latest_model = _sw_models[-1] if _sw_models else ''
    _sw_name = ''
    for _cand in reversed(list(mdb.jobs.keys())):
        if _sw_latest_model and _sw_job_model(mdb.jobs[_cand]) == _sw_latest_model:
            _sw_name = str(_cand)
            break
    if not _sw_name:
        _sw_names = list(mdb.jobs.keys())
        if _sw_names:
            _sw_name = str(_sw_names[-1])
    if _sw_name:
        try:
            mdb.jobs[_sw_name].setValues(numCpus=4, numDomains=4)
        except Exception:
            pass
        mdb.jobs[_sw_name].writeInput(consistencyChecking=OFF)
        print('SOUTHWEITION_JOB ' + str(_sw_name))
        print('SOUTHWEITION_MODEL ' + str(_sw_latest_model))
        print('SOUTHWEITION_INP ' + _sw_os.path.abspath(str(_sw_name) + '.inp'))
    print('SOUTHWEITION_JOB_PREPARED')
except Exception:
    print('SOUTHWEITION_JOB_PREPARED')