# -*- coding: utf-8 -*-
# Abaqus/CAE Python, units: mm, N, MPa
# Builds, meshes, creates Job and writes INP only.

from abaqus import mdb
from abaqusConstants import *
import math
import mesh
import regionToolset

# ------------------------------------------------------------
# Model and job
# ------------------------------------------------------------
BASE = 'SingleBranchRepairFixed'
MODEL = BASE
k = 1
while MODEL in mdb.models or MODEL + '_Job' in mdb.jobs:
    MODEL = BASE + '_%d' % k
    k += 1

JOB = MODEL + '_Job'
m = mdb.Model(name=MODEL)
a = m.rootAssembly
a.DatumCsysByDefault(CARTESIAN)

# ------------------------------------------------------------
# Parameters
# ------------------------------------------------------------
T = 6.0
RO, RI = 20.0, 17.0
RF, RB = 42.0, 32.0
RS, RH, RHEAD = 2.5, 3.0, 5.0
HEAD, NECK, STUB = 4.0, 14.0, 24.0
R, NB = 85.0, 8

PRELOAD = 10000.0
FORCE = 2500.0
MU = 0.10

BODY_SEED = 5.0
BOLT_SEED = 2.0
ELASTIC_SLIP = 0.01
EPS = 1.0e-3

ZT = T + NECK
ZOUT = ZT + R

assert RI < RO
assert RH > RS
assert RO < RB - RHEAD
assert RB + RHEAD < RF
assert R > RO
assert NB >= 4
assert ELASTIC_SLIP > 0.0

bolt_xy = []
for i in range(NB):
    ang = 2.0 * math.pi * float(i) / float(NB)
    bolt_xy.append((RB * math.cos(ang), RB * math.sin(ang)))

# ------------------------------------------------------------
# Helpers
# ------------------------------------------------------------
def make_annulus(name, ro, ri, depth, holes=False):
    s = m.ConstrainedSketch(name=name + '_sk', sheetSize=400.0)

    for rad in (ro, ri):
        s.CircleByCenterPerimeter(
            center=(0.0, 0.0),
            point1=(rad, 0.0))

    if holes:
        for x, y in bolt_xy:
            s.CircleByCenterPerimeter(
                center=(x, y),
                point1=(x + RH, y))

    p = m.Part(
        name=name,
        dimensionality=THREE_D,
        type=DEFORMABLE_BODY)

    p.BaseSolidExtrude(sketch=s, depth=depth)
    del m.sketches[s.name]
    return p


def add_instance(part, name, xyz):
    inst = a.Instance(name=name, part=part, dependent=ON)
    if xyz != (0.0, 0.0, 0.0):
        a.translate(instanceList=(name,), vector=xyz)
    return inst


def boolean_merge(name, items):
    return a.InstanceFromBooleanMerge(
        name=name,
        instances=tuple(items),
        keepIntersections=ON,
        originalInstances=DELETE,
        domain=GEOMETRY)


def face_box(inst, axis, value):
    b = [-1000.0, -1000.0, -1000.0,
         1000.0, 1000.0, 1000.0]

    b[axis] = value - EPS
    b[axis + 3] = value + EPS

    faces = inst.faces.getByBoundingBox(
        xMin=b[0], yMin=b[1], zMin=b[2],
        xMax=b[3], yMax=b[4], zMax=b[5])

    if len(faces) == 0:
        raise RuntimeError(
            'No face found on %s at coordinate %.6g'
            % (inst.name, value))

    return faces


# ------------------------------------------------------------
# Flange, neck, lower tube
# ------------------------------------------------------------
flange_p = make_annulus(
    'FlangeBlank', RF, RI, T, holes=True)

neck_p = make_annulus(
    'NeckBlank', RO, RI, NECK)

stub_p = make_annulus(
    'StubBlank', RO, RI, STUB)

# ------------------------------------------------------------
# Hollow 90-degree elbow
# ------------------------------------------------------------
s = m.ConstrainedSketch(
    name='Elbow_sk',
    sheetSize=400.0)

s.ConstructionLine(
    point1=(0.0, -160.0),
    point2=(0.0, 160.0))

for rad in (RO, RI):
    s.CircleByCenterPerimeter(
        center=(R, 0.0),
        point1=(R + rad, 0.0))

elbow_p = m.Part(
    name='ElbowBlank',
    dimensionality=THREE_D,
    type=DEFORMABLE_BODY)

elbow_p.BaseSolidRevolve(
    sketch=s,
    angle=90.0,
    flipRevolveDirection=OFF)

del m.sketches[s.name]

elbow = a.Instance(
    name='ElbowBuild',
    part=elbow_p,
    dependent=ON)

a.regenerate()

