# -*- coding: utf-8 -*-
# 在 Abaqus/CAE 中通过 File -> Run Script 执行。
# 单位：m、kg、s、Pa。
# 原文：Smojver & Ivancevic, CST 71 (2011), 489-498.
# DOI: 10.1016/j.compscitech.2010.12.024
# 原文对象为内侧襟翼，不是完整飞机机翼。
# 保留：4.38 m展长、1.05 m弦长、1.81 kg、100 m/s、55度、
# 938 kg/m3、USUP(1480,0,0)、23层/2.875 mm撞击区蒙皮。
# 近似：翼型、肋位、支承、局部厚度、夹芯等效、断裂能。
# 原文缺少完整CAD与损伤演化数据，不能视为严格定量复现。
# 原文鸟体长径描述存在歧义；这里取圆柱段长度=2R，含端总长=4R。
# 默认约75000个欧拉单元；细小飞溅需要加密并做收敛分析。
from abaqus import mdb, session
from abaqusConstants import *
import mesh
import regionToolset
import section
import math
import os

BASE = 'SouthweitionAI_CEL'
name = BASE
counter = 1
while name in mdb.models or name in mdb.jobs or os.path.exists(name + '.inp'):
    name = BASE + '_%02d' % counter
    counter += 1
model = mdb.Model(name=name, modelType=STANDARD_EXPLICIT)
assembly = model.rootAssembly
assembly.DatumCsysByDefault(CARTESIAN)

CHORD, SPAN = 1.05, 4.38
NX, NZ = 70, 146
HE = 0.020
MASS, RHO, SPEED = 1.81, 938.0, 100.0
END = 0.006
angle = math.radians(55.0)
direction = (math.cos(angle), math.sin(angle), 0.0)
radius = (MASS / (RHO * math.pi * (10.0 / 3.0))) ** (1.0 / 3.0)

# 原文表1；未给出的横向剪切模量取近似值。
cfrp = model.Material(name='CFRP')
cfrp.Density(table=((1600.0,),))
cfrp.Elastic(type=LAMINA,
    table=((181.e9, 10.3e9, 0.28, 7.17e9, 7.17e9, 3.8e9),))
cfrp.HashinDamageInitiation(
    table=((1500.e6, 1500.e6, 40.e6, 246.e6, 68.e6, 68.e6),))
# 以下断裂能是演示假设，原文未提供，单位J/m2。
cfrp.hashinDamageInitiation.DamageEvolution(
    type=ENERGY, softening=LINEAR,
    table=((15000.0, 15000.0, 500.0, 1500.0),))

for matname, rho, young, yieldstress in (
        ('Al2024', 2768.0, 73.77e9, 263.e6),
        ('Al7075', 2796.0, 71.e9, 434.e6)):
    mat = model.Material(name=matname)
    mat.Density(table=((rho,),))
    mat.Elastic(table=((young, 0.33),))
    # 简化为理想弹塑性，不虚构原文未给出的完整硬化曲线。
    mat.Plastic(table=((yieldstress, 0.0),))
    model.HomogeneousShellSection(
        name=matname, material=matname, thickness=0.002,
        integrationRule=SIMPSON, numIntPts=5)

bird = model.Material(name='Bird')
bird.Density(table=((RHO,),))
bird.Eos(type=USUP, table=((1480.0, 0.0, 0.0),))
model.EulerianSection(name='EulerSection', data={'Bird': 'Bird'})

# 23层：以一个90度中间层实现对称，而不是误写成24层。
half = [45, 0, -45, 90, -45, 0, 45, 45, 0, 0, -45]
angles = half + [90] + list(reversed(half))
layers = tuple(section.SectionLayer(
    material='CFRP', thickness=0.000125,
    orientAngle=float(a), numIntPts=3) for a in angles)
model.CompositeShellSection(
    name='Laminate23', layup=layers,
    preIntegrate=OFF, idealization=NO_IDEALIZATION,
    integrationRule=SIMPSON, thicknessType=UNIFORM)

# 共节点孤立网格：蒙皮、梁、肋、长桁直接连接。
wing = model.Part(name='Flap', dimensionality=THREE_D,
                  type=DEFORMABLE_BODY)
