python script to get element-based variables(e.g. stress) from nodal points

I already checked this topic out in Simulia community, but haven't had a clear clue in regarding how to change my script to get stress from nodes, not from the integration points. 

Here is my script. It saves "element ID, integration point (x,y,z) and Principal Max, Min and Med" to csv file for postprocessing. 

This is what I was trying to:

1. When try with (position=ELEMENT_NODAL), it looks like it shows stresses at 8 nodes per each element (I use C3D8 element), which exceeds the total number of nodes in my model. This maybe due to that the stress values are not averaged at each nodes. 

2. When try with (position=UNIQUE_NODAL), it shows nothing. Not sure why. 

3. I'd like to get averaged, unique stress at each node, which needs to add averaging algorithm. or is there any way to do so, not using averaging technique? 

Can anyone give me a good idea and show how to change my script? 

Thanks,

-------------------------------------------------------------------------------------------------

My edited script to get averaged stresses at nodes. I tested. works good.

-------------------------------------------------------------------------------------------------

##import abaqus module
from abaqus import *
from abaqusConstants import *

##open myodb file
from odbAccess import openOdb
myodb=openOdb('odbpath')

##define stress variable from fieldOutputs repository (Store stress values at all elements from the last frame of the step 'Step-1') 
##use getSubset function to extract stress values at the integration points
S=myodb.steps['Step-1'].frames[-1].fieldOutputs['S'].getSubset(position=NODAL)

resultsFile=open('csvpath','w')
resultsFile.write('%s, %s, %s, %s, %s, %s, %s\n' % ('NodeID', 'coord[x]', 'coord[y]', 'coord[z]', 'Smax(MPa)', 'Smin(MPa)', 'Smid(MPa)'))        

##restrict the stress and coordinate info to a particular element set
myinstance= myodb.rootAssembly.instances['PART-4-1']    
elset= myinstance.nodeSets['SET-1']    
stress=S.getSubset(region=elset)

##print total number of nodes and elements
numnodes=len(myinstance.nodes)
numelements=len(myinstance.elements)
print numnodes, numelements

##print node number and coordinates
for node in myinstance.nodes:    
        resultsFile.write('%s, %s, %s, %s\n' % (node.label, node.coordinates[0], node.coordinates[1], node.coordinates[2]))        
        
##print stresses at nodes
for stressvalue in stress.values:    
      smax=stressvalue.minPrincipal  
      smin=stressvalue.maxPrincipal
      smid=stressvalue.midPrincipal    
      resultsFile.write('%10.4E, %10.4E, %10.4E\n' % (smax/1.0E6, smin/1.0E6, smid/1.0E6))    
resultsFile.close()