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 math
import os
import re
base = 'REC_N_PLA_BudgetFix'
model_name = base
suffix = 1
while model_name in mdb.models or 'Job_' + model_name in mdb.jobs:
    model_name = '%s_%02d' % (base, suffix)
    suffix += 1
model = mdb.Model(name=model_name)
model.setValues(description='REC-N 3x3 reconstruction from Polymers 2025,17,1547. Units mm,N,MPa,tonne. Density from Materials 2023,16,5177. Design assumptions: 28 mm common chord,90/135 degree arcs; PLA E=2000 MPa,nu=0.35 and assumed plastic hardening. Same PLA for core and fillers. Adaptive depth mesh budget.')
WIDTH = 135.0
CORE_HEIGHT = 85.76
FILL = 20.0
DEPTH = 15.0
WALL = 1.76
CELL_HEIGHT = 28.0
CELL_LENGTH = 32.0
PITCH = 45.0
NX = 3
NY = 3
HEIGHT = CORE_HEIGHT + 2.0 * FILL
COMPRESSION = 2.0
GLOBAL_SIZE = 3.0
LOCAL_SIZE = 0.8
DEPTH_TRIALS = (4, 3, 2)
MAX_ELEMENTS = 180000
TOL = 1e-05
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)
components = []

def extrude_profile(sketch):
    name = 'Solid_%03d' % (len(components) + 1)
    solid = model.Part(name=name, dimensionality=THREE_D, type=DEFORMABLE_BODY)
    solid.BaseSolidExtrude(sketch=sketch, depth=DEPTH)
    del model.sketches[sketch.name]
    instance_name = name + '-1'
    assembly.Instance(name=instance_name, part=solid, dependent=ON)
    components.append(instance_name)

def add_bar(x0, y0, x1, y1):
    sketch = model.ConstrainedSketch(name='ConstructionProfile', sheetSize=400.0)
    sketch.rectangle(point1=(x0, y0), point2=(x1, y1))
    extrude_profile(sketch)

def add_arc(endpoint_x, center_y, angle, left):
    half_angle = math.radians(angle * 0.5)
    radius = CELL_HEIGHT / (2.0 * math.sin(half_angle))
    offset = radius * math.cos(half_angle)
    center_x = endpoint_x - offset if left else endpoint_x + offset
    direction = 0.0 if left else math.pi
    angle0 = direction - half_angle
    angle1 = direction + half_angle

    def xy(r, a):
        return (center_x + r * math.cos(a), center_y + r * math.sin(a))
    outer0 = xy(radius + WALL * 0.5, angle0)
    outer1 = xy(radius + WALL * 0.5, angle1)
    inner0 = xy(radius - WALL * 0.5, angle0)
    inner1 = xy(radius - WALL * 0.5, angle1)
    sketch = model.ConstrainedSketch(name='ConstructionProfile', sheetSize=400.0)
    sketch.ArcByCenterEnds(center=(center_x, center_y), point1=outer0, point2=outer1, direction=COUNTERCLOCKWISE)
    sketch.Line(point1=outer1, point2=inner1)
    sketch.ArcByCenterEnds(center=(center_x, center_y), point1=inner1, point2=inner0, direction=CLOCKWISE)
    sketch.Line(point1=inner0, point2=outer0)
    extrude_profile(sketch)
add_bar(0.0, 0.0, WIDTH, FILL)
add_bar(0.0, FILL + CORE_HEIGHT, WIDTH, HEIGHT)
joint0 = FILL + WALL * 0.5
column_edges = []
for column in range(NX):
    center = (column + 0.5) * PITCH
    column_edges.append((center - CELL_LENGTH * 0.5, center + CELL_LENGTH * 0.5))
for row in range(NY + 1):
    joint_y = joint0 + row * CELL_HEIGHT
    for (xleft, xright) in column_edges:
        add_bar(xleft, joint_y - WALL * 0.5, xright, joint_y + WALL * 0.5)
