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_FG_Clean'
name = base
number = 1
while name in mdb.models or 'Job_' + name in mdb.jobs:
    name = '%s_%02d' % (base, number)
    number += 1
model = mdb.Model(name=name)
model.setValues(description='REC-N 3x3 four-layer reconstruction. Units mm,N,MPa,tonne. Dimensions and PLA fractions follow the supplied paper. Arc radii follow common 28 mm chord and 90/135 degree angles. Horizontal end trimming and joint outline are reconstruction assumptions. Elastic and plastic parameters are design assumptions, not fitted data. One geometry union with removed intersection boundaries; four non-overlapping depth partitions with shared mesh nodes.')
WIDTH = 135.0
FILL = 20.0
CORE = 85.76
DEPTH = 15.0
HEIGHT = 2.0 * FILL + CORE
WALL = 1.76
CELL_H = 28.0
CELL_L = 32.0
PITCH = 45.0
NX = 3
NY = 3
LAYER_T = DEPTH / 4.0
PLA_FRACTIONS = (1.0, 0.9, 0.8, 0.7)
LOCAL_SIZE = 0.85
GLOBAL_SIZE = 3.5
MAX_ELEMENTS = 180000
TOL = 1e-05
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)
instance_names = []
source_parts = []
joint0 = FILL + WALL / 2.0
trim_half_height = CELL_H / 2.0 - WALL / 2.0

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

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

def radius_data(angle):
    half = math.radians(angle / 2.0)
    radius = CELL_H / (2.0 * math.sin(half))
    return (radius, radius * math.cos(half))

def add_trimmed_wall(endpoint_x, center_y, angle, left):
    (radius, offset) = radius_data(angle)
    rinner = radius - WALL / 2.0
    router = radius + WALL / 2.0
    center_x = endpoint_x - offset if left else endpoint_x + offset
    side = 1.0 if left else -1.0
    dy = trim_half_height
    if dy >= rinner:
        raise RuntimeError('Invalid arc end trimming dimensions.')
    xouter = center_x + side * math.sqrt(router * router - dy * dy)
    xinner = center_x + side * math.sqrt(rinner * rinner - dy * dy)
    outer_bottom = (xouter, center_y - dy)
    outer_top = (xouter, center_y + dy)
    inner_bottom = (xinner, center_y - dy)
    inner_top = (xinner, center_y + dy)
    sketch = model.ConstrainedSketch(name='BuildProfile', sheetSize=400.0)
    outer_direction = COUNTERCLOCKWISE if left else CLOCKWISE
    inner_direction = CLOCKWISE if left else COUNTERCLOCKWISE
    sketch.ArcByCenterEnds(center=(center_x, center_y), point1=outer_bottom, point2=outer_top, direction=outer_direction)
    sketch.Line(point1=outer_top, point2=inner_top)
    sketch.ArcByCenterEnds(center=(center_x, center_y), point1=inner_top, point2=inner_bottom, direction=inner_direction)
    sketch.Line(point1=inner_bottom, point2=outer_bottom)
    make_solid(sketch)
joint_extension = 0.0
for angle in (90.0, 135.0):
    (radius, offset) = radius_data(angle)
    inner = radius - WALL / 2.0
    extension = offset - math.sqrt(inner * inner - trim_half_height ** 2)
    joint_extension = max(joint_extension, extension)
columns = []
for column in range(NX):
    xcenter = (column + 0.5) * PITCH
    columns.append((xcenter - CELL_L / 2.0, xcenter + CELL_L / 2.0))
add_bar(0.0, 0.0, WIDTH, FILL)
add_bar(0.0, FILL + CORE, WIDTH, HEIGHT)
for row in range(NY + 1):
    yjoint = joint0 + row * CELL_H
    for (xleft, xright) in columns:
        add_bar(xleft - joint_extension, yjoint - WALL / 2.0, xright + joint_extension, yjoint + WALL / 2.0)
sag90 = CELL_H / 2.0 * (math.sqrt(2.0) - 1.0)
for row in range(NY):
    ycenter = joint0 + (row + 0.5) * CELL_H
    for (xleft, xright) in columns:
        for angle in (90.0, 135.0):
            add_trimmed_wall(xleft, ycenter, angle, True)
            add_trimmed_wall(xright, ycenter, angle, False)
    bridges = [(0.0, columns[0][0] + sag90)]
    for column in range(NX - 1):
        bridges.append((columns[column][1] - sag90, columns[column + 1][0] + sag90))
    bridges.append((columns[-1][1] - sag90, WIDTH))
    for (x0, x1) in bridges:
        add_bar(x0, ycenter - WALL / 2.0, x1, ycenter + WALL / 2.0)
merged = assembly.InstanceFromBooleanMerge(name='CleanSpecimen', instances=tuple((assembly.instances[n] for n in instance_names)), keepIntersections=OFF, originalInstances=DELETE, domain=GEOMETRY)
part = model.parts[merged.partName]
merged_name = merged.name
for source_name in source_parts:
    del model.parts[source_name]
for z in (LAYER_T, 2.0 * LAYER_T, 3.0 * LAYER_T):
    datum = part.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=z)
    part.PartitionCellByDatumPlane(datumPlane=part.datums[datum.id], cells=part.cells[:])
