Prendiamo il seguente programma per allocare un vettore di “ilim” elementi e vediamo cosa succederebbe se provassimo ad utilizzare un numero maggiore di elementi rispetto a quelli dichiarati:
c234567==============================================================72!
program pocciavec
implicit none
c allocazione di un vettore di ilim elementi
real vec
integer i,ilim
parameter(ilim=10)
dimension vec(ilim)
c noto ora che invertendo l'ordine delle due righe sopra gfortran (né ifort) non compilano più
c----- fine dichiarazioni, inizio istruzioni ---------------------------
c inizializzo la variabile i al valore zero
i=0
c inserisco un ciclo infinito che legge i valori casuali del vettore
10 i=i+1
c procedo a riempimento vettore
write (*,*) i,vec(i)
goto 10
stop
end
Nel programma si dichiara un vettore vec di ilim=10 elementi tuttavia, attraverso l’istruzione go to, si incrementa la variabile di conteggio i un numero illimitato di volte leggendo di volta in volta valori casuali del vettore. Attivando il programma da terminale, superato il numero di elementi dichiarato, il FORTRAN comincerà a stampare numeri presi dalla memoria ad esso dedicata, pari a 4 GB, fino al suo esaurimento e terminerà con la stampa dell'errore: “Segmentation fault” che indica che il programma sta accedendo a una memoria non assegnata. A differenza dell’output che è gestito dal programma, l’errore è invece gestito dal sistema operativo (nello specifico è il kernel del sistema operativo che interviene), che blocca tutto quando il programma tenta l'accesso a una memoria superiore a quella ad esso assegnata.
Nel caso in cui:
10 i=i+1
c procedo a riempimento vettore
vec(i)=0.01
write (*,*) i,vec(i)
goto 10
(è stata inserita la riga vec(i)=0.01).
Da terminale sarebbero letti solamente 11 elementi del vettore, tutti pari a 9.999999E-3 (=0.01).
Se invece di vec(i)=0.01 avessimo scritto vec(i)=0.00, da terminale il vettore sarebbe stato letto ciclicamente all'infinito ripetendo le righe dalla 1 alla 11 con elementi tutti nulli; questo probabilmente perché con questa istruzione vado ad azzerare anche il contatore i perciò di volta in volta riscrivo e rileggo sempre gli stessi elementi del vettore (sempre gli stessi spazi di memoria).
Quindi è molto rischioso lavorare fuori dai limiti di un vettore perché potrei generare dei loop pericolosi fino a provocare un blocco del programma oppure un ciclo infinito di operazioni.
Per evitare questi problemi devo aggiungere un opzione di controllo sui limiti di allocazione di un vettore ovvero:
-fbounds-check
Questa opzione( che compila controllando i limiti delle matrici o dei vettori e avverte qualora ci siano questi tipi di problemi) deve essere aggiunta quando il programma viene richiamato da terminale, nel modo seguente:
gfortran -fbounds-check- pocciavec.for && ./a.out
Il codice si interrompe dopo il decimo elemento del vettore (così come dovrebbe essere) senza provocare alcun tipo di problema.
Usando la SUBROUTINE per costruire lo stesso programma:
c234567==============================================================72!
program pocciavec
implicit none
c allocazione di un vettore di ilim elementi
real vec
integer i,ilim,imax
c
parameter(ilim=10)
dimension vec(ilim)
c----- fine dichiarazioni, inizio istruzioni ---------------------------
c definisco una variabile imax<=ilim
do imax=1,ilim
c procedo a riempimento vettore
call riempi(imax,vec)
enddo
stop
end
subroutine riempi(n,vec)
implicit none
integer i,n
real vec
dimension vec(n)
do i=1,n
vec(i)=2.0*i
enddo
write(*,*) (vec(i), i=1,n)
return
end
Il programma riempie il vettore vec di elementi 2.0*i tramite subroutine. Nella chiamata del programma da terminale, senza BOUNDS CHECK, non ho nessun problema perchè il vettore è stato specificato di 10 elementi in fase di dichiarazione e così risulta essere anche durante la chiamata della subroutine (n=ilim, il ciclo do legge 10 elementi).
Modificando il ciclo do del programma prima della chiamata della subroutine con
do imax=1,ilim+100
Ottenendo:
do imax=1,ilim+100
c procedo a riempimento vettore
call riempi(imax,vec)
enddo
e avviando il programma da terminale, sempre senza BOUNDS CHECK, mi da errore perché finisco in memorie in cui non posso scrivere perché non sono state assegnate al vettore durante la dichiarazione. Infatti chiedo di leggere ilim+100 (10+100) elementi quando il vettore è stato dichiarato di soli ilim=10 elementi.
Se provo a dare il BOUNDS CHECK durante la chiamata del programma da terminale, il vettore dovrebbe terminare dopo il 10° elemento ma così non è perché la funzione bounds check risolve il problema solamente se le matrici o i vettori sono dichiarati all'interno di uno stesso blocco, nel nostro caso il vettore è dichiarato nella subroutine e richiamato nel programma perciò in due blocchi diversi.
In fortran è possibile costruire matrici con più di 2 dimensioni, fino ad un limite di 7. Prendiamo come esempio il seguente caso
c234567==============================================================72!
program leaddim
implicit none
integer l,i,j,lmax,imax,jmax
real mat
parameter(lmax=2,imax=3,jmax=4)
c allocazione un (iper)array a 3 dimensioni
dimension mat(lmax,imax,jmax)
Il codice sopra indicato rappresenta la parte di dichiarazioni di un programma per la costruzione di una matrice in 3 dimensioni dove con lmax si indica il numero di LAYER, con imax il numero di RIGHE e con jmax il numero di COLONNE. Possiamo notare che la matrice può anche essere dichiarata specificandone la dimensione in un secondo momento tramite il comando DIMENSION come sopra indicato.
Per riempire la matrice posso usare il seguente ciclo “do” nidificato attraverso una subroutine
subroutine riempimat(lmax,imax,jmax,mat)
implicit none
real mat
integer l,i,j,lmax,imax,jmax
dimension mat(lmax,imax,jmax)
do l=1,lmax
do i=1,imax
do j=1,jmax
mat(l,i,j)= 100.*l + 1.*i + 0.01*j
enddo
enddo
enddo
return
end
che richiamerò all'interno del mio programma nel modo seguente dopo la parte dichiarazioni precedentemente indicata
call riempimat(lmax,imax,jmax,mat)
Infine per visualizzare i risultati su terminale una volta che il programma verrà chiamato dovtò usare un secondo ciclo “do” nidificato che scorre su colonne, righe ed infine layer
do l=1,lmax
do i=1,imax
write(*,*) 'layer',l,'row',i,(mat(l,i,j),j=1,jmax)
enddo
enddo
Inoltre è possibile, durante la chiamata della subroutine nel programma, non rispettare le dimensioni esatte della matrice come specificato nella subroutine stessa purché sia rispettato il numero totale di elementi da inserire nella matrice. Per esempio avrei potuto chiamare la subroutine nel modo seguente:
call riempimat(1,lmax*imax*jmax,1,mat)
In questo modo non avrò più una matrice 2x3x4 come dichiarato da programma ma otterrò una matrice 1x24x1 che, come è possibile notare, avrà lo stesso numero di elementi di quella precedente.
Un'ultima cosa che posso fare è di cambiare anche la forma della matrice, cioè di trasformarla per esempio in un vettore, perciò riduco le dimensioni da 3 ad 1. Per eseguire questo procedimento posso operare con la seguente subroutine:
subroutine riempivec(n,vec)
implicit none
real vec
integer i,n
dimension vec(n)
do i=1,n
vec(i)= 1.* i
enddo
return
end
Questa subroutine prende in ingresso due variabili, una reale(vec) ed una intera(n) e riempie il vettore vec di elementi 1.*i tramite il ciclo “do” indicato. Per richiamare questa subroutine all'interno del mio programma scrivo:
call riempivec(lmax*imax*jmax,mat)
Il quale prende in ingresso la mia matrice 2x3x4 di 24 elementi e la trasforma in un vettore di 24 elementi. Questo è possibile sempre se il numero totale di elementi prima e dopo combacia.
Il linguaggio di programmazione FORTRAN consente di suddividere le operazioni tramite sottoprogrammi che ritornano un risultato al programma chiamante, sviluppando un approccio top-down. L'utilizzo dei sottoprogrammi è molto utile quando si deve ripetere una stessa operazione più volte, evitando così tutte le volte di copiare e incollare nel programma principale la stessa operazione. Questi sottoprogrammi possono essere del tipo FUNZIONI oppure SUBROUTINE. Funzioni e subroutine devono avere nome univoco e per il resto sono indipendenti dal programma che le ha richiamate. Il costrutto FUNZIONE restituisce un singolo valore scalare, motivo per cui spesso si ricorre al costrutto SUBROUTINE. A differenza dei programmi che terminano con “ stop end”, i sottoprogrammi terminano con “return end”. Nel sottoprogramma le etichette sono autonome perciò per dichiarare le variabili locali, si possono utilizzare etichette aventi stesso nome e tipologia di quelle già dichiarate nel programma principale (MAIN). Programma e sottoprogramma comunicano inviandosi indirizzi di memoria perciò è fondamentale che comunichino con la stessa tipologia di elementi: in particolare il programma invia i parametri passati al sottoprogramma che restituisce un valore.
In questo paragrafo viene trattata solo la FUNCTION dato che la spiegazione della subroutine è già contenuta in lezioni precedenti.
la FUNCTION può essere considerata una subroutine “speciale” che restituisce un singolo valore scalare. E' indispensabile dichiarare la natura delle FUNZIONI in base al risultato che riportano al programma chiamante e consistente in una delle 6 tipologie utilizzate per dichiarare le variabili nel programma principale quindi : integer, real, double precision, complex, logical e character. La dichiarazione del tipo di funzione è seguita da un nome assegnato e dai parametri di comunicazione col programma principale.
La struttura di una generica FUNCTION è del tipo seguente:
<tipo> <nome funzione> (<parametri passati>) <comandi> return end
Per richiamarla all'interno del programma principale usiamo lo seguente struttura:
c Prima è necessario definirla in fase di dichiarazione
real <nome funzione>
c Dopodiché è possibile usarla come una qualsiasi funzione implicita a fortran, per calcolare un valore scalare
<nome funzione> (numero o parametro)
Prendiamo il seguente caso come esempio: Definiamo una Function Sinxx
real function sinxx(var)
implicit none
real var
if (var .ne. 0.0) then
sinxx=sin(var)/var
else
sinxx=1.0
endif
return
end
La function va inserita sempre all'esterno del program, esattamente come una subroutine. In questo caso la funzione ha il compito di calcolare:
$$ \frac{sin(var)}{var}\ $$
Restituisce il valore esatto nel caso in cui “var” sia diversa da 0, altrimenti restituisce il valore 1 perchè in var=0 la funzione non è definite perciò ne fa il limite per var–>0
IMPORTANTE!!: Quando arrivo al return la variabile omonomia alla function (in tal caso sinxx) deve avere qualcosa di assegnato (in questo caso sinxx=sin(var)/var), se questo non accade la function non funziona.
Per richiamare la FUNCTION all'interno del programma:
program prova
real sinxx
write(*,*) sinxx(2.5)
stop
end
Vedo che è necessario ridichiarare all'interno del programma una variabile con lo stesso nome della function (sinxx), che poi usandola semplicemente come un qualsiasi comando implicito a fortran mi calcola il valore della funzione in 2.5.
Il programma completo quindi risulta essere:
c234567==============================================================72!
c dichiarazione di una function
real function sinxx(var)
implicit none
real var
if (var .ne. 0.0) then
sinxx=sin(var)/var
else
sinxx=1.0
endif
return
end
program prova
real sinxx
write(*,*) sinxx(2.5)
stop
end
Iniziamo ad implementare un codice ad elementi finiti in fortran con elementi triangolari a 3 nodi. In questo modo potremo avere a disposizione un programma opensource per eseguire verifiche resistenziali su oggetti.
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)
| nodo 1 | nodo 2 | nodo 3 | nodo … | |
| x | x1 | x2 | x3 | |
| y | y1 | y2 | y3 |
esempio:
| elemento 1 | elemento 2 | elemento 3 | elemento … | |
| i | 1 | 2 | 5 | |
| j | 5 | 4 | 4 | |
| k | 2 | 8 | 7 |
Esempio di IFXY: Forza nodale applicata al nodo 8
| nodo | … | 8 | … |
Esempio di FXY:
| x | … | 1.5 | … |
| y | … | -2.5 | … |
Esempio di ICNXY: Vincolo con carrello ad asse verticale applicato al nodo 13
| nodo | … | 13 | … |
| tipo di vincolo | … | 2 | … |
nella riga del tipo di vincolo: 1 è vincolo ad asse verticale, 2 è vincolo ad asse orizzontale, 3 è vincolo in entrambe le direzioni (cerniera)
Esempio di CNXY:
| U | … | ? | … |
| V | … | 1.0 | … |
spostamento imposto di 1 mm al nodo 13 lungo y; lungo x è stato inserito “?” perchè può essere un valore qualsiasi dato che in teoria non dovrei leggerlo
K*$\delta$=F coincide con STRUTK*UV=FORCE
$$ \cfrac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \cfrac{1-\nu}{2} \end{bmatrix} $$
$$ \frac{1}{2A} \begin{bmatrix} y_j - y_k & 0 & y_k - y_i & 0 & y_i - y_j & 0 \\ 0 & x_k - x_j & 0 & x_i - x_k & 0 & x_j - x_i \\ x_k - x_j & y_j - y_k & x_i - x_k & y_k - y_i & x_j - x_i & y_i - y_j \end{bmatrix} $$
$$ \begin{bmatrix} k11 & k12 & k13 & k14 & k15 & k16 \\ k21 & k22 & k23 & k24 & k25 & k26 \\ k31 & k32 & k33 & k34 & k35 & k36 \\ k41 & k42 & k43 & k44 & k45 & k46 \\ k51 & k52 & k53 & k54 & k55 & k56 \\ k61 & k62 & k63 & k64 & k65 & k66 \end{bmatrix} $$
While we are on the subject of IF's, you should be aware of one archaic form called the arithmetic IF. I introduce this only because some of you will be asked to improve an existing program. You should not use this in any new programming. There is a good chance that it will be dropped from some future Fortran standard. The arithmetic IF can be recognized by the fact that the expression within the parentheses results in a numerical result, not a logical (true, false) result. The other give-away is the string of three integers (statement label references) after the right parenthesis.
if ( b**2-4*a*c) 100,200,300
100 print *, 'Roots are Complex'
go to 400
200 print *, 'Single Real Root'
go to 400
300 print *, 'Two Real Roots'
400 continue
If the arithmetic expression produces a value less than zero, execution branches to the first label in the list. If it is zero, control branches to the second label, and if the value is greater than zero the program branches to the third label. Repetition of label numbers is permitted in the IF's label list. Note that I must include extra “GO TO” statements to prevent execution of some of the PRINT statements after the initial branch to another PRINT.1)
c234567==============================================================72!
program leaddim
implicit none
integer l,i,j,lmax,imax,jmax
real mat
parameter(lmax=2,imax=3,jmax=4)
c allocazione un (iper)array a 3 dimensioni
dimension mat(lmax,imax,jmax)
c----- fine dichiarazioni, inizio istruzioni ---------------------------
c call riempimat(lmax,imax,jmax,mat)
c call riempimat(1,lmax*imax*jmax,1,mat)
call riempivec(lmax*imax*jmax,mat)
cc riempio un elemento a piacere
c mat(2,1,3) = 1.234567
do l=1,lmax
do i=1,imax
write(*,*) 'layer',l,'row',i,(mat(l,i,j),j=1,jmax)
enddo
enddo
stop
end
c l'opzione -fbounds-check inserisce un controllo sui limiti di allocazione di un vettore
c se non inizializzo entro il codice, con -finit-real=nan forzo un'inizializzazione da compilatore post allocazione
subroutine riempimat(lmax,imax,jmax,mat)
implicit none
real mat
integer l,i,j,lmax,imax,jmax
dimension mat(lmax,imax,jmax)
do l=1,lmax
do i=1,imax
do j=1,jmax
mat(l,i,j)= 100.*l + 1.*i + 0.01*j
enddo
enddo
enddo
return
end
subroutine riempivec(n,vec)
implicit none
real vec
integer i,n
dimension vec(n)
do i=1,n
vec(i)= 1.* i
enddo
return
end
c234567==============================================================72!
c dichiarazione di una function
real function sinxx(var)
implicit none
real var
if (var .ne. 0.0) then
sinxx=sin(var)/var
else
sinxx=1.0
endif
return
end
program prova
real sinxx
external sinxx
write(*,*) sinxx(2.5)
stop
end
c234567==============================================================72!
program pocciavec
implicit none
c allocazione di un vettore di ilim elementi
real vec
integer i,ilim,imax
c
parameter(ilim=10)
dimension vec(ilim)
c noto ora che invertendo l'ordine delle due righe sopra gfortran (né ifort) non compilano più
c----- fine dichiarazioni, inizio istruzioni ---------------------------
c definisco una variabile imax<=ilim
do imax=1,ilim+100
c procedo a riempimento vettore
call riempi(imax,vec)
enddo
stop
end
subroutine riempi(n,vec)
implicit none
integer i,n
real vec
dimension vec(n)
do i=1,n
vec(i)=2.0*i
enddo
write(*,*) (vec(i), i=1,n)
return
end