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 math
import os
try:
    import builtins as py
except ImportError:
    import __builtin__ as py
MODEL_ROOT = 'TypeC_Fast_Repaired'
CIRC_DIV = 48
HEAD_DIV = 10
TARGET_SIZE = 0.8
INNER_R = 6.7
CYL_HALF = 13.3
CYL_T = 0.05
HEAD_T = 0.025
RADIUS = INNER_R + CYL_T / 2.0
SUPPORT_X = (-10.0, 10.0)
WOOD_WIDTH = 0.5
WOOD_DEPTH = 0.3
WEB_T = 0.025
RING_T = 0.02
RING_DEPTH = 0.35
SUPPORT_BREADTH = 16.95
SUPPORT_LENGTH = 1.8
SUPPORT_HEIGHT = 7.4
VAPOR_PRESSURE = 1500000.0
LIQUID_DENSITY = 1100.0
GRAVITY = 9.81
CASE_NAMES = ('LC1', 'LC2', 'LC3', 'LC4', 'LC5')
model_name = MODEL_ROOT
index = 1
while model_name in mdb.models or model_name + 'Job' in mdb.jobs:
    model_name = MODEL_ROOT + '_%d' % index
    index += 1
job_name = model_name + 'Job'
input_path = os.path.abspath(model_name + '_mesh.inp')
(coordinates, connectivity, node_cache) = ({}, {}, {})
(element_groups, temperature) = ({}, {})
foundation = py.set()

def add_node(x, y, z):
    key = (py.round(x, 9), py.round(y, 9), py.round(z, 9))
    if key not in node_cache:
        label = py.len(coordinates) + 1
        node_cache[key] = label
        coordinates[label] = (x, y, z)
        temperature[label] = 20.0
    return node_cache[key]

def add_element(element_type, labels, group):
    label = py.len(connectivity) + 1
    connectivity[label] = (element_type, py.tuple(labels))
    element_groups.setdefault(group, []).append(label)
    return label

def ring_node(x, r, j):
    angle = 2.0 * math.pi * j / CIRC_DIV
    return add_node(x, r * math.sin(angle), -r * math.cos(angle))

def dot(a, b):
    return py.sum([a[k] * b[k] for k in py.range(3)])
assert CIRC_DIV % 6 == 0
stations = [-CYL_HALF, 0.0, CYL_HALF]
for x in SUPPORT_X:
    stations.extend([x - WOOD_WIDTH / 2.0, x, x + WOOD_WIDTH / 2.0])
stations = py.sorted(py.set(stations))
axial = [stations[0]]
for i in py.range(py.len(stations) - 1):
    (a, b) = (stations[i], stations[i + 1])
    count = py.max(1, py.int(math.ceil((b - a) / TARGET_SIZE)))
    axial.extend([a + (b - a) * k / count for k in py.range(1, count + 1)])
sections = []
for k in py.range(1, HEAD_DIV):
    angle = math.pi * k / (2.0 * HEAD_DIV)
    sections.append((-CYL_HALF - RADIUS * math.cos(angle), RADIUS * math.sin(angle)))
sections.extend([(x, RADIUS) for x in axial])
for k in py.range(1, HEAD_DIV):
    angle = math.pi * k / (2.0 * HEAD_DIV)
    sections.append((CYL_HALF + RADIUS * math.sin(angle), RADIUS * math.cos(angle)))
shell_rings = []
for (x, r) in sections:
    shell_rings.append([ring_node(x, r, j) for j in py.range(CIRC_DIV)])
pole_left = add_node(-CYL_HALF - RADIUS, 0.0, 0.0)
pole_right = add_node(CYL_HALF + RADIUS, 0.0, 0.0)
for j in py.range(CIRC_DIV):
    q = (j + 1) % CIRC_DIV
    add_element('S3', [pole_left, shell_rings[0][q], shell_rings[0][j]], 'HEAD')
    add_element('S3', [shell_rings[-1][j], shell_rings[-1][q], pole_right], 'HEAD')
