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
pi = math.pi
(sin, cos, sqrt) = (math.sin, math.cos, math.sqrt)
model_name = 'SpacerGenerator'
serial = 1
while model_name in mdb.models or model_name + '_Job' in mdb.jobs or os.path.exists(model_name + '_Job.inp'):
    model_name = 'SpacerGenerator_%d' % serial
    serial += 1
job_name = model_name + '_Job'
inp_path = os.path.abspath(job_name + '.inp')
pivot = 121.5
wire_radius = 319.0
arm_length = wire_radius - sqrt(2.0) * pivot
load_value = 15170.0

def rotated(point, angle):
    (x, y) = point
    return (x * cos(angle) - y * sin(angle), x * sin(angle) + y * cos(angle))

def circular_points(cx, r, a0, a1, n):
    return py.list(((cx + r * cos(a0 + (a1 - a0) * i / n), r * sin(a0 + (a1 - a0) * i / n)) for i in range(n + 1)))

def rectangular_boundary(values):
    (lo, hi) = (values[0], values[-1])
    return [(hi, t) for t in values[:-1]] + [(t, hi) for t in values[:0:-1]] + [(lo, t) for t in values[:0:-1]] + [(t, lo) for t in values[:-1]] + [(hi, lo)]

def rounded(point):
    return py.tuple((round(value, 7) for value in point))

def edge_id(a, b):
    return py.tuple(py.sorted([rounded(a), rounded(b)]))

class HexMesh:

    def __init__(self):
        self.nodes = []
        self.elements = []
        self.node_ids = {}
        self.surfaces = {}

    def node(self, point):
        key = rounded(point)
        if key not in self.node_ids:
            self.nodes.append(py.tuple(point))
            self.node_ids[key] = len(self.nodes)
        return self.node_ids[key]

    def surface_entry(self, name, side, label):
        if name:
            self.surfaces.setdefault(name, {}).setdefault(side, []).append(label)

    def extrude_quad(self, quad, zlevels, edges, bottom=None, top=None):
        area = py.sum((quad[i][0] * quad[(i + 1) % 4][1] - quad[(i + 1) % 4][0] * quad[i][1] for i in range(4)))
        if abs(area) < 1e-09:
            raise ValueError('Degenerate mesh quadrilateral')
        if area < 0.0:
            quad = [quad[0], quad[3], quad[2], quad[1]]
        for k in range(len(zlevels) - 1):
            connectivity = py.tuple((self.node((point[0], point[1], z)) for z in (zlevels[k], zlevels[k + 1]) for point in quad))
            self.elements.append(connectivity)
            label = len(self.elements)
            if k == 0:
                self.surface_entry(bottom, 1, label)
            if k == len(zlevels) - 2:
                self.surface_entry(top, 2, label)
            for side in range(4):
                group = edges.get(edge_id(quad[side], quad[(side + 1) % 4]))
                self.surface_entry(group, side + 3, label)

    def band(self, inner, outer, nrad, zlevels, tags=None):
        (inner, outer) = (py.list(inner), py.list(outer))
        tags = tags or {}
        if len(inner) != len(outer):
            raise ValueError('Boundary point counts differ')
        for i in range(len(inner) - 1):

            def point(h, j):
                t = float(j) / nrad
                return py.tuple(((1 - t) * inner[h][v] + t * outer[h][v] for v in range(2)))
            for j in range(nrad):
                quad = [point(i, j), point(i + 1, j), point(i + 1, j + 1), point(i, j + 1)]
                edges = {}
                conditions = [(j == 0, 0, 1, 'inner'), (j == nrad - 1, 2, 3, 'outer'), (i == 0, 3, 0, 'start'), (i == len(inner) - 2, 1, 2, 'end')]
                for (active, a, b, key) in conditions:
                    if active and tags.get(key):
                        edges[edge_id(quad[a], quad[b])] = tags[key]
                self.extrude_quad(quad, zlevels, edges, tags.get('bottom'), tags.get('top'))

    def cylinder(self, center, r, zlevels, outer=None, bottom=None):
        a = 0.4 * r
        grid = py.list((-a + 2 * a * i / 4.0 for i in range(5)))
        square = rectangular_boundary(grid)
        inner = py.list(((center[0] + x, center[1] + y) for (x, y) in square))
        circle = py.list(((center[0] + r * x / sqrt(x * x + y * y), center[1] + r * y / sqrt(x * x + y * y)) for (x, y) in square))
        self.band(inner, circle, 2, zlevels, {'outer': outer, 'bottom': bottom})
        for i in range(4):
            for j in range(4):
                quad = py.list(((center[0] + grid[u], center[1] + grid[v]) for (u, v) in [(i, j), (i + 1, j), (i + 1, j + 1), (i, j + 1)]))
                self.extrude_quad(quad, zlevels, {}, bottom)

    def export(self, stream, name, material, element_type):
        stream.write('*Part, name=%s\n*Node\n' % name)
        stream.writelines(('%d, %.10g, %.10g, %.10g\n' % (label, p[0], p[1], p[2]) for (label, p) in enumerate(self.nodes, 1)))
        stream.write('*Element, type=%s\n' % element_type)
        stream.writelines(('%d, %s\n' % (label, ', '.join((str(v) for v in conn))) for (label, conn) in enumerate(self.elements, 1)))
        stream.write('*Elset, elset=ALL, generate\n1, %d, 1\n' % len(self.elements))
        stream.write('*Solid Section, elset=ALL, material=%s\n,\n' % material)
        for (group, faces) in self.surfaces.items():
            for (face, labels) in faces.items():
                stream.write('*Elset, elset=E_%s_%d\n' % (group, face))
                stream.writelines((', '.join((str(v) for v in labels[j:j + 16])) + '\n' for j in range(0, len(labels), 16)))
            stream.write('*Surface, type=ELEMENT, name=%s\n' % group)
            stream.writelines(('E_%s_%d, S%d\n' % (group, face, face) for face in faces))
        stream.write('*End Part\n')
