# -*- coding: utf-8 -*-
from abaqus import mdb, session
import abaqusConstants as ac
import mesh
import math

# Units: mm, N, s, tonne, MPa.
# SINGLE tube, ZERO pre-twist.
# Hole long axes are axial or circumferential on the cylindrical surface.
# Circumferential holes follow cylinder curvature, not a global straight line.
# This is an untwisted comparison model, NOT the paper's pre-twisted specimen.
# A symmetric model need not develop a nonzero compression-induced twist.
# Dimensions and material coefficients remain provisional.

MODEL_BASE = 'Single_Straight_Holes'
JOB_BASE = 'Single_Straight_Compression'
CELL_CIRC, CELL_AXIAL = 8, 4
CELL_SIZE = 10.0
PITCH = CELL_SIZE / 2.0
NC, NZ = 2 * CELL_CIRC, 2 * CELL_AXIAL
HOLE_LONG, HOLE_SHORT = 2.10, 0.60
WALL, COLLAR = 0.50, 1.50

# Reduced mesh: 16 segments per hole, 2 radial subdivisions.
# Keep 2 elements through thickness.
NANG, NRAD, NTHICK = 16, 2, 2
STRAIN = 0.10
STEP_TIME = 0.05
DENSITY = 1.10e-9
# Illustrative Ogden data; not the paper's material fit.
OGDEN_COEFFS = (1.5, 2.0, 0.3, 4.0, 0.2, 6.0, 0.04, 0.0, 0.0)

CIRC = CELL_CIRC * CELL_SIZE
ACTIVE_HEIGHT = CELL_AXIAL * CELL_SIZE
HEIGHT = ACTIVE_HEIGHT + 2.0 * COLLAR
RADIUS = CIRC / (2.0 * math.pi)

model_name, job_name = MODEL_BASE, JOB_BASE
index = 1
while model_name in mdb.models or job_name in mdb.jobs:
    model_name = MODEL_BASE + '_%02d' % index
    job_name = JOB_BASE + '_%02d' % index
    index += 1

model = mdb.Model(name=model_name, modelType=ac.STANDARD_EXPLICIT)
mat = model.Material(name='TPU_Provisional')
mat.Density(table=((DENSITY,),))
mat.Hyperelastic(
    materialType=ac.ISOTROPIC, type=ac.OGDEN, n=3,
    testData=ac.OFF, volumetricResponse=ac.VOLUMETRIC_DATA,
    table=(OGDEN_COEFFS,))
model.HomogeneousSolidSection(
    name='WallSection', material=mat.name, thickness=None)

# Unwrapped coordinates: s = circumferential, z = axial.
uv, quads, lookup = [], [], {}

def add_uv(s, z):
    s = s % CIRC
    if abs(s) < 1.0e-7 or abs(s - CIRC) < 1.0e-7:
        s = 0.0
    key = (round(s, 7), round(z, 7))
    if key not in lookup:
        lookup[key] = len(uv)
        uv.append((s, z))
    return lookup[key]

# Alternating vertical/horizontal ellipse holes.
for iz in range(NZ):
    for ic in range(NC):
        cs = (ic + 0.5) * PITCH
        cz = COLLAR + (iz + 0.5) * PITCH
        ax, az = HOLE_SHORT, HOLE_LONG
        if (ic + iz) % 2:
            ax, az = HOLE_LONG, HOLE_SHORT
        rings = []
        for j in range(NRAD + 1):
            fraction = float(j) / NRAD
            ring = []
            for k in range(NANG):
                angle = 2.0 * math.pi * k / NANG
                co, si = math.cos(angle), math.sin(angle)
                inner = 1.0 / math.sqrt((co / ax)**2 + (si / az)**2)
                outer = 0.5 * PITCH / max(abs(co), abs(si))
                distance = inner + fraction * (outer - inner)
                ring.append(add_uv(cs + distance * co,
                                   cz + distance * si))
            rings.append(ring)
        for j in range(NRAD):
            for k in range(NANG):
                kn = (k + 1) % NANG
                quads.append((rings[j][k], rings[j + 1][k],
                              rings[j + 1][kn], rings[j][kn]))

# Conforming end collars.
for interface_z, end_z in ((COLLAR, 0.0),
                           (COLLAR + ACTIVE_HEIGHT, HEIGHT)):
    edge = sorted((point[0], i) for i, point in enumerate(uv)
                  if abs(point[1] - interface_z) < 1.0e-7)
    end_nodes = [add_uv(s, end_z) for s, i in edge]
    for k in range(len(edge)):
        kn = (k + 1) % len(edge)
        if end_z < interface_z:
            quads.append((end_nodes[k], end_nodes[kn],
                          edge[kn][1], edge[k][1]))
        else:
            quads.append((edge[k][1], edge[kn][1],
                          end_nodes[kn], end_nodes[k]))

part = model.Part(name='SingleTube',
                  dimensionality=ac.THREE_D, type=ac.DEFORMABLE_BODY)
nplane = len(uv)

for ir in range(NTHICK + 1):
    radius = RADIUS - WALL / 2.0 + WALL * float(ir) / NTHICK
    for i, (s, z) in enumerate(uv):
        # IMPORTANT: angle is independent of z. No geometric pre-twist.
        phi = 2.0 * math.pi * s / CIRC
        part.Node(
            coordinates=(radius * math.cos(phi),
                         radius * math.sin(phi), z),
            label=ir * nplane + i + 1)

