from abaqus import *
from abaqusConstants import *
from caeModules import *
import regionToolset
import visualization
from odbAccess import *
import numpy as np
import re
import math
StudentV= False
## Condizioni generali
Simmetria_laminato = False      # Attivare se si vuole laminato simmetrico rispetto al piano medio 
Open_Hole = False               # Attivare se apertura lungo la piastra, genera SOLO UN QUARTO della geometria
Open_Hole_completo = False            # Attivare se apertura lungo la piastra, genera la geometria COMPLETA

# Dimensioni modello
Modello_3D = False             # Attivare per analisi tridimensionale; Si caratterizza nella riga 180-ish
Modello_2D = True              # Attivare per analisi tridimensionale; Si caratterizza alla riga 556-ish

# Tipologia analisi. Attivarne solo una alla volta
Statica = False                  # Attivare per analisi statica 
Frequenze_Naturali = False       # Attivare per analisi delle frequenze naturali 
Buckling_a_compressione = True  # Attivare per analisi di buckling; nella riga 688-ish si possono cambiare i lati ai quali sono applicate le BC
Trazione = False                  # nella riga 123 si possono cambiare i lati ai quali sono applicate le BC

# Parametrizzazione nomi 
if Modello_3D:
    if Open_Hole or Open_Hole_completo:
        part_name='Strati_Fin'   # Nome della parte
        Inp_name='Job-Fin'       # Nome del file dal quale verranno lette le informazioni
        Job_name = 'Job-Fin3D'   # Nome del file lanciato in analisi
    else:
        part_name='square_plate_part_3D'
        Inp_name='Job-3D'
        Job_name = 'Job-Plate3D'

elif Modello_2D:
    if Open_Hole or Open_Hole_completo:
        part_name='Open_Hole_2D'
        Inp_name='Job_Fin-2D'
        Job_name = 'Job-Fin2D'
    else:
        part_name='square_plate_part_2D'
        Inp_name='Job-2D'

# PARAMETRI Geometrici
L1 = 0.3 #0.3048   # 0.6096/2   # dimesione x della lamina 
L2 = 0.5 # 0.9144   # 1.8288/2   # dimensione y della lamina
L3 = 0.05 #0.0009525 = 2 strati    #0.00047625 = 4 strati    0.000238125 = 8 strati # spessore singolo layer    
Ry = 0.2286   # Raggio semiasse x
Rx = 0.1397   # Raggio semiasse y
              # se si vuole un cerchio porre Rx=Ry

Material_name = 'Graphite-Epoxy'
# Material_name = 'Toray'
# costanti ingegneristiche, in Metri e Pascal; 
# E11 = 1.48e+12 
# E22 = 8.48e+09
# E33 = 8.48e+09
# G12 = 3.94e+09
# G13 = 3.94e+09
# G23 = 3.235e+09
# v12 = 0.329
E11 = 137.9e+09     # Modulo di Young
E22 = 8.96e+09
E33 = 8.96e+09
G12 = 7.1e+09       # Modulo di taglio
G13 = 7.1e+09
G23 = 6.21e+09
v12 = 0.3           # Coeff. Poisson
v13 = 0.3
# v23 = 0.32
v23 = 0.49
Density=1540
pressure = 10000             # valore positivo perche' applicata alla top surface e Abaqus in automatico la considera in compressione
Carico_Trazione = -100.0   # Netwon

# PARAMETRI MESH
n = 40 
nz = 16      # Numero di elementi finiti per direzione y
nc= 14         # Numero di elementi finiti ai bordi dell'apertura
only_S3 = False # Tipi di elementi  in analisi bidimensionale con lamina che presenta apertura

# VALORI THETA           #IMPORTANTE METTERE GLI ANGOLI CON VALORI DECIMALI
theta = [0.0, 45.0, 90.0, 45.0] # 90.0, 15.0, 60.0, 0.0, 65.0, 10.0, 60.0, 0.0, 65.0, 10.0, 60.0, 0.0] #, 90.0, 75.0, 30.0, 60.0 ]
# Esempio: theta = [ 0.0, 0.0, 45.0, 45.0, -45.0, -45.0, 90.0, 90.0 ] 
#theta_values = [theta0_first_ply, theta1_first_ply, theta0_second_ply, theta1_second_ply, theta0_third_ply, theta1_third_ply, ecc,  ]
phi=0

Mdb()  # crea un nuovo database
 
#Impostazioni Grafiche per sfondo
session.graphicsOptions.setValues(backgroundStyle=SOLID, backgroundColor='#FFFFFF')

if Statica:
    session.viewports['Viewport: 1'].viewportAnnotationOptions.setValues(title=OFF, state=OFF, compass=OFF)
else:
    session.viewports['Viewport: 1'].viewportAnnotationOptions.setValues(title=OFF, compass=OFF)         # lascio lo state per vedere i modi

# impostazioni iniziali, il sheetSize is utile nel sketch mode
session.viewports['Viewport: 1'].setValues(displayedObject=None)
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=100.0)

# creazione GEOMETRIA per coordinate parametrizzate
if not Open_Hole:
    P1 = 0.0, 0.0       # Punto
    P2 = L1, 0.0
    P3 = L1, L2
    P4 = 0.0, L2

    g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
    s.setPrimaryObject(option=STANDALONE)
    # prima linea
    s.Line(point1=P1, point2=P2)
    s.HorizontalConstraint(entity=g[2], addUndoState=False)
    # seconda linea
    s.Line(point1=P2, point2=P3)
    s.VerticalConstraint(entity=g[3], addUndoState=False)
    s.PerpendicularConstraint(entity1=g[2], entity2=g[3], addUndoState=False)
    # terza linea
    s.Line(point1=P3, point2=P4)
    s.HorizontalConstraint(entity=g[4], addUndoState=False)
    s.PerpendicularConstraint(entity1=g[3], entity2=g[4], addUndoState=False)
    # quarta linea
    s.Line(point1=P4, point2=P1)
    s.VerticalConstraint(entity=g[5], addUndoState=False)
    s.PerpendicularConstraint(entity1=g[4], entity2=g[5], addUndoState=False)
    if Open_Hole_completo:
        Cx = L1/2      # Coordinata X del centro dell' apertura
        Cy = L2/2      # Coordinata Y del centro dell' apertura
        s.EllipseByCenterPerimeter(center=(Cx, Cy), axisPoint1=(Cx+Rx, Cy), axisPoint2=(Cx,Cy+ Ry))
        print("Entro qui")    