# Detect generated quarter-torus direction.
zlo = 0.0
zhi = 0.0
for v in elbow.vertices:
    zv = float(v.pointOn[0][2])
    zlo = min(zlo, zv)
    zhi = max(zhi, zv)

if zhi > RO:
    rot_axis = (0.0, 0.0, 1.0)
elif zlo < -RO:
    rot_axis = (0.0, 1.0, 0.0)
else:
    raise RuntimeError('Cannot determine elbow orientation.')

a.rotate(
    instanceList=(elbow.name,),
    axisPoint=(0.0, 0.0, 0.0),
    axisDirection=rot_axis,
    angle=180.0)

a.translate(
    instanceList=(elbow.name,),
    vector=(R, 0.0, ZT))

a.regenerate()

# ------------------------------------------------------------
# Merge aluminium bodies
# ------------------------------------------------------------
upper_items = [
    add_instance(
        flange_p, 'UpperFlangeBuild',
        (0.0, 0.0, 0.0)),
    add_instance(
        neck_p, 'NeckBuild',
        (0.0, 0.0, T)),
    elbow
]

lower_items = [
    add_instance(
        flange_p, 'LowerFlangeBuild',
        (0.0, 0.0, -T)),
    add_instance(
        stub_p, 'StubBuild',
        (0.0, 0.0, -T - STUB))
]

upper = boolean_merge('UpperBody', upper_items)
lower = boolean_merge('LowerBody', lower_items)
a.regenerate()

# ------------------------------------------------------------
# Equivalent solid bolt
# ------------------------------------------------------------
s = m.ConstrainedSketch(
    name='Bolt_sk',
    sheetSize=100.0)

s.ConstructionLine(
    point1=(0.0, -50.0),
    point2=(0.0, 50.0))

profile = [
    (0.0, -T - HEAD),
    (RHEAD, -T - HEAD),
    (RHEAD, -T),
    (RS, -T),
    (RS, T),
    (RHEAD, T),
    (RHEAD, T + HEAD),
    (0.0, T + HEAD)
]

for i in range(len(profile)):
    s.Line(
        point1=profile[i],
        point2=profile[(i + 1) % len(profile)])

bolt_p = m.Part(
    name='Fastener',
    dimensionality=THREE_D,
    type=DEFORMABLE_BODY)

bolt_p.BaseSolidRevolve(
    sketch=s,
    angle=360.0,
    flipRevolveDirection=OFF)

del m.sketches[s.name]

dp = bolt_p.DatumPlaneByPrincipalPlane(
    principalPlane=XZPLANE,
    offset=0.0)

bolt_p.PartitionCellByDatumPlane(
    cells=bolt_p.cells[:],
    datumPlane=bolt_p.datums[dp.id])

pretension_face = bolt_p.faces.findAt(
    ((0.5 * RS, 0.0, 0.0),))

bolt_p.Surface(
    name='PretensionCut',
    side1Faces=pretension_face)

# ------------------------------------------------------------
# Materials and sections
# ------------------------------------------------------------
al = m.Material(name='Aluminium')
al.Elastic(table=((70000.0, 0.30),))

steel = m.Material(name='BoltSteel')
steel.Elastic(table=((210000.0, 0.30),))
steel.Plastic(table=((800.0, 0.0),))

for section_name, material_name in (
        ('BodySection', 'Aluminium'),
        ('BoltSection', 'BoltSteel')):
    m.HomogeneousSolidSection(
        name=section_name,
        material=material_name,
        thickness=None)

# ------------------------------------------------------------
# Section assignments and mesh
# ------------------------------------------------------------
mesh_specs = (
    ('UpperBody', 'BodySection', BODY_SEED),
    ('LowerBody', 'BodySection', BODY_SEED),
    ('Fastener', 'BoltSection', BOLT_SEED)
)

elem_type = mesh.ElemType(
    elemCode=C3D10M,
    elemLibrary=STANDARD)

for part_name, section_name, size in mesh_specs:
    p = m.parts[part_name]

    p.SectionAssignment(
        region=regionToolset.Region(cells=p.cells[:]),
        sectionName=section_name)

    p.setMeshControls(
        regions=p.cells[:],
        elemShape=TET,
        technique=FREE)

    p.setElementType(
        regions=(p.cells[:],),
        elemTypes=(elem_type,))

    p.seedPart(
        size=size,
        deviationFactor=0.10,
        minSizeFactor=0.15)

    p.generateMesh()

    if len(p.elements) == 0:
        raise RuntimeError('Mesh failed: ' + part_name)

# ------------------------------------------------------------
# Bolt instances
# ------------------------------------------------------------
bolts = []