node_map = {}
groups = dict((k, []) for k in
              ('Upper', 'Lower', 'Stringers', 'Ribs', 'Spars'))

def thickness_y(x):
    q = max(0.0, min(1.0, x / CHORD))
    return 5.0 * 0.15 * CHORD * (
        0.2969 * math.sqrt(q) - 0.1260*q - 0.3516*q*q
        + 0.2843*q**3 - 0.1036*q**4)

def node_at(point):
    key = tuple(round(v, 10) for v in point)
    if key not in node_map:
        label = len(node_map) + 1
        wing.Node(coordinates=point, label=label)
        node_map[key] = label
    return wing.nodes[node_map[key] - 1]

def face(points, group):
    nodes = []
    for pt in points:
        nd = node_at(pt)
        if not any(old.label == nd.label for old in nodes):
            nodes.append(nd)
    if len(nodes) < 3:
        return
    label = len(wing.elements) + 1
    wing.Element(nodes=tuple(nodes),
                 elemShape=QUAD4 if len(nodes) == 4 else TRI3,
                 label=label)
    groups[group].append(label)

xs = [CHORD*i/NX for i in range(NX+1)]
zs = [SPAN*j/NZ for j in range(NZ+1)]
for side, group in ((1.0, 'Upper'), (-1.0, 'Lower')):
    for i in range(NX):
        for j in range(NZ):
            pts = [(xs[k], side*thickness_y(xs[k]), zs[l])
                   for k, l in ((i,j), (i+1,j), (i+1,j+1), (i,j+1))]
            face(pts if side > 0 else list(reversed(pts)), group)

# 前后梁间距0.525 m；绝对位置是几何近似。
for x in (0.210, 0.735):
    h = thickness_y(x)
    for j in range(NZ):
        for k in range(6):
            y0, y1 = -h+2*h*k/6.0, -h+2*h*(k+1)/6.0
            face([(x,y0,zs[j]), (x,y1,zs[j]),
                  (x,y1,zs[j+1]), (x,y0,zs[j+1])], 'Spars')

# 近似7道实腹肋，未虚构原文复杂开孔和连接耳片。
for j in (0, 24, 49, 73, 97, 122, 146):
    for i in range(NX):
        x0, x1 = xs[i], xs[i+1]
        h0, h1 = thickness_y(x0), thickness_y(x1)
        for k in range(6):
            f0, f1 = -1.0+k/3.0, -1.0+(k+1)/3.0
            face([(x0,h0*f0,zs[j]), (x1,h1*f0,zs[j]),
                  (x1,h1*f1,zs[j]), (x0,h0*f1,zs[j])], 'Ribs')

# 上6、下4长桁；间距分别为75 mm和105 mm。
for side, locations in (
        (1.0, [0.285+0.075*k for k in range(6)]),
        (-1.0, [0.315+0.105*k for k in range(4)])):
    for x in locations:
        y0 = side*thickness_y(x)
        y1 = y0-side*0.025
        for j in range(NZ):
            face([(x,y0,zs[j]), (x,y1,zs[j]),
                  (x,y1,zs[j+1]), (x,y0,zs[j+1])], 'Stringers')

for key, labels in groups.items():
    region = wing.Set(name=key,
        elements=wing.elements.sequenceFromLabels(labels=tuple(labels)))
    sec = 'Al2024' if key == 'Ribs' else (
          'Al7075' if key == 'Spars' else 'Laminate23')
    wing.SectionAssignment(region=region, sectionName=sec)

wing.Set(name='All', elements=wing.elements[:])
wing.setElementType(regions=(wing.elements[:],), elemTypes=(
    mesh.ElemType(elemCode=S4R, elemLibrary=EXPLICIT,
                  hourglassControl=ENHANCED, elemDeletion=ON,
                  maxDegradation=0.95),
    mesh.ElemType(elemCode=S3R, elemLibrary=EXPLICIT)))
datum = wing.DatumCsysByThreePoints(
    name='SpanDirection', coordSysType=CARTESIAN,
    origin=(0,0,0), point1=(0,0,1), point2=(1,0,0))
for key in ('Upper', 'Lower', 'Stringers'):
    wing.MaterialOrientation(
        region=wing.sets[key], orientationType=SYSTEM,
        axis=AXIS_3, localCsys=wing.datums[datum.id],
        stackDirection=STACK_3)
