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 mesh
import regionToolset
model_name = 'XFEM_Fast_Demo'
index = 1
while model_name in mdb.models.keys() or model_name in mdb.jobs.keys():
    model_name = 'XFEM_Fast_Demo_%d' % index
    index += 1
model = mdb.Model(name=model_name)
width = 0.2
height = 0.2
thickness = 0.01
seed_length = 0.04
crack_y = 0.1013
opening = 0.0002
tol = 1e-07
profile = model.ConstrainedSketch(name='Specimen_Profile', sheetSize=0.5)
profile.rectangle(point1=(0.0, 0.0), point2=(width, height))
specimen = model.Part(name='Concrete_Specimen', dimensionality=THREE_D, type=DEFORMABLE_BODY)
specimen.BaseSolidExtrude(sketch=profile, depth=thickness)
material = model.Material(name='Concrete_XFEM_Demo')
material.Density(table=((2400.0,),))
material.Elastic(table=((30000000000.0, 0.2),))
material.MaxpsDamageInitiation(table=((2000000.0,),))
material.maxpsDamageInitiation.DamageEvolution(type=ENERGY, softening=LINEAR, table=((100.0,),))
material.maxpsDamageInitiation.DamageStabilizationCohesive(cohesiveCoeff=1e-05)
model.HomogeneousSolidSection(name='Concrete_Section', material='Concrete_XFEM_Demo')
specimen.Set(name='ALL_CELLS', cells=specimen.cells[:])
specimen.SectionAssignment(region=specimen.sets['ALL_CELLS'], sectionName='Concrete_Section')
specimen.seedPart(size=0.005, deviationFactor=0.1, minSizeFactor=0.1)
specimen.setMeshControls(regions=specimen.cells[:], elemShape=HEX, technique=STRUCTURED)
specimen.setElementType(regions=(specimen.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D8, elemLibrary=STANDARD),))
specimen.generateMesh()
if len(specimen.elements) == 0:
    raise RuntimeError('Specimen mesh is empty.')
if specimen.getUnmeshedRegions() is not None:
    raise RuntimeError('Specimen mesh is incomplete.')
for element in specimen.elements:
    if len(element.connectivity) != 8:
        raise RuntimeError('Non-hexahedral specimen element detected.')
root = model.rootAssembly
root.DatumCsysByDefault(CARTESIAN)
instance = root.Instance(name='Specimen-1', part=specimen, dependent=ON)
crack_profile = model.ConstrainedSketch(name='Crack_Profile', sheetSize=0.3)
crack_profile.rectangle(point1=(0.0, 0.0), point2=(seed_length, thickness))
crack_part = model.Part(name='Initial_Crack_Surface', dimensionality=THREE_D, type=DEFORMABLE_BODY)
crack_part.BaseShell(sketch=crack_profile)
crack_instance = root.Instance(name='Initial_Crack-1', part=crack_part, dependent=ON)
root.rotate(instanceList=('Initial_Crack-1',), axisPoint=(0.0, 0.0, 0.0), axisDirection=(1.0, 0.0, 0.0), angle=90.0)
root.translate(instanceList=('Initial_Crack-1',), vector=(0.0, crack_y, 0.0))
root.engineeringFeatures.XFEMCrack(name='Edge_Crack', crackDomain=regionToolset.Region(cells=instance.cells[:]), crackLocation=regionToolset.Region(faces=crack_instance.faces[:]), allowCrackGrowth=ON)
bottom = instance.faces.getByBoundingBox(xMin=-tol, xMax=width + tol, yMin=-tol, yMax=tol, zMin=-tol, zMax=thickness + tol)
top = instance.faces.getByBoundingBox(xMin=-tol, xMax=width + tol, yMin=height - tol, yMax=height + tol, zMin=-tol, zMax=thickness + tol)
if len(bottom) == 0 or len(top) == 0:
    raise RuntimeError('Loading faces were not selected.')
