Differenze

Queste sono le differenze tra la revisione selezionata e la versione attuale della pagina.

--- wikipaom2015:lez14 [2015/04/24 15:56] – [MAIN] 165198
+++ wikipaom2015:lez14 [2015/04/27 20:20] (versione attuale) – 166078
@@ Linea 1: / Linea 1: @@
+======Programma agli elementi finiti triangolari======
+=====Dimensionamento=====
+La prima parte del programma si occupa della dichiarazione delle variabili necessarie per definire la struttura;
+i parametri NODES e NELEMS identificano una mesh di 9 nodi e 8 elementi.
+      PARAMETER (NODES=9)
+      PARAMETER (NELEMS=8)
+      PARAMETER (IFMAX=3,ICNMAX=3)
+      DIMENSION XY(2,NODES),NVERT(3,NELEMS),IFXY(IFMAX),FXY(2,IFMAX),
+      ICNXY(2,ICNMAX),CNXY(2,ICNMAX),STRUTK(2*NODES,2*NODES),
+      FORCE(2*NODES),UV(2*NODES)
+**IFMAX**: numero di nodi caricati\\
+**ICNMAX**: numero di nodi vincolati
+La stringa "Dimension" permette al programma di allocare le matrici necessarie per i calcoli, dedicandogli un appropriato spazio in memoria.\\
+**XY**: matrice che contiene le coordinate locali\\
+**NVERT**: matrice che contiene il numero dei nodi ai tre vertici di ogni elemento\\
+**IFXY**: vettore che contiene i nomi (i numeri) dei nodi caricati\\
+**FXY**: matrice che contiene le componenti lungo x ed y delle forze applicate\\
+**ICNXY**: matrice che contiene i nomi dei nodi vincolati e gli associa il tipo di vincolamento.
+           - vincolato solo lungo X
+           - vincolato solo lungo Y
+           - vincolato lungo X ed Y
+**CNXY**: matrice che contiene gli spostamenti imposti lungo X ed Y\\
+**STRUTK**: matrice di rigidezza dell'intera struttura\\
+**TK**: vettore dei termini noti\\
+**FORCE**: vettore delle forze imposte\\
+**UV**: vettore spostamenti nodali
+      DIMENSION D(3,3),B(3,6),BT(6,3),DB(3,6),ELK(6,6),IPOINT(6),
+      DELTAEL(6),SIGMAEL(3)
+Dopo il dimensionamento degli elementi variabili si termina l'allocazione di matrici e vettori defininendo le grandezze dipendenti dalla tipologia di risoluzione scelta (tria 3 nodi).\\
+**D**: matrice di legame costitutivo (legge di Hooke), lavorando in 3D diventa una matrice 6x6\\
+**B**: matrice che lega spostamenti e deformazione\\
+**BT**: trasposta di B\\
+**DB**: prodotto tra matrice D e B (matrice temporanea)\\
+**ELK**: matrice di rigidezza degli elementi (aggiornata per ogni elemento, evitando così di avere un array in tre dimensioni 6x6xNELEMS)\\
+**IPOINT**: mappatura dei gdl globali e locali per ogni elemento\\
+**DELTAEL**: vettore di 6 spostamenti locali dei 3 vertici dell'elemento\\
+**SIGMAEL**: vettore tre elementi che contiene la σ<sub>x</sub>, la σ<sub>y</sub> e la τ<sub>xy</sub>
+=====MAIN=====
+      OPEN(1,FILE='mesh_tria_comp.dat',STATUS='OLD')
+      OPEN(2,FILE='mesh_tria_comp.3.txt',STATUS='UNKNOWN')
+      CALL READIN(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,FXY,
+      ICNMAX,ICNXY,CNXY)
+      CALL ECHO(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,FXY,ICNMAX,
+      ICNXY,CNXY)
+      CALL CLEAR(2*NODES,2*NODES,STRUTK)
+      CALL CLEAR(2*NODES,1,FORCE)
+  C     CALL CLEAR(3,6,B)
+  C     CALL CLEAR(6,6,ELK)
+      CALL DMAT(YM,PR,ITDP,D)
+====Subroutine Readin====
+Lettura del Modulo di Young e del Coefficiente di Poisson e scelta dello stato di Tensione o Deformazione Piana.\\
+Lettura delle coordinate nodali e dei vertici degli elementi triangolari.\\ Lettura delle forze nodali e delle condizioni di vincolamento.
+     SUBROUTINE READIN(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,
+      FXY,ICNMAX,ICNXY,CNXY)
+      DIMENSION XY(2,NODES),NVERT(3,NELEMS),IFXY(IFMAX),FXY(2,IFMAX),
+      ICNXY(2,ICNMAX),CNXY(2,ICNMAX)
+  C     Lettura di YM, PR e ITDP:
+      READ(1,*)YM,PR,ITDP
+  C     Lettura delle coordinate nodali (x,y):
+      DO 10,I10=1,NODES
+      READ(1,*)NODE,XY(1,I10),XY(2,I10)
+      IF(I10.NE.NODE)THEN
+      write(*,*) 'errore nella numerazione dei nodi sul file di input'
+      stop
+      END IF
+CONTINUE
+  C     Lettura dei vertici degli elementi triangolari:
+      DO 20,I20=1,NELEMS
+      READ(1,*)NELEM,NVERT(1,I20),NVERT(2,I20),NVERT(3,I20)
+      IF(I20.NE.NELEM)THEN
+      write(*,*) 'errore nella numerazione degli elementi'
+      stop
+      END IF
+CONTINUE
+  C     Lettura delle forze nodali:
+      DO 30,I30=1,IFMAX
+      READ(1,*)IFXY(I30),FXY(1,I30),FXY(2,I30)
+CONTINUE
+  C     Lettura delle condizioni di vincolamento:
+      DO 40,I40=1,ICNMAX
+      READ(1,*)ICNXY(1,I40),ICNXY(2,I40),CNXY(1,I40),CNXY(2,I40)
+CONTINUE
+      RETURN
+      END
+Per ogni ciclo DO presente nella subroutine il programma esegue un controllo di sequenzialità dei parametri (coordianate nodali, spostamenti e forze applicate) dei file di input.
+====Subroutine Echo====
+Subroutine che controlla la correttezza di lettura dei dati di input.
+     SUBROUTINE ECHO(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,FXY,
+      ICNMAX,ICNXY,CNXY)
+      DIMENSION XY(2,NODES),NVERT(3,NELEMS),IFXY(IFMAX),FXY(2,IFMAX),
+      ICNXY(2,ICNMAX),CNXY(2,ICNMAX)
+  C     Controllo di YM e PR:
+      WRITE(2,*)'MODULO DI YOUNG: ',YM
+      WRITE(2,*)'COEFFICIENTE DI POISSON: ',PR
+  C     Controllo di ITDP:
+      IF(ITDP.EQ.0)THEN
+      WRITE(2,*)'STATO DI TENSIONE PIANA'
+      END IF
+      IF(ITDP.EQ.1)THEN
+      WRITE(2,*)'STATO DI DEFORMAZIONE PIANA'
+      END IF
+  C     Controllo delle coordinate nodali:
+      WRITE(2,*)'NODO: X: Y: '
+      DO 10,I10=1,NODES
+      WRITE(2,*)I10,XY(1,I10),XY(2,I10)
+CONTINUE
+  C     Controllo dei vertici degli elementi trangolari:
+      WRITE(2,*)'ELEMENTO: Vi: Vj: Vk: '
+      DO 20,I20=1,NELEMS
+      WRITE(2,*)I20,NVERT(1,I20),NVERT(2,I20),NVERT(3,I20)
+CONTINUE
+  C     Controllo delle forze:
+      WRITE(2,*)'NODO: FX: FY: '
+      DO 30,I30=1,IFMAX
+      WRITE(2,*)IFXY(I30),FXY(1,I30),FXY(2,I30)
+CONTINUE
+  C     Controllo delle condizioni di vincolamento:
+      WRITE(2,*)'NODO: TIPO VINCOLAMENTO: DELTAX: DELTAY: '
+      DO 40,I40=1,ICNMAX
+      WRITE(2,*)ICNXY(1,I40),ICNXY(2,I40),CNXY(1,I40),CNXY(2,I40)
+CONTINUE
+      RETURN
+      END
+====Subroutine Clear====
+Inizializzazione a zero degli elementi di una matrice.
+     SUBROUTINE CLEAR(K1,K2,A)
+      DIMENSION A(K1,K2)
+      DO 10,I10=1,K1
+      DO 20,I20=1,K2
+      A(I10,I20)=0.
+CONTINUE
+CONTINUE
+      RETURN
+      END
+====Subroutine Dmat====
+Costruzione della matrice D di applicazione del legame Hookeano tra tensioni e deformazioni: **σ** = D * **ε**.
+     SUBROUTINE DMAT(YM,PR,ITDP,D)
+      DIMENSION D(3,3)
+      IF(ITDP.EQ.0)GO TO 100
+      IF(ITDP.EQ.1)GO TO 200
+  C     Per Stato di Tensione Piana:
+CONTINUE
+      CTP=YM/(1-PR**2)
+      G=YM/(2*(1+PR))
+  C     Struttura della matrice D per Stato di Tensione Piana:
+  C                     | 1*CTP    PR*CTP   0 |
+  C                D =  | PR*CTP   1*CTP    0 |
+  C                     | 0          0      G |
+      D(1,1)=1.
+      D(1,2)=PR
+      D(2,1)=PR
+      D(2,2)=1.
+      DO 10,I10=1,2
+      DO 20,I20=1,2
+      D(I10,I20)=D(I10,I20)*CTP
+CONTINUE
+CONTINUE
+      D(1,3)=0.
+      D(2,3)=0.
+      D(3,1)=0.
+      D(3,2)=0.
+      D(3,3)=G
+      GOTO 300
+  C     Per Stato di Deformazione Piana:
+CONTINUE
+      CDP=YM*(1-PR)/((1+PR)*(1-2*PR))
+      CPR=PR/(1-PR)
+      G=YM/((1+PR)*2)
+  C     Struttura della matrice D per Stato di Deformazione Piana:
+  C                     | 1*CDP     CDP*CPR   0 |
+  C                D =  | CDP*CPR    1*CDP    0 |
+  C                     | 0            0      G |
+      D(1,1)=1.
+      D(1,2)=CPR
+      D(2,1)=CPR
+      D(2,2)=1.
+      DO 30,I30=1,2
+      DO 40,I40=1,2
+      D(I30,I40)=D(I30,I40)*CDP
+CONTINUE
+CONTINUE
+      D(1,3)=0.
+      D(2,3)=0.
+      D(3,1)=0.
+      D(3,2)=0.
+      D(3,3)=G
+      GOTO 300
+CONTINUE
+      RETURN
+      END
+Per decidere il tipo di legame costitutivo da applicare per la generazione della matrice D, la subroutine legge il valore di ITDP attraverso due costrutti IF: questo genere di approccio potrebbe portare ad errore nel caso in cui ITDP fosse uguale a qualunque altro numero che non sia zero o uno, in quanto la sub continuerebbe con l'esecuzione delle istruzioni nell'ordine in cui sono fornite; assocerebbe quindi ad uno stato tensionale sconosciuto (ITDP diverso da 0 e 1) lo stato di tensione piana (poiché è il primo che legge).
+=====MAIN=====
+   DO 10,I10=1,NELEMS
+      XI=XY(1,NVERT(1,I10))
+      YI=XY(2,NVERT(1,I10))                     !CALL AREA, KELMAT
+      XJ=XY(1,NVERT(2,I10))                     !POINTVT eASSEMBL
+      YJ=XY(2,NVERT(2,I10))                     !sono dentro lo stesso
+      XK=XY(1,NVERT(3,I10))                     !ciclo DO 10, cos
+      YK=XY(2,NVERT(3,I10))                     !vengono calcolati per
+      CALL AREA(XI,YI,XJ,YJ,XK,YK,I10,TWODEL)   !ogni elemento.
+      CALL KELMAT(XI,YI,XJ,YJ,XK,YK,TWODEL,I10,D,ELK)
+  C     *****Assemblaggio***********
+      CALL POINTVT(NELEMS,I10,NVERT,IPOINT)
+      CALL ASSEMBL(ELK,NODES,IPOINT,STRUTK)
+CONTINUE
+DO principale per la costruzione della matrice di rigidezza dell'elemento generico I10.
+Scorro tutti gli elementi (I10=1,NELEMS); per ogni elemento prendo il nome dei nodi che gli appartengono dalla matrice NVERT. Scelto il nodo, prendo le coordinate x e y all'interno della matrice XY (scegliendo tra la riga 1 o la riga 2).
+====Subroutine Area====
+Calcolo della superficie DEL di ciascun elemento finito triangolare a partire dalle coordinate nodali dei vertici. Controllo di dimensioni e segno dell'area calcolata. Le coordinate nodali vengono prese, come di consueto, in verso antiorario. Con TWODEL si identifica il doppio della superficie.
+        SUBROUTINE AREA(XI,YI,XJ,YJ,XK,YK,NELEM,TWODEL)
+        TWODEL=XJ*YK+XK*YI+XI*YJ-(XK*YJ+XI*YK+XJ*YI)
+        AREAMIN=1.E-6
+  C     Controllo sulle dimensioni dell'area calcolata:
+        IF(ABS((TWODEL/2)).LE.AREAMIN)THEN
+        WRITE(*,*)'ATTENZIONE! AREA DELL*ELEMENTO TROPPO PICCOLA'
+        WRITE(*,*)NELEM
+        WRITE(*,*)'AREA=',TWODEL/2
+        END IF
+  C     Controllo sul segno dell'area calcolata:
+        IF((TWODEL/2).LT.-AREAMIN)THEN
+        WRITE(*,*)'ATTENZIONE! AREA DELL*ELEMENTO NEGATIVA'
+        WRITE(*,*)NELEM
+        WRITE(*,*)'AREA=',TWODEL/2
+        END IF
+        RETURN
+        END
+TWODEL è il determinante della matrice delle coordinate locali. La subroutine effettua un controllo sulle dimensioni minime dell'area per permettere una corretta esecuzione dei calcoli successivi.\\
+**NB**: La mantissa di un numero di macchina dedica 7 bit per la codifica dell'esponente, perciò sommando due numeri con eccessiva differenza tra ordini di grandezza (es. 1 + 1e-8 = 1) potrebbe venire a perdersi il minore.
+====Subroutine Kelmat====
+Costruzione della matrice di rigidezza ELK di ciascun elemento finito triangolare attraverso la formula: F = ELK * DELTA => ELK = A * (BT * D * B). Le coordinate nodali vengono prese, come di consueto, in verso antiorario.
+      SUBROUTINE KELMAT(XI,YI,XJ,YJ,XK,YK,TWODEL,NELEM,D,ELK)
+      DIMENSION B(3,6),BT(6,3),D(3,3),DB(3,6),ELK(6,6)
+      CALL CLEAR(3,6,B)
+      CALL BMAT(XI,YI,XJ,YJ,XK,YK,TWODEL,B)
+      CALL PRODMAT(3,3,6,D,B,DB)
+      CALL TRSPMAT(3,6,B,BT)
+      CALL PRODMAT(6,3,6,BT,DB,ELK)
+      DEL=TWODEL/2
+      DO 10,I10=1,6
+      DO 20,I20=1,6
+      ELK(I10,I20)=DEL*ELK(I10,I20)
+CONTINUE
+CONTINUE
+      RETURN
+      END
+La matrice D, passata dal main, viene allocata in memoria al lancio della subroutine, mentre la matrice B (locale) richiede un allocamento all'inizio della subroutine, con relativa chiamata della subroutine CLEAR per inizializzazione della matrice stessa.
+====Subroutine Pointvt====
+Costruzione del vettore puntatore IPOINT necessario ad effettuare l'assemblaggio della matrice ELK di rigidezza di ciascun elemento finito triangolare (NELEM) nella matrice STRUTK di rigidezza globale della struttura.
+      SUBROUTINE POINTVT(NELEMS,NELEM,NVERT,IPOINT)
+      DIMENSION NVERT(3,NELEMS),IPOINT(6)
+      IPOINT(1)=NVERT(1,NELEM)*2-1
+      IPOINT(2)=NVERT(1,NELEM)*2
+      IPOINT(3)=NVERT(2,NELEM)*2-1
+      IPOINT(4)=NVERT(2,NELEM)*2
+      IPOINT(5)=NVERT(3,NELEM)*2-1
+      IPOINT(6)=NVERT(3,NELEM)*2
+      RETURN
+      END
+====Subroutine Assembl====
+Assemblaggio della matrice ELK di rigidezza di ciascun elemento finito triangolare nella matrice STRUTK di rigidezza globale della struttura.\\
+La matrice STRUTK viene inzializzata a zero nel MAIN PROGRAM, pertanto si ritiene superflua la sua inizializzazione in questa subroutine.
+      SUBROUTINE ASSEMBL(ELK,NODES,IPOINT,STRUTK)
+      DIMENSION ELK(6,6),IPOINT(6),STRUTK(2*NODES,2*NODES)
+      DO 10,I10=1,6
+      DO 20,I20=1,6
+      IND10=IPOINT(I10)
+      IND20=IPOINT(I20)
+      STRUTK(IND10,IND20)=STRUTK(IND10,IND20)+ELK(I10,I20)
+CONTINUE
+CONTINUE
+      RETURN
+      END
+Nota: il vettore NODES viene utilizzato unicamente per il dimensionamento della matrice STRUTK.
+La riga fondamentale della sub è quella che associa ad un valore della matrice di rigidezza globale della struttura il suo vecchio valore + il valore da aggiungere dovuto al nodo considerato, dato che al valore del termine STRUTK(I10,I20) della matrice di rigidezza globale contribuiscono in genere più elementi finiti.
+=====MAIN=====
+  C  Per il Caricamento si richiama questa sub
+     CALL FORCES(IFMAX,IFXY,FXY,NODES,FORCE)
+====Subroutine Forces====
+Riempimento del vettore globale delle forze nodali FORCE.
+      SUBROUTINE FORCES(IFMAX,IFXY,FXY,NODES,FORCE)
+      DIMENSION IFXY(IFMAX),FXY(2,IFMAX),FORCE(2*NODES)
+      DO 10,I10=1,IFMAX
+      INDX=IFXY(I10)*2-1
+      INDY=IFXY(I10)*2
+      FORCE(INDX)=FXY(1,I10)
+      FORCE(INDY)=FXY(2,I10)
+CONTINUE
+      RETURN
+      END
+Dove la sub importa dal main il numero di nodi caricati (IFMAX), il vettore contenente il nome dei nodi caricati (IFXY) e la matrice contenente i valori delle forze lungo x e y per ogni singolo nodo. Il ciclo DO legge i nodi caricati e salva i valori delle forze nel vettore dei termini noti.
+=====MAIN=====
+  C Vincolamento
+    CALL CNSTNG(ICNMAX,ICNXY,CNXY,NODES,FORCE,STRUTK)
+====Subroutine Cnstng====
+Imposizione delle condizioni di vincolamento specifiche alla matrice STRUTK di rigidezza globale della struttura ed al vettore globale delle forze nodali FORCE.
+        SUBROUTINE CNSTNG(ICNMAX,ICNXY,CNXY,NODES,FORCE,STRUTK)
+        DIMENSION ICNXY(2,ICNMAX),CNXY(2,ICNMAX),FORCE(2*NODES),
+       $STRUTK(2*NODES,2*NODES)
+        DO 10,I10=1,ICNMAX
+        IF(ICNXY(2,I10).EQ.1)GO TO 100
+        IF(ICNXY(2,I10).EQ.2)GO TO 200
+        IF(ICNXY(2,I10).EQ.3)GO TO 300
+        WRITE(*,*)'ERRORE NEL VINCOLAMENTO DEL NODO',I10
+        WRITE(*,*)'L*INDICE DI VINCOLAMENTO DEVE ESSERE 1, 2 o 3'
+        STOP
+  C     Per vincolamento di x:
+CONTINUE
+        INDX=ICNXY(1,I10)*2-1
+        DO 110,I110=1,2*NODES
+        STRUTK(INDX,I110)=0.
+        FORCE(I110)=FORCE(I110)-STRUTK(I110,INDX)*CNXY(1,I10)
+        STRUTK(I110,INDX)=0.
+CONTINUE
+        STRUTK(INDX,INDX)=1.
+        FORCE(INDX)=CNXY(1,I10)
+        GO TO 10
+  C   Per vincolamento di y:
+CONTINUE
+        INDY=ICNXY(1,I10)*2
+        DO 210,I210=1,2*NODES
+        STRUTK(INDY,I210)=0.
+        FORCE(I210)=FORCE(I210)-STRUTK(I210,INDY)*CNXY(2,I10)
+        STRUTK(I210,INDY)=0.
+CONTINUE
+        STRUTK(INDY,INDY)=1.
+        FORCE(INDY)=CNXY(2,I10)
+        GO TO 10
+  C     Per vincolamento di x e y:
+CONTINUE
+        INDX=ICNXY(1,I10)*2-1
+        INDY=ICNXY(1,I10)*2
+        DO 310,I310=1,2*NODES
+        STRUTK(INDX,I310)=0.
+        STRUTK(INDY,I310)=0.
+        FORCE(I310)=FORCE(I310)-STRUTK(I310,INDX)*CNXY(1,I10)-
+       $STRUTK(I310,INDY)*CNXY(2,I10)
+        STRUTK(I310,INDX)=0.
+        STRUTK(I310,INDY)=0.
+CONTINUE
+        STRUTK(INDX,INDX)=1.
+        STRUTK(INDY,INDY)=1.
+        FORCE(INDX)=CNXY(1,I10)
+        FORCE(INDY)=CNXY(2,I10)
+        GO TO 10
+CONTINUE
+        RETURN
+        END
+La subroutine importa il numero di nodi vincolati con il tipo di vincolo ad essi assegnato dal main. Le operazioni che permettono la modifica della matrice di rigidezza mantenendo la sua simmetria sono:\\
+_annullare tutti i coefficienti della matrice di rigidezza corrispondente alla riga ed alla colonna del grado di libertà vincolato, tranne il termine diagonale che verrà posto uguale ad 1.\\
+_si impone il temine noto corrispondente al grado di libertà vincolato pari al valore dello spostamento imposto.\\
+_si sottrae da tutti gli altri termini noti il contributo dovuto al prodotto tra i termini della colonna annullati ed il valore imposto dal nostro grado di libertà.\\
+La subroutine considera in maniera distinta i 3 tipi di problema: vincolamento lungo x, vincolamento lungo y, vincolamento lungo x ed y.\\ Per ogni tipo di vincolamento la sub rimanda alle operazioni di modifica della matrice STRUTK attraverso 3 diversi //GOTO//: questo modo di operare rende ridondante il codice, rischiando di moltiplicare eventuali errori tra le diverse operazioni, sarebbe stato più consono un costrutto del tipo IF-ELSE IF- ELSE- ENDIF.
+=====MAIN=====
+  C Soluzione in termini di spostamenti nodali
+      CALL GAUSS(2*NODES,STRUTK,FORCE,UV)
+      WRITE(2,*)'inodo Ux Uy'
+      DO 20,I20=1,NODES
+      WRITE(2,*)I20,UV(I20*2-1),UV(I20*2)
+CONTINUE
+====Subroutine Gauss====
+Risoluzione di un sistema algebrico lineare di N equazioni in N incognite attraverso il metodo di Gauss (triangolarizzazione della matrice dei coefficienti).\\
+La matrice dei coefficienti viene denominata C; il vettore dei termini noti B ed il vettore delle incognite V.
+        SUBROUTINE GAUSS(N,C,B,V)
+        DIMENSION C(N,N),B(N),V(N)
+  C     Processo di triangolarizzazione della matrice dei coefficienti:
+        DO 10,IMAIN=1,N-1
+        IBETTER=IMAIN
+        DO 101,ITRYME=IMAIN+1,N         !DO 101 scorre tutte le righe
+        VALBETTER=ABS(C(IBETTER,IMAIN)) !dalla IMAIN fino all'ultima.
+        VALTRYME=ABS(C(ITRYME,IMAIN))
+        IF(VALBETTER.LT.VALTRYME)THEN   !Se VALBETTER<VALTRYME, allora
+        IBETTER=ITRYME                  !pongo che ITRYME divenga la
+        END IF                          !riga "migliore" (con valore
+CONTINUE                        !assoluto maggiore).
+  C
+  C                                     !Scambio delle righe: si porta la
+        IF(IBETTER.NE.IMAIN)THEN        !riga con valore assoluto maggiore
+  C                                     !alla riga IMAIN.
+        DO 102,ISWAP=IMAIN,N            !DO 102 permette di scorrere
+        TEMP=C(IMAIN,ISWAP)             !tutti gli elementi della riga
+        C(IMAIN,ISWAP)=C(IBETTER,ISWAP) !(per spostare la riga si deve
+        C(IBETTER,ISWAP)=TEMP           !infatti spostare un elemento per
+CONTINUE                        !volta).
+        TEMP=B(IMAIN)
+        B(IMAIN)=B(IBETTER)             !Allo stesso modo si scambiano i
+        B(IBETTER)=TEMP                 !valori dei termini noti.
+        END IF
+  C
+        DO 103,ITRIANGLEME=IMAIN+1,N
+        RATIO=-C(ITRIANGLEME,IMAIN)/C(IMAIN,IMAIN)
+        C(ITRIANGLEME,IMAIN)=0.
+        DO 104,ISCROLL=IMAIN+1,N
+        C(ITRIANGLEME,ISCROLL)=C(ITRIANGLEME,ISCROLL)+
+       $RATIO*C(IMAIN,ISCROLL)
+CONTINUE
+        B(ITRIANGLEME)=B(ITRIANGLEME)+RATIO*B(IMAIN)
+CONTINUE
+CONTINUE
+  C     A questo punto la matrice dei coefficienti dovrebbe essere stata
+  C     triangolarizzata.
+  C
+  C     Risoluzione del sistema di equazioni:
+  C     Risoluzione dell'ultima equazione del sistema:
+        V(N)=B(N)/C(N,N)
+  C     Risoluzione delle altre equazioni:
+        DO 20,IEQNSCROLL=N-1,1,-1
+        SOMMA=0.
+        DO 202,ISOMMA=IEQNSCROLL+1,N
+        SOMMA=SOMMA+C(IEQNSCROLL,ISOMMA)*V(ISOMMA)
+CONTINUE
+        V(IEQNSCROLL)=(B(IEQNSCROLL)-SOMMA)/(C(IEQNSCROLL,IEQNSCROLL))
+CONTINUE
+        RETURN
+        END
+=====MAIN=====
+  C  Postprocessor per calcolare le tensioni
+      WRITE(2,*)'ielem SigmaX SigmaY Txy EqVonMises'
+      DO 30,I30=1,NELEMS
+      XI=XY(1,NVERT(1,I30))
+      YI=XY(2,NVERT(1,I30))
+      XJ=XY(1,NVERT(2,I30))
+      YJ=XY(2,NVERT(2,I30))
+      XK=XY(1,NVERT(3,I30))
+      YK=XY(2,NVERT(3,I30))
+      CALL BMAT(XI,YI,XJ,YJ,XK,YK,TWODEL,B)
+      CALL PRODMAT(3,3,6,D,B,DB)            !Prodotto tra matrice D e B
+      DELTAEL(1)=UV(2*NVERT(1,I30)-1)
+      DELTAEL(2)=UV(2*NVERT(1,I30))       !Riempie il vettore DELTAEL,
+      DELTAEL(3)=UV(2*NVERT(2,I30)-1)     !che contiene gli spostamenti
+      DELTAEL(4)=UV(2*NVERT(2,I30))       !nodali (x,y) dei 3 nodi di
+      DELTAEL(5)=UV(2*NVERT(3,I30)-1)     !ogni elemento.
+      DELTAEL(6)=UV(2*NVERT(3,I30))
+      CALL PRODMAT(3,6,1,DB,DELTAEL,SIGMAEL)
+====Subroutine Bmat====
+Costruzione della matrice B di collegamento tra gli spostamenti nodali e le deformazioni di ciascun elemento finito triangolare attraverso la formula: EPSILON = B * DELTA.\\
+Le coordinate nodali vengono prese, come di consueto, in verso antiorario.\\
+Infine, poiché i termini nulli, sempre nelle stesse posizioni, vengono gi inizializzati a zero inizialmente, se ne riportano le rispettive assegnazioni in commento.
+      SUBROUTINE BMAT(XI,YI,XJ,YJ,XK,YK,TWODEL,B)
+      DIMENSION B(3,6)
+      AI=XJ*YK-XK*YJ   !
+      AJ=XK*YI-XI*YK   ! Questi valori servono per il calcolo dell'area.
+      AK=XI*YJ-XJ*YI   !
+      BI=YJ-YK
+      BJ=YK-YI
+      BK=YI-YJ
+      CI=XK-XJ
+      CJ=XI-XK
+      CK=XJ-XI
+      B(1,1)=BI
+      B(1,2)=0.
+      B(1,3)=BJ
+      B(1,4)=0.
+      B(1,5)=BK
+      B(1,6)=0.
+      B(2,1)=0.
+      B(2,2)=CI
+      B(2,3)=0.
+      B(2,4)=CJ
+      B(2,5)=0.
+      B(2,6)=CK
+      B(3,1)=CI
+      B(3,2)=BI
+      B(3,3)=CJ
+      B(3,4)=BJ
+      B(3,5)=CK
+      B(3,6)=BK
+      TWODEL=AI+AJ+AK
+      DO 10,I10=1,3
+      DO 20,I20=1,6
+      B(I10,I20)=B(I10,I20)/TWODEL
+CONTINUE
+CONTINUE
+      RETURN
+      END
+====Subroutine Prodmat====
+Esecuzione del prodotto AB di due matrici A e B.
+      SUBROUTINE PRODMAT(K1,K2,K3,A,B,AB)
+      DIMENSION A(K1,K2),B(K2,K3),AB(K1,K3)
+      DO 10,I10=1,K1
+      DO 30,I30=1,K3
+      AB(I10,I30)=0.   !Inizializzazione a zero della matrice prodotto.
+      DO 20,I20=1,K2
+      AB(I10,I30)=AB(I10,I30)+A(I10,I20)*B(I20,I30)   !Prodotto.
+CONTINUE
+CONTINUE
+CONTINUE
+      RETURN
+      END
+Nota: date due matrici quadrate di ordine "n", la complessità computazionale del calcolo del loro prodotto risulta di n^3 (oneroso per il calcolatore).
+=====MAIN=====
+  C  Calcolo della tensione ideale secondo von Mises (solo in pstress)
+      SIGID1=SQRT(SIGMAEL(1)**2+SIGMAEL(2)**2-SIGMAEL(1)*SIGMAEL(2)+
+      $3.*SIGMAEL(3)**2)
+      WRITE(2,*)I30,'    ',SIGMAEL(1),'  ',SIGMAEL(2),
+      $'  ',SIGMAEL(3),'  ',SIGID1
+CONTINUE
+      STOP
+      END