Costruzione della topologia della matrice di rigidezza

Il metodo degli elementi finiti applicato alla soluzione di un'equazione differenziale lineare e stazionaria, cioè indipendente dal tempo, genera un sistema lineare la cui soluzione rappresenta il vettore dei valori che la funzione incognita assume sui nodi con i quali è stato discretizzato il dominio di integrazione. La matrice di tale sistema è generalmente chiamata matrice di rigidezza e nel seguito verrà indicata con il simbolo $H$. Ad esempio, nella soluzione del problema dell'equilibrio elastico di una struttura gli elementi di $H$ corrispondono alle rigidezze associate ai nodi della griglia di calcolo e risultano funzione del modulo di Young e del rapporto di Poisson del materiale che costituisce il corpo in esame. Analogamente, nella soluzione di un problema di flusso $H$ contiene i termini di trasmissività del corpo drenante, funzione della permeabilità e degli spessori associati a ciascun nodo della griglia.

La topologia della matrice di rigidezza, vale a dire l'ubicazione dei coefficienti non nulli in $H$, descritta nella memorizzazione compatta CRS dai vettori interi JA e IA, è determinata univocamente dalla mesh di calcolo con cui si è discretizzato il dominio di interesse. In particolare, è la rete dei contatti nodali nella maglia computazionale a stabilire la posizione degli elementi non nulli di $H$. Gli elementi non nulli della riga $i$ corrispondono infatti agli indici di colonna $j$ dei nodi con cui $i$ risulta in contatto.

Per chiarire il concetto, ci si riferisca al patch di elementi riportato in Figura 1. Ad esempio, nella riga 7, relativa al valore assunto dalla funzione incognita nel nodo 7, gli elementi di $H$ non nulli saranno:

$$h_{7,3} \;\;\; h_{7,6} \;\;\; h_{7,7} \;\;\; h_{7,8} \;\;\; h_{7,11} \;\;\; h_{7,12}$$

Figura 1: Patch di elementi triangolari. In corsivo sono indicati gli indici dei nodi ed in grassetto gli indici degli elementi.

Ciò è naturale conseguenza del fatto che i polinomi che interpolano la soluzione approssimata in ciascun elemento finito hanno supporto locale. È ovvio che se il nodo $i$ è in contatto con il nodo $j$, anche $j$ sarà in contatto con $i$, e quindi la matrice $H$ avrà una topologia simmetrica. Si osservi che, tuttavia, la simmetria della topologia non implica necessariamente la simmetria di $H$ che, invece, nei casi di nostro interesse è assicurata dalla forma dei contributi locali, come si vedrà meglio in seguito. Risulta pertanto sufficiente memorizzare le posizioni $h_{i,j}$ per cui $j \geq i$.

Figura 2: Patch di due elementi triangolari.

A partire dalla tabella dei contatti nodali, definita dalla maglia computazionale, è dunque possibile costruire i vettori JA e IA che individuano la topologia di $H$. Vediamo nel dettaglio l'implementazione di questa procedura sull'esempio di Figura 2. La mesh di calcolo viene descritta dalla successione dei nodi e dalla topologia degli elementi. Avremo pertanto in input le seguenti strutture dati:

• coordinate nodali (es. file coord generato da MeshMaker):

  	4		: numero di nodi $n$
  1	$x_1$	$y_1$
  2	$x_2$	$y_2$
  3	$x_3$	$y_3$
  4	$x_4$	$y_4$

• successione dei nodi in senso antiorario e indice di zona per ciascun elemento (es. file topol generato da MeshMaker):

  	2				: numero di elementi $ne$
  1	2	1	4	1
  2	4	3	2	1

Dai contatti nodali, si deduce che la matrice di rigidezza, di dimensione $n \times n$, avrà la seguente struttura:

$$\left[ \begin{array}{cccc} h_{1,1} & h_{1,2} & 0 & h_{1,4} \\ h_... ..._{3,3} & h_{3,4} \\ h_{4,1} & h_{4,2} & h_{4,3} & h_{4,4} \end{array} \right]$$

