Capstone B W13 logbook

Small gap between updates. Many things were done, such as:

• Developing Python code to create the input files.
• Developing Python code that could run all the jobs in a batched way
• Developing python code to read the ODBs and get key information
• Use that information in MATLAB to then get more data.
• Develop Python code to get the contour images for me (Why do something repetitive for 2 hours if I can just spend a mere 6 hours automating it)
• Actual analysis of results

I ended up creating 434 jobs, which took over ~10 hours to run

I’m simply going to embed some files here and some in the pdf I’ll submit, such as the matlab code.

PatchAnalysisPy

from odbAccess import *
from abaqus import *
import os
# from loadinginfo import createLoadingData
from typing import Any
import numpy as np
import datetime
import traceback
 
mdb: Mdb = openMdb("../patch.cae")
model = mdb.models["Model-1"]
# set loading
 
 
def jobWithLoading(mdb, model, dx, dz, dr, with_patch):
    # suppressing/enabling
    def setFeatureStatus(parent, objects, enable):
        if enable:
            return parent.resumeFeatures(objects)
        else:
            return parent.suppressFeatures(objects)
    def setStatus(parent, objects, enable):
        for object in objects:
            if enable: parent[object].resume()
            else: parent[object].suppress()
    setFeatureStatus(model.rootAssembly, ("Adhesive-Layer-1", "Patch-1"), with_patch)
    setStatus(model.constraints, ("Adhesive-Patch-1", "Surface-Adhesive"), with_patch)
    setStatus(model.fieldOutputRequests, ("Adhesive", "Composite"), with_patch)
    model.ExpressionField("RotationY", "Z*sin(%f) + Y*(cos(%f)-1)" % (dr, dr))
    model.ExpressionField(
        "RotationZ", "Z * (cos(%f)-1) + Y*sin(%f) + %f" % (dr, dr, dz)
    )
    model.boundaryConditions["StretchR"].setValues(u1=dx, ur1=-dr) # must be negative rotation
    patchText = "" if with_patch else "less"
    name = ("Patch%s%+04.2fX%+04.2fZ%+05.2fR" % (patchText, dx, dz, dr * 180 / pi)).replace(".", "_")
    print(name)
    return mdb.Job(name, model), name
 
 
def writeInputsInRange(mdb, model, dxs, dzs, drs, patch_statuses):
    jobs = {}
    dxmax = max(1e-6,max(dxs))
    dzmax = max(1e-6,max(dzs))
    drmax = max(1e-6,max(drs))
    print(dxmax, dzmax, drmax)
    for with_patch in patch_statuses:
        for dx in dxs:
            for dz in dzs:
                for dr in drs:
                    if dx != dxmax and dz != dzmax and dr != drmax:
                        continue
                    job, name = jobWithLoading(mdb, model, dx, dz, dr, with_patch)
                    job.writeInput()
                    job = mdb.JobFromInputFile(name, name + ".inp")
                    jobs[name] = (job, dx, dz, dr, with_patch)
                    print(jobs[name])
    return jobs
 
 
 
def fieldComponent(fieldOutputs, fieldOutput, region, getComponent):
    return np.array(np.nan) if not region else np.array([
        getComponent(value)
        for value in fieldOutputs[fieldOutput].getSubset(region=region).values
    ])
 
 
def volumeMean(fieldOutputs, fieldOutput, region, getComponent):
    fo = fieldComponent(fieldOutputs, fieldOutput, region, getComponent)
    return np.mean(fo)
    # print("FO:" + str(fo.shape))
    # vol = fieldComponent(fieldOutputs, "EVOL", region, lambda v: v.data)
    # print("VOL:" + str(vol.shape))
    # return np.average(fo, weights=vol)
 
 