# Creazione dell'apertura utilizzando le coordinate parametriche
if Open_Hole:
    P1 = Rx, 0.0       # Punto
    P2 = L1, 0.0
    P3 = L1, L2
    P4 = 0.0, L2
    P5 = 0.0, Ry

    g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
    s.setPrimaryObject(option=STANDALONE)
    # prima linea
    s.Line(point1=P1, point2=P2)
    s.HorizontalConstraint(entity=g[2], addUndoState=False)
    # seconda linea
    s.Line(point1=P2, point2=P3)
    s.VerticalConstraint(entity=g[3], addUndoState=False)
    s.PerpendicularConstraint(entity1=g[2], entity2=g[3], addUndoState=False)
    # terza linea
    s.Line(point1=P3, point2=P4)
    s.HorizontalConstraint(entity=g[4], addUndoState=False)
    s.PerpendicularConstraint(entity1=g[3], entity2=g[4], addUndoState=False)
    # quarta linea
    s.Line(point1=P4, point2=P5)
    s.VerticalConstraint(entity=g[5], addUndoState=False)
    s.PerpendicularConstraint(entity1=g[4], entity2=g[5], addUndoState=False)
    # s.CircleByCenterPerimeter(center=(Cx, Cy), point1=(Cx - Raggio, Cy))
    s.EllipseByCenterPerimeter(center=(0.0, 0.0), axisPoint1=(Rx, 0.0), axisPoint2=(0.0, Ry))
    s.autoTrimCurve(curve1=g[6], point1=(-Rx, 0.0))        # Calcella la parte di Ellisse in eccesso

# Creazione materiale
mdb.models['Model-1'].Material(name=Material_name)
mdb.models['Model-1'].materials[Material_name].Elastic(
    type=ENGINEERING_CONSTANTS, table=((E11, E22, E33, v12, v13, v23, G12, G13, G23), ))
if Frequenze_Naturali:
    mdb.models['Model-1'].materials[Material_name].Density(table=((Density, ), ))

if Modello_2D and Modello_3D:
    Modello_2D=False
    Modello_3D=False
    print('Both Model 2D and 3D are selected')

if Simmetria_laminato:      # Specchia i valori del vettore teta
    if len(theta)%2==0:
        switched=[]
        for i in range(0, len(theta), 2)[::-1]:
            switched.append(theta[i])
            switched.append(theta[i+1])
        theta=theta+switched

kk =len(theta) // 2         # Numero di ply
hh = kk * L3                # Altezza del solido in direzione Z
print('Numero strati',kk)

