from abaqus import mdb
from abaqusConstants import OFF
import os as _sw_os
_sw_before_jobs = {str(n): mdb.jobs[n] for n in mdb.jobs.keys()}
_sw_written_jobs = []
def _sw_record_write(_job, *args, **kwargs):
    _answer = _job.writeInput(*args, **kwargs)
    _name = str(getattr(_job, 'name', ''))
    if not _name or _name not in mdb.jobs.keys():
        _names = [str(n) for n in mdb.jobs.keys() if mdb.jobs[n] is _job]
        if len(_names) != 1:
            raise RuntimeError('Cannot identify the job which wrote the INP')
        _name = _names[0]
    _sw_written_jobs.append((_name, _sw_os.path.abspath(_name + '.inp')))
    return _answer

from abaqus import *
from abaqusConstants import *
import regionToolset
source = 'Buried_Tunnel_XFEM_Stable'
base = 'Buried_Tunnel_XFEM_Dynamic_Repair'
if source not in mdb.models.keys():
    raise RuntimeError('Source model missing: ' + source)
name = base
index = 1
while name in mdb.models.keys() or name in mdb.jobs.keys():
    name = '%s_%d' % (base, index)
    index += 1
model = mdb.Model(name=name, objectToCopy=mdb.models[source])
a = model.rootAssembly
for inst_name in ('Lining-1', 'Soil-1'):
    if inst_name not in a.instances.keys():
        raise RuntimeError('Missing instance: ' + inst_name)
    if len(a.instances[inst_name].elements) == 0:
        raise RuntimeError('Empty mesh: ' + inst_name)
lining = a.instances['Lining-1']
soil = a.instances['Soil-1']
for elem in lining.elements:
    if len(elem.connectivity) != 8:
        raise RuntimeError('Lining requires HEX8 elements.')
if 'Soil_Lining_Tie' not in model.constraints.keys():
    raise RuntimeError('Missing soil-lining Tie.')
if 'Crown_XFEM' not in a.engineeringFeatures.cracks.keys():
    raise RuntimeError('Missing XFEM crack.')
for material_name in ('Concrete_XFEM', 'Elastic_Soil'):
    if material_name not in model.materials.keys():
        raise RuntimeError('Missing material: ' + material_name)
for section_name in ('Lining_Section', 'Soil_Section'):
    if section_name not in model.sections.keys():
        raise RuntimeError('Missing section: ' + section_name)
for bc_name in ('Left_Arch_Opening', 'Right_Arch_Opening', 'Left_Support', 'Right_Support'):
    if bc_name in model.boundaryConditions.keys():
        del model.boundaryConditions[bc_name]
for bc_name in ('Soil_Bottom', 'Soil_Lateral', 'Axial_Stabilizer'):
    if bc_name not in model.boundaryConditions.keys():
        raise RuntimeError('Missing soil BC: ' + bc_name)
concrete = model.materials['Concrete_XFEM']
concrete.maxpsDamageInitiation.setValues(table=((2000000.0,),))
concrete.maxpsDamageInitiation.DamageEvolution(type=ENERGY, softening=LINEAR, table=((100.0,),))
concrete.maxpsDamageInitiation.DamageStabilizationCohesive(cohesiveCoeff=0.0001)
initial_step = 'Initial_Equilibrium'
preload_step = 'Surface_Load_Crack_Growth'
fracture_step = 'Fracture_QuasiStatic'
for step_name in (initial_step, preload_step):
    if step_name not in model.steps.keys():
        raise RuntimeError('Missing step: ' + step_name)
model.steps[initial_step].setValues(timePeriod=1.0, initialInc=0.25, minInc=1e-07, maxInc=0.5, maxNumInc=30, nlgeom=OFF)
if 'Load_Ramp_Stable' not in model.amplitudes.keys():
    model.SmoothStepAmplitude(name='Load_Ramp_Stable', timeSpan=STEP, data=((0.0, 0.0), (1.0, 1.0)))
preload_ground = 20000.0
preload_tunnel = 40000.0
for (load_name, magnitude) in (('Surface_Surcharge', preload_ground), ('Tunnel_Internal_Pressure', preload_tunnel)):
    if load_name not in model.loads.keys():
        raise RuntimeError('Missing load: ' + load_name)
    model.loads[load_name].setValues(magnitude=magnitude, amplitude='Load_Ramp_Stable')