label = 1
for ir in range(NTHICK):
    for quad in quads:
        ids = [ir * nplane + i for i in quad]
        ids += [(ir + 1) * nplane + i for i in quad]
        part.Element(nodes=tuple(part.nodes[i] for i in ids),
                     elemShape=ac.HEX8, label=label)
        label += 1

part.Set(name='ALL', elements=part.elements[:])
part.setElementType(
    regions=(part.elements[:],),
    elemTypes=(mesh.ElemType(
        elemCode=ac.C3D8R, elemLibrary=ac.EXPLICIT,
        hourglassControl=ac.ENHANCED),))
part.SectionAssignment(region=part.sets['ALL'],
                       sectionName='WallSection')

for name, target in (('BOTTOM', 0.0), ('TOP', HEIGHT)):
    labels = tuple(ir * nplane + i + 1
                   for ir in range(NTHICK + 1)
                   for i, point in enumerate(uv)
                   if abs(point[1] - target) < 1.0e-7)
    part.Set(name=name,
             nodes=part.nodes.sequenceFromLabels(labels=labels))

# A tracked rim node helps inspect rotation.
top_labels = [i + 1 for i, point in enumerate(uv)
              if abs(point[1] - HEIGHT) < 1.0e-7]
part.Set(name='TOP_MARKER',
         nodes=part.nodes.sequenceFromLabels(labels=(top_labels[0],)))

assembly = model.rootAssembly
assembly.DatumCsysByDefault(ac.CARTESIAN)
instance = assembly.Instance(name='Tube-1', part=part, dependent=ac.ON)

# Ideal end fixtures represented by reference points and couplings.
for name, z in (('BOTTOM', 0.0), ('TOP', HEIGHT)):
    feature = assembly.ReferencePoint(point=(0.0, 0.0, z))
    assembly.Set(
        name='RP_' + name,
        referencePoints=(assembly.referencePoints[feature.id],))
    model.Coupling(
        name='End_' + name,
        controlPoint=assembly.sets['RP_' + name],
        surface=instance.sets[name],
        influenceRadius=ac.WHOLE_SURFACE,
        couplingType=ac.KINEMATIC, localCsys=None,
        u1=ac.ON, u2=ac.ON, u3=ac.ON,
        ur1=ac.OFF, ur2=ac.OFF, ur3=ac.OFF)

model.ExplicitDynamicsStep(
    name='Compress', previous='Initial',
    timePeriod=STEP_TIME, nlgeom=ac.ON)
model.SmoothStepAmplitude(
    name='Ramp', timeSpan=ac.STEP,
    data=((0.0, 0.0), (STEP_TIME, 1.0)))
model.EncastreBC(
    name='BottomFixed', createStepName='Initial',
    region=assembly.sets['RP_BOTTOM'])
model.DisplacementBC(
    name='TopGuide', createStepName='Initial',
    region=assembly.sets['RP_TOP'],
    u1=0.0, u2=0.0, u3=ac.UNSET,
    ur1=0.0, ur2=0.0, ur3=ac.UNSET)
model.DisplacementBC(
    name='CompressionLoad', createStepName='Compress',
    region=assembly.sets['RP_TOP'],
    u3=-STRAIN * ACTIVE_HEIGHT, amplitude='Ramp')

model.ContactProperty('Frictionless')
model.interactionProperties['Frictionless'].NormalBehavior(
    pressureOverclosure=ac.HARD, allowSeparation=ac.ON)
contact = model.ContactExp(name='SelfContact', createStepName='Initial')
contact.includedPairs.setValuesInStep(
    stepName='Initial', useAllstar=ac.ON)
contact.contactPropertyAssignments.appendInStep(
    stepName='Initial',
    assignments=((ac.GLOBAL, ac.SELF, 'Frictionless'),))

model.fieldOutputRequests['F-Output-1'].setValues(
    variables=('S', 'U'), numIntervals=15)
for name in list(model.historyOutputRequests.keys()):
    del model.historyOutputRequests[name]
for name in ('BOTTOM', 'TOP'):
    model.HistoryOutputRequest(
        name='History_' + name, createStepName='Compress',
        region=assembly.sets['RP_' + name],
        variables=('U3', 'UR3', 'RF3', 'RM3'), numIntervals=150)
model.HistoryOutputRequest(
    name='Marker', createStepName='Compress',
    region=instance.sets['TOP_MARKER'],
    variables=('U1', 'U2', 'U3'), numIntervals=150)
model.HistoryOutputRequest(
    name='Energy', createStepName='Compress',
    variables=('ALLIE', 'ALLKE', 'ALLAE', 'ALLVD', 'ALLWK', 'ETOTAL'),
    numIntervals=150)

assembly.regenerate()
job = mdb.Job(
    name=job_name, model=model_name, type=ac.ANALYSIS,
    numCpus=4, numDomains=4,
    parallelizationMethodExplicit=ac.DOMAIN,
    explicitPrecision=ac.DOUBLE,
    memory=80, memoryUnits=ac.PERCENTAGE)
job.writeInput(consistencyChecking=ac.ON)

if len(session.viewports):
    viewport = session.viewports[list(session.viewports.keys())[0]]
    viewport.setValues(displayedObject=assembly)
    viewport.view.fitView()

print('Model: ' + model_name)
print('Job Manager job: ' + job_name)
print('Input file: ' + job_name + '.inp')
print('One tube; zero geometric pre-twist.')
print('Solid elements: %d' % len(part.elements))
print('Top rotation in degrees: History_TOP.UR3 * 180/pi.')
print('Check energy ratios and mesh/time convergence.')
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')