root.Set(name='TOP_FACE', faces=top)
root.Set(name='BOTTOM_FACE', faces=bottom)
model.DisplacementBC(name='Bottom_Fixed', createStepName='Initial', region=root.sets['BOTTOM_FACE'], u1=0.0, u2=0.0, u3=0.0)
top_nodes = instance.nodes.getByBoundingBox(xMin=-tol, xMax=width + tol, yMin=height - tol, yMax=height + tol, zMin=-tol, zMax=thickness + tol)
if len(top_nodes) == 0:
    raise RuntimeError('Top loading nodes are missing.')
root.Set(name='TOP_NODES', nodes=top_nodes)
corner_nodes = instance.nodes.getByBoundingBox(xMin=-tol, xMax=tol, yMin=height - tol, yMax=height + tol, zMin=-tol, zMax=tol)
if len(corner_nodes) != 1:
    raise RuntimeError('Upper lateral restraint node is not unique.')
root.Set(name='UPPER_ANCHOR', nodes=corner_nodes)
model.DisplacementBC(name='Upper_Lateral_Anchor', createStepName='Initial', region=root.sets['UPPER_ANCHOR'], u1=0.0, u3=0.0)
model.StaticStep(name='Crack_Opening', previous='Initial', timePeriod=1.0, initialInc=0.0001, minInc=1e-10, maxInc=0.005, maxNumInc=10000, nlgeom=OFF, matrixStorage=UNSYMMETRIC, stabilizationMethod=DISSIPATED_ENERGY_FRACTION, stabilizationMagnitude=1e-05)
model.SmoothStepAmplitude(name='Smooth_Opening', timeSpan=STEP, data=((0.0, 0.0), (1.0, 1.0)))
model.DisplacementBC(name='Top_Opening', createStepName='Crack_Opening', region=root.sets['TOP_FACE'], u2=opening, amplitude='Smooth_Opening')
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='XFEM_Fields', createStepName='Crack_Opening', variables=('S', 'U', 'RF', 'PHILSM', 'PSILSM', 'STATUSXFEM'), numIntervals=39, timeMarks=ON)
model.HistoryOutputRequest(name='Top_Reactions', createStepName='Crack_Opening', variables=('U2', 'RF2'), region=root.sets['TOP_NODES'], frequency=1)
model.HistoryOutputRequest(name='Energy_Check', createStepName='Crack_Opening', variables=('ALLIE', 'ALLSE', 'ALLWK', 'ALLSD', 'ALLVD'), frequency=1)
root.regenerate()
model.setValues(description='Independent SI-unit XFEM demonstration, not the buried tunnel. Concrete coupon 0.20 x 0.20 x 0.010 m. Initial left-edge crack length 0.040 m at y=0.1013 m, through the coupon thickness. Top opening 0.00020 m; MAXPS=2 MPa; Gf=100 N/m. Viscosity=1e-5; static stabilization energy fraction=1e-5. Validate crack growth and artificial energy after solving.')
job = mdb.Job(name=model_name, model=model_name, type=ANALYSIS, numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE)
_sw_record_write(job, consistencyChecking=ON)
viewport = session.viewports[session.currentViewportName]
viewport.setValues(displayedObject=root)
viewport.assemblyDisplay.setValues(mesh=ON, loads=OFF, bcs=OFF, interactions=OFF)
viewport.view.setValues(session.views['Front'])
viewport.view.fitView()
print('Input written: ' + job.name)
print('Elements: %d' % len(specimen.elements))
print('Initial crack: x=0..0.040, y=0.1013, z=0..0.010.')
print('Expected propagation direction: increasing x.')
print('Requested field frames: initial plus 39 intervals.')
print('Check PHILSM and STATUSXFEM against initial crack extent.')
print('Sum TOP_NODES RF2 for the reaction-displacement curve.')
print('Check ALLSD and ALLVD relative to physical energy.')
print('No solver submission or result verification was performed.')
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')