model.steps[preload_step].setValues(timePeriod=1.0, initialInc=0.02, minInc=1e-08, maxInc=0.1, maxNumInc=200, nlgeom=OFF, matrixStorage=UNSYMMETRIC, stabilizationMethod=DISSIPATED_ENERGY_FRACTION, stabilizationMagnitude=0.0002)
growth = model.XFEMCrackGrowth(name='XFEM_Growth_Control', createStepName=initial_step, crackName='Crown_XFEM', allowGrowth=OFF)
model.ImplicitDynamicsStep(name=fracture_step, previous=preload_step, timePeriod=0.5, application=QUASI_STATIC, timeIncrementationMethod=AUTOMATIC, initialInc=0.001, minInc=1e-07, maxInc=0.01, maxNumInc=600, nlgeom=OFF, matrixStorage=UNSYMMETRIC)
growth.setValuesInStep(stepName=fracture_step, allowGrowth=ON)
model.SmoothStepAmplitude(name='Ground_Fracture_Ramp', timeSpan=STEP, data=((0.0, 1.0), (0.25, 1.25), (0.5, 1.75), (0.75, 2.4), (1.0, 3.0)))
model.SmoothStepAmplitude(name='Tunnel_Fracture_Ramp', timeSpan=STEP, data=((0.0, 1.0), (0.25, 1.4), (0.5, 2.3), (0.75, 3.6), (1.0, 5.0)))
model.loads['Surface_Surcharge'].setValuesInStep(stepName=fracture_step, magnitude=preload_ground, amplitude='Ground_Fracture_Ramp')
model.loads['Tunnel_Internal_Pressure'].setValuesInStep(stepName=fracture_step, magnitude=preload_tunnel, amplitude='Tunnel_Fracture_Ramp')
for key in list(model.fieldOutputRequests.keys()):
    del model.fieldOutputRequests[key]
for key in list(model.historyOutputRequests.keys()):
    del model.historyOutputRequests[key]
model.FieldOutputRequest(name='Initial_Output', createStepName=initial_step, variables=('S', 'U', 'RF'), numIntervals=2, timeMarks=ON)
model.fieldOutputRequests['Initial_Output'].deactivate(preload_step)
model.FieldOutputRequest(name='Preload_Output', createStepName=preload_step, variables=('S', 'U', 'RF'), numIntervals=7, timeMarks=ON)
model.fieldOutputRequests['Preload_Output'].deactivate(fracture_step)
model.FieldOutputRequest(name='Fracture_XFEM_Output', createStepName=fracture_step, variables=('S', 'U', 'RF', 'PHILSM', 'STATUSXFEM'), numIntervals=28, timeMarks=ON)
model.HistoryOutputRequest(name='Fracture_Energy_Check', createStepName=fracture_step, variables=('ALLIE', 'ALLSE', 'ALLKE', 'ALLWK', 'ALLVD'), frequency=1)
model.HistoryOutputRequest(name='Static_Energy_Check', createStepName=preload_step, variables=('ALLIE', 'ALLSE', 'ALLWK', 'ALLSD', 'ALLVD'), frequency=1)
a.regenerate()
model.setValues(description='3D buried tunnel XFEM convergence-repair demo. Original soil, lining, hex mesh, tie and crack retained. No conflicting BC on lining Tie secondary nodes. XFEM growth disabled during static preload. XFEM growth activated during implicit dynamic step. Initial ground pressure 20 kPa and inner pressure 40 kPa. Final ground pressure 60 kPa and inner pressure 200 kPa. Concrete E=30 GPa, tensile strength=2 MPa, fracture energy=100 N/m. No gravity or geostatic initial stress. This is a demonstration, not experimental validation. Check kinetic energy, viscous energy and crack growth.')
job = mdb.Job(name=name, model=name, type=ANALYSIS, numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE)
_sw_record_write(job, consistencyChecking=ON)
try:
    vp = session.viewports[session.currentViewportName]
    vp.setValues(displayedObject=a)
    vp.assemblyDisplay.setValues(mesh=ON, loads=OFF, bcs=OFF, interactions=OFF)
    vp.view.setValues(session.views['Iso'])
    vp.view.fitView()
except Exception as exc:
    print('Viewport warning: %s' % str(exc))
print('MODEL: %s' % name)
print('JOB: %s' % job.name)
print('SOIL ELEMENTS: %d' % len(soil.elements))
print('LINING ELEMENTS: %d' % len(lining.elements))
print('STATIC PRELOAD: XFEM GROWTH OFF')
print('DYNAMIC FRACTURE: XFEM GROWTH ON')
print('FINAL GROUND PRESSURE: 60000 Pa')
print('FINAL TUNNEL PRESSURE: 200000 Pa')
print('REQUESTED OUTPUT: APPROXIMATELY 40 FRAMES')
print('XFEM FIELDS: PHILSM, STATUSXFEM')
print('INPUT WRITTEN; JOB NOT SUBMITTED')
print('VERIFY REAL CRACK GROWTH AND ENERGY BALANCE')
print('SOUTHWEITION_JOB_PREPARED')

if not _sw_written_jobs:
    _sw_candidates = [n for n in mdb.jobs.keys() if n not in _sw_before_jobs or mdb.jobs[n] is not _sw_before_jobs[n]]
    if len(_sw_candidates) != 1:
        raise RuntimeError('Cannot identify this script job: call writeInput() on exactly the intended job')
    _sw_record_write(mdb.jobs[_sw_candidates[0]], consistencyChecking=OFF)
_sw_name, _sw_inp = _sw_written_jobs[-1]
if not _sw_os.path.isfile(_sw_inp):
    raise RuntimeError('This script did not produce its INP: ' + _sw_inp)
print('SOUTHWEITION_CWD ' + _sw_os.path.dirname(_sw_inp))
print('SOUTHWEITION_JOB ' + _sw_name)
print('SOUTHWEITION_MODEL ' + str(getattr(mdb.jobs[_sw_name], 'model', '')))
print('SOUTHWEITION_INP ' + _sw_inp)
print('SOUTHWEITION_JOB_PREPARED')