from abaqus import *
import os
os.chdir('D:\\SeekSook\\Data\\bridge-data\\output\\runs\\job_no_submit\\1d0d15bfe940')

from abaqusConstants import *

import mesh

import os



base_name = 'Grand_CableStayed_Bridge_Fixed'

model_name = base_name

index = 1

while model_name in mdb.models:

    model_name = base_name + '_%02d' % index

    index += 1



model = mdb.Model(name=model_name)

part = model.Part(

    name='BridgeStructure',

    dimensionality=THREE_D,

    type=DEFORMABLE_BODY

)



for name, density, modulus, poisson in [

    ('StructuralSteel', 7.85e-9, 206000.0, 0.3),

    ('DeckConcrete', 2.50e-9, 34500.0, 0.2),

    ('CableSteel', 7.85e-9, 195000.0, 0.3)

]:

    material = model.Material(name=name)

    material.Density(table=((density,),))

    material.Elastic(table=((modulus, poisson),))



for name, a, b, t in [

    ('MainBox', 1800.0, 2400.0, 35.0),

    ('CrossBox', 700.0, 1100.0, 22.0),

    ('TowerBox', 2200.0, 2600.0, 55.0),

    ('BraceBox', 650.0, 650.0, 20.0),

    ('RailBox', 120.0, 120.0, 8.0)

]:

    model.BoxProfile(

        name=name, a=a, b=b,

        uniformThickness=ON, t1=t

    )

    model.BeamSection(

        name=name + 'Section',

        integration=DURING_ANALYSIS,

        poissonRatio=0.3,

        profile=name,

        material='StructuralSteel',

        temperatureVar=LINEAR,

        consistentMassMatrix=False

    )



model.HomogeneousShellSection(

    name='DeckSection',

    material='DeckConcrete',

    thicknessType=UNIFORM,

    thickness=250.0,

    integrationRule=SIMPSON,

    numIntPts=5

)

model.TrussSection(

    name='CableSection', material='CableSteel', area=6500.0

)



node_map = {}

groups = {}

element_count = [0]

member_keys = set()



def add_node(point):

    key = tuple(round(float(value), 6) for value in point)

    if key not in node_map:

        label = len(node_map) + 1

        part.Node(coordinates=key, label=label)

        node_map[key] = label

    return node_map[key]



def add_element(labels, shape, group):

    element_count[0] += 1

    label = element_count[0]

    node_objects = []

    for node_label in labels:

        selected = part.nodes.sequenceFromLabels(labels=(node_label,))

        if len(selected) != 1:

            raise RuntimeError('Invalid node label: %d' % node_label)

        node_objects.append(selected[0])

    part.Element(

        nodes=tuple(node_objects),

        elemShape=shape,

        label=label

    )

    groups.setdefault(group, []).append(label)



def member(start, end, group, divisions=1):

    labels = []

    for i in range(divisions + 1):

        fraction = float(i) / divisions

        point = tuple(

            start[j] + fraction * (end[j] - start[j])

            for j in range(3)

        )

        labels.append(add_node(point))

    for i in range(divisions):

        pair = tuple(sorted((labels[i], labels[i + 1])))

        if pair not in member_keys:

            member_keys.add(pair)

            add_element(labels[i:i + 2], LINE2, group)



xs = [float(-120000 + i * 10000) for i in range(25)]

ys = [-9000.0, -4500.0, 0.0, 4500.0, 9000.0]

deck_nodes = {}



for i, x in enumerate(xs):

    for j, y in enumerate(ys):

        deck_nodes[(i, j)] = add_node((x, y, 0.0))



for i in range(len(xs) - 1):

    for j in range(len(ys) - 1):

        add_element(

            (deck_nodes[(i, j)],

             deck_nodes[(i + 1, j)],

             deck_nodes[(i + 1, j + 1)],

             deck_nodes[(i, j + 1)]),

            QUAD4, 'DECK'

        )



for y in ys:

    for i in range(len(xs) - 1):

        member((xs[i], y, 0.0), (xs[i + 1], y, 0.0), 'MAIN')



for x in xs:

    for j in range(len(ys) - 1):

        member((x, ys[j], 0.0), (x, ys[j + 1], 0.0), 'CROSS')



for i in range(len(xs) - 1):

    for side in (-1.0, 1.0):

        outer = side * 9000.0

        inner = side * 4500.0

        member(

            (xs[i], outer, 0.0),

            (xs[i + 1], inner, 0.0), 'BRACE'

        )

        member(

            (xs[i], inner, 0.0),

            (xs[i + 1], outer, 0.0), 'BRACE'

        )



anchor_levels = [40500.0 + k * 1500.0 for k in range(6)]

tower_levels = sorted(set(

    [-20000.0, -10000.0, 0.0, 12000.0, 24000.0, 36000.0]

    + anchor_levels

))

base_nodes = []