for i, xy in enumerate(bolt_xy):
    name = 'Bolt_%02d' % (i + 1)

    inst = a.Instance(
        name=name,
        part=bolt_p,
        dependent=ON)

    a.rotate(
        instanceList=(name,),
        axisPoint=(0.0, 0.0, 0.0),
        axisDirection=(1.0, 0.0, 0.0),
        angle=90.0)

    a.translate(
        instanceList=(name,),
        vector=(xy[0], xy[1], 0.0))

    bolts.append(inst)

a.regenerate()

# ------------------------------------------------------------
# Analysis steps
# ------------------------------------------------------------
m.StaticStep(
    name='Preload',
    previous='Initial',
    nlgeom=ON,
    initialInc=0.05,
    minInc=1.0e-7,
    maxInc=0.10,
    maxNumInc=300)

m.StaticStep(
    name='Service',
    previous='Preload',
    nlgeom=ON,
    initialInc=0.05,
    minInc=1.0e-7,
    maxInc=0.10,
    maxNumInc=400)

# ------------------------------------------------------------
# Contact
# Repair for the reported error:
# use ABSOLUTE_DISTANCE and do not pass fraction.
# ------------------------------------------------------------
cp = m.ContactProperty(name='Friction')

cp.NormalBehavior(
    pressureOverclosure=HARD,
    allowSeparation=ON)

cp.TangentialBehavior(
    formulation=PENALTY,
    directionality=ISOTROPIC,
    table=((MU,),),
    maximumElasticSlip=ABSOLUTE_DISTANCE,
    absoluteDistance=ELASTIC_SLIP)

gc = m.ContactStd(
    name='AssemblyContact',
    createStepName='Initial')

gc.includedPairs.setValuesInStep(
    stepName='Initial',
    useAllstar=ON)

gc.contactPropertyAssignments.appendInStep(
    stepName='Initial',
    assignments=((GLOBAL, SELF, 'Friction'),))

# ------------------------------------------------------------
# Boundary condition
# ------------------------------------------------------------
a.Set(
    name='FixedEnd',
    faces=face_box(
        lower,
        2,
        -T - STUB))

m.EncastreBC(
    name='FixLowerTube',
    createStepName='Initial',
    region=a.sets['FixedEnd'])

# ------------------------------------------------------------
# Outlet reference point, coupling and service load
# ------------------------------------------------------------
a.Surface(
    name='Outlet',
    side1Faces=face_box(
        upper,
        0,
        R))

rp = a.ReferencePoint(
    point=(R, 0.0, ZOUT))

a.Set(
    name='LoadPoint',
    referencePoints=(
        a.referencePoints[rp.id],
    ))

m.Coupling(
    name='OutletCoupling',
    controlPoint=a.sets['LoadPoint'],
    surface=a.surfaces['Outlet'],
    influenceRadius=WHOLE_SURFACE,
    couplingType=DISTRIBUTING,
    weightingMethod=UNIFORM,
    u1=ON, u2=ON, u3=ON,
    ur1=ON, ur2=ON, ur3=ON)

m.ConcentratedForce(
    name='ServiceForce',
    createStepName='Service',
    region=a.sets['LoadPoint'],
    cf1=-FORCE)

# ------------------------------------------------------------
# Bolt pretension
# ------------------------------------------------------------
for i, inst in enumerate(bolts):
    load = m.BoltLoad(
        name='Tighten_%02d' % (i + 1),
        createStepName='Preload',
        region=inst.surfaces['PretensionCut'],
        magnitude=PRELOAD,
        boltMethod=APPLY_FORCE)

    load.setValuesInStep(
        stepName='Service',
        boltMethod=FIX_LENGTH)

# ------------------------------------------------------------
# Output
# ------------------------------------------------------------
m.fieldOutputRequests['F-Output-1'].setValues(
    variables=(
        'S', 'U', 'RF', 'PEEQ',
        'CSTRESS', 'CDISP'),
    frequency=1)

m.HistoryOutputRequest(
    name='OutletHistory',
    createStepName='Preload',
    region=a.sets['LoadPoint'],
    variables=('U1', 'U2', 'U3'),
    frequency=1)

m.HistoryOutputRequest(
    name='EnergyHistory',
    createStepName='Preload',
    variables=('ALLSE', 'ALLPD', 'ALLFD', 'ALLWK'),
    frequency=1)

a.regenerate()

# ------------------------------------------------------------
# Create Abaqus job and write CAE / INP only
# ------------------------------------------------------------
job = mdb.Job(
    name=JOB,
    model=MODEL,
    type=ANALYSIS,
    numCpus=4,
    numDomains=4,
    memory=80,
    memoryUnits=PERCENTAGE)

mdb.saveAs(pathName=JOB + '.cae')
job.writeInput()

print('Model created: ' + MODEL)
print('Job created: ' + JOB)
print('CAE written: ' + JOB + '.cae')
print('INP written: ' + JOB + '.inp')
print('No analysis was submitted.')
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')