def getStatsForFrame(frame, assembly, stepLoading):
    hasPatch = stepLoading[3] != 0 # means we don't have to worry about weird float issues
    frameLoading = stepLoading[0:3] * frame.frameValue
    fo = frame.fieldOutputs
    # instances
    instances = assembly.instances
    crack_proximity = instances["SURFACE-TOP"].elementSets["CRACK-ANALYSIS"]
    crack = instances["SURFACE-TOP"].elementSets["SURFACE-WEAK"]
    surfaceElem = instances["SURFACE-TOP"].elementSets["PATCH-NO-CRACK"]
    surfaceNode = instances["SURFACE-TOP"].nodeSets["PATCH-NO-CRACK"]
    adhesive = None if ~hasPatch else instances["ADHESIVE-LAYER-1"]
    patch = None if ~hasPatch else instances["PATCH-1"]
    stretch = assembly.nodeSets["STRETCH"]
    # surface strains n rotation
    sE11 =  volumeMean(fo, "LE", surfaceElem, lambda v: v.data[0])
    sE12 =  volumeMean(fo, "LE", surfaceElem, lambda v: v.data[3])
    sUR =  np.mean(fieldComponent(fo, "UR", surfaceNode, lambda v: v.data[0]))
    patchStressMax = np.max(fieldComponent(fo, "S", patch, lambda v: v.mises))
    surfaceStressMax = np.max(fieldComponent(fo, "S", surfaceElem, lambda v: v.mises))
    crackStressMax = np.max(fieldComponent(fo, "S", crack_proximity, lambda v: v.mises))
    crackStrainE11 = np.max(fieldComponent(fo, "LE", crack, lambda v: v.data[0]))
    crackStrainEMax = np.max(fieldComponent(fo, "LE", crack, lambda v: v.maxInPlanePrincipal))
    crackStrainE22 = np.max(fieldComponent(fo, "LE", crack, lambda v: v.data[1]))
    crackStrainE12 = np.max(fieldComponent(fo, "LE", crack, lambda v: v.data[3]))
    adhesiveDegradationMax = np.max(fieldComponent(fo, "SDEG", adhesive, lambda v: v.data))
    adhesiveDegradationAvg= np.mean(fieldComponent(fo, "SDEG", adhesive, lambda v: v.data))
    sRFX = fieldComponent(fo, "RF", stretch, lambda v: v.data[0])
    sRFY = fieldComponent(fo, "RF", stretch, lambda v: v.data[1])
    sRFZ = fieldComponent(fo, "RF", stretch, lambda v: v.data[2])
    sPY =  fieldComponent(fo, "COORD", stretch, lambda v: v.data[1]) + fieldComponent(fo, "U", stretch, lambda v: v.data[1])
    sPZ =  fieldComponent(fo, "COORD", stretch, lambda v: v.data[2]) + fieldComponent(fo, "U", stretch, lambda v: v.data[2])
    stretchRMX = sPZ * sRFY - sPY * sRFZ
    # hashin it
    hashinFC = np.max(fieldComponent(fo, "HSNFCCRT", patch, lambda v: v.data))
    hashinFT = np.max(fieldComponent(fo, "HSNFTCRT", patch, lambda v: v.data))
    hashinMC = np.max(fieldComponent(fo, "HSNMCCRT", patch, lambda v: v.data))
    hashinMT = np.max(fieldComponent(fo, "HSNMTCRT", patch, lambda v: v.data))
    hashinModes = [hashinFC, hashinFT, hashinMC, hashinMT]
    hashinMax = max(hashinModes)
    return {
        "t"  : frame.frameValue,
        "hasPatch": hasPatch,
        "dx" : frameLoading[0],
        "dz" : frameLoading[1],
        "dr" : frameLoading[2], 
        "E11": sE11,
        "E12": sE12,
        "UR1": sUR,
        "RFX": np.sum(sRFX),
        "RFY": np.sum(sRFY),
        "RFZ": np.sum(sRFZ),
        "RMX": np.sum(stretchRMX),
        "AdhDMax": adhesiveDegradationMax,
        "AdhDAvg": adhesiveDegradationAvg,
        "PtchSMax": patchStressMax,
        "SurfSMax": surfaceStressMax,
        "CrckSMax": crackStressMax,
        "CrckE11": crackStrainE11,
        "CrckE12": crackStrainE12,
        "CrckE22": crackStrainE22,
        "CrckEMax": crackStrainEMax,
        "HshnMax": hashinMax,
        "HshnType": hashinModes.index(hashinMax),
    }
 
 
def getStatsForStep(step, assembly, stepLoading):
    return[getStatsForFrame(frame, assembly, stepLoading) for frame in step.frames]
 
 
def dictToCsvString(data: list[dict[str,Any]]):
    return "\n".join([", ".join([f"{value}" for value in data[0].keys()])] + \
              [", ".join([f"{value:+.8e}" for value in row.values()]) for row in data])
 
 