for k in py.range(py.len(shell_rings) - 1):
    middle = 0.5 * (sections[k][0] + sections[k + 1][0])
    group = 'CYLINDER' if py.abs(middle) < CYL_HALF else 'HEAD'
    for j in py.range(CIRC_DIV):
        q = (j + 1) % CIRC_DIV
        add_element('S4R', [shell_rings[k][j], shell_rings[k][q], shell_rings[k + 1][q], shell_rings[k + 1][j]], group)
pressure_elements = element_groups['CYLINDER'] + element_groups['HEAD']
cold_nodes = py.set()
for e in pressure_elements:
    cold_nodes.update(connectivity[e][1])
for x in (SUPPORT_X[0], 0.0, SUPPORT_X[1]):
    for j in py.range(CIRC_DIV):
        q = (j + 1) % CIRC_DIV
        labels = [ring_node(x, RADIUS, j), ring_node(x, RADIUS, q), ring_node(x, RADIUS - RING_DEPTH, q), ring_node(x, RADIUS - RING_DEPTH, j)]
        add_element('S4R', labels, 'STIFFENER')
        cold_nodes.update(labels)
for n in cold_nodes:
    temperature[n] = -35.0
half_arc = CIRC_DIV // 6
arc_indices = py.list(py.range(-half_arc, half_arc + 1))
inner_pad_r = RADIUS + CYL_T / 2.0
outer_pad_r = inner_pad_r + WOOD_DEPTH
base_z = -0.5 * outer_pad_r - SUPPORT_HEIGHT
contact_elements = [[], []]
for (support, center_x) in py.enumerate(SUPPORT_X):
    faces_x = [center_x - WOOD_WIDTH / 2.0, center_x + WOOD_WIDTH / 2.0]
    pad = {}
    for layer in py.range(3):
        radius = inner_pad_r + WOOD_DEPTH * layer / 2.0
        for side in py.range(2):
            for j in arc_indices:
                n = ring_node(faces_x[side], radius, j)
                pad[side, layer, j] = n
                temperature[n] = -35.0 + 27.5 * layer
    for layer in py.range(2):
        for j in arc_indices[:-1]:
            labels = [pad[0, layer, j], pad[0, layer, j + 1], pad[1, layer, j + 1], pad[1, layer, j], pad[0, layer + 1, j], pad[0, layer + 1, j + 1], pad[1, layer + 1, j + 1], pad[1, layer + 1, j]]
            e = add_element('C3D8R', labels, 'WOOD')
            if layer == 0:
                contact_elements[support].append(e)
    for j in arc_indices[:-1]:
        add_element('S4R', [pad[0, 2, j], pad[0, 2, j + 1], pad[1, 2, j + 1], pad[1, 2, j]], 'SADDLE')
    web = {}
    for side in py.range(2):
        for j in arc_indices:
            angle = 2.0 * math.pi * j / CIRC_DIV
            (top_y, top_z) = (outer_pad_r * math.sin(angle), -outer_pad_r * math.cos(angle))
            bottom_y = SUPPORT_BREADTH * j / (2.0 * half_arc)
            for k in py.range(9):
                fraction = k / 8.0
                n = add_node(faces_x[side], top_y * (1.0 - fraction) + bottom_y * fraction, top_z * (1.0 - fraction) + base_z * fraction)
                web[side, j, k] = n
                if k == 8:
                    foundation.add(n)
        for j in arc_indices[:-1]:
            for k in py.range(8):
                add_element('S4R', [web[side, j, k], web[side, j + 1, k], web[side, j + 1, k + 1], web[side, j, k + 1]], 'SADDLE')
    for j in arc_indices[::4]:
        for k in py.range(8):
            add_element('S4R', [web[0, j, k], web[1, j, k], web[1, j, k + 1], web[0, j, k + 1]], 'SADDLE')
    foot_x = [center_x - SUPPORT_LENGTH / 2.0, faces_x[0], faces_x[1], center_x + SUPPORT_LENGTH / 2.0]
    for i in py.range(3):
        for j in arc_indices[:-1]:
            ya = SUPPORT_BREADTH * j / (2.0 * half_arc)
            yb = SUPPORT_BREADTH * (j + 1) / (2.0 * half_arc)
            labels = [add_node(foot_x[i], ya, base_z), add_node(foot_x[i + 1], ya, base_z), add_node(foot_x[i + 1], yb, base_z), add_node(foot_x[i], yb, base_z)]
            add_element('S4R', labels, 'SADDLE')
            foundation.update(labels)
