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 part, material, section, assembly, step, interaction, load, mesh, job
import math, os, time
NAME = 'WingAssemblyMerged'
i = 1
while NAME in mdb.models.keys() or NAME in mdb.jobs.keys():
    NAME = 'WingAssemblyMerged_%02d' % i
    i += 1
OUT = os.path.abspath(NAME + '_' + time.strftime('%Y%m%d_%H%M%S'))
if not os.path.isdir(OUT):
    os.makedirs(OUT)
os.chdir(OUT)
(L, FULL_HALF, SUPPORT, COLLAR) = (304.9, 457.2, 292.1, 12.7)
(SIZE, DIV, NZ, PLY) = (6.0, 2, 4, 0.122)
BASE_Y = [0.0, 0.04, 0.125, 0.25, 0.4, 0.55, 0.7, 0.85, 1.0]
BASE_H = [0.14, 0.73, 1.0, 0.87, 0.7, 0.54, 0.37, 0.21, 0.045]
PLIES = {'FU': [0] * 4 + [-45] + [0] * 4 + [45, 90, -45] + [0] * 4 + [45] + [0] * 4 + [-45] + [0] * 4 + [45] + [0] * 3 + [-45, 90, 45] + [0] * 3, 'FL': [0] * 4 + [-45] + [0] * 4 + [45] + [0] * 3 + [-45, 90, 45] + [0] * 3, 'CU': [45, 0, -45, 90, 45] + [0] * 3, 'CL': [45, -45, 90, 45] + [0] * 3, 'LE': [-45] * 2 + [0] * 2 + [45] * 2 + [0] * 2 + [-45] * 2 + [90] + [45] * 2 + [0] * 3, 'W1': [-45, 0, 45, 0, 45, 0, 45, 0, -45] + [-45] * 2 + [45] * 2, 'W2': [45, -45, 90, 45] + [0] * 3, 'SK': [0, -45, 0, 45, -45]}
(xs, ys, hs) = ([], [], [])
cuts = [0.0, COLLAR, SUPPORT, L]
for (left, right) in zip(cuts[:-1], cuts[1:]):
    n = int(math.ceil((right - left) / SIZE))
    xs.extend((left + (right - left) * k / n for k in range(n)))
xs.append(L)
for j in range(len(BASE_Y) - 1):
    for k in range(DIV):
        t = float(k) / DIV
        ys.append(BASE_Y[j] + t * (BASE_Y[j + 1] - BASE_Y[j]))
        hs.append(BASE_H[j] + t * (BASE_H[j + 1] - BASE_H[j]))
ys.append(1.0)
hs.append(BASE_H[-1])
(nx, ny, nz) = (len(xs), len(ys), NZ + 1)
webs = [min(range(ny), key=lambda j: abs(ys[j] - v)) for v in (0.25, 0.7)]

def node(i, j, k):
    return 1 + (i * ny + j) * nz + k
xyz = {}
for (i, x) in enumerate(xs):
    (c, h) = (101.6 - 50.8 * x / FULL_HALF, 25.4 - 12.7 * x / FULL_HALF)
    for (j, y) in enumerate(ys):
        for k in range(nz):
            xyz[node(i, j, k)] = (x, c * (y - 0.25), 0.5 * h * hs[j] * (-1.0 + 2.0 * k / NZ))
(shell, solid, groups) = ([], [], {})
foam_faces = {1: [], 2: [], 3: [], 5: []}
(skin_ids, web_ids, bottom_ids) = ([], [], [])

def shell_element(conn, key, side):
    label = len(shell) + 1
    shell.append(conn)
    groups.setdefault(key + '_' + side, []).append(label)
    (web_ids if key in ('W1', 'W2') else skin_ids).append(label)
    if side == 'B':
        bottom_ids.append(label)