# jobs = writeInputsInRange(
#     mdb, model, np.linspace(0, 0.7, 8), np.linspace(0, 1, 8), np.linspace(0, pi/18, 8), [True, False]
# )
jobs = writeInputsInRange(
    mdb, model, np.linspace(0, 0.5, 9), np.linspace(0, 0.9, 9), np.linspace(0, pi/18, 9), [True, False]
)
 
 
# uniaxial cases:
if(False):
    jobs = {}
    jobs.update(writeInputsInRange(mdb, model, [1], [0], [0], [False, True]))
    jobs.update(writeInputsInRange(mdb, model, [0], [1], [0], [False, True]))
    jobs.update(writeInputsInRange(mdb, model, [0], [0], [pi/12], [False, True]))
 
 
def writeOdbStats(jobName, job, batchFilename):
    count_errored = 0
    odb = None
    try:
        odb = openOdb("./" + jobName + ".odb", readOnly=True)
        step = odb.steps["Step-1"]
        loading = np.array(jobs[jobName][1:5])
        csvstring = dictToCsvString(getStatsForStep(step, odb.rootAssembly, loading))
    except Exception:
        count_errored+=1
        extxt = traceback.format_exc()
        print(extxt)
        if odb: odb.close()
        return 1
        raise
    with open("../csvs2/" + jobName + "-data.csv", "w") as file:
        file.write(csvstring)
    if batchFilename:
        with open("../csvs2/" + batchFilename + ".csv", "a") as file:
            file.write("\n"+csvstring)
    with open("../csvs2/collated-data.csv", "a") as file:
        file.write("\n"+csvstring)
    if odb: odb.close()
    return 0
 
 
n = 6
count_total = len(jobs)
count_completed = 0
count_errored = 0
pending_jobs = list(jobs.keys())
stats = {}
batchFilename = "Data-" + datetime.datetime.now().strftime("%Y-%m-%dT%H_%M_%S")
while pending_jobs:
    jobsToComplete: list[tuple[str,Any]] = []
    jobsAdded = 0
    while pending_jobs and jobsAdded < 6:
        jobName = pending_jobs.pop()
        if os.path.exists("../csvs2/" + jobName + "-data.csv"):
            print(f"Skipped submitting {jobName}. data exists")
            continue # don't redo existing jobs
        # todo rewrite
        job = jobs[jobName][0]
        print(f"submitting {jobName}")
        job.submit()
        jobsToComplete.append((jobName,job))
        jobsAdded+=1
    for jobName, job in jobsToComplete:
        job.waitForCompletion()
        count_completed+=1
        odb = None
        count_errored += writeOdbStats(jobName, jobs[jobName], batchFilename)
    print(f"Jobs completed {count_completed}/{count_total} ({count_errored} errors)")
 

Link to original

AbaqusMacrosPy

# -*- coding: mbcs -*-
# Do not delete the following import lines
from abaqus import *
from abaqusConstants import *
import __main__
 
FRAME_FOR_LOADING = {
    "+0_50X+0_00Z+0_00R": 8,
    "+0_00X+0_90Z+0_00R": 9,
    "+0_00X+0_00Z+10_00R": 17,
    "+0_12X+0_90Z+2_50R": 8,
    "+0_50X+0_00Z+10_00R": 9
    }
 
def createNextViewport():
    maxNum = max(map(lambda vp: int(vp.split(" ")[1]), session.viewports.keys()))
    session.Viewport(f"Viewport: {maxNum+1}")
 
def manageViewportCount(desiredCount):
    print(f"Want {desiredCount} viewports and have {len(session.viewports)}")
    sortedKeys = sorted(session.viewports.keys(), key = lambda vp: int(vp.split(" ")[1]))
    if desiredCount < len(session.viewports) and False: # actually don't ever delete viewports
        session.viewports[sortedKeys[0]].makeCurrent()
        print("Removing excess viewports")
        for i, key in enumerate(sortedKeys):
            if i >= desiredCount:
                del session.viewports[key]
    elif desiredCount > len(session.viewports):
        print("Adding needed viewports")
        for i in range(len(session.viewports), desiredCount):
            createNextViewport()
    # re-sort
    sortedKeys = sorted(session.viewports.keys(), key = lambda vp: int(vp.split(" ")[1]))
    vps = {i:session.viewports[i] for i in sortedKeys}
    return vps
 