wi = assembly.Instance(name='Flap-1', part=wing, dependent=ON)

# 原文欧拉域1x1x0.6 m；20 mm规则网格约为原文单元数的1/8。
sketch = model.ConstrainedSketch(name='EulerBox', sheetSize=3.0)
sketch.rectangle(point1=(-0.2,-0.5), point2=(0.8,0.5))
ep = model.Part(name='EulerDomain', dimensionality=THREE_D, type=EULERIAN)
ep.BaseSolidExtrude(sketch=sketch, depth=0.6)
ep.Set(name='All', cells=ep.cells[:])
ep.SectionAssignment(region=ep.sets['All'], sectionName='EulerSection')
ep.setMeshControls(regions=ep.cells[:], elemShape=HEX, technique=STRUCTURED)
ep.setElementType(regions=(ep.cells[:],), elemTypes=(
    mesh.ElemType(elemCode=EC3D8R, elemLibrary=EXPLICIT),))
ep.seedPart(size=HE, deviationFactor=0.1, minSizeFactor=0.1)
ep.generateMesh()
ei = assembly.Instance(name='Euler-1', part=ep, dependent=ON)
assembly.translate(instanceList=('Euler-1',), vector=(0,0,2.6))
assembly.regenerate()

# 鸟从下方斜撞蒙皮，留初始间隙；x为弦向，z为展向。
hit = (0.38, -thickness_y(0.38), 2.9)
distance = (radius*(1.0+direction[1])+0.012)/direction[1]
center = tuple(hit[k]-distance*direction[k] for k in range(3))
bins = {}
cell_volume = None
offsets = [(i+0.5)/4.0 for i in range(4)]
for element in ep.elements:
    pts = [n.coordinates for n in element.getNodes()]
    lo = [min(p[k] for p in pts) for k in range(3)]
    hi = [max(p[k] for p in pts) for k in range(3)]
    lo[2] += 2.6
    hi[2] += 2.6
    mid = [(lo[k]+hi[k])*0.5 for k in range(3)]
    if any(abs(mid[k]-center[k]) > 2*radius+HE for k in range(3)):
        continue
    count = 0
    for u in offsets:
        for v in offsets:
            for w in offsets:
                q = [lo[k]+(hi[k]-lo[k])*f-center[k]
                     for k, f in enumerate((u,v,w))]
                axial = sum(q[k]*direction[k] for k in range(3))
                nearest = max(-radius, min(radius, axial))
                d2 = sum((q[k]-nearest*direction[k])**2 for k in range(3))
                count += int(d2 <= radius*radius)
    if count:
        bins.setdefault(count, []).append(element.label)
        cell_volume = math.prod(hi[k]-lo[k] for k in range(3)) if hasattr(
            math, 'prod') else (hi[0]-lo[0])*(hi[1]-lo[1])*(hi[2]-lo[2])

assert bins and cell_volume, 'Bird initialization is empty.'
target = MASS/RHO/cell_volume
filled = sum(len(labels)*k/64.0 for k, labels in bins.items())
assert abs(filled-target)/target < 0.05, 'Initial bird mesh is too coarse.'
delta = target-filled
capacity = sum(len(labels)*((1-k/64.0) if delta >= 0 else k/64.0)
               for k, labels in bins.items() if k < 64)
assignments = []
for k, labels in sorted(bins.items()):
    fraction = k/64.0
    if k < 64:
        fraction += delta*((1-fraction) if delta >= 0 else fraction)/capacity
    assert 0.0 <= fraction <= 1.0
    region = assembly.Set(name='EVF_%02d' % k,
        elements=ei.elements.sequenceFromLabels(labels=tuple(labels)))
    assignments.append((region, (fraction,)))
model.MaterialAssignment(name='InitialBird', instanceList=(ei,),
                         useFields=OFF, assignmentList=tuple(assignments))
velocity_region = assembly.Set(name='EulerNodes', nodes=ei.nodes[:])
model.Velocity(name='BirdVelocity', region=velocity_region,
    velocity1=SPEED*direction[0], velocity2=SPEED*direction[1], velocity3=0.0)