for i in range(nx - 1):
    for j in range(ny - 1):
        for k in range(NZ):
            solid.append((node(i, j, k), node(i + 1, j, k), node(i + 1, j + 1, k), node(i, j + 1, k), node(i, j, k + 1), node(i + 1, j, k + 1), node(i + 1, j + 1, k + 1), node(i, j + 1, k + 1)))
            e = len(solid)
            for (face, flag) in ((1, k == 0), (2, k == NZ - 1), (3, j == 0), (5, j == ny - 2 or j in [w - 1 for w in webs])):
                if flag:
                    foam_faces[face].append(e)
        mid = 0.5 * (ys[j] + ys[j + 1])
        for (k, side) in ((NZ, 'T'), (0, 'B')):
            key = 'LE' if mid < 0.125 else ('FU' if side == 'T' else 'FL') if mid < 0.25 else ('CU' if side == 'T' else 'CL') if mid < 0.7 else 'SK'
            conn = (node(i, j, k), node(i + 1, j, k), node(i + 1, j + 1, k), node(i, j + 1, k))
            shell_element(conn if side == 'T' else conn[::-1], key, side)
    for (j, key, side) in ((0, 'LE', 'V'), (ny - 1, 'SK', 'R'), (webs[0], 'W1', 'V'), (webs[1], 'W2', 'V')):
        for k in range(NZ):
            conn = (node(i, j, k), node(i + 1, j, k), node(i + 1, j, k + 1), node(i, j, k + 1))
            shell_element(conn if side == 'V' else conn[::-1], key, side)
m = mdb.Model(name=NAME)
for (name, data) in (('Foam', (92.0, 0.3)), ('Steel', (210000.0, 0.3))):
    mat = m.Material(name=name)
    mat.Elastic(table=(data,))
    m.HomogeneousSolidSection(name=name, material=name)
mat = m.Material(name='USN125B')
mat.Elastic(type=LAMINA, table=((140000.0, 7100.0, 0.3, 2500.0, 2500.0, 2500.0),))
mat.HashinDamageInitiation(table=((1384.7, 710.6, 22.3, 95.0, 62.52, 62.52),))
(unique_layups, section_map) = ({}, {})
for key in sorted(PLIES):
    signature = ('USN125B', PLY, tuple(PLIES[key]), 3)
    unique_layups.setdefault(signature, []).append(key)
for (signature, keys) in sorted(unique_layups.items()):
    (material_name, ply_thickness, angles, npoints) = signature
    section_name = 'CF_' + '_'.join(keys)
    layers = tuple((section.SectionLayer(material=material_name, thickness=ply_thickness, orientAngle=float(v), numIntPts=npoints, plyName=section_name + '_P%02d' % (j + 1)) for (j, v) in enumerate(angles)))
    m.CompositeShellSection(name=section_name, layup=layers, preIntegrate=OFF, symmetric=OFF, integrationRule=SIMPSON)
    for key in keys:
        section_map[key] = section_name
assert section_map['CL'] == section_map['W2']

def make_part(name, points, conn, code):
    p = m.Part(name=name, dimensionality=THREE_D, type=DEFORMABLE_BODY)
    used = sorted(set((n for element in conn for n in element)))
    for n in used:
        p.Node(label=n, coordinates=points[n])
    index = dict(((n, j) for (j, n) in enumerate(used)))
    for (j, element) in enumerate(conn):
        p.Element(label=j + 1, elemShape=QUAD4 if code == S4R else HEX8, nodes=tuple((p.nodes[index[n]] for n in element)))
    p.Set(name='ALL', elements=p.elements[:])
    p.setElementType(regions=p.sets['ALL'], elemTypes=(mesh.ElemType(elemCode=code, elemLibrary=STANDARD),))
    return p
wing = make_part('WingShellAndSpars', xyz, shell, S4R)
foam = make_part('FoamCore', xyz, solid, C3D8R)
foam.SectionAssignment(region=foam.sets['ALL'], sectionName='Foam')
systems = {}
for (side, v) in (('T', (0.0, 1.0, 0.0)), ('B', (0.0, -1.0, 0.0)), ('V', (0.0, 0.0, 1.0)), ('R', (0.0, 0.0, -1.0))):
    d = wing.DatumCsysByThreePoints(name='CS_' + side, coordSysType=CARTESIAN, origin=(0.0, 0.0, 0.0), point1=(1.0, 0.0, 0.0), point2=v)
    systems[side] = wing.datums[d.id]
(section_regions, orientation_regions) = ({}, {})
for (name, labels) in sorted(groups.items()):
    (key, side) = name.split('_')
    section_regions.setdefault(section_map[key], []).extend(labels)
    orientation_regions.setdefault(side, []).extend(labels)