sag = CELL_HEIGHT * 0.5 * (math.sqrt(2.0) - 1.0)
for row in range(NY):
    center_y = joint0 + (row + 0.5) * CELL_HEIGHT
    for (xleft, xright) in column_edges:
        for angle in (90.0, 135.0):
            add_arc(xleft, center_y, angle, True)
            add_arc(xright, center_y, angle, False)
    bridges = [(0.0, column_edges[0][0] + sag)]
    for column in range(NX - 1):
        bridges.append((column_edges[column][1] - sag, column_edges[column + 1][0] + sag))
    bridges.append((column_edges[-1][1] - sag, WIDTH))
    for (x0, x1) in bridges:
        add_bar(x0, center_y - WALL * 0.5, x1, center_y + WALL * 0.5)
merged = assembly.InstanceFromBooleanMerge(name='Specimen', instances=tuple((assembly.instances[n] for n in components)), keepIntersections=ON, originalInstances=SUPPRESS, domain=GEOMETRY)
instance_name = merged.name
part = model.parts[merged.partName]
part.Set(name='ALL_SOLID', cells=part.cells[:])
material = model.Material(name='PLA_Design')
material.Density(table=((1.24e-09,),))
material.Elastic(table=((2000.0, 0.35),))
material.Plastic(table=((30.0, 0.0), (38.0, 0.02), (45.0, 0.06)))
model.HomogeneousSolidSection(name='PLA_Section', material='PLA_Design', thickness=None)
part.SectionAssignment(region=part.sets['ALL_SOLID'], sectionName='PLA_Section')
part.setMeshControls(regions=part.cells[:], elemShape=HEX, technique=SWEEP, algorithm=MEDIAL_AXIS)
brick = mesh.ElemType(elemCode=C3D8R, elemLibrary=STANDARD, kinematicSplit=AVERAGE_STRAIN, hourglassControl=ENHANCED)
part.setElementType(regions=(part.cells[:],), elemTypes=(brick,))
accepted_depth = None
for depth_divisions in DEPTH_TRIALS:
    if len(part.elements):
        part.deleteMesh(regions=part.cells[:])
    part.seedPart(size=GLOBAL_SIZE, deviationFactor=0.08, minSizeFactor=0.1)
    for edge in part.edges:
        points = [part.vertices[i].pointOn[0] for i in edge.getVertices()]
        if not points:
            continue
        (xs, ys, zs) = zip(*points)
        depth_edge = max(xs) - min(xs) < TOL and max(ys) - min(ys) < TOL and (max(zs) - min(zs) > TOL)
        if depth_edge:
            part.seedEdgeByNumber(edges=(edge,), number=depth_divisions, constraint=FIXED)
        elif max(ys) >= FILL - TOL and min(ys) <= FILL + CORE_HEIGHT + TOL:
            if edge.getSize(printResults=False) <= WALL * 1.15:
                part.seedEdgeByNumber(edges=(edge,), number=3, constraint=FINER)
            else:
                part.seedEdgeBySize(edges=(edge,), size=LOCAL_SIZE, deviationFactor=0.08, minSizeFactor=0.1, constraint=FINER)
    part.generateMesh()
    count = len(part.elements)
    print('Depth divisions=%d; elements=%d; budget=%d' % (depth_divisions, count, MAX_ELEMENTS))
    if not count:
        raise RuntimeError('Mesh generation produced no elements.')
    if part.getUnmeshedRegions() is not None:
        raise RuntimeError('Some geometric regions remain unmeshed.')
    if any((len(element.connectivity) != 8 for element in part.elements)):
        raise RuntimeError('Non-eight-node topology found in hex mesh.')
    if count <= MAX_ELEMENTS:
        accepted_depth = depth_divisions
        break
if accepted_depth is None:
    raise RuntimeError('Depth mesh trials exhausted; final element count=%d.' % len(part.elements))
part.setElementType(regions=(part.cells[:],), elemTypes=(brick,))
reported_types = {}
for element in part.elements:
    key = str(element.type).strip().upper()
    reported_types[key] = reported_types.get(key, 0) + 1
