##############################################################
## Get nodal stress from a deformed state
## https://anning.me/get-nodal-stress/
##############################################################

jobName = '''Name of Job'''
stepName = '''Name of Step'''
outputSetName = '''Name of Output Set'''
#
from odbAccess import*
from abaqusConstants import*
import string
import numpy as np
import os
#
odb = openOdb(path = jobName+'.odb')

outfile = open(jobName + '.csv', 'w')
outfile.write('Node No.' + ',' + 'NStress' + '\n')

#
fm = -1
timeFrame = odb.steps[stepName].frames[fm]
readNode = odb.rootAssembly.nodeSets[outputSetName]
Stress = timeFrame.fieldOutputs['S']  # Remember to set field outputs manually
readNodeStress = Stress.getSubset(position=ELEMENT_NODAL, region=readNode)
readNodeStressValue = readNodeStress.bulkDataBlocks[0].data
readNodeLabels = readNodeStress.bulkDataBlocks[0].nodeLabels

NodeLabels_unique, unq_idx = np.unique(readNodeLabels, return_inverse=True)
Values_Averaged=np.zeros((NodeLabels_unique.size,readNodeStressValue.shape[1]))
unq_counts = np.bincount(unq_idx)

for i in xrange(0, readNodeStressValue.shape[1]):
    ValuesTemp = [item[i] for item in readNodeStressValue]
    unq_sum = np.bincount(unq_idx, weights=ValuesTemp)
    Values_Averaged[:,i] = unq_sum / unq_counts

BoundaryStress = np.zeros(Values_Averaged.shape[0])
for i in range(0, Values_Averaged.shape[0]):
    BoundaryStress[i] = Values_Averaged[i, 0]
    outfile.write(str(i) + ',' + str(BoundaryStress[i]) + ',' + '\n')
outfile.close()
#

SDeviation = np.std(BoundaryStress)
SMax = np.amax(BoundaryStress)
f = open(jobName + '_output.dat', 'w')
f.write(str(SMax))
f.write('\n')
f.write(str(SDeviation))
f.close()


odb.close()