for (name, labels) in sorted(section_regions.items()):
    region = wing.Set(name='REG_' + name, elements=wing.elements.sequenceFromLabels(tuple(sorted(labels))))
    wing.SectionAssignment(region=region, sectionName=name, offset=0.0, offsetType=MIDDLE_SURFACE, thicknessAssignment=FROM_SECTION)
for (side, labels) in sorted(orientation_regions.items()):
    region = wing.Set(name='DIR_' + side, elements=wing.elements.sequenceFromLabels(tuple(sorted(labels))))
    wing.MaterialOrientation(region=region, orientationType=SYSTEM, axis=AXIS_3, localCsys=systems[side], stackDirection=STACK_3)
wing.Surface(name='BOND', side2Elements=wing.elements.sequenceFromLabels(tuple(skin_ids)), side1Elements=wing.elements.sequenceFromLabels(tuple(web_ids)))
wing.Surface(name='LOWER', side1Elements=wing.elements.sequenceFromLabels(tuple(bottom_ids)))
foam.Surface(name='BOND', **dict((('face%dElements' % f, foam.elements.sequenceFromLabels(tuple(labels))) for (f, labels) in foam_faces.items())))
collar_nodes = tuple((n.label for n in wing.nodes if n.coordinates[0] <= COLLAR + 1e-08))
wing.Set(name='COLLAR', nodes=wing.nodes.sequenceFromLabels(collar_nodes))
free_root = tuple((node(0, j, k) for j in range(1, ny - 1) if j not in webs for k in range(1, NZ)))
foam.Set(name='ROOT_INTERIOR', nodes=foam.nodes.sequenceFromLabels(free_root))
a = m.rootAssembly
a.DatumCsysByDefault(CARTESIAN)
wi = a.Instance(name='Wing-1', part=wing, dependent=ON)
fi = a.Instance(name='Foam-1', part=foam, dependent=ON)

def pair_call(method, main, secondary, **kwargs):
    try:
        return method(main=main, secondary=secondary, **kwargs)
    except TypeError:
        return method(master=main, slave=secondary, **kwargs)
pair_call(m.Tie, wi.surfaces['BOND'], fi.surfaces['BOND'], name='FoamBond', adjust=OFF, thickness=OFF, tieRotations=OFF, positionToleranceMethod=SPECIFIED, positionTolerance=0.001, constraintEnforcement=NODE_TO_SURFACE)
m.DisplacementBC(name='CoreSymmetry', createStepName='Initial', region=fi.sets['ROOT_INTERIOR'], u1=0.0)

def fixture(name, boxes, rp_point):
    (points, conn) = ({}, [])
    for (x0, x1, y0, y1, z0, z1) in boxes:
        start = len(points) + 1
        vertices = [(x0, y0, z0), (x1, y0, z0), (x1, y1, z0), (x0, y1, z0), (x0, y0, z1), (x1, y0, z1), (x1, y1, z1), (x0, y1, z1)]
        for (j, v) in enumerate(vertices):
            points[start + j] = v
        conn.append(tuple(range(start, start + 8)))
    p = make_part(name, points, conn, C3D8R)
    p.SectionAssignment(region=p.sets['ALL'], sectionName='Steel')
    p.Surface(name='TOP', face2Elements=p.elements[:])
    instance = a.Instance(name=name + '-1', part=p, dependent=ON)
    rp = a.ReferencePoint(point=rp_point)
    region = a.Set(name=name + '_RP', referencePoints=(a.referencePoints[rp.id],))
    m.RigidBody(name=name + '_Rigid', refPointRegion=region, bodyRegion=instance.sets['ALL'])
    return (instance, region)