for (layer, fraction) in enumerate(PLA_FRACTIONS):
    material_name = 'FG_PLA%03d' % int(round(100.0 * fraction))
    section_name = 'Layer_%d_Section' % (layer + 1)
    material = model.Material(name=material_name)
    density = fraction * 1.24e-09 + (1.0 - fraction) * 1.04e-09
    modulus = fraction * 2000.0 + (1.0 - fraction) * 1600.0
    poisson = fraction * 0.35 + (1.0 - fraction) * 0.37
    yield_stress = fraction * 30.0 + (1.0 - fraction) * 25.0
    material.Density(table=((density,),))
    material.Elastic(table=((modulus, poisson),))
    material.Plastic(table=((yield_stress, 0.0), (yield_stress + 8.0, 0.02), (yield_stress + 15.0, 0.06)))
    model.HomogeneousSolidSection(name=section_name, material=material_name, thickness=None)
    cells = part.cells.getByBoundingBox(xMin=-TOL, xMax=WIDTH + TOL, yMin=-TOL, yMax=HEIGHT + TOL, zMin=layer * LAYER_T - TOL, zMax=(layer + 1) * LAYER_T + TOL)
    if not len(cells):
        raise RuntimeError('Layer cell selection failed: %d' % (layer + 1))
    region = part.Set(name='LAYER_%d' % (layer + 1), cells=cells)
    part.SectionAssignment(region=region, sectionName=section_name)
part.setMeshControls(regions=part.cells[:], elemShape=HEX, technique=SWEEP, algorithm=MEDIAL_AXIS)
brick = mesh.ElemType(elemCode=C3D8R, elemLibrary=STANDARD, hourglassControl=ENHANCED)
part.setElementType(regions=(part.cells[:],), elemTypes=(brick,))
accepted = False
for local_size in (LOCAL_SIZE, 1.0, 1.2):
    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[v].pointOn[0] for v in edge.getVertices()]
        if not points:
            continue
        (xs, ys, zs) = zip(*points)
        axial = max(xs) - min(xs) < TOL and max(ys) - min(ys) < TOL and (max(zs) - min(zs) > TOL)
        if axial:
            part.seedEdgeByNumber(edges=(edge,), number=1, constraint=FIXED)
        elif max(ys) >= FILL - TOL and min(ys) <= FILL + CORE + TOL:
            length = edge.getSize(printResults=False)
            if length <= 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('Local seed=%g; element count=%d' % (local_size, count))
    if count == 0 or part.getUnmeshedRegions() is not None:
        raise RuntimeError('Incomplete sweep mesh; inspect clean source faces.')
    if any((len(e.connectivity) != 8 for e in part.elements)):
        raise RuntimeError('Unexpected non-eight-node mesh topology.')
    if count <= MAX_ELEMENTS:
        accepted = True
        break
if not accepted:
    raise RuntimeError('Mesh budget exceeded after automatic reseeding.')
part.setElementType(regions=(part.cells[:],), elemTypes=(brick,))
assembly.regenerate()
instance = assembly.instances[merged_name]
bottom = instance.faces.getByBoundingBox(xMin=-TOL, xMax=WIDTH + TOL, yMin=-TOL, yMax=TOL, zMin=-TOL, zMax=DEPTH + TOL)
top = instance.faces.getByBoundingBox(xMin=-TOL, xMax=WIDTH + TOL, yMin=HEIGHT - TOL, yMax=HEIGHT + TOL, zMin=-TOL, zMax=DEPTH + TOL)
if not len(bottom) or not len(top):
    raise RuntimeError('End-face selection failed.')
assembly.Set(name='BOTTOM', faces=bottom)
assembly.Surface(name='TOP', side1Faces=top)
rp = assembly.ReferencePoint(point=(WIDTH / 2.0, HEIGHT, DEPTH / 2.0))
assembly.Set(name='TOP_RP', referencePoints=(assembly.referencePoints[rp.id],))
model.Coupling(name='TopVertical', controlPoint=assembly.sets['TOP_RP'], surface=assembly.surfaces['TOP'], influenceRadius=WHOLE_SURFACE, couplingType=KINEMATIC, localCsys=None, u1=OFF, u2=ON, u3=OFF, ur1=OFF, ur2=OFF, ur3=OFF)
model.EncastreBC(name='FixBottom', createStepName='Initial', region=assembly.sets['BOTTOM'])
model.DisplacementBC(name='UnusedRP', createStepName='Initial', region=assembly.sets['TOP_RP'], u1=0.0, u3=0.0, ur1=0.0, ur2=0.0, ur3=0.0)
model.StaticStep(name='Compression', previous='Initial', nlgeom=ON, initialInc=0.02, minInc=1e-07, maxInc=0.05, maxNumInc=500)
model.DisplacementBC(name='Compress2mm', createStepName='Compression', region=assembly.sets['TOP_RP'], u2=-2.0)
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 = mdb.Job(name='Job_' + name, model=name, numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE, description='Clean REC-N FG geometry; assumed material calibration')
_sw_record_write(job, consistencyChecking=ON)
input_path = os.path.abspath(job.name + '.inp')
element_types = set()
with open(input_path, 'r') as stream:
    for line in stream:
        text = line.strip().upper()
        if re.match('^\\*ELEMENT(?:\\s*,|\\s*$)', text):
            match = re.search('(?:^|,)\\s*TYPE\\s*=\\s*([^,\\s]+)', text)
            if match:
                element_types.add(match.group(1))
if element_types != set(('C3D8R',)):
    raise RuntimeError('Unexpected INP element types: %s' % sorted(element_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('One specimen instance; four shared-node material layers.')
print('Elements=%d; nodes=%d' % (len(part.elements), len(part.nodes)))
print('INP verified: %s' % input_path)
print('Review geometry and mesh, then submit 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')