cdm.ing.unimore.it

Questa pagina è in sola lettura. Puoi visualizzare il sorgente, ma non puoi modificarlo. Contatta l'amministratore se pensi che ci sia un errore.
====== FORTRAN ======

====== Controllo sui limiti di allocazione di matrici o vettori ======

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:

<code>
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
</code>

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:

<code>
   10 i=i+1
c     procedo a riempimento vettore
      vec(i)=0.01
      write (*,*) i,vec(i)
      goto 10
</code>

(è 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:

<code>
-fbounds-check
</code>

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:

<code>
gfortran -fbounds-check- pocciavec.for && ./a.out
</code>

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:

<code>
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

</code>

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

<code>
do imax=1,ilim+100
</code>

Ottenendo:

<code>
      do imax=1,ilim+100
c     procedo a riempimento vettore
      call riempi(imax,vec)
      enddo
</code>

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.

====== Allocazione Matrici con 3 dimensioni ======

In fortran è possibile costruire matrici con più di 2 dimensioni, fino ad un limite di 7. Prendiamo come esempio il seguente caso

<code>
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)
</code>

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

<code>
      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
</code>

che richiamerò all'interno del mio programma nel modo seguente dopo la parte dichiarazioni precedentemente indicata

<code>
      call riempimat(lmax,imax,jmax,mat)
</code>

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

<code>
      do l=1,lmax
       do i=1,imax
        write(*,*) 'layer',l,'row',i,(mat(l,i,j),j=1,jmax)
       enddo
      enddo
</code>

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:

<code>
       call riempimat(1,lmax*imax*jmax,1,mat)
</code>

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.

==== Trasformazione di una matrice in un vettore ====

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:

<code>
      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
</code>

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:

<code>
       call riempivec(lmax*imax*jmax,mat)
</code>

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.

======Generazione di una FUNCTION======

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.

====Struttura della Function====

La struttura di una generica FUNCTION è del tipo seguente:

<code>
<tipo> <nome funzione> (<parametri passati>)

<comandi>

return
end
</code>

Per richiamarla all'interno del programma principale usiamo lo seguente struttura:

<code>
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)       
</code>

==== Esempio di programma ====

Prendiamo il seguente caso come esempio:
Definiamo una Function Sinxx

<code>
      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
</code>

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:

<code>
      program prova
      real sinxx
      write(*,*) sinxx(2.5)
      stop
      end
</code>

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:

<code>
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

</code>

====== Codice FEM Tria3 in Fortran ======

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.

==== Fase di dichiarazione ====

<code>
      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)
</code>

  * NODES=9 : 9 nodi di struttura totali
  * NELEMS=8 : 8 elementi triangolari 3 nodi in cui la struttura è divisa
  * IFMAX=3 : 3 forze esterne nodali 
  * ICNMAX=3 : 3 nodi a cui è applicato un vincolo
  * XY(2,NODES) : Matrice che contiene le coordinate x,y dei nodi

|   |nodo 1|nodo 2|nodo 3|nodo ...|
| x |  x1  |  x2  |  x3  |        | 
| y |  y1  |  y2  |  y3  |        |

  * NVERT(3,NODES) : Matrice che contiene i label dei nodi che definiscono l'elemento (i,j,k), lega numerazione          locale a globale

esempio: 

|   |elemento 1|elemento 2|elemento 3|elemento ...|
| i |     1    |    2     |    5     |            | 
| j |     5    |    4     |    4     |            |
| k |     2    |    8     |    7     |            |

  * Matrici e vettori che vengono sempre usati a coppie: IFXY (vettore che indica a quale nodo è applicata una forza esterna),FXY (matrice 2xIFMAX che indica la componente lungo x ed y della forza applicata al nodo)

Esempio di IFXY: Forza nodale applicata al nodo 8

|nodo| ... | 8 | ... |

Esempio di FXY:

|  x  | ... | 1.5 | ... |
|  y  | ... |-2.5 | ... |

  * Matrici che vengono usate a coppie: ICNXY(matrice 2xICNMAX che indica a quale nodo è applicato il vincolo e che tipo di vincolo), CNXY(matrice 2xICNMAX che indica l'entità dello spostamento imposto al nodo vincolato)

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

  * STRUTK(2*NODES,2*NODES) : Matrice rigidezza di struttura di 2*NODES x 2*NODES elementi perché ho 2 gradi di libertà per ogni nodo
  * FORCE(2*NODES) : vettore dei termini noti del sistema
  * UV(2*NODES) : vettore degli spostamenti nodali cioè vettore delle incognite 

K*$\delta$=F coincide con STRUTK*UV=FORCE

==== Dimensionamenti fissi ====

  * D(3,3) : matrice di legame elastico

$$
   \cfrac{E}{1-\nu^2}
   \begin{bmatrix} 1 & \nu & 0 \\
                   \nu & 1 & 0 \\
                   0 & 0 & \cfrac{1-\nu}{2} \end{bmatrix}                   
$$

  * B(3,6) : Matrice di legame spostamenti-deformazioni che è diversa da elemento a elemento perciò sarà allocata di volta in volta per ogni elemento

$$
\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}
$$

  * BT(6,3) : Trasposta della matrice B
  * DB(3,6) : Matrice prodotto tra matrice D e matrice B
  * ELK(6,6) : Matrice rigidezza di elemento

$$
    \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}
$$  
  * IPOINT(6) : Vettore mappatura tra GDL locali (dell'elemento considerato) e globali
  * DELTAEL(6) : Vettore che contiene gli spostamenti nodali di elemento
  * SIGMAEL(3) : Vettore che contiene le tensioni del singolo elemento $\sigma_{x}$,$\sigma_{y}$,$\tau_{xy}$
===== Note fortran: costrutti obsoleti ======
 
===== ARITHMETIC IF =====

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.
<code>
      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
</code>

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.((da http://www.personal.psu.edu/jhm/f90/fortran.html))

===== codice odierno =====
<code>
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

</code>

<code>
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

</code>

<code>
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

</code>

===== esempio codice fortran tria3  =====

{{:wikipaom2015:main_tria_base3.for| codice tria3}}
{{:wikipaom2015:main_iso4_base3.for| codice iso4}}