if Modello_3D:  
    # Creazione parte 3D
    p = mdb.models['Model-1'].Part(name='Part-1', dimensionality=THREE_D, type=DEFORMABLE_BODY)
    p = mdb.models['Model-1'].parts['Part-1']
    p.BaseSolidExtrude(sketch=s, depth=L3)               #differenza col metodo 2D, p.BaseShell(sketch=s)    
    s.unsetPrimaryObject()
    del mdb.models['Model-1'].sketches['__profile__']
    
    # Creazione dell'assembly
    a = mdb.models['Model-1'].rootAssembly
    a.DatumCsysByDefault(CARTESIAN)
    p = mdb.models['Model-1'].parts['Part-1']
    a.Instance(name='Part-1-1', part=p, dependent=ON)
    # print(math.pi)
    # Creazione di un pattern lineare:
    a = mdb.models['Model-1'].rootAssembly
    a.LinearInstancePattern(instanceList=('Part-1-1', ), direction1=(
        0.0, 0.0, 1.0), direction2=(0.0, 1.0, 0.0), number1=kk, number2=1, 
        spacing1=L3, spacing2=0.1)    
    
    # Creazione della lista delle istanze
    instance_list = ['Part-1-1']
    for i in range(2, kk+1):
        instance_name = 'Part-1-1-lin-' + str(i) + '-1'
        instance_list.append(instance_name)

    # Unione delle istanze usando operazione booleana
    a.InstanceFromBooleanMerge(name=part_name, 
                            instances=[a.instances[inst] for inst in instance_list], 
                            keepIntersections=ON, 
                            originalInstances=DELETE, 
                            domain=GEOMETRY)  

    # Create step
    mdb.models['Model-1'].StaticStep(name='Step-1', previous='Initial')
    session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='Step-1')
    mdb.models['Model-1'].fieldOutputRequests['F-Output-1'].setValues(
        variables=('S', 'MISES', 'E', 'PE', 'PEEQ', 'PEMAG', 'U', 'RF'))     # da qui si possono modificare gli output dell'analisi
    
    #Crea section per cella
    mdb.models['Model-1'].HomogeneousSolidSection(name='Section-1', 
        material=Material_name, thickness=None)

    # Definizione della maschera di base
    mask_base = '[#{} ]'

    # Creazione e assegnazione della sezione per ogni set
    for i in range(kk):
        # Creazione della maschera per il set corrente
        mask_value = hex(1 << i)[2:]  # Converti il numero in esadecimale per creare la maschera corretta
        mask = mask_base.format(mask_value)
        
        # Selezione delle celle
        p = mdb.models['Model-1'].parts[part_name]
        c = p.cells
        cells = c.getSequenceFromMask(mask=(mask, ), )
        
        # Creazione del set
        set_name = 'Set-'+str(i+1)
        region = p.Set(cells=cells, name=set_name)
        
        # Assegnazione della sezione
        p.SectionAssignment(region=region, sectionName='Section-1', offset=0.0, 
                            offsetType=MIDDLE_SURFACE, offsetField='', 
                            thicknessAssignment=FROM_SECTION)
    
    # Crea Load
    if Open_Hole or Open_Hole_completo:
        if Statica:
            a = mdb.models['Model-1'].rootAssembly
            s1 = a.instances['{}-{}'.format(part_name, 1)].faces
            side1Faces1 = s1.getSequenceFromMask(mask=('[#20 ]', ), )
            region = a.Surface(side1Faces=side1Faces1, name='Surf-1')
            mdb.models['Model-1'].Pressure(name='Load-1', createStepName='Step-1', 
                region=region, distributionType=UNIFORM, field='', magnitude=pressure)  # viene applicata pressione su faccia superiore
            print('Carico creato correttamente')
        if Trazione:
            a = mdb.models['Model-1'].rootAssembly
            s1 = a.instances['{}-{}'.format(part_name, 1)].faces
            num_load = kk #* 2
            load_count = 0 
            for i in range(29):
                mask_value = hex(1 << i)[2:] 
                if Open_Hole:
                    if ((i + 1) % 6 == 3 ): # or ((i + 1) % 6 == 1):  
                        side1Faces1 = s1.getSequenceFromMask(mask=('[#' + mask_value + ' ]', ), )
                        region = regionToolset.Region(side1Faces=side1Faces1)
                        mdb.models['Model-1'].Pressure(name='Load-' + str(load_count + 1), createStepName='Step-1', region=region, distributionType=UNIFORM, field='', magnitude=-pressure, 
                                amplitude=UNSET)
                    else:
                        continue
               
                load_count += 1  # Incrementa il contatore
                if load_count >= num_load:  # Ferma il ciclo se raggiunge il numero di BC
                    break
    else:
        a = mdb.models['Model-1'].rootAssembly
        s1 = a.instances['{}-{}'.format(part_name, 1)].faces
        side1Faces1 = s1.getSequenceFromMask(mask=('[#10 ]', ), )
        region = regionToolset.Region(side1Faces=side1Faces1)
        mdb.models['Model-1'].Pressure(name='Load-1', createStepName='Step-1', 
            region=region, distributionType=UNIFORM, field='', magnitude=pressure, 
            amplitude=UNSET)

    # Crea BC
    if Statica:      
        a = mdb.models['Model-1'].rootAssembly
        f1 = a.instances['{}-{}'.format(part_name, 1)].faces

        # Numero totale di BC da creare
        num_bc = kk * 4  # Imposta il numero di BC

        # Prima parte: Maschere fino a #80000000
        bc_count = 0  # Contatore delle BC create
        for i in range(29):
            mask_value = hex(1 << i)[2:]  # Genera la maschera in esadecimale
            if Open_Hole_completo:
                if ((i + 1) % 6 == 5 ) or ((i + 1) % 6 == 0):  # Salta ogni quinta o sesta maschera
                    continue
            else:
                if (i + 1) % 5 == 0:  # Salta ogni quinta maschera
                    continue
            faces1 = f1.getSequenceFromMask(mask=('[#' + mask_value + ' ]', ), )
            region = regionToolset.Region(faces=faces1)
            mdb.models['Model-1'].EncastreBC(name='BC-' + str(bc_count + 1), createStepName='Step-1', region=region, localCsys=None)
            
            bc_count += 1  # Incrementa il contatore
            if bc_count >= num_bc:  # Ferma il ciclo se raggiunge il numero di BC
                break

        # Seconda parte: Imporre #0 #1 e #0 #2 subito dopo #80000000
        high_bit_masks = ['[#40000000 ]', '[#80000000 ]', '[#0 #1 ]', '[#0 #2 ]']
        for idx, mask_str in enumerate(high_bit_masks):
            if bc_count >= num_bc:  # Ferma il ciclo se raggiunge il numero di BC
                break
            faces1 = f1.getSequenceFromMask(mask=(mask_str, ), )
            region = regionToolset.Region(faces=faces1)
            mdb.models['Model-1'].EncastreBC(name='BC-' + str(bc_count + 1), createStepName='Step-1', region=region, localCsys=None)
            bc_count += 1  # Incrementa il contatore

        # Terza parte: Ripartire da #0 #8 e saltare ogni quinta maschera dopo #0 #8
        high_bit = True  # Flag per ripartire da #0 #8
        for i in range(3, 31):  # Inizia da #0 #8 (corrisponde a i=3 nel nuovo contesto)
            
            if bc_count >= num_bc:  # Ferma il ciclo se raggiunge il numero di BC
                break
            
            if high_bit:
                mask_value = hex(1 << i)[2:]  # Calcola la maschera dopo #0 #2
                mask_str = "[#0 #" + mask_value + " ]"
                
                if Open_Hole_completo:
                    if ((i + 1) % 6 == 5 ) or ((i + 1) % 6 == 0):  # Salta ogni quinta o sesta maschera
                        continue
                else:
                    if (i + 1) % 5 == 0:  # Salta ogni quinta maschera
                        continue

                faces1 = f1.getSequenceFromMask(mask=(mask_str, ), )
                region = regionToolset.Region(faces=faces1)
                mdb.models['Model-1'].EncastreBC(name='BC-' + str(bc_count + 1), createStepName='Step-1', region=region, localCsys=None)
            
                bc_count += 1  # Incrementa il contatore
    if Trazione:
        a = mdb.models['Model-1'].rootAssembly
        f1 = a.instances['{}-{}'.format(part_name, 1)].faces
        num_bc = kk #* 2        
        bc_count = 0  
        for i in range(29):
            mask_value = hex(1 << i)[2:]  # Genera la maschera in esadecimale
            if Open_Hole:
                if ((i + 1) % 6 == 1 ):  # or ((i + 1) % 6 == 4): 
                    faces1 = f1.getSequenceFromMask(mask=('[#' + mask_value + ' ]', ), )
                    region = regionToolset.Region(faces=faces1)
                    mdb.models['Model-1'].DisplacementBC(name='BC-' + str(bc_count + 1), createStepName='Step-1', 
                        region=region, u1=0.0, u2=UNSET, u3=0.0, ur1=0.0, ur2=UNSET, ur3=0.0, 
                        amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
                        localCsys=None)
                else:
                    continue    
                bc_count += 1  # Incrementa il contatore
                if bc_count >= num_bc:  # Ferma il ciclo se raggiunge il numero di BC
                    break 
    
        a = mdb.models['Model-1'].rootAssembly
        f1 = a.instances['{}-{}'.format(part_name, 1)].faces
        num_bc = kk       
        bc_count = 0  
        for i in range(29):
            mask_value = hex(1 << i)[2:]  # Genera la maschera in esadecimale
            if Open_Hole:
                if ((i + 1) % 6 == 4 ):   
                    faces1 = f1.getSequenceFromMask(mask=('[#' + mask_value + ' ]', ), )
                    region = regionToolset.Region(faces=faces1)
                    mdb.models['Model-1'].DisplacementBC(name='BC-' + str(bc_count + kk+1), createStepName='Step-1', 
                        region=region, u1=UNSET, u2=0.0, u3=0.0, ur1=UNSET, ur2=0.0, ur3=0.0, 
                        amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
                        localCsys=None)
                else:
                    continue    
                bc_count += 1  # Incrementa il contatore
                if bc_count >= num_bc:  # Ferma il ciclo se raggiunge il numero di BC
                    break            

    # ## Crea Mesh     
    if Open_Hole or Open_Hole_completo:
        m = [1, 3, 7, 0x1F, 0x3F, 0x7F, 0xFF, 0x1FF, 0x3FF]  # per un massimo di 10 strati
        k = kk - 1  
        mask_value = m[k]
        #print(hex(mask_value)[2:])
        p = mdb.models['Model-1'].parts[part_name]
        c = p.cells
        pickedRegions = c.getSequenceFromMask(mask=('[#' + hex(mask_value)[2:].upper() + ' ]',), ) 
        p.setMeshControls(regions=pickedRegions, technique=SWEEP, algorithm=ADVANCING_FRONT)
        # pickedRegions =(cells, )
        p = mdb.models['Model-1'].parts[part_name]
        p.seedPart(size=L1/n, deviationFactor=0.1, minSizeFactor=0.1)
        elemType1 = mesh.ElemType(elemCode=C3D8R, elemLibrary=STANDARD, 
            kinematicSplit=AVERAGE_STRAIN, secondOrderAccuracy=OFF, 
            hourglassControl=DEFAULT, distortionControl=DEFAULT)
        elemType2 = mesh.ElemType(elemCode=C3D6, elemLibrary=STANDARD)
        elemType3 = mesh.ElemType(elemCode=C3D4, elemLibrary=STANDARD, 
            secondOrderAccuracy=OFF, distortionControl=DEFAULT)
        p = mdb.models['Model-1'].parts[part_name]
        c = p.cells
        cells = c.getSequenceFromMask(mask=('[#' + str(mask_value) + ' ]',), )
        pickedRegions =(cells, )   
        p.setElementType(regions=pickedRegions, elemTypes=(elemType1, elemType2, elemType3))    
        p = mdb.models['Model-1'].parts[part_name]
        p.generateMesh()
    elif StudentV:    
        p = mdb.models['Model-1'].parts[part_name]
        p.seedPart(size=L2/n, deviationFactor=0.1, minSizeFactor=0.1)
        p = mdb.models['Model-1'].parts[part_name]
        p.generateMesh()
    else:

        p = mdb.models['Model-1'].parts[part_name]
        e = p.edges
        pickedEdges = e.getSequenceFromMask(mask=('[#4 ]', ), )
        p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER)
        p = mdb.models['Model-1'].parts[part_name]
        e = p.edges
        pickedEdges = e.getSequenceFromMask(mask=('[#1 ]', ), )
        p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER)
        p = mdb.models['Model-1'].parts[part_name]
        e = p.edges
        pickedEdges = e.getSequenceFromMask(mask=('[#2000 ]', ), )
        p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER)
        p = mdb.models['Model-1'].parts[part_name]
        e = p.edges
        pickedEdges = e.getSequenceFromMask(mask=('[#40 ]', ), )
        p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER)
        p = mdb.models['Model-1'].parts[part_name]
        e = p.edges
        pickedEdges = e.getSequenceFromMask(mask=('[#10 ]', ), )
        p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER)
        p = mdb.models['Model-1'].parts[part_name]
        e = p.edges
        pickedEdges = e.getSequenceFromMask(mask=('[#10000 ]', ), )
        p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER)
        
        p = mdb.models['Model-1'].parts[part_name]
        e = p.edges
        pickedEdges = e.getSequenceFromMask(mask=('[#2 ]', ), )
        p.seedEdgeByNumber(edges=pickedEdges, number=nz, constraint=FINER)
        p = mdb.models['Model-1'].parts[part_name]
        e = p.edges
        pickedEdges = e.getSequenceFromMask(mask=('[#1000 ]', ), )
        p.seedEdgeByNumber(edges=pickedEdges, number=nz, constraint=FINER)
        
        p = mdb.models['Model-1'].parts[part_name]
        p.generateMesh()
        

    # Crea .inp dal quale leggere le informazioni
    mdb.Job(name=Inp_name, model='Model-1', description='', type=ANALYSIS, 
        atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, 
        memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, 
        explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, 
        modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', 
        scratch='', resultsFormat=ODB)
    mdb.jobs[Inp_name].writeInput(consistencyChecking=OFF)

    def get_data_from_inp(filename, dataname, cast=str):
        data = []
        with open(filename, 'r') as file:
            str_file = file.read()
            lines = str_file.split('\n')
            
            index = next((i for i, line in enumerate(lines) if '*'+dataname in line), None) + 1
            while '*' not in lines[index][0]:
                data.append([cast(line.strip()) for line in lines[index].split(',')])
                index += 1
        return data

    matrice_nodi=(get_data_from_inp("{}.inp".format(Inp_name),'Node',lambda x: float(x+'0' if '.' in x else x)))
    print('Nodi',len(matrice_nodi))
    matrice_Assembly= (get_data_from_inp("{}.inp".format(Inp_name),'Element',int))
    # print(matrice_Assembly)
    print('Numero elementi:',len(matrice_Assembly))
    Numero_elementi= len(matrice_Assembly)
    nodi_ele=len(matrice_Assembly[0])-1
    print('Nodi per elemento:',nodi_ele)

    X = [riga[1] for riga in matrice_nodi]
    Y = [riga[2] for riga in matrice_nodi]
    Z = [riga[3] for riga in matrice_nodi]
    # print ("Vettore x:", X)
    # print(Elementi)

    xc = [0] * len(matrice_Assembly)
    yc = [0] * len(matrice_Assembly)
    zc = [0] * len(matrice_Assembly)

    for element_num in range(len(matrice_Assembly)):
        for vertix_num in matrice_Assembly[element_num][1:]:
                 
            xc[element_num] += X[vertix_num-1]
            yc[element_num] += Y[vertix_num-1]
            zc[element_num] += Z[vertix_num-1]

        xc[element_num] /= len(matrice_Assembly[0][1:])
        yc[element_num] /= len(matrice_Assembly[0][1:])
        zc[element_num] /= len(matrice_Assembly[0][1:])
        # print(element_num+1, xc[element_num]) #, il centroide calcolato correttamente
     
    # Imposto centro di riferimento funzione al centro della piastra
    x_coord = [0] * len(matrice_Assembly)
    y_coord = [0] * len(matrice_Assembly)

    for element_num in range(len(matrice_Assembly)):
        if Open_Hole:
            x_coord [element_num]= xc[element_num]
            y_coord [element_num]= yc[element_num]
        else:
            x_coord [element_num]= L1/2-xc[element_num]
            y_coord [element_num]= L2/2-yc[element_num]
        # print(x_coord [element_num])

    def VAT (x, y, phi, T0, T1, L):
        if not Open_Hole:
            return phi + T0 + (T1-T0)/L * np.abs((x * np.cos(np.deg2rad(phi))) + (y * np.sin(np.deg2rad(phi))))          
        else:
            return phi + T0 + (T1 - T0) / L * ((x * np.cos(np.deg2rad(phi))) + (y * np.sin(np.deg2rad(phi))))           # per piastra simmetrica
  
    
    thetaC=[]
    
    for element_num in range(len(matrice_Assembly)):
        LS=int(element_num/(len(matrice_Assembly)/kk))
        # print(element_num, LS)
        if Open_Hole:
                thetaC.append(VAT(x_coord[element_num],yc[element_num],phi,theta[LS*2],theta[(LS*2)+1],L1)) 
        else:
                thetaC.append(VAT(x_coord[element_num],yc[element_num],phi,theta[LS*2],theta[(LS*2)+1],L1/2)) 
        # print("Futuro Field", element_num+1, thetaC[element_num])

    # Ora creo il field su ABAQUS   
    model_name = 'Model-1'
    field_type = SCALAR
    location = ELEMENTS
    data_width = 1
    default_values = (0.0,)  
    
    field_names = []
    elementi =[]
    values = []
   
    for element_num in range(len(matrice_Assembly)):
        elementi.append(int(element_num+1))
        values.append(float(thetaC[element_num]))
    # print("Valori",values)
    # print("Elementi",elementi)
    
    for i in range(kk):
        field_name = 'DiscField-' + str(i + 1)
        print('Creato',field_name)
        ele_s=len(matrice_Assembly)/kk  
        field_val= values[ele_s*i:ele_s*(i+1)]          # (ultimovalorenoncompreso)
        field_ele= elementi[ele_s*i:ele_s*(i+1)]

        mdb.models[model_name].DiscreteField(
                name=field_name,
                description='',
                location=location,
                fieldType=field_type,
                dataWidth=data_width,
                defaultValues=default_values,
                # data=(('', 1, elementi, values),))
                data=(('', 1, field_ele, field_val),))
        field_names.append(field_name)

    # assegnare orientamento alle diverse celle
    # Parametri
    p = mdb.models['Model-1'].parts[part_name]

    # Ciclo per iterare su ciascuno dei kk strati
    for i in range(kk):
        # Genera i nomi dei set e dei campi discreti dinamicamente
        set_name = 'Set-' + str(i + 1)                  # Set-1, Set-2, ..., Set-kk
        disc_field_name = 'DiscField-' + str(i + 1)     # DiscField-1, DiscField-2, ..., DiscField-kk

        # Seleziona la regione corrispondente al set corrente
        region = p.sets[set_name]
        
        # Assegna l'orientamento del materiale
        mdb.models['Model-1'].parts[part_name].MaterialOrientation(
            region=region, 
            orientationType=SYSTEM, 
            axis=AXIS_3, 
            localCsys=None, 
            fieldName='', 
            additionalRotationType=ROTATION_FIELD, 
            angle=0.0, 
            additionalRotationField=disc_field_name,  # Nome del campo discreto dinamico
            stackDirection=STACK_3
        )
        print('Associato orientamento allo al i-esimo strato')

    if Statica:
        mdb.Job(name=part_name, model='Model-1', description='', type=ANALYSIS, 
            atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, 
            memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, 
            explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, 
            modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', 
            scratch='', resultsFormat=ODB)
        mdb.jobs[part_name].submit(consistencyChecking=OFF)

    if Trazione:
        mdb.Job(name=part_name, model='Model-1', description='', type=ANALYSIS, 
            atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, 
            memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, 
            explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, 
            modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', 
            scratch='', resultsFormat=ODB)
        mdb.jobs[part_name].submit(consistencyChecking=OFF)


    