print('CAE element type counts: %s' % reported_types)
assembly.regenerate()
instance = assembly.instances[instance_name]
bottom_faces = instance.faces.getByBoundingBox(xMin=-TOL, xMax=WIDTH + TOL, yMin=-TOL, yMax=TOL, zMin=-TOL, zMax=DEPTH + TOL)
top_faces = instance.faces.getByBoundingBox(xMin=-TOL, xMax=WIDTH + TOL, yMin=HEIGHT - TOL, yMax=HEIGHT + TOL, zMin=-TOL, zMax=DEPTH + TOL)
if not len(bottom_faces) or not len(top_faces):
    raise RuntimeError('Specimen end-face selection failed.')
assembly.Set(name='BOTTOM', faces=bottom_faces)
assembly.Surface(name='TOP_SURFACE', side1Faces=top_faces)
reference = assembly.ReferencePoint(point=(WIDTH * 0.5, HEIGHT, DEPTH * 0.5))
assembly.Set(name='TOP_RP', referencePoints=(assembly.referencePoints[reference.id],))
model.Coupling(name='TopVerticalCoupling', controlPoint=assembly.sets['TOP_RP'], surface=assembly.surfaces['TOP_SURFACE'], influenceRadius=WHOLE_SURFACE, couplingType=KINEMATIC, localCsys=None, u1=OFF, u2=ON, u3=OFF, ur1=OFF, ur2=OFF, ur3=OFF)
model.EncastreBC(name='BottomFixed', createStepName='Initial', region=assembly.sets['BOTTOM'])
model.DisplacementBC(name='UnusedRP_DOF', createStepName='Initial', region=assembly.sets['TOP_RP'], u1=0.0, u2=UNSET, u3=0.0, ur1=0.0, ur2=0.0, ur3=0.0)
model.StaticStep(name='Compression', previous='Initial', nlgeom=ON, timePeriod=1.0, initialInc=0.02, minInc=1e-07, maxInc=0.05, maxNumInc=500)
model.DisplacementBC(name='CompressionLoading', createStepName='Compression', region=assembly.sets['TOP_RP'], u2=-COMPRESSION)
model.fieldOutputRequests['F-Output-1'].setValues(variables=('S', 'U', 'RF', 'LE', 'PEEQ'), numIntervals=20)
model.HistoryOutputRequest(name='ForceDisplacement', createStepName='Compression', variables=('U2', 'RF2'), region=assembly.sets['TOP_RP'], frequency=1)
job_name = 'Job_' + model_name
job = mdb.Job(name=job_name, model=model_name, type=ANALYSIS, description='REC-N 3x3; budget-controlled C3D8R; design PLA', numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE)
_sw_record_write(job, consistencyChecking=ON)
input_path = os.path.abspath(job_name + '.inp')
input_types = set()
with open(input_path, 'r') as stream:
    for line in stream:
        keyword = line.strip().upper()
        if re.match('^\\*ELEMENT(?:\\s*,|\\s*$)', keyword):
            match = re.search('(?:^|,)\\s*TYPE\\s*=\\s*([^,\\s]+)', keyword)
            if match:
                input_types.add(match.group(1))
if input_types != set(('C3D8R',)):
    raise RuntimeError('Unexpected INP element types: %s.' % sorted(input_types))
if 'Viewport: 1' in session.viewports:
    viewport = session.viewports['Viewport: 1']
    viewport.setValues(displayedObject=assembly)
    viewport.view.setValues(session.views['Iso'])
    viewport.view.fitView()
print('Model: %s' % model_name)
print('Accepted depth divisions: %d' % accepted_depth)
print('Nodes=%d; elements=%d' % (len(part.nodes), len(part.elements)))
print('Verified INP element types: %s' % sorted(input_types))
print('Input written: %s' % input_path)
print('Material and arc construction retain declared design assumptions.')
print('Review mesh and submit the existing job through Job Manager.')
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')