frame = HexMesh()
frame_z = [-8.0, -4.0, 0.0, 4.0, 8.0]
for quadrant in range(4):
    angle = quadrant * pi / 2
    for center_x in (-50.0, 50.0):
        (inner, outer) = ([], [])
        for (side, count) in enumerate((6, 20, 6, 20)):
            for i in range(count):
                t = float(i) / count
                if side == 0:
                    a = -pi / 2 + pi * t
                    (q, p) = ((34 + 6 * cos(a), 6 * sin(a)), (50, -18 + 36 * t))
                elif side == 1:
                    (q, p) = ((34 - 68 * t, 6), (50 - 100 * t, 18))
                elif side == 2:
                    a = pi / 2 + pi * t
                    (q, p) = ((-34 + 6 * cos(a), 6 * sin(a)), (-50, 18 - 36 * t))
                else:
                    (q, p) = ((-34 + 68 * t, -6), (-50 + 100 * t, -18))
                inner.append(rotated((q[0] + center_x, q[1] + 125), angle))
                outer.append(rotated((p[0] + center_x, p[1] + 125), angle))
        frame.band(inner + [inner[0]], outer + [outer[0]], 3, frame_z)
    boundary = rectangular_boundary([100, 107, 113, 119, 125, 131, 137, 143])
    hole = []
    for (x, y) in boundary:
        distance = sqrt((x - pivot) ** 2 + (y - pivot) ** 2)
        hole.append((pivot + 9 * (x - pivot) / distance, pivot + 9 * (y - pivot) / distance))
    frame.band((rotated(p, angle) for p in hole), (rotated(p, angle) for p in boundary), 5, frame_z, {'inner': 'HOLE%d' % quadrant})
(body, cover, pad, pin, locks) = py.tuple((HexMesh() for i in range(5)))
arm_z = [8.0, 12.0, 16.0, 20.0, 24.0]
L = arm_length
body.band(circular_points(0, 9, 0, 2 * pi, 32), circular_points(0, 18, 0, 2 * pi, 32), 3, arm_z, {'inner': 'EYE'})
body.band(circular_points(L, 21, pi / 2, 3 * pi / 2, 16), circular_points(L, 30, pi / 2, 3 * pi / 2, 16), 3, arm_z, {'inner': 'LINER', 'start': 'SPLIT', 'end': 'SPLIT', 'top': 'TOP'})
body.band(circular_points(0, 18, -pi / 4, pi / 4, 8), circular_points(L, 30, 5 * pi / 4, 3 * pi / 4, 8), 28, arm_z)
cover.band(circular_points(L, 21, -pi / 2, pi / 2, 16), circular_points(L, 30, -pi / 2, pi / 2, 16), 3, arm_z, {'inner': 'LINER', 'start': 'SPLIT', 'end': 'SPLIT', 'top': 'TOP'})
pad.band(circular_points(L, 17, 0, 2 * pi, 32), circular_points(L, 21, 0, 2 * pi, 32), 2, arm_z, {'inner': 'BORE', 'outer': 'OUTER'})
pin.cylinder((0, 0), 9, frame_z, outer='LOWER')
pin.cylinder((0, 0), 9, arm_z, outer='UPPER')
for y in (-25.5, 25.5):
    locks.cylinder((L, y), 2.0, [24.0, 26.0, 28.0], bottom='BOTTOM')
specs = [('FRAME', frame, 'ZL104', 'C3D8R'), ('BODY', body, 'ZL104', 'C3D8R'), ('COVER', cover, 'ZL104', 'C3D8R'), ('PAD', pad, 'PAD_EQ', 'C3D8RH'), ('PIN', pin, 'STEEL_EQ', 'C3D8R'), ('LOCKS', locks, 'STEEL_EQ', 'C3D8R')]
with py.open(inp_path, 'w') as stream:
    stream.write('*Heading\nApproximate spacer with generator-based mesh\n')
    for (part_name, data, material, code) in specs:
        data.export(stream, part_name, material, code)
    stream.write('*Assembly, name=ASSEMBLY\n')
    stream.writelines(('*Instance, name=%s_SEED, part=%s\n*End Instance\n' % (item[0], item[0]) for item in specs))
    stream.write('*End Assembly\n')
    for (mat, e, nu, rho) in [('ZL104', 69000.0, 0.34, 2.65e-09), ('STEEL_EQ', 210000.0, 0.3, 7.85e-09), ('PAD_EQ', 100.0, 0.49, 1.1e-09)]:
        stream.write('*Material, name=%s\n*Elastic\n%.10g, %.10g\n*Density\n%.10g\n' % (mat, e, nu, rho))
