# -*- coding: mbcs -*-
from abaqus import *
import os
os.chdir('C:\\Users\\zengweicheng\\AppData\\Roaming\\yuanxinplatform-desktop\\bridge-data\\output\\runs\\job_no_submit\\337620b0e223')
from abaqusConstants import *
import regionToolset
import mesh

# ------------------------------------------------------------
# Model initialization
# ------------------------------------------------------------
modelName = 'CompositeBuckling3D'
jobName = 'CompositeBuckling3D'

if modelName in mdb.models.keys():
    del mdb.models[modelName]
model = mdb.Model(name=modelName)

# Units: mm, N, MPa
height = 160.0
width = 80.0
thickness = 1.89
plyThickness = 0.105

# Paper stacking sequence:
# [0_3/45/-45/0/-45/45/-45/45/-45/45/90/45/-45/90_3]T
plyAngles = [
    0.0, 0.0, 0.0,
    45.0, -45.0, 0.0,
    -45.0, 45.0, -45.0,
    45.0, -45.0, 45.0,
    90.0, 45.0, -45.0,
    90.0, 90.0, 90.0
]

# ------------------------------------------------------------
# Three-dimensional specimen
# Outer size: 80 x 160 x 1.89 mm
# Central rounded slot: 20 x 100 mm
# ------------------------------------------------------------
sketch = model.ConstrainedSketch(
    name='SpecimenSketch',
    sheetSize=300.0
)
sketch.rectangle(
    point1=(-width / 2.0, -height / 2.0),
    point2=(width / 2.0, height / 2.0)
)

# Rounded slot: overall length 100 mm, width 20 mm
slotRadius = 10.0
slotCenterOffset = 40.0

sketch.Line(
    point1=(-slotRadius, -slotCenterOffset),
    point2=(-slotRadius, slotCenterOffset)
)
sketch.Line(
    point1=(slotRadius, slotCenterOffset),
    point2=(slotRadius, -slotCenterOffset)
)
sketch.ArcByCenterEnds(
    center=(0.0, slotCenterOffset),
    point1=(-slotRadius, slotCenterOffset),
    point2=(slotRadius, slotCenterOffset),
    direction=CLOCKWISE
)
sketch.ArcByCenterEnds(
    center=(0.0, -slotCenterOffset),
    point1=(slotRadius, -slotCenterOffset),
    point2=(-slotRadius, -slotCenterOffset),
    direction=CLOCKWISE
)

part = model.Part(
    name='CompositePlate',
    dimensionality=THREE_D,
    type=DEFORMABLE_BODY
)
part.BaseSolidExtrude(sketch=sketch, depth=thickness)
del model.sketches['SpecimenSketch']

# ------------------------------------------------------------
# Sets and surfaces
# ------------------------------------------------------------
allCells = part.cells[:]
part.Set(name='ALL_CELLS', cells=allCells)

bottomFaces = part.faces.getByBoundingBox(
    xMin=-40.1, xMax=40.1,
    yMin=-80.1, yMax=-79.9,
    zMin=-0.01, zMax=thickness + 0.01
)
topFaces = part.faces.getByBoundingBox(
    xMin=-40.1, xMax=40.1,
    yMin=79.9, yMax=80.1,
    zMin=-0.01, zMax=thickness + 0.01
)
frontFaces = part.faces.getByBoundingBox(
    xMin=-40.1, xMax=40.1,
    yMin=-80.1, yMax=80.1,
    zMin=thickness - 0.01, zMax=thickness + 0.01
)

part.Set(name='BOTTOM_END', faces=bottomFaces)
part.Set(name='TOP_END', faces=topFaces)
part.Surface(name='FRONT_SURFACE', side1Faces=frontFaces)

# Material orientation: local 1-axis parallel to specimen height
datum = part.DatumCsysByThreePoints(
    name='MaterialCSYS',
    coordSysType=CARTESIAN,
    origin=(0.0, 0.0, 0.0),
    point1=(0.0, 1.0, 0.0),
    point2=(1.0, 0.0, 0.0)
)
materialCsys = part.datums[datum.id]

# ------------------------------------------------------------
# Carbon/epoxy properties from the paper
# ------------------------------------------------------------
material = model.Material(name='CarbonEpoxy')

material.Elastic(
    type=ENGINEERING_CONSTANTS,
    table=((
        143528.5, 5826.3, 5826.3,
        0.36, 0.36, 0.40,
        3845.5, 3845.5, 2080.8
    ),)
)

# Representative density for carbon/epoxy in tonne/mm^3
material.Density(table=((1.55e-9,),))

# Hashin strengths from Tables 1 and 2
hashin = material.HashinDamageInitiation(
    table=((
        2220.7, 641.0,
        49.0, 114.0,
        83.5, 57.0
    ),),
    alpha=1.0
)

# Paper fracture energies: Gft, Gfc, Gmt, Gmc
hashin.DamageEvolution(
    type=ENERGY,
    table=((133.0, 10.0, 0.5, 1.6),)
)

hashin.DamageStabilization(
    fiberTensileCoeff=0.0005,
    fiberCompressiveCoeff=0.0005,
    matrixTensileCoeff=0.0005,
    matrixCompressiveCoeff=0.0005
)

# ------------------------------------------------------------
# Eighteen-ply continuum-shell layup
# ------------------------------------------------------------
layup = part.CompositeLayup(
    name='PaperLayup18',
    description='18-ply layup reproduced from Composite Structures 273',
    elementType=CONTINUUM_SHELL,
    symmetric=OFF
)
layup.Section(
    preIntegrate=OFF,
    integrationRule=SIMPSON,
    thicknessType=UNIFORM,
    poissonDefinition=DEFAULT,
    temperature=GRADIENT,
    useDensity=OFF
)

