# -*- coding: ascii -*-
from __future__ import print_function
from abaqus import mdb
from abaqusConstants import *
import mesh
import math
import os
import time

# Units: mm, N, s, tonne, MPa
HEIGHT = 100.0
BOTTOM_RADIUS = 28.0
TOP_RADIUS = 40.0
WALL = 3.0
BASE = 4.0
SEED = 5.0

DROP_CLEARANCE = 3.0 * HEIGHT
CUP_SPACING = 120.0
TILT = 25.0
SPEED = 3500.0
DURATION = 0.12

E_MODULUS = 70000.0
POISSON = 0.22
DENSITY = 2.50e-9
TENSILE_STRENGTH = 45.0
FRACTURE_ENERGY = 0.010
FAIL_OPENING = 2.0 * FRACTURE_ENERGY / TENSILE_STRENGTH

tag = time.strftime('%Y%m%d_%H%M%S')
model_name = 'TwoCupDrop_' + tag
job_name = 'TwoCupDropJob_' + tag

index = 1
while (model_name in mdb.models or
       job_name in mdb.jobs or
       os.path.exists(job_name + '.inp')):
    model_name = 'TwoCupDrop_' + tag + '_' + str(index)
    job_name = 'TwoCupDropJob_' + tag + '_' + str(index)
    index += 1

model = mdb.Model(
    name=model_name,
    modelType=STANDARD_EXPLICIT
)

# ------------------------------------------------------------
# Cup geometry
# ------------------------------------------------------------
slope = (TOP_RADIUS - BOTTOM_RADIUS) / HEIGHT
radial_wall = WALL * math.sqrt(1.0 + slope * slope)
inner_base_radius = BOTTOM_RADIUS + slope * BASE - radial_wall

profile = (
    (0.0, 0.0),
    (BOTTOM_RADIUS, 0.0),
    (TOP_RADIUS, HEIGHT),
    (TOP_RADIUS - radial_wall, HEIGHT),
    (inner_base_radius, BASE),
    (0.0, BASE)
)

sketch = model.ConstrainedSketch(
    name='CupProfile',
    sheetSize=300.0
)
sketch.ConstructionLine(
    point1=(0.0, -20.0),
    point2=(0.0, 140.0)
)

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

cup = model.Part(
    name='Cup',
    dimensionality=THREE_D,
    type=DEFORMABLE_BODY
)
cup.BaseSolidRevolve(
    sketch=sketch,
    angle=360.0,
    flipRevolveDirection=OFF
)
del model.sketches['CupProfile']

# ------------------------------------------------------------
# Glass material
# ------------------------------------------------------------
glass = model.Material(name='Glass')
glass.Density(table=((DENSITY,),))
glass.Elastic(table=((E_MODULUS, POISSON),))

glass.BrittleCracking(
    type=DISPLACEMENT,
    table=(
        (TENSILE_STRENGTH, 0.0),
        (0.0, FAIL_OPENING)
    )
)

shear_strain = FAIL_OPENING / SEED
glass.brittleCracking.BrittleShear(
    type=RETENTION_FACTOR,
    table=(
        (1.0, 0.0),
        (0.10, shear_strain),
        (0.01, 10.0 * shear_strain)
    )
)

glass.brittleCracking.BrittleFailure(
    failureCriteria=UNIDIRECTIONAL,
    table=((FAIL_OPENING,),)
)

model.HomogeneousSolidSection(
    name='GlassSection',
    material='Glass',
    thickness=None
)

cup.Set(
    name='GlassCells',
    cells=cup.cells[:]
)
cup.SectionAssignment(
    region=cup.sets['GlassCells'],
    sectionName='GlassSection'
)

# ------------------------------------------------------------
# Cup mesh
# ------------------------------------------------------------
cup.setMeshControls(
    regions=cup.cells[:],
    elemShape=TET,
    technique=FREE
)

cup.setElementType(
    regions=(cup.cells[:],),
    elemTypes=(
        mesh.ElemType(
            elemCode=C3D4,
            elemLibrary=EXPLICIT,
            elemDeletion=ON
        ),
    )
)

cup.seedPart(
    size=SEED,
    deviationFactor=0.15,
    minSizeFactor=0.30
)
cup.generateMesh()

if len(cup.elements) == 0:
    raise RuntimeError('Cup mesh contains no elements.')

cup.Set(
    name='ALL_ELEMENTS',
    elements=cup.elements[:]
)
cup.Set(
    name='ALL_NODES',
    nodes=cup.nodes[:]
)

# ------------------------------------------------------------
# Rigid floor
# ------------------------------------------------------------
floor_sketch = model.ConstrainedSketch(
    name='FloorProfile',
    sheetSize=1000.0
)
floor_sketch.rectangle(
    point1=(-300.0, -300.0),
    point2=(300.0, 300.0)
)

floor = model.Part(
    name='Floor',
    dimensionality=THREE_D,
    type=DISCRETE_RIGID_SURFACE
)
floor.BaseShell(sketch=floor_sketch)
del model.sketches['FloorProfile']