def ShowFieldOutput(viewportSettings):
    
    vps = manageViewportCount(len(viewportSettings))
    usedViewports = [];
                     
    for settings, viewport in zip(viewportSettings,vps.values()):
        viewport.view.setProjection(projection=PARALLEL)
        viewport.odbDisplay.contourOptions.setValues(maxAutoCompute=ON, minAutoCompute=ON)
        viewport.makeCurrent()
        viewport.odbDisplay.contourOptions.setValues(
        colorByMatchingPlies=ON if settings[3] else OFF)
        if settings[1]:
            viewport.odbDisplay.setPrimaryVariable(
                variableLabel=settings[0],
                outputPosition=INTEGRATION_POINT,
                refinement=settings[1]
                )
        else:
            viewport.odbDisplay.setPrimaryVariable(
                variableLabel=settings[0],
                outputPosition=INTEGRATION_POINT
                )
        
        viewport.odbDisplay.basicOptions.setValues(
            sectionResults=USE_ENVELOPE,
            envelopeCriteria=settings[2]
            ) 
        viewport.odbDisplay.display.setValues(plotState=(CONTOURS_ON_UNDEF,))
        usedViewports.append(viewport)
        
    
    # get mins and maxes
    legendGroups = set([settings[4] for settings in viewportSettings])
    for lg in legendGroups:
        if lg is None: continue
        connectedVPsIx = [i for i,settings in enumerate(viewportSettings) if settings[4] == lg]
        connectedVPs = [usedViewports[i] for i in connectedVPsIx]
        lgMin = 1e9
        lgMax = -1e9
        for vp in connectedVPs:
            lgMax = max(lgMax, vp.odbDisplay.contourOptions.autoMaxValue)
            lgMin = min(lgMax, vp.odbDisplay.contourOptions.autoMinValue)
        # now that we have min/max set to that
        for vp in connectedVPs:
            vp.odbDisplay.contourOptions.setValues(
                maxAutoCompute=OFF, maxValue=lgMax, minAutoCompute=OFF, minValue=lgMin)
            
    
    return usedViewports
        
 
def FOStress():
    viewportSettings = [
        ("S", (INVARIANT, "Mises"), MAX_VALUE, False, 0),
        ("S", (INVARIANT, "Mises"), MIN_VALUE, False, 0),
        ("S", (COMPONENT, "S11"), MAX_VALUE, False, 1),
        ("S", (COMPONENT, "S11"), MIN_VALUE, False, 1),
        ("S", (COMPONENT, "S22"), MAX_VALUE, False, 2),
        ("S", (COMPONENT, "S22"), MIN_VALUE, False, 2),
        ("S", (COMPONENT, "S12"), MAX_VALUE, False, 3),
        ("S", (COMPONENT, "S12"), MIN_VALUE, False, 3)
        ]
    return ShowFieldOutput(viewportSettings)
 
    
def FOHashin():
    viewportSettings = [
        ("HSNFCCRT", None, MAX_VALUE, False, None),
        ("HSNFCCRT", None, MIN_VALUE, True, None),
        ("HSNFTCRT", None, MAX_VALUE, False, None),
        ("HSNFTCRT", None, MIN_VALUE, True, None),
        ("HSNMCCRT", None, MAX_VALUE, False, None),
        ("HSNMCCRT", None, MIN_VALUE, True, None),
        ("HSNMTCRT", None, MAX_VALUE, False, None),
        ("HSNMTCRT", None, MIN_VALUE, True, None)
        ]
    return ShowFieldOutput(viewportSettings)
    
    
def FOAdhesive():
    viewportSettings = [
        ("SDEG", None, MAX_VALUE, False, None),
        ("S", (INVARIANT, "Mises"), MAX_VALUE, False, None)
        ]
    return ShowFieldOutput(viewportSettings)
    
    
def ViewSettings():
    import connectorBehavior
    viewport = session.viewports.values()[0]
    viewport.viewportAnnotationOptions.setValues(title=OFF, 
        state=OFF, compass=OFF)
    viewport.viewportAnnotationOptions.setValues(
        triadPosition=(95, 90))
    viewport.viewportAnnotationOptions.setValues(
        legendFont='-*-verdana-medium-r-normal-*-*-70-*-*-p-*-*-*')
 
 
def StressSurface():
    session.viewports[session.currentViewportName].odbDisplay.contourOptions.setValues(
        maxAutoCompute=OFF, maxValue=345, minAutoCompute=OFF, minValue=-345)
        
        
def StressPatch():
    session.viewports[session.currentViewportName].odbDisplay.contourOptions.setValues(
        maxAutoCompute=OFF, maxValue=400, minAutoCompute=OFF, minValue=-400)
 
 