layup.ReferenceOrientation(
    orientationType=SYSTEM,
    localCsys=materialCsys,
    fieldName='',
    additionalRotationType=ROTATION_NONE,
    angle=0.0,
    axis=AXIS_3,
    stackDirection=STACK_3
)

cellRegion = regionToolset.Region(cells=allCells)

for index, plyAngle in enumerate(plyAngles):
    layup.CompositePly(
        suppressed=False,
        plyName='Ply-%02d' % (index + 1),
        region=cellRegion,
        material='CarbonEpoxy',
        thicknessType=SPECIFY_THICKNESS,
        thickness=plyThickness,
        orientationType=SPECIFY_ORIENT,
        orientationValue=plyAngle,
        additionalRotationType=ROTATION_NONE,
        additionalRotationField='',
        axis=AXIS_3,
        angle=0.0,
        numIntPoints=3
    )

# ------------------------------------------------------------
# Mesh: reduced model size for fast calculation
# SC8R is the continuum-shell element used in the paper
# ------------------------------------------------------------
elemType = mesh.ElemType(
    elemCode=SC8R,
    elemLibrary=STANDARD,
    secondOrderAccuracy=OFF,
    hourglassControl=DEFAULT
)
part.setElementType(
    regions=(allCells,),
    elemTypes=(elemType,)
)
part.seedPart(
    size=5.0,
    deviationFactor=0.1,
    minSizeFactor=0.2
)
part.generateMesh()

# ------------------------------------------------------------
# Assembly
# ------------------------------------------------------------
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)
instance = assembly.Instance(
    name='CompositePlate-1',
    part=part,
    dependent=ON
)
assembly.regenerate()

# ------------------------------------------------------------
# Nonlinear compression step
# Paper test cross-head speed: 2 mm/min
# Total imposed shortening selected as 2 mm
# ------------------------------------------------------------
model.StaticStep(
    name='AxialCompression',
    previous='Initial',
    description='Geometrically nonlinear axial compression',
    nlgeom=ON,
    timePeriod=1.0,
    initialInc=0.002,
    minInc=1.0e-7,
    maxInc=0.02,
    maxNumInc=1000,
    stabilizationMagnitude=2.0e-4,
    stabilizationMethod=DISSIPATED_ENERGY_FRACTION,
    continueDampingFactors=False,
    adaptiveDampingRatio=0.05
)

# Articulated-test approximation:
# bottom fixed in translations; top moves axially and is laterally restrained
model.DisplacementBC(
    name='BottomSupport',
    createStepName='Initial',
    region=instance.sets['BOTTOM_END'],
    u1=0.0,
    u2=0.0,
    u3=0.0,
    amplitude=UNSET,
    distributionType=UNIFORM,
    fieldName='',
    localCsys=None
)

model.DisplacementBC(
    name='TopGuide',
    createStepName='Initial',
    region=instance.sets['TOP_END'],
    u1=0.0,
    u3=0.0,
    amplitude=UNSET,
    distributionType=UNIFORM,
    fieldName='',
    localCsys=None
)

model.DisplacementBC(
    name='Compression',
    createStepName='AxialCompression',
    region=instance.sets['TOP_END'],
    u2=-2.0,
    amplitude=UNSET,
    fixed=OFF,
    distributionType=UNIFORM,
    fieldName='',
    localCsys=None
)

# Small transverse pressure triggers the experimentally observed
# single-half-wave buckling mode without a separate eigenvalue job.
model.Pressure(
    name='BucklingTrigger',
    createStepName='AxialCompression',
    region=instance.surfaces['FRONT_SURFACE'],
    magnitude=1.0e-4,
    distributionType=UNIFORM,
    field=''
)

# ------------------------------------------------------------
# Outputs
# ------------------------------------------------------------
model.fieldOutputRequests['F-Output-1'].setValues(
    variables=(
        'S', 'E', 'LE', 'U', 'RF',
        'HSNFTCRT', 'HSNFCCRT',
        'HSNMTCRT', 'HSNMCCRT',
        'DAMAGEFT', 'DAMAGEFC',
        'DAMAGEMT', 'DAMAGEMC',
        'DAMAGESHR', 'STATUS'
    ),
    frequency=5
)

model.HistoryOutputRequest(
    name='TopReaction',
    createStepName='AxialCompression',
    variables=('U2', 'RF2'),
    region=instance.sets['TOP_END'],
    sectionPoints=DEFAULT,
    rebar=EXCLUDE,
    frequency=1
)

# ------------------------------------------------------------
# Job creation and INP generation
# ------------------------------------------------------------
if jobName in mdb.jobs.keys():
    del mdb.jobs[jobName]

job = mdb.Job(
    name=jobName,
    model=modelName,
    description='3D composite plate buckling based on the supplied paper',
    type=ANALYSIS,
    memory=90,
    memoryUnits=PERCENTAGE,
    explicitPrecision=SINGLE,
    nodalOutputPrecision=SINGLE,
    echoPrint=OFF,
    modelPrint=OFF,
    contactPrint=OFF,
    historyPrint=OFF,
    userSubroutine='',
    scratch='',
    resultsFormat=ODB,
    multiprocessingMode=DEFAULT,
    numCpus=4,
    numDomains=4,
    numGPUs=0
)

job.writeInput(consistencyChecking=OFF)

# SEEKOOK_CAE_SAVE
mdb.saveAs(pathName='C:\\Users\\zengweicheng\\AppData\\Roaming\\yuanxinplatform-desktop\\bridge-data\\output\\runs\\job_no_submit\\337620b0e223\\abaqus_model.cae')