for x in (-60000.0, 60000.0):

    for side in (-1.0, 1.0):

        y = side * 12000.0

        base_nodes.append(add_node((x, y, -20000.0)))

        for k in range(len(tower_levels) - 1):

            member(

                (x, y, tower_levels[k]),

                (x, y, tower_levels[k + 1]), 'TOWER'

            )

        member((x, side * 9000.0, 0.0), (x, y, 0.0), 'TOWER')



    for z in (0.0, 24000.0, 48000.0):

        member(

            (x, -12000.0, z), (x, 12000.0, z),

            'TOWER', divisions=4

        )



    for lower, upper in ((0.0, 24000.0), (24000.0, 48000.0)):

        middle = (x, 0.0, 0.5 * (lower + upper))

        for side in (-1.0, 1.0):

            member((x, side * 12000.0, lower), middle, 'BRACE')

            member(middle, (x, -side * 12000.0, upper), 'BRACE')



    for direction in (-1.0, 1.0):

        for k in range(1, 7):

            for side in (-1.0, 1.0):

                member(

                    (x, side * 12000.0, anchor_levels[k - 1]),

                    (x + direction * k * 10000.0, side * 9000.0, 0.0),

                    'CABLE'

                )



for side in (-1.0, 1.0):

    y = side * 9000.0

    for x in xs:

        member((x, y, 0.0), (x, y, 1400.0), 'RAIL')

    for i in range(len(xs) - 1):

        member(

            (xs[i], y, 1400.0),

            (xs[i + 1], y, 1400.0), 'RAIL'

        )



section_names = {

    'MAIN': 'MainBoxSection',

    'CROSS': 'CrossBoxSection',

    'TOWER': 'TowerBoxSection',

    'BRACE': 'BraceBoxSection',

    'RAIL': 'RailBoxSection',

    'DECK': 'DeckSection',

    'CABLE': 'CableSection'

}



for group, labels in groups.items():

    elements = part.elements.sequenceFromLabels(labels=tuple(labels))

    region = part.Set(name=group, elements=elements)

    part.SectionAssignment(region=region, sectionName=section_names[group])

    if group == 'DECK':

        code = S4R

    elif group == 'CABLE':

        code = T3D2

    else:

        code = B31

        orientation = (0.0, 1.0, 0.0)

        if group == 'CROSS':

            orientation = (0.0, 0.0, 1.0)

        elif group in ('TOWER', 'BRACE'):

            orientation = (1.0, 0.0, 0.0)

        part.assignBeamSectionOrientation(

            region=region, method=N1_COSINES, n1=orientation

        )

    part.setElementType(

        regions=(elements,),

        elemTypes=(mesh.ElemType(elemCode=code, elemLibrary=STANDARD),)

    )



part.Set(

    name='TOWER_BASES',

    nodes=part.nodes.sequenceFromLabels(labels=tuple(base_nodes))

)

for name, i in [('LEFT_END', 0), ('RIGHT_END', len(xs) - 1)]:

    labels = tuple(deck_nodes[(i, j)] for j in range(len(ys)))

    part.Set(

        name=name,

        nodes=part.nodes.sequenceFromLabels(labels=labels)

    )



assembly = model.rootAssembly

assembly.DatumCsysByDefault(CARTESIAN)

instance = assembly.Instance(name='Bridge-1', part=part, dependent=ON)

assembly.regenerate()

assembly.Surface(

    name='ROAD_SURFACE',

    side1Elements=instance.elements.sequenceFromLabels(

        labels=tuple(groups['DECK'])

    )

)



model.EncastreBC(

    name='TowerFoundations', createStepName='Initial',

    region=instance.sets['TOWER_BASES']

)

model.DisplacementBC(

    name='LeftBearings', createStepName='Initial',

    region=instance.sets['LEFT_END'], u1=0.0, u2=0.0, u3=0.0

)

model.DisplacementBC(

    name='RightBearings', createStepName='Initial',

    region=instance.sets['RIGHT_END'], u2=0.0, u3=0.0

)

for name, previous in [('SelfWeight', 'Initial'), ('TrafficLoad', 'SelfWeight')]:

    model.StaticStep(

        name=name, previous=previous, nlgeom=ON,

        initialInc=0.01, minInc=1.0e-8,

        maxInc=0.05, maxNumInc=1000

    )

model.Gravity(

    name='Gravity', createStepName='SelfWeight', comp3=-9810.0

)

model.Pressure(

    name='RoadLoad', createStepName='TrafficLoad',

    region=assembly.surfaces['ROAD_SURFACE'], magnitude=0.003

)

model.fieldOutputRequests['F-Output-1'].setValues(

    variables=('S', 'LE', 'U', 'RF'), numIntervals=40

)



model.keywordBlock.synchVersions(storeNodesAndElements=False)

active = False

for i, block in enumerate(model.keywordBlock.sieBlocks):

    text = block.strip().upper().replace(' ', '')

    if text.startswith('*MATERIAL,'):

        active = 'NAME=CABLESTEEL' in text.splitlines()[0]

    elif active and text.startswith('*ELASTIC'):

        model.keywordBlock.insert(i, '\n*NO COMPRESSION\n')

        break

else:

    raise RuntimeError('Cable elastic keyword not found.')



job = mdb.Job(

    name=model_name + '_Job', model=model_name,

    type=ANALYSIS, numCpus=4, numDomains=4

)

job.writeInput(consistencyChecking=ON)



if len(session.viewports):

    viewport = session.viewports[list(session.viewports.keys())[0]]

    viewport.setValues(displayedObject=assembly)

    viewport.view.setValues(session.views['Iso'])

    viewport.view.fitView()



print('Nodes: %d' % len(part.nodes))

print('Elements: %d' % len(part.elements))

print(os.path.abspath(job.name + '.inp'))

# SEEKOOK_CAE_SAVE
mdb.saveAs(pathName='D:\\SeekSook\\Data\\bridge-data\\output\\runs\\job_no_submit\\1d0d15bfe940\\abaqus_model.cae')