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
import math
import random
base = 'Tunnel_3D_MultiCrack_XFEM'
name = base
n = 1
while name in mdb.models.keys() or name in mdb.jobs.keys():
    name = '%s_%d' % (base, n)
    n += 1
m = mdb.Model(name=name)
(ri, ro) = (2.0, 2.18)
L = 3.6
(half_width, ground_y) = (6.0, 7.0)
tol = 1e-05
rng = random.Random(20261023)
centers = (0.68, 1.81, 2.94)
zones = ((0.22, 1.12), (1.37, 2.24), (2.5, 3.38))
s = m.ConstrainedSketch(name='Lining_Profile', sheetSize=16.0)
s.ArcByCenterEnds(center=(0.0, 0.0), point1=(ro, 0.0), point2=(-ro, 0.0), direction=COUNTERCLOCKWISE)
s.Line(point1=(-ro, 0.0), point2=(-ri, 0.0))
s.ArcByCenterEnds(center=(0.0, 0.0), point1=(-ri, 0.0), point2=(ri, 0.0), direction=CLOCKWISE)
s.Line(point1=(ri, 0.0), point2=(ro, 0.0))
p = m.Part(name='Tunnel_Lining', dimensionality=THREE_D, type=DEFORMABLE_BODY)
p.BaseSolidExtrude(sketch=s, depth=L)
ss = m.ConstrainedSketch(name='Soil_Profile', sheetSize=20.0)
outline = [(-half_width, 0.0), (-half_width, ground_y), (half_width, ground_y), (half_width, 0.0), (ro, 0.0)]
for j in range(len(outline) - 1):
    ss.Line(point1=outline[j], point2=outline[j + 1])
ss.ArcByCenterEnds(center=(0.0, 0.0), point1=(ro, 0.0), point2=(-ro, 0.0), direction=COUNTERCLOCKWISE)
ss.Line(point1=(-ro, 0.0), point2=(-half_width, 0.0))
soil = m.Part(name='Soil_Block', dimensionality=THREE_D, type=DEFORMABLE_BODY)
soil.BaseSolidExtrude(sketch=ss, depth=L)
con = m.Material(name='XFEM_Concrete')
con.Density(table=((2400.0,),))
con.Elastic(table=((30000000000.0, 0.2),))
con.MaxpsDamageInitiation(table=((2000000.0,),))
con.maxpsDamageInitiation.DamageEvolution(type=ENERGY, softening=LINEAR, table=((100.0,),))
con.maxpsDamageInitiation.DamageStabilizationCohesive(cohesiveCoeff=5e-05)
sm = m.Material(name='Elastic_Soil')
sm.Density(table=((1900.0,),))
sm.Elastic(table=((50000000.0, 0.3),))
for (part, mat_name, sec) in ((p, 'XFEM_Concrete', 'Lining_Section'), (soil, 'Elastic_Soil', 'Soil_Section')):
    m.HomogeneousSolidSection(name=sec, material=mat_name)
    part.Set(name='ALL_CELLS', cells=part.cells[:])
    part.SectionAssignment(region=part.sets['ALL_CELLS'], sectionName=sec)
for z in (0.22, 1.12, 1.37, 2.24, 2.5, 3.38):
    d = p.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=z)
    p.PartitionCellByDatumPlane(cells=p.cells[:], datumPlane=p.datums[d.id])
for theta in (1.16, 1.98):
    d0 = p.DatumPointByCoordinate(coords=(0.0, 0.0, 0.0))
    d1 = p.DatumPointByCoordinate(coords=(math.cos(theta), math.sin(theta), 0.0))
    d2 = p.DatumPointByCoordinate(coords=(0.0, 0.0, L))
    plane = p.DatumPlaneByThreePoints(point1=p.datums[d0.id], point2=p.datums[d1.id], point3=p.datums[d2.id])
    p.PartitionCellByDatumPlane(cells=p.cells[:], datumPlane=p.datums[plane.id])
p.seedPart(size=0.24, deviationFactor=0.1, minSizeFactor=0.1)
radial = p.edges[0:0]
refined = p.edges[0:0]
for (j, edge) in enumerate(p.edges):
    vids = edge.getVertices()
    if len(vids) != 2:
        continue
    q1 = p.vertices[vids[0]].pointOn[0]
    q2 = p.vertices[vids[1]].pointOn[0]
    r1 = math.hypot(q1[0], q1[1])
    r2 = math.hypot(q2[0], q2[1])
    dz = abs(q1[2] - q2[2])
    dx = abs(q1[0] - q2[0])
    dy = abs(q1[1] - q2[1])
    zm = 0.5 * (q1[2] + q2[2])
    crown = min(q1[1], q2[1]) > 1.7
    if abs(r1 - r2) > 0.12 and dz < tol:
        radial = radial + p.edges[j:j + 1]
    elif crown and dz < tol and (abs(r1 - r2) < tol):
        refined = refined + p.edges[j:j + 1]
    elif crown and dz > 0.1 and (dx < tol) and (dy < tol):
        if any((lo < zm < hi for (lo, hi) in zones)):
            refined = refined + p.edges[j:j + 1]