volume_samples = []
total_length = 2.0 * (CYL_HALF + INNER_R)
for i in py.range(80):
    x = -total_length / 2.0 + total_length * (i + 0.5) / 80.0
    distance = py.max(0.0, py.abs(x) - CYL_HALF)
    radius = math.sqrt(py.max(0.0, INNER_R ** 2 - distance ** 2))
    weight = math.pi * radius ** 2 * total_length / (80.0 * 8.0 * 24.0)
    for k in py.range(8):
        radial = radius * math.sqrt((k + 0.5) / 8.0)
        for j in py.range(24):
            angle = 2.0 * math.pi * (j + 0.5) / 24.0
            volume_samples.append(((x, radial * math.sin(angle), radial * math.cos(angle)), weight))

def free_surface(acceleration, full):
    magnitude = math.sqrt(dot(acceleration, acceleration))
    if full:
        return -py.abs(acceleration[0]) * CYL_HALF - magnitude * INNER_R
    values = py.sorted([(dot(acceleration, p), w) for (p, w) in volume_samples])
    limit = 0.1 * py.sum([item[1] for item in values])
    accumulated = 0.0
    for (potential, weight) in values:
        accumulated += weight
        if accumulated >= limit:
            return potential
    return values[-1][0]
deck = ['*HEADING', model_name, '*PREPRINT, ECHO=NO, MODEL=NO, HISTORY=NO', '*NODE']
for n in py.sorted(coordinates.keys()):
    deck.append('%d, %.10g, %.10g, %.10g' % ((n,) + coordinates[n]))
for element_type in ('S4R', 'S3', 'C3D8R'):
    deck.append('*ELEMENT, TYPE=' + element_type)
    for e in py.sorted(connectivity.keys()):
        (kind, labels) = connectivity[e]
        if kind == element_type:
            deck.append(', '.join([py.str(v) for v in (e,) + labels]))

def write_set(kind, name, labels):
    deck.append('*%s, %s=%s' % (kind, kind, name))
    values = py.sorted(py.set(py.list(labels)))
    for start in py.range(0, py.len(values), 16):
        deck.append(', '.join([py.str(v) for v in values[start:start + 16]]))
for group in py.sorted(element_groups.keys()):
    write_set('ELSET', group, element_groups[group])
write_set('ELSET', 'TANK_SHELL', pressure_elements)
write_set('ELSET', 'ALL_ELEMENTS', py.list(connectivity.keys()))
write_set('NSET', 'ALL_NODES', py.list(coordinates.keys()))
write_set('NSET', 'FOUNDATION', py.list(foundation))
for (material, elastic, density, expansion) in [('TANK_STEEL', 200000000000.0, 7850, 1.1e-05), ('SADDLE_STEEL', 200000000000.0, 7850, 1.1e-05), ('SUPPORT_WOOD', 8500000000.0, 1350, 9.2e-06)]:
    deck.extend(['*MATERIAL, NAME=' + material, '*ELASTIC', '%.9g, 0.3' % elastic, '*DENSITY', '%.9g' % density, '*EXPANSION', '%.9g' % expansion])
for (group, material, thickness) in [('CYLINDER', 'TANK_STEEL', CYL_T), ('HEAD', 'TANK_STEEL', HEAD_T), ('STIFFENER', 'TANK_STEEL', RING_T), ('SADDLE', 'SADDLE_STEEL', WEB_T)]:
    deck.extend(['*SHELL SECTION, ELSET=%s, MATERIAL=%s' % (group, material), '%.9g, 5' % thickness])
