# -*- coding: mbcs -*-
# python script for model analysis related to "Incorporating pathological gait into patient‑specific finite element models of the haemophilic ankle"
# Used and tested with Abaqus 2017.
# Mengoni - 2024
# ===================================
# Script working process
    # 1. scan the content of folder ODBFilesFolderPath for all existing .odb files (the baseline models)
    # 2. for each loading step in each odb, extract contact outputs on the cartilage-to-cartilage interaction and cvon Mises stress in the bones
    # 3. writes outputs into csv file for each odb (total contact area, total contact force, max contact pressure, mean and max stress in tibia and talus)
# ===================================

# ==== USER DEFINED INFO ===
ODBFilesFolderPath = os.getcwd() # path to folder with baseline input files
# ==========================


from abaqus import *
from abaqusConstants import *

import os, fnmatch

for root, dirnames, filenames in os.walk(ODBFilesFolderPath):
    for jobName in fnmatch.filter(filenames, '*.odb'):
        csvName = jobName.split('.')[0] + '_outputs.csv'
        with open(csvName, 'w') as masterFile:
            masterFile.write('step name,area,force,max contact pressure,tibia mean stress,tibia max stress,talus mean stress,talus max stress\n')        
        thisOdb = session.openOdb(name=jobName)
        for stepName in thisOdb.steps.keys():
            if stepName.startswith('move') or stepName.startswith('offload'): continue
            wholeModel = thisOdb.steps[stepName].historyRegions['NodeSet  Z000001']
            contactArea = wholeModel.historyOutputs['CAREA    ASSEMBLY_TALARCARTILAGE/ASSEMBLY_TIBIALCARTILAGE'].data[-1][1]
            contactForce = wholeModel.historyOutputs['CFNM     ASSEMBLY_TALARCARTILAGE/ASSEMBLY_TIBIALCARTILAGE'].data[-1][1]
            contactPressure = thisOdb.steps[stepName].frames[-1].fieldOutputs['CPRESS   ASSEMBLY_TALARCARTILAGE/ASSEMBLY_TIBIALCARTILAGE-1'].values

            tibia = thisOdb.rootAssembly.instances['PART-1-1'].elementSets['PT_TIBIA']
            talus = thisOdb.rootAssembly.instances['PART-1-1'].elementSets['PT_TALUS']
            maxPress = max(value.data for value in contactPressure)
            
            allStress = thisOdb.steps[stepName].frames[-1].fieldOutputs['S'].getScalarField(invariant=MISES)
            tibiaStress = [ptValue.data for ptValue in allStress.getSubset(region=tibia).values]
            talusStress = [ptValue.data for ptValue in allStress.getSubset(region=talus).values]
           
            toWrite = "%s,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f\n"%(stepName,contactArea,contactForce,maxPress,sum(tibiaStress)/len(tibiaStress),max(tibiaStress),sum(talusStress)/len(talusStress),max(talusStress))
            with open(csvName, 'a') as masterFile:
                masterFile.write(toWrite)