collar_boxes = [(0.0, COLLAR, -35.0, 86.0, 16.0, 24.0), (0.0, COLLAR, -35.0, 86.0, -24.0, -16.0), (0.0, COLLAR, -35.0, -27.0, -16.0, 16.0), (0.0, COLLAR, 78.0, 86.0, -16.0, 16.0)]
(ci, cr) = fixture('LoadCollar', collar_boxes, (0.0, 25.4, 0.0))
lower_t = max((len(PLIES[k]) * PLY for k in ('LE', 'FL', 'CL', 'SK')))
support_z = min((v[2] for v in xyz.values() if SUPPORT - 12.7 - 1e-08 <= v[0] <= SUPPORT + 12.7 + 1e-08)) - 0.5 * lower_t - 0.01
(si, sr) = fixture('JigSupport', [(SUPPORT - 12.7, SUPPORT + 12.7, -40.0, 100.0, support_z - 12.0, support_z)], (SUPPORT, 30.0, support_z - 6.0))
m.Coupling(name='CollarGrip', controlPoint=cr, surface=wi.sets['COLLAR'], couplingType=KINEMATIC, influenceRadius=WHOLE_SURFACE, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)
m.DisplacementBC(name='LoadGuide', createStepName='Initial', region=cr, u1=0.0, u2=0.0, u3=0.0, ur1=0.0, ur2=0.0, ur3=0.0)
m.EncastreBC(name='FixedJig', createStepName='Initial', region=sr)
prop = m.ContactProperty('FrictionlessSupport')
prop.NormalBehavior(pressureOverclosure=HARD, allowSeparation=ON, constraintEnforcementMethod=DEFAULT)
pair_call(m.SurfaceToSurfaceContactStd, si.surfaces['TOP'], wi.surfaces['LOWER'], name='SupportContact', createStepName='Initial', sliding=FINITE, interactionProperty='FrictionlessSupport', thickness=ON)
m.StaticStep(name='Bending', previous='Initial', nlgeom=ON, timePeriod=1.0, initialInc=0.001, maxInc=0.02, minInc=1e-08, maxNumInc=1000)
m.boundaryConditions['LoadGuide'].setValuesInStep(stepName='Bending', u3=-50.8)
for repo in (m.fieldOutputRequests, m.historyOutputRequests):
    for key in list(repo.keys()):
        del repo[key]
m.FieldOutputRequest(name='Motion', createStepName='Bending', variables=('U', 'RF'), frequency=1)
m.FieldOutputRequest(name='Ply', createStepName='Bending', region=wi.sets['ALL'], variables=('S', 'LE', 'HSNFTCRT', 'HSNFCCRT', 'HSNMTCRT', 'HSNMCCRT'), frequency=1)
m.FieldOutputRequest(name='Core', createStepName='Bending', region=fi.sets['ALL'], variables=('S', 'LE'), frequency=1)
m.FieldOutputRequest(name='Contact', createStepName='Bending', variables=('CSTRESS', 'CDISP'), frequency=1)
for (name, region) in (('Load', cr), ('Support', sr)):
    m.HistoryOutputRequest(name=name, createStepName='Bending', region=region, variables=('U3', 'RF3'), frequency=1)
m.HistoryOutputRequest(name='Energy', createStepName='Bending', variables=('ALLSE', 'ALLAE', 'ALLWK'), frequency=1)
a.regenerate()
m.keywordBlock.synchVersions(storeNodesAndElements=False)
for (i, block) in enumerate(m.keywordBlock.sieBlocks):
    lines = block.splitlines()
    changed = False
    for (j, line) in enumerate(lines):
        if line.strip().upper().startswith('*ELEMENT OUTPUT'):
            if 'ALLSECTIONPTS' not in line.upper():
                lines[j] = line.rstrip().rstrip(',') + ', ALLSECTIONPTS'
                changed = True
    if changed:
        m.keywordBlock.replace(i, '\n'.join(lines))
assert len(m.materials) == 3 and len(m.sections) == 9
job = mdb.Job(name=NAME, model=NAME, type=ANALYSIS, numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE, description='Assembly with deduplicated layups')
_sw_record_write(job, consistencyChecking=ON)
mdb.saveAs(pathName=os.path.join(OUT, NAME + '.cae'))
if len(session.viewports.keys()):
    vp = session.viewports[list(session.viewports.keys())[0]]
    vp.setValues(displayedObject=a)
    vp.view.setValues(cameraPosition=(420.0, -500.0, 320.0), cameraTarget=(152.0, 25.0, 0.0), cameraUpVector=(0.0, 0.0, 1.0))
    vp.view.fitView()
print('PARTS: WingShellAndSpars, FoamCore, LoadCollar, JigSupport')
print('MATERIALS: 3; COMPOSITE SECTIONS: 7; SOLID SECTIONS: 2')
print('SHARED SECTION: CL + W2 -> ' + section_map['CL'])
print('FULL_WING_FORCE = 2 * abs(LoadCollar_RP.RF3)')
print('OUTPUT: ' + OUT)
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')