# -*- 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。
# 保留原截面、连续橡胶外环、材料、接触和C3D8H网格。
# Z向载荷作用于内孔刚性耦合参考点，外圆固定。
# 这是悬置轴向相对运动，橡胶主要发生剪切，并非端面拉伸。
# 12000 N为试算幅值，不是论文载荷或已验证的承载极限。
# 原外环厚1.2 mm、轴向厚30 mm及空腔样条近似继续保留。
# 只建立模型、保存CAE、写INP；不提交计算、不读取ODB。

OUT_ROOT = os.path.join(tempfile.gettempdir(), 'AbaqusRubberRuns')
if not os.path.isdir(OUT_ROOT):
    os.makedirs(OUT_ROOT)
OUTDIR = tempfile.mkdtemp(prefix='AxialZ_', dir=OUT_ROOT)
BASE = 'RubberZ_' + 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 = 12000.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 model: ' + BASE)
    report('Output directory: ' + OUTDIR)
    if not 0.0 < RING_T < 2.0:
        raise ValueError('Ring thickness must be between 0 and 2 mm.')
    if FMAX <= 0.0:
        raise ValueError('FMAX must be positive.')
    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.')

    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 section.')
    report('2/7 Generating original quadrilateral section 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:
        enodes = py.list(element.getNodes())
        if py.len(enodes) != 4:
            raise RuntimeError('Non-quadrilateral element detected.')
        for node in enodes:
            parent.setdefault(node.label, node.label)
        first = root(enodes[0].label)
        for node in enodes[1:]:
            parent[root(node.label)] = first
    if py.len(py.set([root(n) for n in py.list(parent.keys())])) != 1:
        raise RuntimeError('Disconnected rubber mesh.')

    report('3/7 Extruding into six layers of 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 unchanged NR45 material and 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 Setting axial guide and original 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],))
    driver = assembly.sets['DRIVER']
    m.Coupling(name='InnerDriver', controlPoint=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)
    # 释放U3，约束U1/U2及转动，避免原脚本将Z向锁死。
    m.DisplacementBC(name='DriverGuide', createStepName='Initial',
        region=driver, u1=0.0, u2=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 Z-force preload and fully reversed cycle.')
    # 自动增量、全牛顿法；不降低默认平衡精度，不加人工稳定阻尼。
    for name, previous, duration, inc, maxinc in (
        ('Preload', 'Initial', 1.0, 0.002, 0.02),
        ('FatigueCycle', 'Preload', PERIOD, PERIOD/400.0, PERIOD/80.0)):
        m.StaticStep(name=name, previous=previous, nlgeom=ON,
            timePeriod=duration, initialInc=inc,
            minInc=duration*1.e-9, maxInc=maxinc, maxNumInc=5000,
            timeIncrementationMethod=AUTOMATIC, solutionTechnique=FULL_NEWTON)

    m.SmoothStepAmplitude(name='PreloadRamp', timeSpan=STEP,
                         data=((0.0, 0.0), (1.0, 1.0)))
    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='AxialForce', createStepName='Preload',
        region=driver, cf3=FMAX, amplitude='PreloadRamp', follower=OFF)
    m.loads['AxialForce'].setValuesInStep(
        stepName='FatigueCycle', amplitude='CycleAmplitude')

    # 输出积分点应力/应变及材料方向，保留40个循环输出区间。
    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=('U3', 'CF3', 'RF3'), region=driver, frequency=1)
    m.HistoryOutputRequest(name='EnergyHistory', createStepName='Preload',
        variables=('ALLSE', 'ALLWK'), region=MODEL, frequency=1)
    # 受力控制的自由U3方向，参考点RF3不等于施加载荷；查看CF3。
    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 four-CPU Job and writing input.')
    job = mdb.Job(name=BASE, model=BASE, type=ANALYSIS,
        numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE)
    job.writeInput(consistencyChecking=ON)
    cae = os.path.join(OUTDIR, BASE + '.cae')
    mdb.saveAs(pathName=cae)
    report('CAE: ' + cae)
    report('INP: ' + os.path.join(OUTDIR, BASE + '.inp'))
    report('JOB: ' + BASE)
    report('HEX elements: %d' % py.len(p.elements))
    report('Z force cycle: +12000 -> 0 -> -12000 -> 0 -> +12000 N')
    report('Submit this existing Job in CAE Job Manager.')
    report('ODB is produced only after the solver runs.')
    report('For rubber fatigue use complete FatigueCycle S/NE history.')
    report('Large-strain material validity and convergence need verification.')

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')