if len(radial) == 0 or len(refined) == 0:
    raise RuntimeError('Local refinement edges not found.')
p.seedEdgeByNumber(edges=radial, number=5, constraint=FINER)
p.seedEdgeBySize(edges=refined, size=0.12, constraint=FINER)
p.setMeshControls(regions=p.cells[:], elemShape=HEX, technique=SWEEP)
p.setElementType(regions=(p.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D8, elemLibrary=STANDARD),))
p.generateMesh()
soil.seedPart(size=0.75, deviationFactor=0.1, minSizeFactor=0.1)
soil.setMeshControls(regions=soil.cells[:], elemShape=TET, technique=FREE)
soil.setElementType(regions=(soil.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D4, elemLibrary=STANDARD),))
soil.generateMesh()
for part in (p, soil):
    if not len(part.elements) or part.getUnmeshedRegions() is not None:
        raise RuntimeError('Incomplete mesh: ' + part.name)
for elem in p.elements:
    if len(elem.connectivity) != 8:
        raise RuntimeError('Lining requires all HEX8 elements.')
a = m.rootAssembly
a.DatumCsysByDefault(CARTESIAN)
li = a.Instance(name='Lining-1', part=p, dependent=ON)
si = a.Instance(name='Soil-1', part=soil, dependent=ON)
crack_info = []
for (k, zc) in enumerate(centers):
    xc = 0.05 if k % 2 == 0 else -0.055
    xc += rng.uniform(-0.012, 0.012)
    axial_len = 0.42 + rng.uniform(-0.035, 0.035)
    inner_y = math.sqrt(ri * ri - xc * xc)
    y0 = inner_y - 0.008
    y1 = inner_y + 0.07
    z0 = zc - axial_len / 2.0
    csk = m.ConstrainedSketch(name='Crack_Sketch_%d' % (k + 1), sheetSize=1.0)
    csk.rectangle(point1=(0.0, 0.0), point2=(axial_len, y1 - y0))
    cp = m.Part(name='Crack_Surface_%d' % (k + 1), dimensionality=THREE_D, type=DEFORMABLE_BODY)
    cp.BaseShell(sketch=csk)
    inst_name = 'Initial_Crack_%d' % (k + 1)
    a.Instance(name=inst_name, part=cp, dependent=ON)
    a.rotate(instanceList=(inst_name,), axisPoint=(0.0, 0.0, 0.0), axisDirection=(0.0, 1.0, 0.0), angle=-90.0)
    a.translate(instanceList=(inst_name,), vector=(xc, y0, z0))
    a.regenerate()
    domain = li.cells.findAt(((0.0, 0.5 * (ri + ro), zc),))
    a.engineeringFeatures.XFEMCrack(name='XFEM_Crack_%d' % (k + 1), crackDomain=regionToolset.Region(cells=domain), crackLocation=regionToolset.Region(faces=a.instances[inst_name].faces[:]), allowCrackGrowth=ON)
    crack_info.append((xc, z0, axial_len, y1 - y0))
outer_faces = li.faces[0:0]
inner_faces = li.faces[0:0]
for (j, face) in enumerate(li.faces):
    q = face.pointOn[0]
    r = math.hypot(q[0], q[1])
    if abs(r - ro) < 0.001 and q[1] > 0.0:
        outer_faces = outer_faces + li.faces[j:j + 1]
    elif abs(r - ri) < 0.001 and q[1] > 0.0:
        inner_faces = inner_faces + li.faces[j:j + 1]
if len(outer_faces) == 0 or len(inner_faces) == 0:
    raise RuntimeError('Curved lining faces not selected.')