def ContourColors():
    session.viewports[session.currentViewportName].odbDisplay.contourOptions.setValues(
        outsideLimitsAboveColor='#FF00F8', outsideLimitsBelowColor='#6C008A')
 
 
def SaveToImage(imageInfo, viewports):
    size = (3840, 2010)# (2560, 1340)
    session.printOptions.setValues(vpDecorations=OFF, reduceColors=False)
    session.pngOptions.setValues(imageSize = size)
    odbName = session.viewports.values()[0].displayedObject.name.split("/")[-1].strip(".odb");
    loading = odbName.strip("Patch").strip("less");
    hasPatch = not "less" in odbName
    fileName = loading + "-" + imageInfo
    if imageInfo == "SurfaceStress":
        fileName +=  "-Patch" if hasPatch else "-Patchless"
    filePatch = 'C:/Users/S3943498/OneDrive - RMIT University/Documents/Capstone personal/Patch/images/' + fileName
    session.printToFile(
        fileName=filePatch, 
        format=PNG, canvasObjects=tuple(viewports))
    print(f"Saving {imageInfo} to ./images/{fileName}")
 
 
def ShowSurface():
    import displayGroupOdbToolset as dgo
    leaf = dgo.LeafFromElementSets(elementSets=('SURFACE-TOP.REMESH', ))
    session.viewports[session.currentViewportName].odbDisplay.displayGroup.replace(leaf=leaf)
    session.viewports[session.currentViewportName].view.setValues(session.views['User-1'])
 
 
 
def ShowPatch():
    import displayGroupOdbToolset as dgo
    leaf = dgo.LeafFromPartInstance(partInstanceName=('PATCH-1', ))
    session.viewports[session.currentViewportName].odbDisplay.displayGroup.replace(leaf=leaf)
    session.viewports[session.currentViewportName].view.setValues(session.views['User-3'])
 
 
def ShowAdhesive():
    import displayGroupOdbToolset as dgo
    leaf = dgo.LeafFromPartInstance(partInstanceName=('ADHESIVE-LAYER-1', ))
    session.viewports[session.currentViewportName].odbDisplay.displayGroup.replace(leaf=leaf)
    session.viewports[session.currentViewportName].view.setValues(session.views['User-2'])
 
    
def APatchStress():
    ShowPatch()
    SaveToImage("PatchStress", FOStress())
    
    
def APatchHashin():
    ShowPatch()
    SaveToImage("PatchHashin", FOHashin())
    
    
def ASurfaceStress():
    ShowSurface()
    SaveToImage("SurfaceStress", FOStress())
    
def AAdhesiveAll():
    ShowAdhesive()
    SaveToImage("AdhesiveAll", FOAdhesive())
 
def AAAGetImagesForOdb(odb):
    odbName = odb.name.split("/")[-1].strip(".odb")
    loading = odbName.strip("Patch").strip("less")
    # get frame for this odb
    if loading in FRAME_FOR_LOADING:
        frameIx = FRAME_FOR_LOADING[loading]
        frame = odb.steps['Step-1'].frames[frameIx]
    else:
        frameIx = len(odb.steps['Step-1'].frames)-1
        
    print(odbName + " : " + loading + ", Frame " + str(frameIx))
    for viewport in session.viewports.values():
        viewport.setValues(displayedObject = odb)
        viewport.odbDisplay.setFrame(step='Step-1', frame=frameIx) 
    manageViewportCount(8)
    if "less" in odbName:
        # patchless, only do surface
        ASurfaceStress()
    else:
        ASurfaceStress()
        APatchHashin()
        APatchStress() # have to do stress before changing element set
        AAdhesiveAll()
    
def AAAGetImagesForViewportOdb():
    odb = session.viewports[session.currentViewportName].displayedObject
    AAAGetImagesForOdb(odb)
 
 
def AAAAGetImagesForLoadedOdbs():
    for odb in session.odbs.values():
        AAAGetImagesForOdb(odb)
 
 
def ZZImports():
    import section
    import regionToolset
    import displayGroupMdbToolset as dgm
    import part
    import material
    import assembly
    import step
    import interaction
    import load
    import mesh
    import optimization
    import job
    import sketch
    import visualization
    import xyPlot
    import displayGroupOdbToolset as dgo
    import connectorBehavior
    pass
 

Link to original