deck.extend(['*SOLID SECTION, ELSET=WOOD, MATERIAL=SUPPORT_WOOD', ',', '*SURFACE, TYPE=ELEMENT, NAME=TANK_OUTSIDE', 'TANK_SHELL, SPOS'])
for i in py.range(2):
    write_set('ELSET', 'PAD_TOP_%d' % i, contact_elements[i])
    deck.extend(['*SURFACE, TYPE=ELEMENT, NAME=PAD_SURFACE_%d' % i, 'PAD_TOP_%d, S1' % i])
deck.extend(['*TIE, NAME=FIXED_SUPPORT, ADJUST=YES, POSITION TOLERANCE=0.06', 'PAD_SURFACE_0, TANK_OUTSIDE', '*SURFACE INTERACTION, NAME=FRICTIONLESS', '*SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE=HARD', '*CONTACT PAIR, INTERACTION=FRICTIONLESS, TYPE=SURFACE TO SURFACE, SMALL SLIDING, ADJUST=0.04', 'PAD_SURFACE_1, TANK_OUTSIDE', '*BOUNDARY', 'FOUNDATION, 1, 3, 0.', '*INITIAL CONDITIONS, TYPE=TEMPERATURE', 'ALL_NODES, 20.'])
outputs = ['*OUTPUT, FIELD, FREQUENCY=999', '*NODE OUTPUT', 'U, RF, NT', '*ELEMENT OUTPUT', 'S, E', '*CONTACT OUTPUT', 'CSTRESS, CDISP', '*OUTPUT, HISTORY, FREQUENCY=1', '*ENERGY OUTPUT', 'ALLIE, ALLSE, ALLWK']
deck.extend(['*STEP, NAME=Thermal, NLGEOM=NO, INC=150', '*STATIC', '0.2, 1.0, 1.e-7, 0.5', '*TEMPERATURE'])
for n in py.sorted(temperature.keys()):
    if temperature[n] != 20.0:
        deck.append('%d, %.9g' % (n, temperature[n]))
deck.extend(outputs + ['*END STEP'])
g = GRAVITY
accelerations = {'LC1': (0.25 * g, 0.0, -g), 'LC2': (0.0, 0.7 * g, -g), 'LC3': (0.0, 0.0, -1.87 * g), 'LC4': (0.0, 0.5 * g, -math.sqrt(3.0) * g / 2.0), 'LC5': (0.5 * g, 0.0, -g)}
for case in CASE_NAMES:
    acceleration = accelerations[case]
    level = free_surface(acceleration, case == 'LC5')
    magnitude = math.sqrt(dot(acceleration, acceleration))
    direction = py.tuple([v / magnitude for v in acceleration])
    deck.extend(['*STEP, NAME=%s, NLGEOM=NO, INC=150' % case, '*STATIC', '0.2, 1.0, 1.e-7, 0.5', '*DLOAD, OP=NEW', 'ALL_ELEMENTS, GRAV, %.9g, %.9g, %.9g, %.9g' % ((magnitude,) + direction)])
    for e in pressure_elements:
        labels = connectivity[e][1]
        center = [py.sum([coordinates[n][k] for n in labels]) / py.len(labels) for k in py.range(3)]
        pressure = VAPOR_PRESSURE + LIQUID_DENSITY * py.max(0.0, dot(acceleration, center) - level)
        deck.append('%d, P, %.10g' % (e, pressure))
    deck.extend(outputs + ['*END STEP'])
with py.open(input_path, 'w') as stream:
    stream.write('\n'.join(deck) + '\n')
model = mdb.ModelFromInputFile(name=model_name, inputFileName=input_path)
job = mdb.Job(name=job_name, model=model_name, type=ANALYSIS, numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE)
_sw_record_write(job, consistencyChecking=ON)
if 'Viewport: 1' in session.viewports:
    session.viewports['Viewport: 1'].setValues(displayedObject=model.rootAssembly)
    session.viewports['Viewport: 1'].view.fitView()
print('Model: ' + model_name)
print('Job: ' + job_name)
print('Nodes: %d; Elements: %d' % (py.len(coordinates), py.len(connectivity)))
print('Input file written. Job has not been 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')