Grazie alla simmetria di $H$, si memorizza la sola parte triangolare alta inclusa la diagonale principale. Si verifica immediatamente che il vettore JA degli indici di colonna di $nt=9$ componenti risulta:

   JA$$= 1 \;\; 2 \;\; 4 \;\; 2 \;\; 3 \;\; 4 \;\; 3 \;\; 4 \;\; 4$$

mentre il vettore IA di $n+1=5$ componenti contenente le posizioni in SYSMAT, e quindi in JA, degli elementi diagonali è:

   IA$$= 1 \;\; 4 \;\; 7 \;\; 9 \;\; 10$$

Applichiamo ora una procedura automatica che consenta di costruire JA dal file topol. Definiamo $n1$ come il massimo numero di contatti nodali ammessi dalla griglia del problema, cioè il massimo numero di elementi non nulli che possiamo trovare su ciascuna riga di $H$. Nel caso in esame, $n1=4$ e di conseguenza il massimo numero di elementi non nulli da memorizzare sarà $nt=n1 \cdot n=16$. Il vettore JA viene generato come una successione di $nt$ componenti intere, inizialmente nulle, che possiamo interpretare come $n$ celle, corrispondenti alle $n$ righe di $H$, ciascuna composta di $n1$ elementi. In esse andranno inseriti gli indici di colonna degli elementi non nulli incontrati in ciascuna riga:

$$\boxed{0 \;\; 0 \;\; 0 \;\; 0}\boxed{0 \;\; 0 \;\; 0 \;\; 0}\boxed{0 \;\; 0 \;\; 0 \;\; 0}\boxed{0 \;\; 0 \;\; 0 \;\; 0}$$

Il primo elemento di ogni cella corrisponde all'elemento diagonale, il cui indice di colonna $j$ è uguale all'indice di riga $i$. Ciò significa che la prima componente di ogni cella deve essere pari all'indice della cella stessa, cioè sarà 1 nella prima cella, 2 nella seconda, e così via:

$$\boxed{1 \;\; 0 \;\; 0 \;\; 0}\boxed{2 \;\; 0 \;\; 0 \;\; 0}\boxed{3 \;\; 0 \;\; 0 \;\; 0}\boxed{4 \;\; 0 \;\; 0 \;\; 0}$$

Per determinare i vari contatti nodali si esegue un ciclo sugli elementi della griglia di calcolo. Si consideri il triangolo 1 definito dalla successione di nodi 2, 1, 4. Poiché nella cella $i$ vanno inseriti solamente i contatti nodali relativi alla parte triangolare alta di $H$, cioè gli indici $j$ tali che $j \geq i$, conviene ordinare gli indici dei nodi in senso crescente:

$$2 \;\; 1 \;\; 4 \;\; \Rightarrow \;\; 1 \;\; 2 \;\; 4$$

Si osservi che tale successione ordinata è utile solamente ai fini della determinazione della topologia di $H$, mentre per il calcolo dei contributi locali va mantenuta la successione originaria in senso antiorario. È consigliabile, pertanto, fare uso di un vettore ausiliario locale I1(3) in cui memorizzare provvisoriamente la sequenza nodale ordinata in senso crescente. Con un ciclo sulle componenti di I1 appare ora chiaro che nella prima cella vanno aggiunti i contatti nodali 2 e 4, nella seconda va aggiunto il contatto 4 e nella quarta non va aggiunto nulla:

$$\boxed{1 \;\; 2 \;\; 4 \;\; 0}\boxed{2 \;\; 4 \;\; 0 \;\; 0}\boxed{3 \;\; 0 \;\; 0 \;\; 0}\boxed{4 \;\; 0 \;\; 0 \;\; 0}$$

Seguiamo la medesima procedura per il triangolo 2. Viene creato il vettore I1 con la sequenza nodale ordinata in senso crescente:

$$4 \;\; 3 \;\; 2 \;\; \Rightarrow \;\; 2 \;\; 3 \;\; 4$$

Nella seconda cella vanno aggiunti i contatti ai nodi 3 e 4. A causa del modo in cui vengono memorizzati gli indici di colonna in JA, l'intero 3 va inserito fra il 2 ed il 4 già presenti, spostando quindi l'ultimo indice avanti di una posizione. Inoltre, poiché il contatto con il nodo 4 è già stato individuato, null'altro va aggiunto alla seconda cella. Infine, nella terza cella va introdotto il contatto con il nodo 4:

$$\boxed{1 \;\; 2 \;\; 4 \;\; 0}\boxed{2 \;\; 3 \;\; 4 \;\; 0}\boxed{3 \;\; 4 \;\; 0 \;\; 0}\boxed{4 \;\; 0 \;\; 0 \;\; 0}$$

Va da sé che nell'ultima cella non andrà mai aggiunto nessun contatto. Si noti ora che il vettore così determinato corrisponde, a meno degli zeri, al vettore JA cercato.

La procedura seguita per la costruzione di JA può quindi essere riassunta nel seguente algoritmo:

001 $nt:=n1 \cdot n$
002 Per $i=1,nt$
003 JA$(i):=0$
004 Fine Per
005 Per $i=1,nt$ con passo $n1$
006 JA $(i):=\left(i+n1-1\right)/n1$
007 Fine Per
008 Per $k=1,ne$
009 ordinamento nodi in senso crescente in I1(3)
010 Per $j=1,2$
011 $m:=\left(\mbox{\tt I1}(j)-1\right) \cdot n1+1$
012 Per $i=j+1,3$
013 $m:=m+1$
014 Se JA$(m)=0$
015 JA$(m):=$ I1$(i)$
016 vai all'istruzione 030
017 Fine Se
018 Se JA$(m)<$ I1$(i)$
019 vai all'istruzione 013
020 Fine Se
021 Se JA$(m)>$ I1$(i)$
022 sposta gli elementi della cella una posizione in avanti
023 controlla che la posizione dell'ultimo elemento sia $< n1 \cdot$   I1$(i)$
024 JA$(m):=$ I1$(i)$
025 vai all'istruzione 030
026 Fine Se
027 Se JA$(m)=$ I1$(i)$
028 vai all'istruzione 030
029 Fine Se
030 continua
031 Fine Per
032 Fine Per
033 Fine Per

La costruzione del vettore IA è ora quasi immediata. La posizione dell'elemento diagonale della riga $i$ corrisponde, infatti, al primo elemento della cella $i$ in JA, vale a dire $(i-1) \cdot n1+1$. Il valore così determinato va quindi depurato del numero di zeri rimasti nelle celle precedenti. Il corrispondente algoritmo può pertanto essere scritto nel modo seguente:

001 IA(1):=1
002 Per $i=1,nt$ con passo $n1$
003 $m:=0$
004 Per $j=1,n1$
005 Se JA $(i+j-1) \neq 0$
006 $m:=m+1$
007 Fine Se
008 Fine Per
009 $k:=(i+n1-1)/n1+1$
010 IA$(k):=$ IA$(k-1)+m$
011 Fine Per

Si noti che mediante l'algoritmo appena riportato la componente $n\!+\!1$-esima di IA viene calcolata automaticamente.

Infine, è necessario compattare il vettore JA eliminando tutti gli zeri in esso contenuti. Per far questo è sufficiente scorrere le componenti di JA e, quando si incontra una componente nulla, spostare tutte le successive indietro di una posizione. Il numero di termini non nulli di JA, che alla fine della procedura sono raggruppati nelle prime posizioni, è il vero $nt$.

Per eseguire in modo più efficiente questa operazione conviene, mano a mano che si scorrono le componenti di JA, contare i termini non nulli e spostare l'$m$-esimo termine diverso da zero in posizione $m$. Tale posizione sarà necessariamente minore o tutt'al più uguale all'indice della componente raggiunta nello scorrimento. Per effettuare la compattazione di JA si può pertanto utilizzare il seguente algoritmo:

001 $m:=0$
002 Per $i=1,nt$
003 Se JA $(i) \neq 0$
004 $m:=m+1$
005 JA$(m):=$ JA$(i)$
006 Fine Se
007 Fine Per
008 $nt:=m$

Poiché sappiamo che IA$(n+1)$ deve essere uguale a $nt+1$, ma i due termini sono stati calcolati indipendentemente uno dall'altro, è sempre consigliabile controllare la consistenza della procedura verificando che effettivamente tale condizione sia rispettata.

Un'implementazione efficiente degli algoritmi per la costruzione della topologia della matrice di rigidezza è proposta nella subroutine TOPOL.F, messa a disposizione dello studente.