# 近似两处前梁安装约束，避免不合理地固支整个襟翼边界。
mounts = [n.label for n in wi.nodes
          if abs(n.coordinates[0]-0.210) < 1.e-7
          and any(abs(n.coordinates[2]-zs[j]) < 1.e-7 for j in (24,122))]
support = assembly.Set(name='Mounts',
    nodes=wi.nodes.sequenceFromLabels(labels=tuple(mounts)))
model.EncastreBC(name='MountBC', createStepName='Initial', region=support)
model.ExplicitDynamicsStep(name='Impact', previous='Initial',
                          timePeriod=END, nlgeom=ON)
model.ContactProperty('Frictionless')
model.interactionProperties['Frictionless'].NormalBehavior(
    pressureOverclosure=HARD, allowSeparation=ON)
contact = model.ContactExp(name='CELContact', createStepName='Initial')
contact.includedPairs.setValuesInStep(stepName='Initial', useAllstar=ON)
contact.contactPropertyAssignments.appendInStep(
    stepName='Initial', assignments=((GLOBAL, SELF, 'Frictionless'),))

for key in list(model.fieldOutputRequests.keys()):
    del model.fieldOutputRequests[key]
model.FieldOutputRequest(name='BirdFields', createStepName='Impact',
    region=ei.sets['All'], variables=('EVF','V','S'), numIntervals=24)
model.FieldOutputRequest(name='FlapFields', createStepName='Impact',
    region=wi.sets['All'], variables=('U','V','S','LE','STATUS',
    'DAMAGEFT','DAMAGEFC','DAMAGEMT','DAMAGEMC','PEEQ'), numIntervals=24)
for key in list(model.historyOutputRequests.keys()):
    del model.historyOutputRequests[key]
model.HistoryOutputRequest(name='Energy', createStepName='Impact',
    variables=('ALLKE','ALLIE','ALLAE','ALLVD','ETOTAL'), numIntervals=120)
# 24帧间隔包含1.5、3、4.5、6 ms；时刻以分析起点计，含初始间隙。
# 求解后在CAE隐藏上蒙皮/上长桁，使用Bird的EVF=0.5等值面观察飞溅。
# 显示剖切不应通过删除实体蒙皮来实现。
job = mdb.Job(name=name, model=name, type=ANALYSIS,
              numCpus=4, numDomains=4,
              parallelizationMethodExplicit=DOMAIN,
              explicitPrecision=DOUBLE, nodalOutputPrecision=SINGLE,
              description='Reduced CEL flap bird impact; assumptions in script')
job.writeInput(consistencyChecking=ON)
print('Model and job created: ' + name)
print('Input written: ' + os.path.abspath(name + '.inp'))
print('Elements: flap=%d, Euler=%d' % (len(wing.elements), len(ep.elements)))
print('Initialized bird mass: %.6f kg' % MASS)
print('Submit this existing job later through the CAE Job Manager.')
print('SOUTHWEITION_JOB_PREPARED')
import os as _sw_os
try:
    print('SOUTHWEITION_CWD ' + _sw_os.getcwd())
    def _sw_job_model(_job):
        _model = getattr(_job, 'model', '')
        try:
            return str(_model.name)
        except Exception:
            return str(_model or '')
    _sw_models = list(mdb.models.keys())
    _sw_latest_model = _sw_models[-1] if _sw_models else ''
    _sw_name = ''
    for _cand in reversed(list(mdb.jobs.keys())):
        if _sw_latest_model and _sw_job_model(mdb.jobs[_cand]) == _sw_latest_model:
            _sw_name = str(_cand)
            break
    if not _sw_name:
        _sw_names = list(mdb.jobs.keys())
        if _sw_names:
            _sw_name = str(_sw_names[-1])
    if _sw_name:
        try:
            mdb.jobs[_sw_name].setValues(numCpus=4, numDomains=4)
        except Exception:
            pass
        mdb.jobs[_sw_name].writeInput(consistencyChecking=OFF)
        print('SOUTHWEITION_JOB ' + str(_sw_name))
        print('SOUTHWEITION_MODEL ' + str(_sw_latest_model))
        print('SOUTHWEITION_INP ' + _sw_os.path.abspath(str(_sw_name) + '.inp'))
    print('SOUTHWEITION_JOB_PREPARED')
except Exception:
    print('SOUTHWEITION_JOB_PREPARED')