Abaqus python reading stresses from a composite shell

Hi everyone,

i have a question regarding reading stresses from an odb. I have a shell model with an composite layUp (e.g. 4 plies) and i try to do the post-procssing withPython . I want to extract the stresses from every ply of the composite laminate. Then i want to calculate a stress interaction criterion and save a this as a new field output variable in the odb. I already achieved it for a whole element with the code below.

Now i have to find a way to get the stresses from the section points of the shell element and then do the same calculation and save it again in the odb. 

Does anybody now, how to do it?

Thank you very much in advance

from odbAccess import *
from abaqusConstants import *
import visualization

# Parameters
name_job = 'Test'         # Name of your job (without .odb extension)
step_name = 'Step-1'      # Step name from your Abaqus analysis
name_set = 'ESET'         # Node set name (ensure this exists in the ODB)
odb_var = 'S'             # Field variable for stress (use 'S' for stress)

R11=2200
R12=120   
    
    

# Open the ODB file
odb = openOdb(path=name_job + '.odb', readOnly=True)

# Access the instance and step in the ODB
instance1 = odb.rootAssembly.instances['PART-1-1']  # Modify 'PART-1-1' if necessary
steps = odb.steps[step_name]
lastFrame = odb.steps[step_name].frames[-1]  # Get the last frame of the step

# Get the region of the node set
region_Auswertung=odb.rootAssembly.elementSets[name_set]
# Check if the requested field exists
if odb_var not in lastFrame.fieldOutputs.keys():
    print "Error: Field output '{}' not found in the last frame.".format(odb_var)


# Extract the stress field output
print "Reading Stress Output for node set '{}'...".format(name_set)
var_output = lastFrame.fieldOutputs[odb_var]
get_elements = var_output.getSubset(region=region_Auswertung)
var_value = get_elements.values

# Check if any stress values are available
if not var_value:
    print "No stress data found for the element set '{}'.".format(name_set)
label_list=[]    
result_list = []
# Iterate over the stress values and print the results
for myStress in var_value:
    element_label = myStress.elementLabel
    print element_label
    stress_components = myStress.data  # This is a tuple of [S11, S22, S33, S12, S13, S23] for 3D
    
    # Print the stress components for each node
    if len(stress_components) == 6:
        S11, S22, S33, S12, S13, S23 = stress_components
    elif len(stress_components) == 4:  # In case of 2D analysis
        S11, S22, S12, S33 = stress_components  # Assuming 2D (modify based on your case)
        print S11, S22, S33, S12
        S13 = S23 = 0  # In 2D, these stress components might be zero
    if S11>0:    
        ya_sun = (((S11/R11)**2)+((abs(S12)/R12)**2))**0.5 
    else: 
        pass
    print 'ya_sun:', ya_sun
    label_list.append(element_label)
    result_list.append((ya_sun,))
    # except KeyError:  
    #     print 'Calculation failed'
    #     pass

try:
    NewField = lastFrame.FieldOutput(
        name='New Var',
        description='Our brand new scalar field',
        type=SCALAR)
    NewField.addData(position=WHOLE_ELEMENT, instance=instance1, labels=label_list, data=result_list)

except Exception as e:
    print 'Failed to create new FieldOutput:', e