if Modello_2D:  
    # Creazione parte 2D
    p = mdb.models['Model-1'].Part(name=part_name, dimensionality=THREE_D, type=DEFORMABLE_BODY)
    p = mdb.models['Model-1'].parts[part_name]
    p.BaseShell(sketch=s)
    s.unsetPrimaryObject()
    del mdb.models['Model-1'].sketches['__profile__']    

    # Crea assembly, utile se il modello dovesse avere more than one part
    a = mdb.models['Model-1'].rootAssembly
    a.DatumCsysByDefault(CARTESIAN)
    a.Instance(name='{}-{}'.format(part_name, 1), part=p, dependent=ON)
    # in assembly inserisco le condizioni al contorno e i carichi

    # Create step statico
    if Statica or Trazione:
        mdb.models['Model-1'].StaticStep(name='Step-1', previous='Initial')

    # Create step frequenza
    if Frequenze_Naturali:   
        mdb.models['Model-1'].FrequencyStep(name='Step-1', previous='Initial', 
            simLinearDynamics=OFF, normalization=MASS, numEigen=10)

    # # Create step buck
    if Buckling_a_compressione:
        mdb.models['Model-1'].BuckleStep(name='Step-1', previous='Initial', numEigen=5, 
            vectors=10, maxIterations=3000)

    ## creazione BOUNDARY CONDITION
    if Statica or Frequenze_Naturali:                # Statica
        # Definizione dei parametri
        instance_name = '{}-{}'.format(part_name, 1)

        if Statica:
            step_name = 'Step-1'
        elif Frequenze_Naturali:
            step_name = 'Initial' 

        if Open_Hole_completo:
            bc_definitions = [              # aggiornando la lista bc_definitions, si modificano i vincoli di incastro
                {'mask': '[#2 ]', 'set_name': 'Set-1', 'bc_name': 'clamped1'},  # lato sinistro, P1-P4
                {'mask': '[#4 ]', 'set_name': 'Set-2', 'bc_name': 'clamped2'},  # lato sotto, P1-P2
                {'mask': '[#8 ]', 'set_name': 'Set-3', 'bc_name': 'clamped3'},  # lato destro, P2-P3
                {'mask': '[#10 ]', 'set_name': 'Set-4', 'bc_name': 'clamped4'}   # lato sopra, P3-P4
            ]
        else:
           bc_definitions = [              
                {'mask': '[#1 ]', 'set_name': 'Set-1', 'bc_name': 'clamped1'},  # lato sinistro, P1-P4
                {'mask': '[#2 ]', 'set_name': 'Set-2', 'bc_name': 'clamped2'},  # lato sotto, P1-P2
                {'mask': '[#4 ]', 'set_name': 'Set-3', 'bc_name': 'clamped3'},  # lato destro, P2-P3
                {'mask': '[#8 ]', 'set_name': 'Set-4', 'bc_name': 'clamped4'}   # lato sopra, P3-P4
            ] 

        # Funzione per creare un vincolo di incastro
        def create_encastre_bc(instance_name, mask, set_name, bc_name, step_name):
            edges = a.instances[instance_name].edges
            edges_selected = edges.getSequenceFromMask(mask=(mask, ), )
            region = a.Set(edges=edges_selected, name=set_name)
            mdb.models['Model-1'].EncastreBC(
                name=bc_name, 
                createStepName=step_name, 
                region=region, 
                localCsys=None
            )

        # Creazione dei vincoli di incastro utilizzando i parametri definiti
        for bc_def in bc_definitions:
            create_encastre_bc(
                instance_name=instance_name,
                mask=bc_def['mask'],
                set_name=bc_def['set_name'],
                bc_name=bc_def['bc_name'],
                step_name=step_name
            )

    if Trazione:
        if Open_Hole:     
            instance_name = '{}-{}'.format(part_name, 1)
            
            a = mdb.models['Model-1'].rootAssembly
            e1 = a.instances[instance_name].edges
            edges1 = e1.getSequenceFromMask(mask=('[#2 ]', ), )
            region = regionToolset.Region(edges=edges1)
            mdb.models['Model-1'].DisplacementBC(name='BC-1', createStepName='Step-1', 
                region=region, u1=UNSET, u2=0.0, u3=0.0, ur1=UNSET, ur2=0.0, ur3=0.0, 
                amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
                localCsys=None)
            
            a = mdb.models['Model-1'].rootAssembly
            e1 = a.instances[instance_name].edges
            edges1 = e1.getSequenceFromMask(mask=('[#10 ]', ), )
            region = regionToolset.Region(edges=edges1)
            mdb.models['Model-1'].DisplacementBC(name='BC-2', createStepName='Step-1', 
                region=region, u1=0.0, u2=UNSET, u3=0.0, ur1=0.0, ur2=UNSET, ur3=0.0, 
                amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
                localCsys=None)
        elif Open_Hole_completo:
            instance_name = '{}-{}'.format(part_name, 1)
            
            a = mdb.models['Model-1'].rootAssembly
            e1 = a.instances[instance_name].edges
            edges1 = e1.getSequenceFromMask(mask=('[#2 ]', ), )
            region = regionToolset.Region(edges=edges1)
            mdb.models['Model-1'].DisplacementBC(name='BC-1', createStepName='Step-1', 
                region=region, u1=0.0, u2=UNSET, u3=0.0, ur1=0.0, ur2=UNSET, ur3=0.0, 
                amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
                localCsys=None)
            
            a = mdb.models['Model-1'].rootAssembly
            e1 = a.instances[instance_name].edges
            edges1 = e1.getSequenceFromMask(mask=('[#8 ]', ), )
            region = regionToolset.Region(edges=edges1)
            mdb.models['Model-1'].DisplacementBC(name='BC-2', createStepName='Step-1', 
                region=region, u1=0.0, u2=UNSET, u3=0.0, ur1=0.0, ur2=UNSET, ur3=0.0, 
                amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
                localCsys=None)

  # # creazione CARICO, 
    if Statica:
        # Definizione dei parametri
        instance_name = '{}-{}'.format(part_name, 1)    
        mask = '[#1 ]'
        surface_name = 'Surf-1'                         # superficie superiore
        pressure_name = 'Pressure'
        step_name = 'Step-1'
        distribution_type = UNIFORM

        # Creazione del carico di pressione
        faces = a.instances[instance_name].faces
        faces1 = faces.getSequenceFromMask(mask=(mask, ), )
        region = a.Surface(side1Faces=faces1, name=surface_name)
        mdb.models['Model-1'].Pressure(
        name=pressure_name, 
        createStepName=step_name, 
        region=region, 
        distributionType=distribution_type, 
        field='', 
        magnitude=pressure              # definita all'inizio del codice 
        )

    if Trazione:
        # Parametri
        instance_name = '{}-{}'.format(part_name, 1)
        step_name = 'Step-1'
        distribution_type = UNIFORM
        if Open_Hole:
            load_definitions = [
                {'mask': '[#4 ]', 'load_name': 'Load-1'},  # 4 e 10 per bordi in y
                # {'mask': '[#10 ]', 'load_name': 'Load-2'}  # 2 e 8 per bordi in x
            ]
        elif Open_Hole_completo:
             load_definitions = [
                {'mask': '[#4 ]', 'load_name': 'Load-1'},  # 4 e 10 per bordi in y
                {'mask': '[#10 ]', 'load_name': 'Load-2'}  # 2 e 8 per bordi in x
            ]
        else:
            load_definitions = [
                {'mask': '[#2 ]', 'load_name': 'Load-1'},  # 2 e 8 per bordi in y
                {'mask': '[#8 ]', 'load_name': 'Load-2'}  # 1 e 4 per bordi in x
            ]

        # Function to create a shell edge load
        def create_shell_edge_load(instance_name, mask, load_name, step_name, magnitude, distribution):
            edges = a.instances[instance_name].edges
            edges_selected = edges.getSequenceFromMask(mask=(mask, ), )
            region = regionToolset.Region(side1Edges=edges_selected)
            mdb.models['Model-1'].ShellEdgeLoad(
                name=load_name,
                createStepName=step_name,
                region=region,
                magnitude=magnitude,
                distributionType=distribution,
                field='',
                localCsys=None
            )

        a = mdb.models['Model-1'].rootAssembly
        for load_def in load_definitions:
            create_shell_edge_load(
                instance_name=instance_name,
                mask=load_def['mask'],
                load_name=load_def['load_name'],
                step_name=step_name,
                magnitude=Carico_Trazione,
                distribution=distribution_type
            )

    # creazione MESH                                                                  
    if Open_Hole:
        if only_S3:
            p = mdb.models['Model-1'].parts[part_name]
            f = p.faces
            pickedRegions = f.getSequenceFromMask(mask=('[#1 ]', ), )
            p.setMeshControls(regions=pickedRegions, elemShape=TRI)
            elemType1 = mesh.ElemType(elemCode=S4R, elemLibrary=STANDARD)
            elemType2 = mesh.ElemType(elemCode=S3, elemLibrary=STANDARD, 
                secondOrderAccuracy=OFF)
            p = mdb.models['Model-1'].parts[part_name]
            f = p.faces
            faces = f.getSequenceFromMask(mask=('[#1 ]', ), )
            pickedRegions =(faces, )
            p.setElementType(regions=pickedRegions, elemTypes=(elemType2, elemType2))
            p = mdb.models['Model-1'].parts[part_name]
            e = p.edges
            pickedEdges = e.getSequenceFromMask(mask=('[#1 ]', ), )
            p.seedEdgeByNumber(edges=pickedEdges, number=nc, constraint=FINER)
            p = mdb.models['Model-1'].parts[part_name]
            e = p.edges
            pickedEdges = e.getSequenceFromMask(mask=('[#2 ]', ), )
            p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER)
        else:
            p = mdb.models['Model-1'].parts[part_name]                                   # Accesso al modello e alla parte
            e = p.edges 
            pickedEdges = e.getSequenceFromMask(mask=('[#1 ]', ), )                     # Seleziona specifici spigoli utilizzando una maschera
            p.seedEdgeByNumber(edges=pickedEdges, number=nc, constraint=FINER)           # Specifica numero elementi sullo spigolo
            pickedEdges = e.getSequenceFromMask(mask=('[#2 ]', ), )
            p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER)          # Definisce il numero di elementi lungo i bordi selezionati
            p = mdb.models['Model-1'].parts[part_name]                                 
            f = p.faces
            pickedRegions = f.getSequenceFromMask(mask=('[#1 ]', ), )
            p.setMeshControls(regions=pickedRegions, technique=STRUCTURED, elemShape=QUAD_DOMINATED)
    elif Open_Hole_completo:
            p = mdb.models['Model-1'].parts[part_name]                                  
            e = p.edges 
            pickedEdges = e.getSequenceFromMask(mask=('[#1 ]', ), )                     
            p.seedEdgeByNumber(edges=pickedEdges, number=nc, constraint=FINER)
            pickedEdges = e.getSequenceFromMask(mask=('[#2 ]', ), )
            p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER) 
            p = mdb.models['Model-1'].parts[part_name]                                 
            f = p.faces
            pickedRegions = f.getSequenceFromMask(mask=('[#1 ]', ), )
            p.setMeshControls(regions=pickedRegions, technique=FREE, elemShape=QUAD_DOMINATED)      #elemShape=TRI se elementi triangolari  
    else:
        p = mdb.models['Model-1'].parts[part_name]                                   # Accesso al modello e alla parte
        e = p.edges 
        pickedEdges = e.getSequenceFromMask(mask=('[#1 ]', ), )                     
        p.seedEdgeByNumber(edges=pickedEdges, number=n, constraint=FINER)  
    p = mdb.models['Model-1'].parts[part_name]
    p.generateMesh()                                                                # Generazione della mesh
                                                            
    # BC buck
    if Buckling_a_compressione:
        a = mdb.models['Model-1'].rootAssembly
        e1 = a.instances['{}-{}'.format(part_name, 1)].edges
        edges1 = e1.getSequenceFromMask(mask=('[#2 ]', ), )  # mettere #2 se si vuole bloccare asse y, #1 se asse x
        region = regionToolset.Region(edges=edges1)
        mdb.models['Model-1'].EncastreBC(name='BC-1', createStepName='Step-1', 
                region=region, localCsys=None, buckleCase=PERTURBATION_AND_BUCKLING)
    
        # buck load
        a = mdb.models['Model-1'].rootAssembly
        s1 = a.instances['{}-{}'.format(part_name, 1)].edges
        side1Edges1 = s1.getSequenceFromMask(mask=('[#8 ]', ), ) # mettere #8 se si vuole caricare asse y, #4 se asse x
        region = regionToolset.Region(side1Edges=side1Edges1)
        mdb.models['Model-1'].ShellEdgeLoad(name='Load-Buck', createStepName='Step-1', 
                region=region, magnitude=1.0, distributionType=UNIFORM, field='', 
                localCsys=None)

    # Crea .inp
    mdb.Job(name=Inp_name, model='Model-1', description='', type=ANALYSIS, 
        atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, 
        memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, 
        explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, 
        modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', 
        scratch='', resultsFormat=ODB)
    mdb.jobs[Inp_name].writeInput(consistencyChecking=OFF)
    print('File inp creato')

    def get_data_from_inp(filename, dataname, cast=str):
        data = []
        with open(filename, 'r') as file:
            str_file = file.read()
            lines = str_file.split('\n')
            
            index = next((i for i, line in enumerate(lines) if '*'+dataname in line), None) + 1
            while '*' not in lines[index][0]:
                data.append([cast(line.strip()) for line in lines[index].split(',')])
                index += 1
        return data

    matrice_nodi=(get_data_from_inp("{}.inp".format(Inp_name),'Node',lambda x: float(x+'0' if '.' in x else x)))
    # print(matrice_nodi)   
    if Open_Hole_completo:
        if only_S3:
            matrice_Assembly= (get_data_from_inp("{}.inp".format(Inp_name),'Element',int))
        else:
            matrice_Assembly_S3= (get_data_from_inp("{}.inp".format(Inp_name),'Element, type=S3',int))
            matrice_Assembly_S4R= (get_data_from_inp("{}.inp".format(Inp_name),'Element, type=S4R',int))
            matrice_Assembly =   matrice_Assembly_S3 + matrice_Assembly_S4R  
    else:
        matrice_Assembly= (get_data_from_inp("{}.inp".format(Inp_name),'Element',int))
    # print(matrice_Assembly)
    print('Numero elementi:',len(matrice_Assembly))
    Numero_elementi= len(matrice_Assembly)
    nodi_ele=len(matrice_Assembly[0])-1
    print('Nodi per elemento:',nodi_ele)

    X = [riga[1] for riga in matrice_nodi]
    Y = [riga[2] for riga in matrice_nodi]
    Z = [riga[3] for riga in matrice_nodi]
    # print ("Vettore x:", X)
    # print(Elementi)

    xc = [0] * len(matrice_Assembly)
    yc = [0] * len(matrice_Assembly)
    zc = [0] * len(matrice_Assembly)

    for element_num in range(len(matrice_Assembly)):
        for vertix_num in matrice_Assembly[element_num][1:]:
                 
            xc[element_num] += X[vertix_num-1]
            yc[element_num] += Y[vertix_num-1]
            zc[element_num] += Z[vertix_num-1]

        xc[element_num] /= len(matrice_Assembly[0][1:])
        yc[element_num] /= len(matrice_Assembly[0][1:])
        zc[element_num] /= len(matrice_Assembly[0][1:])
        # print(element_num+1, xc[element_num]) #, il centroide calcolato correttamente
     
    # Imposto centro di riferimento funzione al centro della piastra
    x_coord = [0] * len(matrice_Assembly)
    y_coord = [0] * len(matrice_Assembly)

    for element_num in range(len(matrice_Assembly)):
        if Open_Hole:
            x_coord [element_num]= xc[element_num]
            y_coord [element_num]= yc[element_num]
        else:
            x_coord [element_num]= L1/2-xc[element_num]
            y_coord [element_num]= L2/2-yc[element_num]
        # print(x_coord [element_num])
    
    def VAT (x, y, phi, T0, T1, L):
        if not Open_Hole:
            return phi + T0 + (T1-T0)/L * np.abs((x * np.cos(np.deg2rad(phi))) + (y * np.sin(np.deg2rad(phi))))          
        else:
            return phi + T0 + (T1 - T0) / L * ((x * np.cos(np.deg2rad(phi))) + (y * np.sin(np.deg2rad(phi))))
  
    # Ora creo il field su ABAQUS   
    for i in range(kk):
        model_name = 'Model-1'
        field_type = SCALAR
        location = ELEMENTS
        data_width = 1
        default_values = (0.0,)  
        
        field_names = []
        elementi =[]
        values = []
        thetaC=[]
        for element_num in range(len(matrice_Assembly)):
            # field_name = 'DiscField-' + str(i + 1) 
            index1 = 2 * i
            index2 = 2 * i + 1
            
            theta0 = theta[index1]
            theta1 = theta[index2]
            
            if Open_Hole:
                thetaC.append(VAT(x_coord[element_num],yc[element_num],phi,theta0,theta1,L1)) 
            else:
                thetaC.append(VAT(x_coord[element_num],yc[element_num],phi,theta0,theta1,L1/2))    
            elementi.append(int(element_num+1))
            values.append(thetaC[element_num])
            # print("Futuro Field", element_num+1, thetaC[element_num])
        field_name = 'DiscField-' + str(i + 1)
        
        mdb.models[model_name].DiscreteField(
                name=field_name,
                description='',
                location=location,
                fieldType=field_type,
                dataWidth=data_width,
                defaultValues=default_values,
                data=(('', 1, elementi, values),))
        field_names.append(field_name)
              
    #Creazione PLY
    # Definizione del numero di ply e dei campi di rotazione
    rotation_fields = ['DiscField-' + str(i+1) for i in range(kk)]  # Generazione dei nomi dei campi di rotazione
    # print(rotation_fields) 
    # Selezione del modello e della parte
    p = mdb.models['Model-1'].parts[part_name]                  # Selezione del modello e la parte
    f = p.faces                                                 # Ottieni le facce della parte
    faces = f.getSequenceFromMask(mask=('[#1 ]', ), )           # Selezione le facce utilizzando una maschera
    region = regionToolset.Region(faces=faces)                  # Definizione della regione

    # Creazione del CompositeLayup
    compositeLayup = p.CompositeLayup(
        name='CompositeLayup-1',                # nome
        description='', 
        elementType=SHELL,                      # tipo di elemento
        offsetType=MIDDLE_SURFACE, 
        symmetric=False,                        # non modificare mai                                     
        thicknessAssignment=FROM_SECTION
    )
    compositeLayup.Section(preIntegrate=OFF, integrationRule=SIMPSON,  thicknessType=UNIFORM,  poissonDefinition=DEFAULT, temperature=GRADIENT, useDensity=OFF)
    compositeLayup.ReferenceOrientation(orientationType=GLOBAL, localCsys=None, fieldName='', additionalRotationType=ROTATION_NONE, angle=0.0, axis=AXIS_3)
    compositeLayup.suppress()
    # prima indico una rotazione nulla per entrambi i ply, poi la si modifica inserendo i rispettivi field
    # Creazione dei ply in un ciclo for
    half_kk = kk // 2

    for i in range(kk):
        ply_name = 'Ply-' + str(i + 1)  
        rotation_field = rotation_fields[i]  
        print(rotation_field)

        compositeLayup.CompositePly(
            suppressed=False, 
            plyName=ply_name, 
            region=region, 
            material=Material_name,                  # richiamo materiale
            thicknessType=SPECIFY_THICKNESS, 
            thickness=L3,                               # Spessore specificato
            orientationType=ANGLE_0, 
            additionalRotationType=ROTATION_FIELD, 
            additionalRotationField=rotation_field,
            axis=AXIS_3, 
            angle=0.0,  
            numIntPoints=3
        )
    compositeLayup.resume()

    # ANALISI
    # Definizione dei parametri variabili
    if Statica:
        job_name = 'Job-Statico_Fibrecurve'                 # nome analisi
        model_name = 'Model-1'                              # modello su cui si svolge l'analisi
        description = 'Analisi statica- fibre curve'        # descrizione opzionale
        results_format = ODB                                # file output

        # Creazione del JOB
        mdb.Job(name=job_name, model=model_name, 
                description=description, type=ANALYSIS, atTime=None, 
                waitMinutes=0, waitHours=0, queue=None, memory=90, 
                memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, 
                explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, 
                modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', 
                scratch='', resultsFormat=results_format)

        # Sottomissione del job
        mdb.jobs[job_name].submit(consistencyChecking=OFF)

        # # Creazione del punto Q per confrontare risultati
        # a = mdb.models['Model-1'].rootAssembly
        # coords=(L1/4,L2/4,0.0)
        # a.DatumPointByCoordinate(coords=coords)

        # nid = ((n - 4) / 2 + 3) ** 2 , per sapere in quale nodo cercare i risultati

    if Frequenze_Naturali:
        job_name = 'Job-Freq_Nat'                                 # nome analisi
        model_name = 'Model-1'                                    # modello su cui si svolge l'analisi
        description = 'Analisi frequenze naturali - fibre curve'  # descrizione opzionale
        results_format = ODB                                      # file output
    
    # Creazione del JOB
        mdb.Job(name=job_name, model=model_name, description=description, 
                type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0,
                queue=None, memory=90, memoryUnits=PERCENTAGE, 
                getMemoryFromAnalysis=True, explicitPrecision=SINGLE,
                nodalOutputPrecision=SINGLE, echoPrint=OFF, 
                modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', 
                scratch='', resultsFormat=results_format)

        # Sottomissione del job
        mdb.jobs[job_name].submit(consistencyChecking=OFF)

    if Buckling_a_compressione:
        job_name = 'Job-Buckling'               # nome analisi
        model_name = 'Model-1'                  # modello su cui si svolge l'analisi
        description = 'Analisi Buckling'        # descrizione opzionale
        results_format = ODB                    # file output   

        mdb.Job(name=job_name, model=model_name, description=description,
                    type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, 
                    memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, 
                    explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, 
                    modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', 
                    scratch='', resultsFormat=results_format)
        mdb.jobs[job_name].submit(consistencyChecking=OFF)

    if Trazione:
        job_name = 'Job-Trazione'               # nome analisi
        model_name = 'Model-1'                  # modello su cui si svolge l'analisi
        description = 'Analisi Open_Hole '        # descrizione opzionale
        results_format = ODB                    # file output   

        mdb.Job(name=job_name, model=model_name, description=description,
                    type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, 
                    memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, 
                    explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, 
                    modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', 
                    scratch='', resultsFormat=results_format)
        mdb.jobs[job_name].submit(consistencyChecking=OFF)

    # # Visualizzazione finale
if Modello_2D or Modello_3D:
    p = mdb.models['Model-1'].parts[part_name]
    session.viewports['Viewport: 1'].setValues(displayedObject=p)