a.Surface(name='LINING_OUTER', side1Faces=outer_faces)
a.Surface(name='LINING_INNER', side1Faces=inner_faces)
soil_arch = si.faces.findAt(((0.0, ro, 1.8),))
a.Surface(name='SOIL_ARCH', side1Faces=soil_arch)
m.Tie(name='Soil_Lining_Tie', main=a.surfaces['SOIL_ARCH'], secondary=a.surfaces['LINING_OUTER'], positionToleranceMethod=COMPUTED, adjust=ON, tieRotations=ON)
bottom = si.faces.findAt(((-4.0, 0.0, L / 2.0),), ((4.0, 0.0, L / 2.0),))
sides = si.faces.findAt(((-half_width, 3.0, L / 2.0),), ((half_width, 3.0, L / 2.0),))
front = si.faces.findAt(((0.0, 6.0, 0.0),))
top = si.faces.findAt(((0.0, ground_y, L / 2.0),))
for (key, faces) in (('SOIL_BOTTOM', bottom), ('SOIL_SIDES', sides), ('SOIL_FRONT', front)):
    a.Set(name=key, faces=faces)
for (key, region, dof) in (('Bottom_U2', 'SOIL_BOTTOM', 2), ('Sides_U1', 'SOIL_SIDES', 1), ('Front_U3', 'SOIL_FRONT', 3)):
    values = {'u1': UNSET, 'u2': UNSET, 'u3': UNSET}
    values['u%d' % dof] = 0.0
    m.DisplacementBC(name=key, createStepName='Initial', region=a.sets[region], **values)
a.Surface(name='GROUND_TOP', side1Faces=top)
m.StaticStep(name='Initial_Check', previous='Initial', timePeriod=1.0, initialInc=0.25, minInc=1e-06, maxInc=0.5, maxNumInc=20, nlgeom=OFF)
m.ImplicitDynamicsStep(name='Crack_Drive', previous='Initial_Check', timePeriod=1.0, application=QUASI_STATIC, timeIncrementationMethod=AUTOMATIC, initialInc=0.002, minInc=1e-07, maxInc=0.03, maxNumInc=800, nlgeom=OFF, matrixStorage=UNSYMMETRIC)
m.SmoothStepAmplitude(name='Pressure_Ramp', timeSpan=STEP, data=((0.0, 0.0), (0.25, 0.08), (0.5, 0.32), (0.75, 0.7), (1.0, 1.0)))
for (key, surface, magnitude) in (('Ground_Surcharge', 'GROUND_TOP', 60000.0), ('Tunnel_Internal_Pressure', 'LINING_INNER', 650000.0)):
    m.Pressure(name=key, createStepName='Crack_Drive', region=a.surfaces[surface], magnitude=magnitude, amplitude='Pressure_Ramp')
for key in list(m.fieldOutputRequests.keys()):
    del m.fieldOutputRequests[key]
for key in list(m.historyOutputRequests.keys()):
    del m.historyOutputRequests[key]
initial_output = m.FieldOutputRequest(name='Equilibrium_Output', createStepName='Initial_Check', variables=('S', 'U', 'RF'), numIntervals=2, timeMarks=ON)
initial_output.deactivate('Crack_Drive')
m.FieldOutputRequest(name='XFEM_Animation', createStepName='Crack_Drive', variables=('S', 'U', 'RF', 'PHILSM', 'STATUSXFEM'), numIntervals=36, timeMarks=ON)
m.HistoryOutputRequest(name='Energy_Monitor', createStepName='Crack_Drive', variables=('ALLIE', 'ALLSE', 'ALLKE', 'ALLWK', 'ALLVD'), frequency=1)
a.regenerate()
m.setValues(description='Lightweight 3D buried tunnel multiple-crack XFEM demo. Soil width 12 m, height 7 m, length 3.6 m. Tunnel Ri 2 m, Ro 2.18 m. Three separated reproducible initial crown cracks. Geometric local refinement, not adaptive remeshing. Concrete E 30 GPa, ft 2 MPa, Gf 100 N/m. Elastic soil, bonded interface, no geostatic stress. Ground pressure 60 kPa, internal pressure 650 kPa. This is a high-load demonstration only. Real crack propagation requires ODB verification.')
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 HEX ELEMENTS: %d' % len(p.elements))
for (k, data) in enumerate(crack_info):
    print('CRACK %d: x=%.4f z0=%.4f length=%.4f depth=%.4f' % ((k + 1,) + data))
print('THREE NONOVERLAPPING XFEM DOMAINS')
print('LOCAL REFINEMENT SIZE: 0.12 m')
print('GROUND PRESSURE: 60000 Pa')
print('TUNNEL PRESSURE: 650000 Pa')
print('TARGET FIELD FRAMES: 40')
print('INP WRITTEN; JOB NOT SUBMITTED')
print('VERIFY PHILSM, STATUSXFEM 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')