# -*- coding: mbcs -*-

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

# ============================================================
# Parameters
# Units: mm, N, MPa
# ============================================================
MODEL_NAME = 'Flange_Model'
JOB_NAME = 'Flange_Job'

OUTER_R = 60.0
INNER_R = 25.0
THICKNESS = 12.0

BOLT_NUM = 8
BOLT_PCD = 90.0
BOLT_HOLE_R = 4.5

E = 210000.0
NU = 0.30

PRESSURE_VALUE = 5.0
MESH_SIZE = 6.0

# ============================================================
# Clean old objects
# ============================================================
if JOB_NAME in mdb.jobs.keys():
    del mdb.jobs[JOB_NAME]

if MODEL_NAME in mdb.models.keys():
    del mdb.models[MODEL_NAME]

model = mdb.Model(name=MODEL_NAME)

# ============================================================
# Sketch
# ============================================================
sketch = model.ConstrainedSketch(
    name='FlangeSketch',
    sheetSize=200.0
)

# Outer boundary
sketch.CircleByCenterPerimeter(
    center=(0.0, 0.0),
    point1=(OUTER_R, 0.0)
)

# Central hole
sketch.CircleByCenterPerimeter(
    center=(0.0, 0.0),
    point1=(INNER_R, 0.0)
)

# Bolt holes
bolt_circle_r = BOLT_PCD / 2.0

for i in range(BOLT_NUM):
    angle = 2.0 * math.pi * float(i) / float(BOLT_NUM)
    x = bolt_circle_r * math.cos(angle)
    y = bolt_circle_r * math.sin(angle)

    sketch.CircleByCenterPerimeter(
        center=(x, y),
        point1=(x + BOLT_HOLE_R, y)
    )

# ============================================================
# Create 3D part
# ============================================================
part = model.Part(
    name='Flange',
    dimensionality=THREE_D,
    type=DEFORMABLE_BODY
)

part.BaseSolidExtrude(
    sketch=sketch,
    depth=THICKNESS
)

del model.sketches['FlangeSketch']

# ============================================================
# Material and section
# ============================================================
steel = model.Material(name='Steel')

steel.Elastic(
    table=((E, NU),)
)

model.HomogeneousSolidSection(
    name='SteelSection',
    material='Steel',
    thickness=None
)

all_cells = part.cells[:]

part.SectionAssignment(
    region=regionToolset.Region(cells=all_cells),
    sectionName='SteelSection'
)

part.Set(
    name='ALL_CELLS',
    cells=all_cells
)

# ============================================================
# Mesh
# ============================================================
part.seedPart(
    size=MESH_SIZE,
    deviationFactor=0.1,
    minSizeFactor=0.1
)

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

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

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

part.generateMesh()

# ============================================================
# Assembly
# ============================================================
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)

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

assembly.regenerate()

# ============================================================
# Static step
# ============================================================
model.StaticStep(
    name='LoadStep',
    previous='Initial',
    nlgeom=OFF,
    initialInc=0.1,
    minInc=1.0e-5,
    maxInc=0.2,
    maxNumInc=100
)

# ============================================================
# Bottom face set for boundary condition
# ============================================================
bottom_faces = instance.faces.getByBoundingBox(
    xMin=-OUTER_R - 1.0,
    xMax=OUTER_R + 1.0,
    yMin=-OUTER_R - 1.0,
    yMax=OUTER_R + 1.0,
    zMin=-0.001,
    zMax=0.001
)

assembly.Set(
    name='BOTTOM_SET',
    faces=bottom_faces
)

model.EncastreBC(
    name='FixBottom',
    createStepName='Initial',
    region=assembly.sets['BOTTOM_SET']
)

# ============================================================
# Top face surface for pressure
# Pressure requires Surface, not Region
# ============================================================
top_faces = instance.faces.getByBoundingBox(
    xMin=-OUTER_R - 1.0,
    xMax=OUTER_R + 1.0,
    yMin=-OUTER_R - 1.0,
    yMax=OUTER_R + 1.0,
    zMin=THICKNESS - 0.001,
    zMax=THICKNESS + 0.001
)

assembly.Surface(
    name='TOP_SURFACE',
    side1Faces=top_faces
)

model.Pressure(
    name='TopPressure',
    createStepName='LoadStep',
    region=assembly.surfaces['TOP_SURFACE'],
    magnitude=PRESSURE_VALUE,
    distributionType=UNIFORM
)

# ============================================================
# Output
# ============================================================
if 'F-Output-1' in model.fieldOutputRequests.keys():
    model.fieldOutputRequests['F-Output-1'].setValues(
        variables=('S', 'E', 'U', 'RF')
    )
else:
    model.FieldOutputRequest(
        name='F-Output-1',
        createStepName='LoadStep',
        variables=('S', 'E', 'U', 'RF')
    )

# ============================================================
# Create job
# ============================================================
job = mdb.Job(
    name=JOB_NAME,
    model=MODEL_NAME,
    description='3D flange static pressure analysis',
    type=ANALYSIS,
    memory=90,
    memoryUnits=PERCENTAGE,
    explicitPrecision=SINGLE,
    nodalOutputPrecision=SINGLE,
    echoPrint=OFF,
    modelPrint=OFF,
    contactPrint=OFF,
    historyPrint=OFF,
    multiprocessingMode=DEFAULT,
    numCpus=4,
    numDomains=4,
    numGPUs=0
)

# Only write INP
job.writeInput()