PROGRAM ELEMENTI FINITI TRIANGOLARI 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) C C ****************************************************************** C Dimensionamenti fissi C ****************************************************************** C DIMENSION D(3,3),B(3,6),BT(6,3),DB(3,6),ELK(6,6),IPOINT(6), $DELTAEL(6),SIGMAEL(3) C C ****************************************************************** C Inizio programma C ****************************************************************** C 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) CALL DMAT(YM,PR,ITDP,D) C C ****************************************************************** C DO principale per la costruzione della matrice di rigidezza C dell'elemento generico I10 C C Scorro tutti gli elementi (I10=1,NELEMS); per ogni elemento prendo C il nome dei nodi che gli appartengono dalla matrice NVERT. C Scelto il nodo, prendo le coordinate x e y all'interno della C matrice XY (scegliendo tra la riga 1 o la riga 2). C ****************************************************************** C 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) 10 CONTINUE C C ****************************************************************** C Caricamento C ****************************************************************** C CALL FORCES(IFMAX,IFXY,FXY,NODES,FORCE) C C ****************************************************************** C Vincolamento C ****************************************************************** C CALL CNSTNG(ICNMAX,ICNXY,CNXY,NODES,FORCE,STRUTK) C C ****************************************************************** C Soluzione in termini di spostamenti nodali C ****************************************************************** C 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) 20 CONTINUE C C ****************************************************************** C Postprocessor per calcolare le tensioni C ****************************************************************** C 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) C C ****************************************************************** C Calcolo della tensione ideale secondo von Mises C ****************************************************************** C 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 30 CONTINUE STOP END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE READIN: Lettura del Modulo di Young e del Coefficiente C di Poisson e scelta dello stato di Tensione o Deformazione Piana. C Lettura delle coordinate nodali e dei vertici degli elementi C triangolari. Lettura delle forze nodali e delle condizioni di C vincolamento. C 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), C l'indice di nodo NODE non viene considerato, C si utilizza una numerazione basata sull'ordine: DO 10,I10=1,NODES READ(1,*)NODE,XY(1,I10),XY(2,I10) 10 CONTINUE C Lettura dei vertici degli elementi triangolari: C l'indice di elemento NELEM non viene considerato, C si utilizza una numerazione basata sull'ordine: DO 20,I20=1,NELEMS READ(1,*)NELEM,NVERT(1,I20),NVERT(2,I20),NVERT(3,I20) 20 CONTINUE C Lettura delle forze nodali: DO 30,I30=1,IFMAX READ(1,*)IFXY(I30),FXY(1,I30),FXY(2,I30) 30 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) 40 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE ECHO: Controllo della correttezza di lettura dei dati C di input. C 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) 10 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) 20 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) 30 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) 40 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE CLEAR: Inizializzazione a zero degli elementi di una C matrice. C SUBROUTINE CLEAR(K1,K2,A) DIMENSION A(K1,K2) DO 10,I10=1,K1 DO 20,I20=1,K2 A(I10,I20)=0. 20 CONTINUE 10 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE DMAT: Costruzione della matrice D di applicazione del C legame Hookeano tra tensioni e deformazioni: SIGMA = D * EPSILON. C 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: 100 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 20 CONTINUE 10 CONTINUE D(1,3)=0. D(2,3)=0. D(3,1)=0. D(3,2)=0. D(3,3)=G RETURN C Per Stato di Deformazione Piana: 200 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 40 CONTINUE 30 CONTINUE D(1,3)=0. D(2,3)=0. D(3,1)=0. D(3,2)=0. D(3,3)=G RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE KELMAT: Costruzione della matrice di rigidezza ELK C di ciascun elemento finito triangolare attraverso la formula: C F = ELK * DELTA ==> ELK = A * (BT * D * B). C Le coordinate nodali vengono prese, come di consueto, in verso C antiorario. C 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) 20 CONTINUE 10 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE BMAT: Costruzione della matrice B di collegamento tra C gli spostamenti nodali e le deformazioni di ciascun elemento C finito triangolare attraverso la formula: EPSILON = B * DELTA. C Le coordinate nodali vengono prese, come di consueto, in verso C antiorario. C Infine, poich i termini nulli, sempre nelle stesse posizioni, C vengono gi inizializzati a zero inizialmente, se ne riportano le C rispettive assegnazioni in commento. C 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 20 CONTINUE 10 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE POINTVT: Costruzione del vettore puntatore IPOINT C necessario ad effettuare l'assemblaggio della matrice ELK di C rigidezza di ciascun elemento finito triangolare nella matrice C STRUTK di rigidezza globale della struttura. C 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 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE FORCES: Riempimento del vettore globale delle forze C nodali FORCE. C 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) 10 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE CNSTNG: Imposizione delle condizioni di vincolamento C specifiche alla matrice STRUTK di rigidezza globale della C struttura ed al vettore globale delle forze nodali FORCE. C 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: 100 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. 110 CONTINUE STRUTK(INDX,INDX)=1. FORCE(INDX)=CNXY(1,I10) GO TO 10 C Per vincolamento di y: 200 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. 210 CONTINUE STRUTK(INDY,INDY)=1. FORCE(INDY)=CNXY(2,I10) GO TO 10 C Per vincolamento di x e y: 300 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. 310 CONTINUE STRUTK(INDX,INDX)=1. STRUTK(INDY,INDY)=1. FORCE(INDX)=CNXY(1,I10) FORCE(INDY)=CNXY(2,I10) GO TO 10 10 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE AREA: Calcolo della superficie DEL di ciascun elemento C finito triangolare a partire dalle coordinate nodali dei vertici. C Controllo di dimensioni e segno dell'area calcolata. C Le coordinate nodali vengono prese, come di consueto, in verso C antiorario. C Con TWODEL si identifica il doppio della superficie. C 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 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE PRODMAT: Esecuzione del prodotto AB di due matrici C A e B. C 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. 20 CONTINUE 30 CONTINUE 10 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE TRSPMAT: Calcolo della matrice trasposta AT di una C matrice A. C SUBROUTINE TRSPMAT(K1,K2,A,AT) DIMENSION A(K1,K2),AT(K2,K1) DO 10,I10=1,K1 DO 20,I20=1,K2 AT(I20,I10)=A(I10,I20) 20 CONTINUE 10 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE ASSEMBL: Assemblaggio della matrice ELK di rigidezza di C ciascun elemento finito triangolare nella matrice STRUTK di C rigidezza globale della struttura. C La matrice STRUTK viene inzializzata a zero nel MAIN PROGRAM, C pertanto si ritiene superflua la sua inizializzazione in questa C subroutine. C 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) 20 CONTINUE 10 CONTINUE RETURN END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SUBROUTINE GAUSS: Risoluzione di un sistema algebrico lineare di C N equazioni in N incognite attraverso il metodo di Gauss C (triangolarizzazione della matrice dei coefficienti). C La matrice dei coefficienti viene denominata C; il vettore dei C termini noti B ed il vettore delle incognite V. C 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