model = mdb.ModelFromInputFile(name=model_name, inputFileName=inp_path)
assembly = model.rootAssembly
for instance_name in py.list(assembly.instances.keys()):
    del assembly.features[instance_name]
parts = py.dict(((key.upper(), model.parts[key]) for key in model.parts.keys()))
assembly.DatumCsysByDefault(CARTESIAN)
fi = assembly.Instance(name='FRAME-1', part=parts['FRAME'], dependent=ON)
model.StaticStep(name='Centripetal', previous='Initial', nlgeom=OFF, timePeriod=1.0, initialInc=1.0, maxInc=1.0)
model.FieldOutputRequest(name='Results', createStepName='Centripetal', variables=('S', 'U', 'RF', 'E'))

def bond(label, main, secondary):
    model.Tie(label, main, secondary, positionToleranceMethod=SPECIFIED, positionTolerance=0.5, adjust=OFF, tieRotations=ON, thickness=ON)
for k in range(4):
    angle = 45.0 + 90.0 * k
    shift = rotated((pivot, pivot), k * pi / 2)
    instances = {}
    for part_name in ('BODY', 'COVER', 'PAD', 'PIN', 'LOCKS'):
        label = '%s-%d' % (part_name, k + 1)
        instances[part_name] = assembly.Instance(name=label, part=parts[part_name], dependent=ON)
        assembly.rotate(instanceList=(label,), axisPoint=(0, 0, 0), axisDirection=(0, 0, 1), angle=angle)
        assembly.translate(instanceList=(label,), vector=(shift[0], shift[1], 0))
    (b, c, r, p, l) = py.tuple((instances[n] for n in ('BODY', 'COVER', 'PAD', 'PIN', 'LOCKS')))
    bond('FramePin%d' % k, p.surfaces['LOWER'], fi.surfaces['HOLE%d' % k])
    bond('BodyPin%d' % k, p.surfaces['UPPER'], b.surfaces['EYE'])
    bond('CoverBody%d' % k, b.surfaces['SPLIT'], c.surfaces['SPLIT'])
    for group in ('LINER', 'TOP'):
        assembly.SurfaceByBoolean(name='%s%d' % (group, k), surfaces=(b.surfaces[group], c.surfaces[group]), operation=UNION)
    bond('PadBond%d' % k, assembly.surfaces['LINER%d' % k], r.surfaces['OUTER'])
    bond('LockBond%d' % k, assembly.surfaces['TOP%d' % k], l.surfaces['BOTTOM'])
    theta = angle * pi / 180.0
    rp = assembly.ReferencePoint(point=(wire_radius * cos(theta), wire_radius * sin(theta), 16.0))
    region = assembly.Set(name='LOAD%d' % k, referencePoints=(assembly.referencePoints[rp.id],))
    model.Coupling(name='Conductor%d' % k, controlPoint=region, surface=r.surfaces['BORE'], influenceRadius=WHOLE_SURFACE, couplingType=KINEMATIC, localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)
    model.ConcentratedForce(name='Force%d' % k, createStepName='Centripetal', region=region, cf1=-load_value * cos(theta), cf2=-load_value * sin(theta))
assembly.regenerate()
for (label, axis) in (('SymX', 0), ('SymY', 1)):
    labels = py.tuple((n.label for n in fi.nodes if abs(n.coordinates[axis]) < 1e-06))
    region = assembly.Set(name=label, nodes=fi.nodes.sequenceFromLabels(labels=labels))
    settings = {'u1': 0.0} if axis == 0 else {'u2': 0.0}
    model.DisplacementBC(name=label, createStepName='Initial', region=region, **settings)
anchor = fi.nodes.getByBoundingSphere(center=(0.0, 125.0, 0.0), radius=0.01)
if len(anchor) != 1:
    raise ValueError('Expected one transverse anchor node')
model.DisplacementBC(name='AnchorZ', createStepName='Initial', region=assembly.Set(name='AnchorZ', nodes=anchor), u3=0.0)
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 len(session.viewports):
    viewport = session.viewports[py.list(session.viewports.keys())[0]]
    viewport.setValues(displayedObject=assembly)
    viewport.assemblyDisplay.setValues(mesh=ON)
    viewport.view.setValues(session.views['Iso'])
    viewport.view.fitView()
print('JOB: ' + job_name)
print('ELEMENTS: %d' % py.sum((len(i.elements) for i in assembly.instances.values())))
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')