rp_feature = floor.ReferencePoint(
    point=(0.0, 0.0, 0.0)
)
floor.Set(
    name='RP',
    referencePoints=(
        floor.referencePoints[rp_feature.id],
    )
)

floor.setMeshControls(
    regions=floor.faces[:],
    elemShape=QUAD,
    technique=STRUCTURED
)

floor.setElementType(
    regions=(floor.faces[:],),
    elemTypes=(
        mesh.ElemType(
            elemCode=R3D4,
            elemLibrary=EXPLICIT
        ),
        mesh.ElemType(
            elemCode=R3D3,
            elemLibrary=EXPLICIT
        )
    )
)

floor.seedPart(
    size=50.0,
    deviationFactor=0.10,
    minSizeFactor=0.10
)
floor.generateMesh()

if len(floor.elements) == 0:
    raise RuntimeError('Floor mesh contains no elements.')

floor.Set(
    name='ALL_ELEMENTS',
    elements=floor.elements[:]
)

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

floor_inst = assembly.Instance(
    name='Floor-1',
    part=floor,
    dependent=ON
)

cup_data = (
    ('Cup-1', -0.5 * CUP_SPACING, -TILT),
    ('Cup-2',  0.5 * CUP_SPACING,  TILT)
)

for name, xpos, tilt in cup_data:
    assembly.Instance(
        name=name,
        part=cup,
        dependent=ON
    )

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

    assembly.rotate(
        instanceList=(name,),
        axisPoint=(0.0, 0.0, 0.0),
        axisDirection=(0.0, 1.0, 0.0),
        angle=tilt
    )

    rim_drop = BOTTOM_RADIUS * math.sin(
        math.radians(abs(tilt))
    )

    assembly.translate(
        instanceList=(name,),
        vector=(
            xpos,
            0.0,
            DROP_CLEARANCE + rim_drop
        )
    )

assembly.regenerate()

# ------------------------------------------------------------
# Floor rigid body and boundary condition
# ------------------------------------------------------------
model.RigidBody(
    name='FixedFloorBody',
    refPointRegion=floor_inst.sets['RP'],
    bodyRegion=floor_inst.sets['ALL_ELEMENTS']
)

model.EncastreBC(
    name='FixedFloor',
    createStepName='Initial',
    region=floor_inst.sets['RP']
)

# ------------------------------------------------------------
# Explicit dynamics
# ------------------------------------------------------------
model.ExplicitDynamicsStep(
    name='Impact',
    previous='Initial',
    timePeriod=DURATION,
    nlgeom=ON
)

for name, xpos, tilt in cup_data:
    model.Velocity(
        name='InitialVelocity_' + name,
        region=assembly.instances[name].sets['ALL_NODES'],
        velocity1=0.0,
        velocity2=0.0,
        velocity3=-SPEED,
        omega=0.0
    )

model.Gravity(
    name='Gravity',
    createStepName='Impact',
    comp3=-9810.0
)

# ------------------------------------------------------------
# General contact
# Positive elastic-slip fraction fixes the previous error.
# ------------------------------------------------------------
contact_prop = model.ContactProperty('GlassContact')

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

contact_prop.TangentialBehavior(
    formulation=PENALTY,
    directionality=ISOTROPIC,
    table=((0.15,),),
    maximumElasticSlip=FRACTION,
    fraction=0.005
)

general_contact = model.ContactExp(
    name='GeneralContact',
    createStepName='Initial'
)

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

general_contact.contactPropertyAssignments.appendInStep(
    stepName='Initial',
    assignments=((GLOBAL, SELF, 'GlassContact'),)
)

# ------------------------------------------------------------
# Output
# ------------------------------------------------------------
for request_name in list(model.fieldOutputRequests.keys()):
    del model.fieldOutputRequests[request_name]

for request_name in list(model.historyOutputRequests.keys()):
    del model.historyOutputRequests[request_name]

model.FieldOutputRequest(
    name='ImpactField',
    createStepName='Impact',
    variables=('U', 'V', 'S', 'LE', 'STATUS'),
    numIntervals=60
)

model.HistoryOutputRequest(
    name='ImpactEnergy',
    createStepName='Impact',
    variables=(
        'ALLKE',
        'ALLIE',
        'ALLAE',
        'ALLVD',
        'ETOTAL'
    ),
    timeInterval=0.0002
)

# ------------------------------------------------------------
# Job: create and write input only
# ------------------------------------------------------------
job = mdb.Job(
    name=job_name,
    model=model_name,
    type=ANALYSIS,
    description='Two brittle glass cups dropping toward rigid floor',
    numCpus=4,
    numDomains=4,
    parallelizationMethodExplicit=DOMAIN,
    explicitPrecision=DOUBLE,
    nodalOutputPrecision=SINGLE
)

job.writeInput()

print('Model created: ' + model_name)
print('Job created: ' + job_name)
print('Cup elements: %d' % len(cup.elements))
print('Drop clearance: %.1f mm' % DROP_CLEARANCE)
print('Downward velocity: %.1f mm/s' % SPEED)
print('Input written: ' + os.path.abspath(job_name + '.inp'))
print('Input generation only; analysis was not 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')