temp

CLAUDE.md

@@ -23,40 +23,59 @@ python simulate.py
 
 # Visualizzare il primo run
 python plot_csv.py
+
+# Animazione della superficie sviluppata e vista 3D isometrica
+python plot_animazione.py
+python plot_animazione_3d.py
 ```
 
 Attivare sempre il venv (`source .venv/bin/activate`) prima di eseguire qualsiasi comando Python.
 
-Non sono configurati test o linter.
+Non sono configurati linter. La validazione del solutore si esegue a mano:
+
+```bash
+python valida_shell.py
+```
+
+Va rieseguita dopo ogni modifica alla fisica o all'assemblaggio: confronta le matrici e il solutore con soluzioni analitiche (raffreddamento convettivo, modi sinusoidali circonferenziali, conservazione dell'energia, potenza della sorgente, periodicità) e stampa gli errori.
 
 ## Architettura
 
 Generatore di dataset per misurazioni termiche pseudo-realistiche di una fascetta (anello cilindrico sottile) riscaldata da sorgenti a induzione in movimento.
 
-**Geometria:** la fascetta ha diametro, spessore e lunghezza configurabili. Il dominio simulato è la sezione rettangolare lunghezza × spessore, con origine (0, 0) nel vertice in alto a sinistra: x = lunghezza (le sorgenti si muovono in direzione -x sul lato esterno), z = spessore (0 = lato esterno, spessore = lato interno). La coordinata circonferenziale y non è risolta: l'offset y delle sorgenti è collassato in un'attenuazione gaussiana del flusso. Il sensore è un pirometro a infrarossi dentro la fascetta che misura la superficie interna in un punto x fisso.
+**Geometria:** la fascetta ha diametro medio, spessore e lunghezza configurabili. Il dominio simulato è l'intera superficie cilindrica, parametrizzata da x (lunghezza, le sorgenti si muovono in direzione -x) e theta (circonferenza, periodica). Lo spessore non è discretizzato: la parete è termicamente sottile, quindi faccia esterna e faccia interna hanno la stessa temperatura. Il sensore è un pirometro a infrarossi dentro la fascetta che misura la superficie interna in un punto x fisso, sulla generatrice theta = 0; l'offset circonferenziale delle sorgenti è misurato rispetto a quella generatrice.
 
 **Flusso dei dati:**
 
 1. `config.py` — tutti i parametri configurabili (dizionari SIMULAZIONE, FASCETTA, ARIA, SORGENTE, SENSORE, RANDOMIZZAZIONE)
 2. `materials.py` — dizionario MATERIALI con proprietà termofisiche ed elettriche per materiale
-3. `simulate.py` — motore principale: genera N run randomizzati, scrive i CSV, scrive `metadata.csv`
-4. `plot_csv.py` — visualizzazione autonoma per un singolo run
+3. `mesh.py` — nodi, connettività Q4 e periodicità circonferenziale della shell cilindrica
+4. `elementi_shell.py` — funzioni di forma bilineari, quadratura di Gauss 2×2, matrici del singolo elemento
+5. `assemblaggio.py` — assemblaggio delle matrici globali e del vettore della sorgente
+6. `sorgente.py` — moto del gruppo di sorgenti e distribuzione gaussiana del flusso
+7. `solutore.py` — Eulero implicito e fattorizzazione LU sparsa
+8. `sensore.py` — interpolazione della temperatura, inerzia, rumore, quantizzazione
+9. `simulate.py` — motore principale: genera N run randomizzati, scrive i CSV, scrive `metadata.csv`
+10. `plot_csv.py`, `plot_animazione.py`, `plot_animazione_3d.py` — visualizzazioni autonome
 
-**Pipeline fisica dentro `simula_singolo()` in [simulate.py](simulate.py):**
+**Pipeline fisica in [simulate.py](simulate.py):**
 
-- La skin depth è calcolata dalla resistività elettrica del materiale e dalla frequenza di induzione (`calcola_skin_depth_m`)
-- Le sorgenti gaussiane in movimento producono un profilo di flusso termico superficiale q(x) sul lato esterno, variabile nel tempo (`profilo_flusso_incidente_W_m2`)
-- Quel flusso è ridistribuito volumetricamente attraverso lo spessore con decadimento esponenziale in z (`profilo_deposizione_z_1_m`): q_vol(x, z) = q(x) · p(z)
-- Uno schema 2D a volumi finiti con Eulero implicito integra l'equazione del calore su `n_nodi_x × n_nodi_z` celle: `prepara_stato_termico` costruisce griglia, coefficienti e matrice sparsa fattorizzata LU una volta per run (`costruisci_solutore_implicito_2d`, che restituisce l'oggetto `splu`), poi `passo_implicito` avanza il campo risolvendo solo il sistema triangolare
-- Le condizioni al contorno sono incorporate nella matrice: convezione su tutti e quattro i lati della sezione, più un termine di conduzione circonferenziale (y) verso il resto della fascetta assunto a temperatura ambiente. Il termine è un'equazione di aletta ricavata sull'intero volume del cilindro: il calore conduce lungo y attraverso l'intero spessore mentre le superfici esterna e interna dell'intero cilindro scambiano per convezione, dando q_y = -(h_esterno + h_interno)/spessore · (T - T_amb), senza parametri di conduzione y configurabili a parte
+- `prepara_stato_termico` costruisce la mesh, assembla C, K_cond e K_conv, fattorizza LU una volta per run e prepara l'assemblatore della sorgente e l'interpolatore del sensore
+- La conduzione assiale e quella circonferenziale nascono dall'assemblaggio: elementi adiacenti condividono i nodi del bordo comune, quindi non serve alcun accoppiamento esplicito tra elementi
+- Convezione su entrambe le facce (h_esterno + h_interno su tutta l'area) e sui due bordi assiali (h_bordi · spessore lungo il perimetro), entrambe incorporate in K_conv e nel termine noto costante `f_ambiente`
+- Le sorgenti gaussiane in movimento producono q''(x, theta, t); la distanza circonferenziale è riavvolta su 2πR, quindi la sorgente scalda anche i punti raggiunti dall'altro verso. Il vettore di carico è ∫ N^T q'' dA, integrato con la stessa quadratura di Gauss delle matrici
+- Il flusso è separabile in un fattore assiale variabile nel tempo e uno circonferenziale costante: `assembla_sorgente` sfrutta la separazione, non è un'approssimazione
+- `passo_termico` avanza il campo risolvendo (C/dt + K) T = C/dt · T_prec + f_ambiente + f_sorgente
 - La temperatura iniziale del campo è la temperatura ambiente (randomizzata per run)
 - L'output del sensore aggiunge inerzia del primo ordine, rumore gaussiano e quantizzazione
 
-`prepara_stato_termico` e `passo_implicito` sono condivisi con `plot_animazione.py`, che riproduce la fisica di run_0001 per animare la sezione: ogni modifica alla fisica va fatta lì, non duplicata.
+La skin depth (`calcola_skin_depth_m`) e il numero di Fourier sullo spessore (`numero_fourier_spessore`) sono calcolati e registrati in `metadata.csv` come diagnostiche: non entrano nella soluzione, servono a verificare che l'ipotesi di parete termicamente sottile regga.
 
-**Randomizzazione per run** (`configurazione_randomizzata`): ogni run perturba velocità, flusso di picco, sigma del punto, offset y, temperatura ambiente e rumore del sensore con estrazioni gaussiane/uniformi da un RNG con seed fisso, garantendo riproducibilità.
+`prepara_stato_termico` e `passo_termico` sono condivisi con `plot_animazione.py`, che riproduce la fisica di run_0001 per animare la superficie, e da lì con `plot_animazione_3d.py`: ogni modifica alla fisica va fatta lì, non duplicata.
 
-**Schema di output** (`dataset/run_XXXX.csv`): serie temporale con colonne `id_run, tempo_s, x_sorgente_m, offset_y_sorgente_m, flusso_termico_sorgente_W_m2, skin_depth_m, T_vera_lato_sensore_C, T_misurata_sensore_C, T_lato_caldo_C, T_ambiente_C, velocita_m_s, sigma_punto_m, flusso_picco_W_m2, materiale`. `metadata.csv` ha una riga per run con tutti i parametri e le temperature di picco.
+**Randomizzazione per run** (`configurazione_randomizzata`): ogni run perturba velocità, flusso di picco, sigma del punto, offset circonferenziale, temperatura ambiente e rumore del sensore con estrazioni gaussiane/uniformi da un RNG con seed fisso, garantendo riproducibilità.
+
+**Schema di output** (`dataset/run_XXXX.csv`): serie temporale con colonne `id_run, tempo_s, x_sorgente_m, offset_y_sorgente_m, flusso_termico_sorgente_W_m2, skin_depth_m, T_vera_lato_sensore_C, T_misurata_sensore_C, T_max_fascetta_C, T_media_fascetta_C, T_ambiente_C, velocita_m_s, sigma_punto_m, flusso_picco_W_m2, materiale`. `metadata.csv` ha una riga per run con tutti i parametri, le dimensioni della mesh, le diagnostiche e le temperature di picco.
 
 ## Convenzioni su `config.py`
 
@@ -64,8 +83,11 @@ Ogni parametro in [config.py](config.py) ha un commento che spiega solo cos'è (
 
 ## Vincoli progettuali chiave
 
-- Il modello è 2D nella sezione (x = lunghezza, z = spessore). La coordinata circonferenziale y non è risolta spazialmente — l'offset y del percorso delle sorgenti è collassato in un'attenuazione gaussiana del flusso, e la conduzione lungo y è un termine di scambio lineare verso la temperatura ambiente; il diametro è registrato solo come geometria del setup.
+- Il modello è una shell con una temperatura per nodo: lo spessore non è discretizzato e la temperatura della faccia esterna coincide con quella della faccia interna. L'ipotesi regge finché `numero_fourier_spessore` (α · σ/|v| / spessore²) resta molto maggiore di 1; sotto ~1 servirebbe una formulazione multistrato. Con i valori di default vale circa 5.4, e il gradiente misurato nello spessore col vecchio modello 2D era di 0.002 °C al picco.
+- La mesh cilindrica è strutturata e uniforme, quindi tutti gli elementi sono lo stesso rettangolo: le matrici elementari si calcolano una volta e si replicano. Se si introducessero mesh non uniformi o spessore variabile, `assemblaggio.py` dovrebbe ciclare sugli elementi.
+- La periodicità circonferenziale è strutturale, non un vincolo imposto: non esiste un nodo separato per theta = 2π, la connettività riavvolge l'ultimo elemento sul primo.
 - Le posizioni di inizio/fine corsa delle sorgenti (`x_inizio_m`, `x_fine_m`) sono distanze dal punto x del sensore lungo il verso di marcia; il segno di `velocita_m_s` determina il verso (negativo = -x).
 - La matrice implicita è costruita e fattorizzata una volta per run (proprietà del materiale costanti, nessun coefficiente dipendente dalla temperatura). Se si aggiungono proprietà dipendenti dalla temperatura, la matrice deve essere ricostruita e rifattorizzata ad ogni passo temporale.
+- La matrice è simmetrica definita positiva: `splu` usa `permc_spec="MMD_AT_PLUS_A"`, che dimezza circa il riempimento rispetto al COLAMD di default.
 - `simulate.py` cancella e ricrea l'intera cartella di output ad ogni esecuzione (`shutil.rmtree`).
 - Aggiungere un nuovo materiale richiede solo una nuova voce nel dizionario `MATERIALI` in [materials.py](materials.py); la chiave del materiale va poi impostata in `FASCETTA["materiale"]` in [config.py](config.py).

README.md

@@ -1,4 +1,4 @@
-# Simulatore Termico 2D — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso
+# Simulatore Termico a Elementi Shell — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso
 
 Questo progetto genera misurazioni CSV pseudo-realistiche della temperatura di una
 fascetta (anello cilindrico sottile) riscaldata da un gruppo di sorgenti a induzione
@@ -9,105 +9,191 @@ per l'addestramento e la validazione di modelli di stima/regressione termica.
 
 La fascetta è un anello cilindrico definito da tre dimensioni:
 
-- **diametro** (default 70 mm) — il diametro del cilindro;
-- **spessore** (default 0.12 mm) — lo spessore della parete;
+- **diametro medio** (default 70 mm) — il diametro della superficie media, su cui
+  giace la mesh;
+- **spessore** (default 0.18 mm) — lo spessore della parete;
 - **lunghezza** (default 100 mm) — l'estensione assiale lungo `x`.
 
-Il dominio simulato è la **sezione rettangolare lunghezza × spessore**. Il sistema di
-coordinate ha l'origine `(0, 0)` nel vertice in alto a sinistra della sezione:
+Il dominio simulato è l'**intera superficie cilindrica**, discretizzata con elementi
+shell quadrilateri. La superficie è parametrizzata da `(x, θ)`:
 
 ```text
-        sorgenti (induttori), in moto verso -x
-              ▼           ▼           ▼
-  (0,0) ─────────────────────────────────────► x
-    │   ┌───────────────────────────────────┐  z = 0        lato ESTERNO
-    │   │  sezione della fascetta           │               (flusso termico)
-    │   └───────────────────────────────────┘  z = spessore lato INTERNO
-    ▼                     ┆
-    z                     ┆  linea di vista
-                          ▲
-                     sensore IR (fisso, x = 50 mm,
-                     a 10 mm dalla parete interna)
+                 sorgenti (induttori), in moto verso -x
+                       ▼           ▼           ▼
+                  ╭────────────────────────────────╮
+     θ = 0  ─ ─ ─ ┤   ●  ← punto osservato dal     ├ ─ ─ ─  generatrice del sensore
+                  │      sensore IR (x = 50 mm)    │
+                  ╰────────────────────────────────╯
+                  x = 0                     x = lunghezza
+
+     superficie sviluppata in piano; il bordo superiore e quello inferiore
+     sono lo stesso punto della fascetta (θ = 0 ≡ θ = 2π)
 ```
 
 - **x** = direzione della lunghezza, da `0` a `lunghezza`. Le sorgenti viaggiano in
-  direzione `-x` sul lato esterno.
-- **z** = direzione dello spessore, da `0` (lato esterno, dove arriva il flusso
-  termico) a `spessore` (lato interno, osservato dal sensore).
-- **y** = coordinata circonferenziale (lungo la circonferenza π·diametro). Non è
-  risolta spazialmente: vedi sotto come viene trattata.
+  direzione `-x` sulla superficie esterna.
+- **θ** = coordinata circonferenziale, **periodica** e risolta spazialmente. La
+  generatrice `θ = 0` è quella del sensore; l'offset circonferenziale delle sorgenti
+  (`offset_y_percorso_m`, misurato lungo l'arco) è riferito ad essa.
+- Lo **spessore non è discretizzato**: la parete è termicamente sottile, quindi la
+  faccia esterna e quella interna hanno la stessa temperatura.
 
 Il **sensore** è un pirometro a infrarossi posto all'interno della fascetta, a una
 distanza configurabile dalla parete interna (default 10 mm). Essendo senza contatto,
 la distanza non influenza la misura: il sensore legge la temperatura della superficie
-interna nel punto `x` configurato (default 50 mm, al centro della lunghezza).
+interna nel punto `x` configurato (default 50 mm, al centro della lunghezza) sulla
+generatrice `θ = 0`.
 
 ## Modello fisico
 
-Non è una simulazione FEM elettromagnetica + termica completa: è un generatore
-pratico di dataset. La catena di approssimazioni è la seguente.
+Il campo elettromagnetico non è simulato: la parte termica è però una FEM vera. La
+catena di modellazione è la seguente.
 
-1. **Sorgenti gaussiane in moto** — ogni sorgente ha un'impronta gaussiana di raggio
-   `sigma_punto_m`. Un gruppo di `numero_sorgenti` sorgenti equidistanti
-   (`distanza_sorgenti_m`) si muove rigidamente a velocità costante. Il profilo di
-   flusso sul lato esterno è la somma dei contributi:
-   `q(x, t) = Σᵢ q_picco · efficienza · exp(-((x - xᵢ(t))² + Δy²) / (2σ²))`.
-   L'offset circonferenziale `Δy` tra il percorso delle sorgenti e il punto osservato
-   dal sensore non è risolto spazialmente: entra come attenuazione gaussiana del flusso.
+1. **Shell termica** — per una parete sottile con temperatura uniforme nello spessore
+   l'equazione del calore sulla superficie cilindrica è
 
-2. **Skin depth** — il riscaldamento a induzione è approssimato come riscaldamento
-   volumetrico che decade esponenzialmente con la profondità `z`:
-   `q_vol(x, z) = q(x) · exp(-z/δ) / (δ·(1 - exp(-spessore/δ)))`, normalizzato in modo
-   da conservare il flusso superficiale. La skin depth `δ = √(2ρₑ/(ωμ))` è calcolata
-   dalla resistività elettrica e dalla permeabilità del materiale alla frequenza di
-   induzione, oppure può essere imposta con `skin_depth_fissa_m`. Per la banda
-   stagnata a 20 kHz risulta ≈ 0.1 mm, confrontabile con lo spessore: la parete è
-   quasi isoterma attraverso lo spessore.
+   ```text
+   ρ cp s ∂T/∂t = k s ( ∂²T/∂x² + (1/R²) ∂²T/∂θ² ) + q'' − q_conv
+   ```
 
-3. **Diffusione 2D del calore** — l'equazione del calore è integrata nella sezione
-   `(x, z)` con volumi finiti ed Eulero implicito (incondizionatamente stabile).
+   La conduzione **assiale** e quella **circonferenziale** sono entrambe risolte.
+   Discretizzando con elementi quadrilateri a 4 nodi e una temperatura per nodo si
+   ottiene `C Ṫ + K T = f(t)`, con
 
-4. **Scambi con l'esterno** — la sezione scambia calore con l'ambiente su tutto il
-   contorno:
-   - convezione sul lato esterno (`h_esterno_W_m2K`), sul lato interno
-     (`h_interno_W_m2K`) e sui due bordi in x (`h_bordi_W_m2K`);
-   - **conduzione circonferenziale**: la sezione cede calore per conduzione lungo `y`
-     al resto della fascetta, assunto a temperatura ambiente. Il termine è
-     un'equazione di aletta ricavata sull'intero volume del cilindro: il calore
-     conduce lungo `y` attraverso l'intero spessore mentre le superfici esterna e
-     interna dell'intero cilindro perdono calore per convezione, dando
-     `q_y = -(h_esterno + h_interno)/spessore · (T - T_amb)` — nessun parametro di
-     conduzione `y` aggiuntivo da configurare.
+   ```text
+   Cₑ      = ∫ ρ cp s Nᵀ N dA          (matrice consistente)
+   Kₑ_cond = ∫ Bᵀ D B s dA             D = k·I  (materiale isotropo)
+   ```
+
+   La conduzione tra elementi adiacenti nasce dall'assemblaggio: due elementi
+   confinanti condividono i nodi del bordo comune, quindi condividono gli stessi
+   gradi di libertà termici. Non serve alcuna legge di scambio tra elementi dello
+   stesso materiale.
+
+2. **Convezione** — nelle shell la faccia esterna e quella interna sono **entrambe**
+   superfici fisiche esposte, quindi ogni elemento scambia su tutta la sua area:
+
+   ```text
+   Kₑ_conv = ∫ (h_est + h_int) Nᵀ N dA        fₑ_conv = ∫ (h_est + h_int) T∞ Nᵀ dA
+   ```
+
+   I due bordi assiali espongono solo lo spessore della lamiera:
+   `Kₑ_bordo = ∫ h_b s Nᵀ N dl`. Con `s = 0.18 mm` su `100 mm` di lunghezza il
+   contributo dei bordi vale ~0.2% dell'area ed è quindi marginale.
+
+3. **Sorgenti gaussiane in moto** — un gruppo di `numero_sorgenti` sorgenti
+   equidistanti (`distanza_sorgenti_m`) si muove rigidamente a velocità costante lungo
+   la generatrice individuata da `offset_y_percorso_m`. Il flusso sulla superficie
+   esterna è
+
+   ```text
+   q''(x, θ, t) = q_picco · efficienza · Σᵢ exp( −((x − xᵢ(t))² + Δy²) / (2σ²) )
+   Δy = R · wrap(θ − θ_s)
+   ```
+
+   La distanza circonferenziale è **riavvolta** su `2πR`: la sorgente scalda anche i
+   punti raggiunti dall'altro verso della circonferenza. Il carico nodale è
+   `fₑ_src = ∫ Nᵀ q'' dA`, integrato con la stessa quadratura di Gauss delle matrici:
+   la posizione circonferenziale dell'induttore è quindi simulata davvero, non
+   rappresentata da un coefficiente di attenuazione.
+
+4. **Skin depth** — `δ = √(2ρₑ/(ωμ))`, calcolata dalla resistività elettrica e dalla
+   permeabilità del materiale, oppure imposta con `skin_depth_fissa_m`. Poiché la
+   shell non discretizza lo spessore, la skin depth **non entra nella soluzione**: è
+   registrata in `metadata.csv` come diagnostica. Per la banda stagnata a 20 kHz vale
+   0.142 mm, confrontabile con lo spessore di 0.18 mm — cioè la deposizione è già
+   quasi uniforme attraverso la parete.
 
 5. **Temperatura iniziale** — il campo parte uniformemente alla temperatura ambiente
    del run (che è randomizzata, quindi varia run per run).
 
-6. **Sensore realistico** — la lettura aggiunge alla temperatura vera della superficie
-   interna: inerzia del primo ordine (`costante_tempo_s`), rumore gaussiano
-   (`rumore_std_C`) e quantizzazione (`quantizzazione_C`).
+6. **Sensore realistico** — la temperatura nel punto osservato è interpolata con le
+   funzioni di forma dell'elemento che lo contiene; la lettura vi aggiunge inerzia del
+   primo ordine (`costante_tempo_s`), rumore gaussiano (`rumore_std_C`) e
+   quantizzazione (`quantizzazione_C`).
+
+### Perché una sola temperatura nello spessore
+
+La parete è termicamente sottile con ampio margine. Il numero di Fourier che confronta
+il tempo di diffusione attraverso lo spessore con il tempo di transito della sorgente,
+
+```text
+Fo = α · (σ/|v|) / s²
+```
+
+vale **5.4** con i parametri di default (ed è riportato per ogni run in
+`metadata.csv`). Il modello 2D `(x, z)` precedente, che risolveva lo spessore con 15
+nodi, dava un salto esterno–interno di **0.002 °C al picco** su una sovratemperatura
+di ~200 °C, e 0.53 °C come massimo su tutto il run. Risolvere lo spessore non aggiunge
+informazione a questi parametri.
+
+Il criterio per sapere quando servirebbe una formulazione multistrato è `Fo ≈ 1`, cioè
+uno spessore critico `s* = √(α σ/|v|)`:
+
+| materiale                | Fo a 0.18 mm | spessore critico |
+|--------------------------|--------------|------------------|
+| alluminio                | 31.3         | 1.01 mm          |
+| ferro                    | 8.4          | 0.52 mm          |
+| banda stagnata / acciaio dolce | 5.4    | 0.42 mm          |
+| inox 304                 | 1.5          | 0.22 mm          |
 
 ## Metodo numerico
 
-- Griglia a volumi finiti `n_nodi_x × n_nodi_z` (default 100 × 15); le incognite sono
-  i centri cella.
-- Eulero implicito con passo `dt_interno_s` (default 0.2 ms), più fine del periodo di
-  campionamento CSV.
-- Tutti i termini (diffusione, convezione, conduzione circonferenziale) sono lineari e
-  costanti nel run: la matrice sparsa viene costruita e **fattorizzata LU una sola
-  volta per run** (`scipy.sparse.linalg.splu`); ogni passo temporale risolve solo i
-  sistemi triangolari. Un run da 30 s simulati richiede ~20 s di calcolo.
+- Mesh strutturata `n_elementi_x × n_elementi_theta` (default 100 × 110, cioè elementi
+  di 1.0 × 2.0 mm) per **11110 nodi**; la periodicità circonferenziale è strutturale,
+  non un vincolo imposto: non esiste un nodo separato per `θ = 2π`.
+- Quadratura di Gauss 2×2 per tutte le matrici elementari e per il vettore della
+  sorgente. La mesh è uniforme, quindi tutti gli elementi sono lo stesso rettangolo:
+  le matrici elementari si calcolano una volta e si replicano.
+- Eulero implicito (incondizionatamente stabile) con passo `dt_interno_s`
+  (default 1 ms), più fine del periodo di campionamento CSV.
+- Tutti i termini sono lineari e costanti nel run: `A = C/dt + K_cond + K_conv` viene
+  **fattorizzata LU una sola volta per run** (`scipy.sparse.linalg.splu`); ogni passo
+  risolve solo i sistemi triangolari. `A` è simmetrica definita positiva, quindi si usa
+  l'ordinamento `MMD_AT_PLUS_A`, che dimezza circa il riempimento rispetto al COLAMD di
+  default. Un run da 30 s simulati richiede ~65 s di calcolo; i run sono generati in
+  parallelo su più processi.
 - Se in futuro si introducessero proprietà dipendenti dalla temperatura, la matrice
   andrebbe ricostruita e rifattorizzata a ogni passo.
 
+### Validazione
+
+`python valida_shell.py` verifica il solutore su casi con soluzione nota. Va rieseguito
+dopo ogni modifica alla fisica o all'assemblaggio. Risultati con la mesh di default:
+
+| Verifica                                        | Esito                                    |
+|-------------------------------------------------|------------------------------------------|
+| capacità totale = `ρ cp s A`                    | esatta a precisione macchina             |
+| somma righe di `K_cond` = 0 (T costante → flusso nullo) | 4·10⁻¹⁸                          |
+| simmetria e semidefinita positività di `K_cond` | verificate                               |
+| raffreddamento convettivo uniforme vs `exp(−t/τ)` | errore 0.011 °C, pari a `dt/(2τ)`      |
+| conservazione dell'energia (adiabatico)         | errore relativo 4·10⁻¹⁴                  |
+| `Σ f_src` vs potenza analitica `2πσ² q_max`     | rapporto 1.0000; **0.5000** con sorgente centrata sul bordo |
+| modi sinusoidali circonferenziali `n = 1, 2, 4` vs `exp(−α n²t/R²)` | errore 0.04% – 0.17%   |
+| convergenza al raffinamento (modo `n = 4`)      | 2.75% → 1.07% → 0.17% → 0.05% da 30 a 220 elementi |
+| simmetria della diffusione attorno a `θ = 0`    | 2·10⁻¹⁴                                  |
+
+L'ultima riga è la prova che la periodicità funziona davvero: un impulso su `θ = 0`
+diffonde in modo identico nei due versi solo se la connettività riavvolge la
+circonferenza.
+
 ## File
 
 ```text
-config.py           tutti i parametri di simulazione
-materials.py        proprietà termofisiche ed elettriche dei materiali
-simulate.py         motore fisico + generazione dei CSV
-plot_csv.py         grafici rapidi (temperature e flusso) del primo run
-plot_animazione.py  animazione della sezione: campo T(x,z), sorgenti, sensore
-dataset/            output generato da simulate.py (ricreato a ogni esecuzione)
+config.py               tutti i parametri di simulazione
+materials.py            proprietà termofisiche ed elettriche dei materiali
+mesh.py                 nodi, connettività Q4, periodicità circonferenziale
+elementi_shell.py       funzioni di forma, quadratura, matrici del singolo elemento
+assemblaggio.py         assemblaggio delle matrici globali e del vettore sorgente
+sorgente.py             moto degli induttori e distribuzione gaussiana del flusso
+solutore.py             integrazione implicita e fattorizzazione sparsa
+sensore.py              interpolazione, inerzia, rumore, quantizzazione
+simulate.py             orchestrazione dei run + generazione dei CSV
+valida_shell.py         verifiche del solutore su casi con soluzione nota
+plot_csv.py             grafici rapidi (temperature e flusso) del primo run
+plot_animazione.py      animazione della superficie sviluppata T(x, θ)
+plot_animazione_3d.py   animazione 3D isometrica del cilindro
+dataset/                output generato da simulate.py (ricreato a ogni esecuzione)
 ```
 
 ## Installazione
@@ -127,19 +213,50 @@ python simulate.py
 # grafici statici del primo run (temperatura e flusso nel tempo)
 python plot_csv.py
 
-# animazione della sezione durante il passaggio delle sorgenti
+# animazione della superficie sviluppata durante il passaggio delle sorgenti
 python plot_animazione.py
+
+# animazione 3D isometrica del cilindro
+python plot_animazione_3d.py
 ```
 
 Gli script di visualizzazione aprono finestre interattive (backend Qt); se il backend
 non è interattivo (es. sessione senza display) salvano automaticamente PNG/GIF in
 `dataset/`.
 
-L'animazione riproduce esattamente la fisica di `run_0001` (stesso seed) e mostra tre
-pannelli allineati: il profilo di flusso `q(x)` con le sorgenti in transito, il campo
-di temperatura nella sezione con il sensore IR, e la temperatura nel punto osservato
-(vera e con inerzia del sensore). Finestra temporale e cadenza dei fotogrammi si
-regolano con le costanti in testa a `plot_animazione.py`.
+Entrambe le animazioni riproducono esattamente la fisica di `run_0001` (stesso seed).
+Quella 2D mostra tre pannelli allineati: il profilo di flusso `q(x)` con le sorgenti in
+transito, il campo di temperatura sulla superficie sviluppata in piano (asse verticale
+= arco circonferenziale, con il sensore a quota 0), e la temperatura nel punto
+osservato (vera e con inerzia del sensore). Quella 3D disegna il campo calcolato sulla
+superficie del cilindro — non è più una ricostruzione, il modello risolve davvero la
+distribuzione circonferenziale. Finestra temporale e cadenza dei fotogrammi si regolano
+con le costanti in testa a `plot_animazione.py`.
+
+La mappa sviluppata usa `pcolormesh` con shading Gouraud, che interpola linearmente tra
+i valori nodali: è esattamente ciò che fanno le funzioni di forma bilineari del Q4,
+quindi si guarda la soluzione dell'elemento finito e non un'interpolazione grafica.
+
+### Vedere gli elementi shell
+
+Entrambe le animazioni hanno una costante `MOSTRA_MESH` che disegna i bordi degli
+elementi. È `False` di default perché con la mesh di produzione gli elementi sono
+1 × 2 mm, più fitti dei pixel disponibili, e il reticolo coprirebbe il campo. Per
+ispezionare la mesh servono due modifiche insieme:
+
+```python
+# in plot_animazione.py e/o plot_animazione_3d.py
+MOSTRA_MESH = True
+
+# in config.py, dentro FASCETTA
+"n_elementi_x": 25,
+"n_elementi_theta": 28,
+```
+
+Così si vedono i singoli quadrilateri, in piano e avvolti sul cilindro. Attenzione: con
+28 elementi sulla circonferenza l'arco vale 7.9 mm contro un `sigma_punto_m` di 12 mm,
+cioè meno di due elementi per sigma — la mesh grossolana serve a **guardare la
+discretizzazione**, non a produrre risultati (vedi la tabella di convergenza sopra).
 
 ## Configurazione
 
@@ -148,8 +265,8 @@ Tutto si modifica in `config.py`. I dizionari principali:
 | Dizionario       | Contenuto                                                                 |
 |------------------|---------------------------------------------------------------------------|
 | `SIMULAZIONE`    | numero di run, durata, campionamento CSV, passo interno, seed, cartella   |
-| `FASCETTA`       | diametro, lunghezza, spessore, griglia, conduzione circonferenziale, materiale |
-| `ARIA`           | temperatura ambiente e coefficienti di convezione dei quattro lati        |
+| `FASCETTA`       | diametro medio, lunghezza, spessore, densità della mesh, materiale        |
+| `ARIA`           | temperatura ambiente e coefficienti di convezione di facce e bordi        |
 | `SORGENTE`       | corsa, velocità (il segno dà il verso), gruppo di sorgenti, gaussiana, flusso, frequenza |
 | `SENSORE`        | posizione (x e distanza dalla parete), inerzia, rumore, quantizzazione    |
 | `RANDOMIZZAZIONE`| entità delle perturbazioni per run                                        |
@@ -191,26 +308,30 @@ Punti da conoscere:
 | `skin_depth_m`                 | skin depth usata (costante nel run)                                |
 | `T_vera_lato_sensore_C`        | temperatura vera della superficie interna nel punto del sensore    |
 | `T_misurata_sensore_C`         | lettura del sensore (inerzia + rumore + quantizzazione)            |
-| `T_lato_caldo_C`               | temperatura della superficie esterna nello stesso punto x          |
+| `T_max_fascetta_C`             | temperatura massima su tutta la superficie in quell'istante        |
+| `T_media_fascetta_C`           | temperatura media su tutta la superficie in quell'istante          |
 | `T_ambiente_C`                 | temperatura ambiente del run                                       |
 | `velocita_m_s`, `sigma_punto_m`, `flusso_picco_W_m2` | parametri randomizzati del run              |
 | `materiale`                    | chiave del materiale                                               |
 
 ### `dataset/metadata.csv` — una riga per run
 
-Contiene tutti i parametri effettivi del run (geometria, griglia, coefficienti di
-scambio, parametri delle sorgenti e del sensore, valori randomizzati) e le temperature
-di picco vera e misurata: utile come ground truth e per filtrare i run.
+Contiene tutti i parametri effettivi del run (geometria, dimensioni della mesh,
+coefficienti di scambio, parametri delle sorgenti e del sensore, valori randomizzati),
+le diagnostiche `skin_depth_m` e `numero_fourier_spessore`, e le temperature di picco:
+utile come ground truth e per filtrare i run.
 
 ## Limitazioni
 
 1. Il campo elettromagnetico non è simulato: l'accoppiamento induttivo è ridotto a
-   impronta gaussiana × efficienza × decadimento esponenziale in z.
-2. La coordinata circonferenziale y non è risolta: offset del percorso e conduzione
-   verso il resto della fascetta sono modelli collassati (attenuazione gaussiana e
-   scambio lineare verso T ambiente).
+   un'impronta gaussiana moltiplicata per un'efficienza costante.
+2. Lo spessore non è discretizzato: il modello non può rappresentare un gradiente tra
+   faccia esterna e faccia interna. Con i parametri di default il gradiente reale è di
+   pochi millesimi di grado (vedi `numero_fourier_spessore`), ma con pareti più spesse
+   di ~0.4 mm, materiali a bassa diffusività o sorgenti molto più veloci servirebbe una
+   shell multistrato.
 3. Le proprietà dei materiali sono costanti con la temperatura; per gli acciai
    ferromagnetici la skin depth reale varia fortemente con temperatura e campo
    (punto di Curie non modellato).
-4. Irraggiamento non modellato: a ~220 °C le perdite radiative non sono del tutto
+4. Irraggiamento non modellato: a ~210 °C le perdite radiative non sono del tutto
    trascurabili rispetto alla convezione.

assemblaggio.py

@@ -0,0 +1,151 @@
+# Assemblaggio delle matrici globali della shell termica.
+#
+# Tutti gli elementi della mesh cilindrica sono lo stesso rettangolo, quindi
+# ogni matrice elementare si calcola una volta e si replica su tutti gli
+# elementi. L'assemblaggio è la somma dei contributi sui gradi di libertà
+# condivisi: due elementi confinanti condividono i due nodi del bordo comune,
+# ed è questo che genera la conduzione tra elementi adiacenti. Non serve
+# nessun accoppiamento aggiuntivo tra elementi dello stesso materiale.
+
+import numpy as np
+import scipy.sparse as sp
+
+import elementi_shell as el
+
+
+def _assembla_matrice(matrice_elemento: np.ndarray, connettivita: np.ndarray, n_nodi: int):
+    # Replica una matrice elementare costante su tutti gli elementi e somma i
+    # contributi sui gradi di libertà condivisi.
+    n_elementi, n_locali = connettivita.shape
+    righe = np.repeat(connettivita, n_locali, axis=1).ravel()
+    colonne = np.tile(connettivita, (1, n_locali)).ravel()
+    dati = np.tile(matrice_elemento.ravel(), n_elementi)
+    return sp.coo_matrix(
+        (dati, (righe, colonne)), shape=(n_nodi, n_nodi)
+    ).tocsr()
+
+
+def _assembla_vettore(vettore_elemento: np.ndarray, connettivita: np.ndarray, n_nodi: int):
+    dati = np.tile(vettore_elemento, connettivita.shape[0])
+    return np.bincount(connettivita.ravel(), weights=dati, minlength=n_nodi)
+
+
+def assembla_capacita(mesh: dict, materiale: dict, spessore_m: float):
+    # C = integrale di rho * cp * spessore * N^T N dA.
+    #
+    # Matrice consistente: con integrazione temporale implicita conserva
+    # meglio l'energia rispetto alla versione concentrata sulla diagonale.
+    capacita_superficiale = (
+        materiale["densita_kg_m3"] * materiale["calore_specifico_J_kgK"] * spessore_m
+    )
+    M_e = el.matrice_massa_elemento(mesh["lato_x_m"], mesh["lato_arco_m"])
+    return _assembla_matrice(
+        capacita_superficiale * M_e, mesh["connettivita"], mesh["n_nodi"]
+    )
+
+
+def assembla_conduzione(mesh: dict, materiale: dict, spessore_m: float):
+    # K_cond = integrale di B^T D B * spessore dA, con D isotropa.
+    K_e = el.matrice_conduzione_elemento(
+        mesh["lato_x_m"],
+        mesh["lato_arco_m"],
+        materiale["conducibilita_termica_W_mK"],
+        spessore_m,
+    )
+    return _assembla_matrice(K_e, mesh["connettivita"], mesh["n_nodi"])
+
+
+def assembla_convezione(mesh: dict, aria: dict, spessore_m: float):
+    # Convezione sulle facce e sui bordi assiali.
+    #
+    # Nelle shell la faccia esterna e quella interna sono entrambe superfici
+    # fisiche esposte all'aria, quindi ogni elemento scambia su tutta la sua
+    # area con il coefficiente combinato h_esterno + h_interno. I bordi
+    # assiali (x = 0 e x = lunghezza) espongono invece solo lo spessore della
+    # lamiera: il loro contributo è proporzionale a h_bordi * spessore ed è
+    # marginale rispetto a quello delle facce.
+    T_ambiente = aria["temperatura_ambiente_C"]
+    h_facce = aria["h_esterno_W_m2K"] + aria["h_interno_W_m2K"]
+
+    M_e = el.matrice_massa_elemento(mesh["lato_x_m"], mesh["lato_arco_m"])
+    f_e = el.vettore_carico_uniforme_elemento(mesh["lato_x_m"], mesh["lato_arco_m"])
+
+    K_conv = _assembla_matrice(h_facce * M_e, mesh["connettivita"], mesh["n_nodi"])
+    f_ambiente = _assembla_vettore(
+        h_facce * T_ambiente * f_e, mesh["connettivita"], mesh["n_nodi"]
+    )
+
+    h_bordi_efficace = aria["h_bordi_W_m2K"] * spessore_m
+    M_bordo = el.matrice_massa_segmento(mesh["lato_arco_m"])
+    f_bordo = el.vettore_carico_segmento(mesh["lato_arco_m"])
+
+    K_conv = K_conv + _assembla_matrice(
+        h_bordi_efficace * M_bordo, mesh["segmenti_bordo"], mesh["n_nodi"]
+    )
+    f_ambiente = f_ambiente + _assembla_vettore(
+        h_bordi_efficace * T_ambiente * f_bordo,
+        mesh["segmenti_bordo"],
+        mesh["n_nodi"],
+    )
+
+    return K_conv, f_ambiente
+
+
+def prepara_assemblatore_sorgente(mesh: dict) -> dict:
+    # Precalcola quanto serve per assemblare il vettore della sorgente a ogni
+    # passo temporale: le funzioni di forma nei punti di Gauss pesate per il
+    # peso di quadratura e lo jacobiano, e le coordinate dei punti di Gauss
+    # separate per direzione.
+    det_j = mesh["lato_x_m"] * mesh["lato_arco_m"] / 4.0
+    pesi = el.pesi_ai_punti_gauss()
+    N_gauss = el.matrice_forme_ai_punti_gauss()
+    offset_x, offset_arco = el.coordinate_locali_punti_gauss(
+        mesh["lato_x_m"], mesh["lato_arco_m"]
+    )
+
+    # Coordinate assolute dei punti di Gauss, per direzione: (n_elementi_x, 2)
+    # e (n_elementi_theta, 2).
+    x_gauss = mesh["x_nodi_m"][: mesh["n_elementi_x"], None] + offset_x[None, :]
+    arco_gauss = mesh["arco_nodi_m"][:, None] + offset_arco[None, :]
+
+    return {
+        "peso_forme": (pesi * det_j)[:, None] * N_gauss,
+        "connettivita": mesh["connettivita"],
+        "n_nodi": mesh["n_nodi"],
+        "n_elementi_x": mesh["n_elementi_x"],
+        "n_elementi_theta": mesh["n_elementi_theta"],
+        "x_gauss_m": x_gauss,
+        "arco_gauss_m": arco_gauss,
+    }
+
+
+def assembla_sorgente(
+    assemblatore: dict, flusso_x_W_m2: np.ndarray, fattore_arco: np.ndarray
+) -> np.ndarray:
+    # f_src = integrale di N^T q'' dA, con q'' fattorizzato come
+    # q''(x, arco) = flusso_x(x) * fattore_arco(arco).
+    #
+    # La fattorizzazione non è un'approssimazione: la gaussiana della
+    # sorgente è separabile in x e in arco, e tutte le sorgenti del gruppo
+    # condividono la stessa posizione circonferenziale, quindi il fattore
+    # circonferenziale è comune e costante nel tempo.
+    #
+    # flusso_x_W_m2 ha forma (n_elementi_x, 2) e fattore_arco
+    # (n_elementi_theta, 2): sono i valori nei punti di Gauss di ciascuna
+    # direzione.
+    #
+    # Per gran parte di un run le sorgenti sono fuori dalla fascetta e il
+    # fattore assiale è identicamente nullo: in quel caso il vettore è nullo
+    # e non serve percorrere la mesh.
+    if not flusso_x_W_m2.any():
+        return np.zeros(assemblatore["n_nodi"])
+
+    q_gauss = (
+        flusso_x_W_m2[:, None, :, None] * fattore_arco[None, :, None, :]
+    ).reshape(-1, 4)
+    contributi = q_gauss @ assemblatore["peso_forme"]
+    return np.bincount(
+        assemblatore["connettivita"].ravel(),
+        weights=contributi.ravel(),
+        minlength=assemblatore["n_nodi"],
+    )

config.py

@@ -1,26 +1,25 @@
-# Configurazione per il simulatore termico 2D della sezione di una fascetta.
+# Configurazione per il simulatore termico a elementi shell di una fascetta.
 #
 # Geometria e modello fisico:
-# - La fascetta è un anello cilindrico con diametro "diametro_mm", spessore
-#   "spessore_mm" e lunghezza "lunghezza_mm".
-# - Il dominio simulato è la sezione rettangolare lunghezza × spessore.
-# - Sistema di coordinate: origine (0, 0) nel vertice in alto a sinistra
-#   della sezione. x = direzione della lunghezza (da 0 a lunghezza),
-#   z = direzione dello spessore (0 = lato esterno, dove agiscono le
-#   sorgenti; spessore = lato interno, osservato dal sensore).
-# - y è la coordinata circonferenziale: non è risolta spazialmente, l'offset
-#   y del percorso delle sorgenti è collassato in un'attenuazione gaussiana
-#   del flusso.
-# - Le sorgenti a induzione si muovono in direzione -x sul lato esterno.
-# - Il riscaldamento a induzione è approssimato come riscaldamento volumetrico
-#   che decade esponenzialmente con la profondità z secondo la skin depth.
-# - La sezione scambia per convezione con l'aria su tutti e quattro i lati
-#   (esterno, interno e i due bordi in x). Scambia inoltre per conduzione
-#   lungo y con il resto della fascetta, assunto a temperatura ambiente:
-#   il calore conduce attraverso l'intero volume dello spessore mentre le
-#   superfici esterna e interna dell'intero cilindro perdono calore per
-#   convezione (equazione dell'aletta), derivato da h_esterno, h_interno e
-#   spessore_mm senza parametri di conduzione y aggiuntivi.
+# - La fascetta è un anello cilindrico con diametro medio "diametro_medio_mm",
+#   spessore "spessore_mm" e lunghezza "lunghezza_mm".
+# - Il dominio simulato è l'intera superficie cilindrica, discretizzata con
+#   elementi shell quadrilateri a 4 nodi e una temperatura per nodo.
+# - Sistema di coordinate: x = direzione della lunghezza (da 0 a lunghezza),
+#   theta = coordinata circonferenziale, periodica (theta = 0 coincide con
+#   theta = 2*pi). Il sensore sta per convenzione su theta = 0.
+# - Lo spessore non è discretizzato: la parete è termicamente sottile, quindi
+#   la faccia esterna e quella interna hanno la stessa temperatura. La
+#   validità dell'ipotesi è riportata a ogni run come numero di Fourier sullo
+#   spessore in metadata.csv.
+# - Le sorgenti a induzione si muovono in direzione -x sulla superficie
+#   esterna, lungo la generatrice individuata da "offset_y_percorso_m".
+# - La conduzione assiale e quella circonferenziale sono entrambe risolte e
+#   nascono dall'assemblaggio: elementi adiacenti condividono i nodi del
+#   bordo comune.
+# - Faccia esterna e faccia interna scambiano entrambe per convezione con
+#   l'aria su tutta l'area; i due bordi assiali scambiano attraverso lo
+#   spessore della lamiera.
 # - La temperatura iniziale della fascetta è pari alla temperatura ambiente.
 # - Il sensore è un pirometro a infrarossi posto all'interno della fascetta,
 #   a distanza "distanza_parete_mm" dalla parete interna: misura senza
@@ -46,7 +45,7 @@ SIMULAZIONE = {
 
     # Passo di integrazione numerica interna.
     # Può essere inferiore al periodo di campionamento CSV.
-    "dt_interno_s": 0.0002,
+    "dt_interno_s": 0.001,
 
     # Seed per la riproducibilità.
     "seed": 42,
@@ -60,8 +59,9 @@ SIMULAZIONE = {
 }
 
 FASCETTA = {
-    # Diametro della fascetta [mm].
-    "diametro_mm": 70.0,
+    # Diametro medio della fascetta [mm], cioè il diametro della superficie
+    # su cui giace la mesh shell.
+    "diametro_medio_mm": 70.0,
 
     # Lunghezza della fascetta lungo x [mm].
     "lunghezza_mm": 100.0,
@@ -69,12 +69,12 @@ FASCETTA = {
     # Spessore della parete [mm].
     "spessore_mm": 0.18,
 
-    # Numero di celle del volume finito lungo x (lunghezza).
-    "n_nodi_x": 100,
+    # Numero di elementi shell lungo x (lunghezza).
+    "n_elementi_x": 100,
 
-    # Numero di celle del volume finito lungo z (spessore).
-    # Più nodi = maggiore risoluzione spaziale, simulazione più lenta.
-    "n_nodi_z": 15,
+    # Numero di elementi shell lungo la circonferenza.
+    # Più elementi = maggiore risoluzione spaziale, simulazione più lenta.
+    "n_elementi_theta": 110,
 
     # Deve corrispondere a una chiave in materials.py.
     "materiale": "banda_stagnata",
@@ -84,13 +84,14 @@ ARIA = {
     # Temperatura dell'aria ambiente.
     "temperatura_ambiente_C": 25.0,
 
-    # Coefficiente di convezione sul lato esterno (z = 0, lato sorgenti).
+    # Coefficiente di convezione sulla faccia esterna (lato sorgenti).
     "h_esterno_W_m2K": 12.0,
 
-    # Coefficiente di convezione sul lato interno (z = spessore, lato sensore).
+    # Coefficiente di convezione sulla faccia interna (lato sensore).
     "h_interno_W_m2K": 8.0,
 
-    # Coefficiente di convezione sui bordi laterali (x = 0 e x = lunghezza).
+    # Coefficiente di convezione sui bordi assiali (x = 0 e x = lunghezza),
+    # che espongono all'aria solo lo spessore della lamiera.
     "h_bordi_W_m2K": 10.0,
 }
 
@@ -103,9 +104,10 @@ SORGENTE = {
     # gruppo più lontana dal sensore (quella che lo supera per ultima).
     "x_fine_m": 5.0,
 
-    # Offset circonferenziale (y) tra il percorso delle sorgenti e il punto
-    # osservato dal sensore. Se diverso da zero, la sorgente passa a lato,
-    # riducendo il picco.
+    # Offset circonferenziale, misurato lungo l'arco, tra il percorso delle
+    # sorgenti e la generatrice osservata dal sensore. Se diverso da zero, la
+    # sorgente passa a lato e il calore deve diffondere lungo la
+    # circonferenza per raggiungere il sensore.
     "offset_y_percorso_m": 0.0,
 
     # Velocità delle sorgenti lungo x. Il segno indica il verso di marcia.
@@ -128,6 +130,9 @@ SORGENTE = {
     "efficienza_riscaldamento": 0.35,
 
     # Frequenza di induzione usata per stimare la skin depth se skin_depth_fissa_m è None.
+    # La shell non discretizza lo spessore, quindi la skin depth non entra
+    # nella soluzione: è calcolata e registrata come diagnostica, per
+    # verificare che resti confrontabile con lo spessore della lamiera.
     "frequenza_hz": 20000.0,
 
     # Override della skin depth. Usare None per calcolarla dalle proprietà elettriche del materiale.
@@ -140,9 +145,11 @@ SORGENTE = {
 
 SENSORE = {
     # Coordinata x del punto della superficie interna osservato dal sensore [mm].
+    # La posizione circonferenziale è theta = 0 per convenzione: è l'offset
+    # delle sorgenti a essere misurato rispetto al sensore.
     "x_mm": 50.0,
 
-    # Distanza del sensore dalla parete interna lungo z [mm].
+    # Distanza del sensore dalla parete interna [mm].
     # Il sensore è a infrarossi: la distanza non influenza la misura,
     # è registrata solo come geometria del setup.
     "distanza_parete_mm": 10.0,

elementi_shell.py

@@ -0,0 +1,135 @@
+# Elemento shell termico quadrilatero a 4 nodi (Q4) con una temperatura per
+# nodo, cioè temperatura uniforme nello spessore.
+#
+# L'elemento vive nel piano (x, arco) della superficie sviluppata. Poiché la
+# mesh cilindrica è strutturata e uniforme, tutti gli elementi sono lo stesso
+# rettangolo lato_x × lato_arco: le matrici elementari si calcolano una volta
+# sola e si replicano in fase di assemblaggio.
+#
+# Le funzioni di forma bilineari sono
+#
+#     N_a(xi, eta) = (1 + xi_a * xi) * (1 + eta_a * eta) / 4
+#
+# con (xi_a, eta_a) i vertici dell'elemento di riferimento. La mappa verso
+# l'elemento fisico è affine, quindi lo jacobiano è costante:
+# det(J) = lato_x * lato_arco / 4.
+#
+# Integrazione con quadratura di Gauss 2x2, esatta per i polinomi che
+# compaiono nelle matrici di massa e conduzione dell'elemento Q4.
+
+import numpy as np
+
+
+# Vertici dell'elemento di riferimento, nello stesso ordine dei nodi locali
+# usati dalla connettività in mesh.py.
+NODI_RIFERIMENTO = np.array([(-1.0, -1.0), (1.0, -1.0), (1.0, 1.0), (-1.0, 1.0)])
+
+# Quadratura di Gauss a 2 punti per direzione.
+PUNTI_GAUSS_1D = np.array([-1.0 / np.sqrt(3.0), 1.0 / np.sqrt(3.0)])
+PESI_GAUSS_1D = np.array([1.0, 1.0])
+
+
+def funzioni_forma(xi: float, eta: float) -> np.ndarray:
+    xi_a, eta_a = NODI_RIFERIMENTO[:, 0], NODI_RIFERIMENTO[:, 1]
+    return 0.25 * (1.0 + xi_a * xi) * (1.0 + eta_a * eta)
+
+
+def derivate_funzioni_forma(xi: float, eta: float) -> tuple[np.ndarray, np.ndarray]:
+    xi_a, eta_a = NODI_RIFERIMENTO[:, 0], NODI_RIFERIMENTO[:, 1]
+    dN_dxi = 0.25 * xi_a * (1.0 + eta_a * eta)
+    dN_deta = 0.25 * eta_a * (1.0 + xi_a * xi)
+    return dN_dxi, dN_deta
+
+
+def _punti_quadratura_2d():
+    # Punti (xi, eta) e pesi della quadratura 2x2, nell'ordine
+    # g = p * 2 + q con xi = punto p (direzione x) ed eta = punto q
+    # (direzione circonferenziale). Lo stesso ordine è atteso dai valori del
+    # flusso passati a vettore_sorgente_elementi().
+    for p, xi in enumerate(PUNTI_GAUSS_1D):
+        for q, eta in enumerate(PUNTI_GAUSS_1D):
+            yield xi, eta, PESI_GAUSS_1D[p] * PESI_GAUSS_1D[q]
+
+
+def matrice_forme_ai_punti_gauss() -> np.ndarray:
+    # Matrice (4 punti di Gauss, 4 nodi) con le funzioni di forma valutate
+    # nei punti di quadratura.
+    return np.array([funzioni_forma(xi, eta) for xi, eta, _ in _punti_quadratura_2d()])
+
+
+def pesi_ai_punti_gauss() -> np.ndarray:
+    # Pesi di quadratura nei 4 punti, nello stesso ordine.
+    return np.array([peso for _, _, peso in _punti_quadratura_2d()])
+
+
+def coordinate_locali_punti_gauss(lato_x_m: float, lato_arco_m: float):
+    # Offset dei punti di Gauss rispetto al vertice di riferimento
+    # dell'elemento, separati per direzione: 2 valori lungo x e 2 lungo
+    # l'arco. La separazione è possibile perché l'elemento è un rettangolo.
+    offset_x = lato_x_m * (1.0 + PUNTI_GAUSS_1D) / 2.0
+    offset_arco = lato_arco_m * (1.0 + PUNTI_GAUSS_1D) / 2.0
+    return offset_x, offset_arco
+
+
+def matrice_massa_elemento(lato_x_m: float, lato_arco_m: float) -> np.ndarray:
+    # Integrale di N^T N sull'area dell'elemento. Matrice consistente (non
+    # concentrata): è la base sia della capacità termica sia della convezione
+    # sulle facce.
+    det_j = lato_x_m * lato_arco_m / 4.0
+    M = np.zeros((4, 4))
+    for xi, eta, peso in _punti_quadratura_2d():
+        N = funzioni_forma(xi, eta)
+        M += peso * det_j * np.outer(N, N)
+    return M
+
+
+def matrice_conduzione_elemento(
+    lato_x_m: float, lato_arco_m: float, conducibilita_W_mK: float, spessore_m: float
+) -> np.ndarray:
+    # Integrale di B^T D B * spessore sull'area, con D = k * identità
+    # (materiale isotropo). Rappresenta la conduzione tangenziale, cioè sia
+    # quella assiale sia quella circonferenziale: sul cilindro sviluppato le
+    # due direzioni sono ortogonali e la derivata circonferenziale rispetto
+    # all'arco equivale a (1/R) d/dtheta.
+    det_j = lato_x_m * lato_arco_m / 4.0
+    K = np.zeros((4, 4))
+    for xi, eta, peso in _punti_quadratura_2d():
+        dN_dxi, dN_deta = derivate_funzioni_forma(xi, eta)
+        # Jacobiano diagonale e costante: la derivata fisica è quella
+        # naturale riscalata dal semilato dell'elemento.
+        B = np.stack([dN_dxi * (2.0 / lato_x_m), dN_deta * (2.0 / lato_arco_m)])
+        K += peso * det_j * conducibilita_W_mK * spessore_m * (B.T @ B)
+    return K
+
+
+def vettore_carico_uniforme_elemento(lato_x_m: float, lato_arco_m: float) -> np.ndarray:
+    # Integrale di N^T sull'area: distribuisce ai nodi un carico superficiale
+    # uniforme unitario.
+    det_j = lato_x_m * lato_arco_m / 4.0
+    f = np.zeros(4)
+    for xi, eta, peso in _punti_quadratura_2d():
+        f += peso * det_j * funzioni_forma(xi, eta)
+    return f
+
+
+def _funzioni_forma_segmento(xi: float) -> np.ndarray:
+    return np.array([0.5 * (1.0 - xi), 0.5 * (1.0 + xi)])
+
+
+def matrice_massa_segmento(lunghezza_m: float) -> np.ndarray:
+    # Integrale di N^T N su un segmento a 2 nodi, per i bordi assiali.
+    det_j = lunghezza_m / 2.0
+    M = np.zeros((2, 2))
+    for xi, peso in zip(PUNTI_GAUSS_1D, PESI_GAUSS_1D):
+        N = _funzioni_forma_segmento(xi)
+        M += peso * det_j * np.outer(N, N)
+    return M
+
+
+def vettore_carico_segmento(lunghezza_m: float) -> np.ndarray:
+    # Integrale di N^T su un segmento a 2 nodi.
+    det_j = lunghezza_m / 2.0
+    f = np.zeros(2)
+    for xi, peso in zip(PUNTI_GAUSS_1D, PESI_GAUSS_1D):
+        f += peso * det_j * _funzioni_forma_segmento(xi)
+    return f

mesh.py

@@ -0,0 +1,119 @@
+# Generazione della mesh a elementi shell quadrilateri sulla superficie
+# cilindrica della fascetta.
+#
+# La superficie è parametrizzata da (x, theta): x è la coordinata assiale,
+# theta quella circonferenziale. Al posto di theta si usa quasi ovunque
+# l'ascissa curvilinea arco = raggio_medio * theta, così la shell diventa un
+# rettangolo lunghezza × circonferenza e gli elementi sono tutti identici.
+#
+# I nodi formano una griglia (n_elementi_x + 1) × n_elementi_theta: in
+# direzione circonferenziale non c'è un nodo finale distinto perché theta = 0
+# e theta = 2*pi sono lo stesso nodo. La periodicità è quindi strutturale,
+# nasce dalla connettività e non da vincoli imposti a posteriori.
+#
+# Ordine dei nodi locali di ogni elemento (antiorario nel piano x-arco):
+#
+#     4 ---- 3      arco
+#     |      |       ^
+#     |      |       |
+#     1 ---- 2       +---> x
+
+import numpy as np
+
+
+def costruisci_mesh_cilindrica(
+    lunghezza_m: float,
+    raggio_medio_m: float,
+    n_elementi_x: int,
+    n_elementi_theta: int,
+) -> dict:
+    # Costruisce nodi e connettività della shell cilindrica.
+    n_nodi_x = n_elementi_x + 1
+    n_nodi = n_nodi_x * n_elementi_theta
+    n_elementi = n_elementi_x * n_elementi_theta
+
+    lato_x_m = lunghezza_m / n_elementi_x
+    circonferenza_m = 2.0 * np.pi * raggio_medio_m
+    lato_arco_m = circonferenza_m / n_elementi_theta
+
+    x_nodi_m = np.arange(n_nodi_x) * lato_x_m
+    theta_nodi_rad = np.arange(n_elementi_theta) * (2.0 * np.pi / n_elementi_theta)
+    arco_nodi_m = theta_nodi_rad * raggio_medio_m
+
+    # Connettività: elemento (i, j) collega i nodi (i, j), (i+1, j),
+    # (i+1, j+1), (i, j+1), con j+1 riavvolto sulla circonferenza.
+    i_elem = np.repeat(np.arange(n_elementi_x), n_elementi_theta)
+    j_elem = np.tile(np.arange(n_elementi_theta), n_elementi_x)
+    j_succ = (j_elem + 1) % n_elementi_theta
+
+    def indice(i, j):
+        return i * n_elementi_theta + j
+
+    connettivita = np.stack(
+        [
+            indice(i_elem, j_elem),
+            indice(i_elem + 1, j_elem),
+            indice(i_elem + 1, j_succ),
+            indice(i_elem, j_succ),
+        ],
+        axis=1,
+    )
+
+    # Segmenti dei due bordi assiali (x = 0 e x = lunghezza): anelli chiusi di
+    # n_elementi_theta segmenti, usati per la convezione sullo spessore.
+    j = np.arange(n_elementi_theta)
+    j_dopo = (j + 1) % n_elementi_theta
+    segmenti_bordo = np.concatenate(
+        [
+            np.stack([indice(0, j), indice(0, j_dopo)], axis=1),
+            np.stack(
+                [indice(n_elementi_x, j), indice(n_elementi_x, j_dopo)], axis=1
+            ),
+        ]
+    )
+
+    return {
+        "n_elementi_x": n_elementi_x,
+        "n_elementi_theta": n_elementi_theta,
+        "n_elementi": n_elementi,
+        "n_nodi_x": n_nodi_x,
+        "n_nodi": n_nodi,
+        "lunghezza_m": lunghezza_m,
+        "raggio_medio_m": raggio_medio_m,
+        "circonferenza_m": circonferenza_m,
+        "lato_x_m": lato_x_m,
+        "lato_arco_m": lato_arco_m,
+        "x_nodi_m": x_nodi_m,
+        "theta_nodi_rad": theta_nodi_rad,
+        "arco_nodi_m": arco_nodi_m,
+        "connettivita": connettivita,
+        "segmenti_bordo": segmenti_bordo,
+        "area_totale_m2": lunghezza_m * circonferenza_m,
+    }
+
+
+def campo_su_griglia(mesh: dict, T: np.ndarray) -> np.ndarray:
+    # Rimappa il vettore nodale sulla griglia (n_nodi_x, n_elementi_theta),
+    # comoda per le mappe sviluppate e per le superfici 3D.
+    return T.reshape(mesh["n_nodi_x"], mesh["n_elementi_theta"])
+
+
+def coordinate_3d(mesh: dict, chiudi_circonferenza: bool = True):
+    # Coordinate cartesiane dei nodi per la vista 3D: x lungo l'asse del
+    # cilindro, y e z sulla sezione circolare.
+    #
+    # Con chiudi_circonferenza si ripete la prima colonna in coda, così la
+    # superficie disegnata non mostra una fessura in theta = 0. È solo una
+    # necessità di disegno: il nodo ripetuto non è un grado di libertà.
+    theta = mesh["theta_nodi_rad"]
+    if chiudi_circonferenza:
+        theta = np.append(theta, 2.0 * np.pi)
+    X, Theta = np.meshgrid(mesh["x_nodi_m"], theta, indexing="ij")
+    R = mesh["raggio_medio_m"]
+    return X, R * np.sin(Theta), R * np.cos(Theta)
+
+
+def chiudi_campo(campo: np.ndarray) -> np.ndarray:
+    # Ripete la prima colonna circonferenziale in coda, in accordo con
+    # coordinate_3d(chiudi_circonferenza=True).
+    return np.concatenate([campo, campo[:, :1]], axis=1)

plot_animazione.py

@@ -1,9 +1,14 @@
-# Animazione della sezione della fascetta durante il passaggio delle sorgenti.
+# Animazione della superficie della fascetta durante il passaggio delle
+# sorgenti.
 #
 # Riproduce la fisica di run_0001 (stesso seed di simulate.py) e mostra:
-# - il profilo di flusso termico q(x) sul lato esterno e le sorgenti in moto;
-# - il campo di temperatura T(x, z) nella sezione lunghezza × spessore;
+# - il profilo di flusso termico q(x) lungo la generatrice delle sorgenti;
+# - il campo di temperatura sulla superficie cilindrica sviluppata in piano,
+#   cioè T(x, arco) con l'arco centrato sulla generatrice del sensore;
 # - il sensore infrarosso e la temperatura nel punto osservato.
+#
+# La mappa sviluppata è periodica in verticale: il bordo superiore e quello
+# inferiore sono lo stesso punto della fascetta.
 
 import random
 from pathlib import Path
@@ -12,13 +17,17 @@ import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib.animation import FuncAnimation, PillowWriter
+from matplotlib.collections import LineCollection
 
+import sensore as sens
+import sorgente as src
 from config import SIMULAZIONE
+from mesh import campo_su_griglia, chiudi_campo
 from simulate import (
+    campo_iniziale,
     configurazione_randomizzata,
-    passo_implicito,
+    passo_termico,
     prepara_stato_termico,
-    profilo_flusso_incidente_W_m2,
 )
 
 # Istante di inizio dei fotogrammi mostrati (la simulazione parte comunque da 0).
@@ -33,27 +42,41 @@ DT_FRAME_S = 0.05
 # Millisecondi tra i fotogrammi in riproduzione.
 INTERVALLO_RIPRODUZIONE_MS = 30
 
+# Disegna i bordi degli elementi shell sopra il campo di temperatura.
+# Con la mesh di default (100 x 110, elementi di 1 x 2 mm) il reticolo è più
+# fitto dei pixel disponibili e copre il campo: per ispezionare la mesh
+# conviene attivarlo e insieme abbassare n_elementi_x e n_elementi_theta in
+# config.py, per esempio a 25 x 28.
+MOSTRA_MESH = False
+
 
 def simula_campi(cfg_run: dict) -> dict:
     # Esegue la simulazione fino a T_FINE_ANIMAZIONE_S salvando, a ogni
     # fotogramma, campo di temperatura, profilo di flusso e stato del sensore.
+    #
+    # I campi sono restituiti come griglie (n_nodi_x, n_elementi_theta) con la
+    # colonna dell'arco ruotata in modo da portare la generatrice del sensore
+    # (theta = 0) al centro: è solo una scelta di visualizzazione.
     fascetta = cfg_run["fascetta"]
     aria = cfg_run["aria"]
     sorgente = cfg_run["sorgente"]
     sensore = cfg_run["sensore"]
 
-    stato = prepara_stato_termico(fascetta, aria, sorgente)
-    n_x = stato["n_x"]
-    n_z = stato["n_z"]
-    x_centri = stato["x_centri_m"]
+    stato = prepara_stato_termico(fascetta, aria, sorgente, sensore)
+    mesh = stato["mesh"]
     dt = stato["dt_s"]
+    interpolatore = stato["interpolatore_sensore"]
 
-    x_sensore = sensore["x_mm"] / 1000.0
-    i_sensore = min(n_x - 1, max(0, int(x_sensore / stato["dx_m"])))
+    # Dopo np.roll(..., rotazione) la colonna theta = 0 finisce all'indice
+    # "rotazione", quindi l'arco con segno della colonna m vale
+    # (m - rotazione) * lato_arco.
+    rotazione = mesh["n_elementi_theta"] // 2
+    arco_centrato_m = (
+        np.arange(mesh["n_elementi_theta"]) - rotazione
+    ) * mesh["lato_arco_m"]
 
-    T = np.full((n_x, n_z), aria["temperatura_ambiente_C"], dtype=float)
-    T_sensore = T[i_sensore, -1]
-    tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
+    T = campo_iniziale(stato)
+    T_letta = sens.temperatura_osservata_C(interpolatore, T)
 
     tempi, campi, flussi, x_riferimenti = [], [], [], []
     T_vere, T_lette = [], []
@@ -61,18 +84,22 @@ def simula_campi(cfg_run: dict) -> dict:
     prossimo_frame_t = 0.0
     t = 0.0
     while t <= T_FINE_ANIMAZIONE_S + 1e-12:
-        x_rif, q_x = profilo_flusso_incidente_W_m2(sorgente, x_sensore, t, x_centri)
-        T = passo_implicito(stato, T, q_x)
-
-        T_sensore += (T[i_sensore, -1] - T_sensore) * dt / tau_sensore
+        T, info = passo_termico(stato, T, t)
+        T_vera = sens.temperatura_osservata_C(interpolatore, T)
+        T_letta = sens.aggiorna_inerzia_C(
+            T_letta, T_vera, dt, sensore["costante_tempo_s"]
+        )
 
         if t + 1e-12 >= prossimo_frame_t:
+            _, q_nodi = src.profilo_flusso_x_W_m2(
+                sorgente, stato["x_sensore_m"], t, mesh["x_nodi_m"]
+            )
             tempi.append(t)
-            campi.append(T.copy())
-            flussi.append(q_x.copy())
-            x_riferimenti.append(x_rif)
-            T_vere.append(T[i_sensore, -1])
-            T_lette.append(T_sensore)
+            campi.append(np.roll(campo_su_griglia(mesh, T), rotazione, axis=1))
+            flussi.append(q_nodi)
+            x_riferimenti.append(info["x_riferimento_m"])
+            T_vere.append(T_vera)
+            T_lette.append(T_letta)
             prossimo_frame_t += DT_FRAME_S
 
         t += dt
@@ -84,11 +111,18 @@ def simula_campi(cfg_run: dict) -> dict:
         "x_riferimenti": np.array(x_riferimenti),
         "T_vere": np.array(T_vere),
         "T_lette": np.array(T_lette),
-        "x_centri_mm": x_centri * 1000.0,
-        "spessore_mm": fascetta["spessore_mm"],
+        "mesh": mesh,
+        "rotazione": rotazione,
+        "x_nodi_mm": mesh["x_nodi_m"] * 1000.0,
+        "arco_centrato_mm": arco_centrato_m * 1000.0,
+        "arco_chiuso_mm": np.append(
+            arco_centrato_m, arco_centrato_m[-1] + mesh["lato_arco_m"]
+        ) * 1000.0,
         "lunghezza_mm": fascetta["lunghezza_mm"],
         "x_sensore_mm": sensore["x_mm"],
-        "sorgente": cfg_run["sorgente"],
+        "offset_sorgente_mm": sorgente["offset_y_percorso_m"] * 1000.0,
+        "sorgente": sorgente,
+        "T_ambiente_C": aria["temperatura_ambiente_C"],
     }
 
 
@@ -102,23 +136,27 @@ def main() -> None:
     n_frame = len(tempi) - indice_inizio
 
     lunghezza_mm = dati["lunghezza_mm"]
-    spessore_mm = dati["spessore_mm"]
     x_vista_mm = (-10.0, lunghezza_mm + 10.0)
+    arco_mm = dati["arco_chiuso_mm"]
+    arco_min_mm, arco_max_mm = arco_mm[0], arco_mm[-1]
 
-    q_max_MW = max(q.max() for q in dati["flussi"]) / 1e6
+    # Il minimo evita un asse degenere se la finestra dell'animazione si
+    # chiude prima che le sorgenti entrino nella fascetta.
+    q_max_MW = max(max(q.max() for q in dati["flussi"]) / 1e6, 1e-3)
     T_max = max(c.max() for c in dati["campi"])
 
     sorgente = dati["sorgente"]
     numero_sorgenti = sorgente.get("numero_sorgenti", 1)
     distanza_mm = sorgente.get("distanza_sorgenti_m", 0.0) * 1000.0
+    offset_mm = dati["offset_sorgente_mm"]
 
-    fig, (ax_flusso, ax_sezione, ax_storia) = plt.subplots(
+    fig, (ax_flusso, ax_superficie, ax_storia) = plt.subplots(
         3, 1, figsize=(10, 8), height_ratios=[1.0, 1.6, 1.2],
         gridspec_kw={"hspace": 0.45},
     )
-    fig.suptitle("Sezione della fascetta: sorgenti in transito e sensore")
+    fig.suptitle("Superficie sviluppata della fascetta: sorgenti in transito e sensore")
 
-    # Pannello 1: profilo di flusso sul lato esterno e posizioni delle sorgenti.
+    # Pannello 1: profilo di flusso lungo x e posizioni delle sorgenti.
     linea_flusso, = ax_flusso.plot([], [], color="tab:red")
     marker_sorgenti, = ax_flusso.plot(
         [], [], "v", color="tab:red", markersize=10, clip_on=False
@@ -135,37 +173,50 @@ def main() -> None:
     ax_flusso.grid(True, alpha=0.3)
     ax_flusso.set_xticklabels([])
 
-    # Pannello 2: campo di temperatura nella sezione (z verso il basso,
-    # origine nel vertice in alto a sinistra come nel modello).
-    immagine = ax_sezione.imshow(
-        dati["campi"][indice_inizio].T,
-        extent=(0.0, lunghezza_mm, spessore_mm, 0.0),
-        aspect="auto",
+    # Pannello 2: temperatura sulla superficie sviluppata. L'asse verticale è
+    # l'arco circonferenziale, con 0 sulla generatrice del sensore.
+    #
+    # pcolormesh con shading "gouraud" interpola linearmente tra i valori
+    # nodali: è esattamente ciò che fanno le funzioni di forma bilineari
+    # dell'elemento Q4, quindi il disegno è la soluzione dell'elemento finito
+    # e non un'interpolazione grafica arbitraria.
+    Xg, Yg = np.meshgrid(dati["x_nodi_mm"], dati["arco_chiuso_mm"], indexing="ij")
+    mappa = ax_superficie.pcolormesh(
+        Xg, Yg,
+        chiudi_campo(dati["campi"][indice_inizio]),
+        shading="gouraud",
         cmap="inferno",
-        vmin=cfg_run["aria"]["temperatura_ambiente_C"],
+        vmin=dati["T_ambiente_C"],
         vmax=T_max,
-        interpolation="bilinear",
     )
-    ax_sezione.set_xlim(*x_vista_mm)
-    ax_sezione.set_ylim(3.2 * spessore_mm, -0.6 * spessore_mm)
-    ax_sezione.set_ylabel("z [mm]")
-    ax_sezione.set_xlabel("x [mm]")
-    # Sensore infrarosso sotto la parete interna (posizione schematica,
-    # non in scala) con linea di vista tratteggiata.
-    x_sens = dati["x_sensore_mm"]
-    ax_sezione.plot([x_sens], [2.4 * spessore_mm], "^", color="tab:blue", markersize=12)
-    ax_sezione.plot(
-        [x_sens, x_sens], [1.1 * spessore_mm, 2.1 * spessore_mm],
-        linestyle="--", color="tab:blue", linewidth=1,
+    if MOSTRA_MESH:
+        # Bordi degli elementi shell: linee lungo x a ogni nodo
+        # circonferenziale e viceversa.
+        segmenti = [
+            [(x, arco_min_mm), (x, arco_max_mm)] for x in dati["x_nodi_mm"]
+        ] + [
+            [(0.0, a), (lunghezza_mm, a)] for a in dati["arco_chiuso_mm"]
+        ]
+        ax_superficie.add_collection(
+            LineCollection(segmenti, colors="white", linewidths=0.5, alpha=0.55)
+        )
+    ax_superficie.set_xlim(*x_vista_mm)
+    ax_superficie.set_ylim(arco_min_mm, arco_max_mm)
+    ax_superficie.set_ylabel("arco circonferenziale [mm]")
+    ax_superficie.set_xlabel("x [mm]")
+    ax_superficie.plot(
+        [dati["x_sensore_mm"]], [0.0], "^", color="tab:blue", markersize=11
     )
-    ax_sezione.text(
-        x_sens + 3, 2.4 * spessore_mm, "sensore IR",
-        color="tab:blue", va="center",
+    ax_superficie.text(
+        dati["x_sensore_mm"] + 3, 0.0, "sensore IR", color="tab:blue", va="center"
+    )
+    ax_superficie.axhline(
+        offset_mm, color="tab:cyan", linewidth=0.8, linestyle="--", alpha=0.7
     )
     # La colorbar è agganciata a tutti i pannelli per non restringere solo
-    # quello della sezione, mantenendo allineati gli assi x.
+    # quello della superficie, mantenendo allineati gli assi x.
     barra = fig.colorbar(
-        immagine, ax=(ax_flusso, ax_sezione, ax_storia), pad=0.02, aspect=35
+        mappa, ax=(ax_flusso, ax_superficie, ax_storia), pad=0.02, aspect=35
     )
     barra.set_label("T [°C]")
 
@@ -186,7 +237,7 @@ def main() -> None:
         k = indice_inizio + frame
         t = tempi[k]
 
-        linea_flusso.set_data(dati["x_centri_mm"], dati["flussi"][k] / 1e6)
+        linea_flusso.set_data(dati["x_nodi_mm"], dati["flussi"][k] / 1e6)
 
         x_sorgenti_mm = (
             dati["x_riferimenti"][k] * 1000.0
@@ -198,7 +249,7 @@ def main() -> None:
             np.full(int(visibili.sum()), q_max_MW * 1.12),
         )
 
-        immagine.set_data(dati["campi"][k].T)
+        mappa.set_array(chiudi_campo(dati["campi"][k]))
 
         linea_vera.set_data(tempi[: k + 1], dati["T_vere"][: k + 1])
         linea_letta.set_data(tempi[: k + 1], dati["T_lette"][: k + 1])
@@ -206,7 +257,7 @@ def main() -> None:
         testo_tempo.set_text(f"t = {t:.3f} s")
 
         return (
-            linea_flusso, marker_sorgenti, immagine,
+            linea_flusso, marker_sorgenti, mappa,
             linea_vera, linea_letta, cursore, testo_tempo,
         )
 
@@ -216,7 +267,7 @@ def main() -> None:
 
     # Se il backend non è interattivo si salva una GIF invece di mostrare la finestra.
     if matplotlib.get_backend().lower() == "agg":
-        percorso = Path("dataset") / "animazione_sezione.gif"
+        percorso = Path("dataset") / "animazione_superficie.gif"
         animazione.save(
             percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS)
         )

plot_animazione_3d.py

@@ -1,26 +1,26 @@
 # Animazione 3D del barattolo in vista isometrica: colormap della temperatura
-# sulla superficie esterna.
+# sulla superficie cilindrica.
 #
-# Il modello risolve solo la sezione x-z (vedi CLAUDE.md): la coordinata
-# circonferenziale y non è simulata, è collassata in un'attenuazione
-# gaussiana del flusso. Per estrudere il campo attorno alla circonferenza si
-# riusa la stessa gaussiana: la sovratemperatura rispetto al piano y=0 (dove
-# si trova il sensore) viene scalata in funzione della distanza
-# circonferenziale dal percorso delle sorgenti. È quindi una ricostruzione
-# visiva, non un calcolo di diffusione in y.
+# Il modello shell risolve la temperatura sull'intera superficie, quindi il
+# campo disegnato qui è quello calcolato: non c'è più nessuna ricostruzione
+# della distribuzione circonferenziale. La faccia esterna e quella interna
+# hanno la stessa temperatura (una temperatura per nodo), quindi la superficie
+# mostrata vale per entrambe.
 
 import random
+from pathlib import Path
 
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib import cm
 from matplotlib.animation import FuncAnimation, PillowWriter
-from pathlib import Path
 
-from config import FASCETTA, SIMULAZIONE
+import sorgente as src
+from config import SIMULAZIONE
+from mesh import chiudi_campo, coordinate_3d
 from plot_animazione import DT_FRAME_S as DT_FRAME_S_BASE
-from plot_animazione import T_FINE_ANIMAZIONE_S, T_INIZIO_ANIMAZIONE_S, simula_campi
+from plot_animazione import T_INIZIO_ANIMAZIONE_S, simula_campi
 from simulate import configurazione_randomizzata
 
 # Tempo simulato tra un fotogramma e il successivo dell'animazione 3D. Più
@@ -31,19 +31,12 @@ DT_FRAME_S = 0.15
 # Millisecondi tra i fotogrammi in riproduzione.
 INTERVALLO_RIPRODUZIONE_MS = 60
 
-# Numero di punti lungo la circonferenza per l'estrusione della superficie.
-N_THETA = 72
-
-
-def attenuazione_circonferenziale(
-    y_m: np.ndarray, offset_y_m: float, sigma_m: float
-) -> np.ndarray:
-    # Fattore che scala la sovratemperatura (T - T_ambiente) del piano y=0
-    # in funzione della distanza circonferenziale y dal piano stesso,
-    # normalizzato in modo da valere 1 in y=0.
-    sigma = max(sigma_m, 1e-9)
-    esponente = -0.5 * ((y_m - offset_y_m) ** 2 - offset_y_m ** 2) / (sigma * sigma)
-    return np.exp(esponente)
+# Disegna i bordi degli elementi shell sulla superficie del cilindro.
+# Con la mesh di default (100 x 110, elementi di 1 x 2 mm) il reticolo satura
+# il disegno e sbiadisce i colori: per ispezionare la mesh conviene attivarlo
+# e insieme abbassare n_elementi_x e n_elementi_theta in config.py, per
+# esempio a 25 x 28.
+MOSTRA_MESH = False
 
 
 def main() -> None:
@@ -56,26 +49,19 @@ def main() -> None:
     passo = max(1, round(DT_FRAME_S / DT_FRAME_S_BASE))
     indici_frame = list(range(indice_inizio, len(tempi), passo))
 
-    T_ambiente = cfg_run["aria"]["temperatura_ambiente_C"]
-    sorgente = cfg_run["sorgente"]
-    sigma_m = sorgente["sigma_punto_m"]
-    offset_y_m = sorgente["offset_y_percorso_m"]
-    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
-    distanza_m = sorgente.get("distanza_sorgenti_m", 0.0)
+    mesh = dati["mesh"]
+    # simula_campi ruota le colonne per centrare il sensore nella mappa
+    # sviluppata: qui serve l'ordine originale della mesh.
+    rotazione = dati["rotazione"]
 
-    raggio_m = (FASCETTA["diametro_mm"] / 1000.0) / 2.0
-    x_centri_m = dati["x_centri_mm"] / 1000.0
-    lunghezza_m = dati["lunghezza_mm"] / 1000.0
+    T_ambiente = dati["T_ambiente_C"]
+    sorgente = dati["sorgente"]
+    raggio_m = mesh["raggio_medio_m"]
+    lunghezza_m = mesh["lunghezza_m"]
 
-    theta = np.linspace(-np.pi, np.pi, N_THETA)
-    y_circ_m = theta * raggio_m
-    attenuazione = attenuazione_circonferenziale(y_circ_m, offset_y_m, sigma_m)
+    Xm, Ym, Zm = coordinate_3d(mesh)
 
-    Xm, Thetam = np.meshgrid(x_centri_m, theta)
-    Ym = raggio_m * np.sin(Thetam)
-    Zm = raggio_m * np.cos(Thetam)
-
-    T_max = max(c[:, 0].max() for c in dati["campi"])
+    T_max = max(c.max() for c in dati["campi"])
     # vmin più basso della temperatura ambiente reale: altrimenti le zone
     # fredde cadrebbero sul nero puro di "inferno" e, essendo lo shading
     # moltiplicativo, nessuna illuminazione basterebbe a renderle visibili.
@@ -101,25 +87,29 @@ def main() -> None:
     barra = fig.colorbar(mappabile, ax=ax, shrink=0.6, pad=0.05)
     barra.set_label("T [°C]")
 
+    theta_sorgente = src.theta_percorso_rad(sorgente, raggio_m)
+
     def disegna_frame(k: int):
         ax.cla()
         ax.view_init(elev=35.264, azim=45)
         ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m))
         ax.set_axis_off()
 
-        T_lato_esterno = dati["campi"][k][:, 0]
-        T_superficie = T_ambiente + (T_lato_esterno[None, :] - T_ambiente) * attenuazione[:, None]
-        colori = cmap(norm(T_superficie))
+        campo = np.roll(dati["campi"][k], -rotazione, axis=1)
+        colori = cmap(norm(chiudi_campo(campo)))
 
+        # Ogni faccia disegnata è un elemento shell della mesh: con
+        # MOSTRA_MESH se ne vedono i bordi.
+        bordi = {"edgecolor": (1.0, 1.0, 1.0, 0.25), "linewidth": 0.15}
         ax.plot_surface(
             Xm, Ym, Zm, facecolors=colori, rstride=1, cstride=1,
-            antialiased=False, shade=True, lightsource=lightsource, linewidth=0,
+            antialiased=False, shade=True, lightsource=lightsource,
+            **(bordi if MOSTRA_MESH else {"linewidth": 0}),
         )
 
-        x_sorgenti_m = dati["x_riferimenti"][k] + np.arange(numero_sorgenti) * distanza_m
+        x_sorgenti_m = src.posizioni_sorgenti_m(sorgente, dati["x_riferimenti"][k])
         visibili = (x_sorgenti_m >= 0.0) & (x_sorgenti_m <= lunghezza_m)
         if visibili.any():
-            theta_sorgente = offset_y_m / raggio_m
             ax.scatter(
                 x_sorgenti_m[visibili],
                 np.full(int(visibili.sum()), raggio_m * np.sin(theta_sorgente) * 1.05),

sensore.py

@@ -0,0 +1,60 @@
+# Pirometro a infrarossi montato all'interno della fascetta.
+#
+# Il sensore osserva senza contatto un punto fisso della superficie interna,
+# individuato dalla coordinata assiale x e dalla posizione circonferenziale
+# theta = 0, che per convenzione è la generatrice di riferimento del modello.
+#
+# Con una shell a una temperatura per nodo la faccia interna e quella esterna
+# hanno la stessa temperatura: il valore osservato è la temperatura della
+# shell nel punto, interpolata con le funzioni di forma dell'elemento che lo
+# contiene. L'approssimazione regge finché la parete è termicamente sottile
+# rispetto al transito della sorgente (vedi numero_fourier_spessore in
+# simulate.py).
+#
+# Alla temperatura vera si sovrappone la catena di misura: inerzia del primo
+# ordine, rumore gaussiano e quantizzazione della lettura.
+
+import numpy as np
+
+import elementi_shell as el
+
+
+def prepara_interpolatore(mesh: dict, x_m: float, theta_rad: float) -> dict:
+    # Individua l'elemento che contiene il punto osservato e i pesi nodali
+    # con cui ricostruirvi la temperatura.
+    n_x = mesh["n_elementi_x"]
+    n_theta = mesh["n_elementi_theta"]
+
+    x_limitato = min(max(x_m, 0.0), mesh["lunghezza_m"])
+    i = min(n_x - 1, max(0, int(x_limitato / mesh["lato_x_m"])))
+    xi = 2.0 * (x_limitato - mesh["x_nodi_m"][i]) / mesh["lato_x_m"] - 1.0
+
+    arco = (theta_rad * mesh["raggio_medio_m"]) % mesh["circonferenza_m"]
+    j = int(arco / mesh["lato_arco_m"]) % n_theta
+    eta = 2.0 * (arco - mesh["arco_nodi_m"][j]) / mesh["lato_arco_m"] - 1.0
+
+    elemento = i * n_theta + j
+    return {
+        "nodi": mesh["connettivita"][elemento],
+        "pesi": el.funzioni_forma(xi, eta),
+        "x_m": x_limitato,
+        "theta_rad": theta_rad,
+    }
+
+
+def temperatura_osservata_C(interpolatore: dict, T: np.ndarray) -> float:
+    return float(T[interpolatore["nodi"]] @ interpolatore["pesi"])
+
+
+def aggiorna_inerzia_C(
+    T_letta_C: float, T_vera_C: float, dt_s: float, costante_tempo_s: float
+) -> float:
+    # Risposta del primo ordine del sensore.
+    tau = max(costante_tempo_s, 1e-9)
+    return T_letta_C + (T_vera_C - T_letta_C) * dt_s / tau
+
+
+def quantizza(valore: float, passo: float) -> float:
+    if passo <= 0.0:
+        return valore
+    return round(valore / passo) * passo

simulate.py

@@ -1,3 +1,11 @@
+# Motore principale: genera N run randomizzati della fascetta riscaldata da
+# sorgenti a induzione in movimento e scrive i CSV del dataset.
+#
+# Il modello termico è una shell a elementi finiti quadrilateri sull'intera
+# superficie cilindrica, con una temperatura per nodo. La pipeline di un run è
+# in prepara_stato_termico() e passo_termico(), condivise con i moduli di
+# animazione: ogni modifica alla fisica va fatta lì, non duplicata.
+
 import csv
 import math
 import os
@@ -8,11 +16,14 @@ from copy import deepcopy
 from pathlib import Path
 
 import numpy as np
-import scipy.sparse as sp
-from scipy.sparse.linalg import splu
 
+import assemblaggio
+import sensore as sens
+import solutore as sol
+import sorgente as src
 from config import ARIA, FASCETTA, RANDOMIZZAZIONE, SENSORE, SIMULAZIONE, SORGENTE
 from materials import MATERIALI
+from mesh import costruisci_mesh_cilindrica
 
 
 MU0 = 4.0 * math.pi * 1e-7
@@ -31,241 +42,120 @@ def calcola_skin_depth_m(materiale: dict, frequenza_hz: float) -> float:
     return math.sqrt(2.0 * rho_e / (omega * mu))
 
 
-def _spread_sorgenti_m(sorgente: dict) -> float:
-    # Distanza lungo x tra la prima e l'ultima sorgente del gruppo.
-    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
-    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
-    return (numero_sorgenti - 1) * distanza
-
-
-def _x_riferimento_iniziale_m(sorgente: dict, x_sensore_m: float) -> float:
-    # Posizione a t=0 della sorgente di indice 0. Le sorgenti i sono a
-    # x_i = riferimento + i * distanza, quindi per v >= 0 l'indice 0 è la
-    # più arretrata nel verso di marcia, per v < 0 è la più avanzata.
-    # x_inizio_m è la distanza dal sensore della sorgente più avanzata
-    # (quella che lo raggiunge per prima).
-    spread = _spread_sorgenti_m(sorgente)
-    x_inizio = sorgente["x_inizio_m"]
-    if sorgente["velocita_m_s"] >= 0:
-        return (x_sensore_m - x_inizio) - spread
-    return x_sensore_m + x_inizio
-
-
-def _x_riferimento_finale_m(sorgente: dict, x_sensore_m: float) -> float:
-    # Posizione di fine corsa della sorgente di indice 0. x_fine_m è la
-    # distanza dal sensore della sorgente più arretrata nel verso di marcia
-    # (quella che lo supera per ultima).
-    spread = _spread_sorgenti_m(sorgente)
-    x_fine = sorgente["x_fine_m"]
-    if sorgente["velocita_m_s"] >= 0:
-        return x_sensore_m + x_fine
-    return (x_sensore_m - x_fine) - spread
-
-
-def _intervallo_attivo(inizio: float, fine: float, v: float, x_m: float) -> bool:
-    if v >= 0:
-        return inizio <= x_m <= fine
-    return fine <= x_m <= inizio
-
-
-def profilo_flusso_incidente_W_m2(
-    sorgente: dict,
-    x_sensore_m: float,
-    t_s: float,
-    x_centri_m: np.ndarray,
-) -> tuple[float, np.ndarray]:
-    # Restituisce x_sorgente_m (posizione della sorgente di riferimento) e
-    # il profilo di flusso termico efficace q(x) [W/m²] sul lato esterno,
-    # somma dei contributi di tutte le sorgenti attive.
+def numero_fourier_spessore(
+    materiale: dict, spessore_m: float, sorgente: dict
+) -> float:
+    # Verifica quantitativa dell'ipotesi di parete termicamente sottile:
     #
-    # Ogni sorgente in moto ha un'impronta gaussiana lungo x, valutata sui
-    # centri cella della sezione. L'offset circonferenziale y non è risolto
-    # spazialmente: entra come attenuazione gaussiana del flusso.
-    x_rif_iniziale = _x_riferimento_iniziale_m(sorgente, x_sensore_m)
-    x_rif_finale = _x_riferimento_finale_m(sorgente, x_sensore_m)
-    x_riferimento = x_rif_iniziale + sorgente["velocita_m_s"] * t_s
-
-    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
-    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
-    v = sorgente["velocita_m_s"]
-    zero_dopo_fine = sorgente.get("zero_dopo_fine", True)
-    dy = sorgente["offset_y_percorso_m"]
-    sigma = sorgente["sigma_punto_m"]
-
-    q_picco = sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"]
-    attenuazione_y = math.exp(-0.5 * (dy * dy) / (sigma * sigma))
-
-    q_x = np.zeros_like(x_centri_m)
-    for i in range(numero_sorgenti):
-        x_i = x_riferimento + i * distanza
-
-        if zero_dopo_fine:
-            inizio_i = x_rif_iniziale + i * distanza
-            fine_i = x_rif_finale + i * distanza
-            if not _intervallo_attivo(inizio_i, fine_i, v, x_i):
-                continue
-
-        dx = x_centri_m - x_i
-        q_x += q_picco * attenuazione_y * np.exp(-0.5 * (dx * dx) / (sigma * sigma))
-
-    return x_riferimento, q_x
-
-
-def profilo_deposizione_z_1_m(
-    z_centri_m: np.ndarray,
-    spessore_m: float,
-    skin_depth_m: float,
-) -> np.ndarray:
-    # Profilo di deposizione volumetrica del flusso superficiale [1/m]:
+    #     Fo = alpha * (sigma / |v|) / spessore^2
     #
-    # p(z) = exp(-z / delta) / (delta * (1 - exp(-spessore / delta)))
-    #
-    # Normalizzato in modo che integrale_0^spessore p(z) dz = 1, così che
-    # q_vol(x, z) = q(x) * p(z) conservi il flusso superficiale.
-    delta = max(skin_depth_m, 1e-9)
-    normalizzazione = delta * (1.0 - math.exp(-spessore_m / delta))
-    return np.exp(-z_centri_m / delta) / normalizzazione
-
-
-def _laplaciano_1d(n: int) -> sp.spmatrix:
-    # Operatore alle differenze -T'' su n celle con bordi adiabatici (Neumann).
-    diagonale = np.full(n, 2.0)
-    diagonale[0] = 1.0
-    diagonale[-1] = 1.0
-    fuori = -np.ones(n - 1)
-    return sp.diags([fuori, diagonale, fuori], [-1, 0, 1])
-
-
-def costruisci_solutore_implicito_2d(
-    n_x: int,
-    n_z: int,
-    dt_s: float,
-    dx_m: float,
-    dz_m: float,
-    spessore_m: float,
-    materiale: dict,
-    h_esterno_W_m2K: float,
-    h_interno_W_m2K: float,
-    h_bordi_W_m2K: float,
-):
-    # Costruisce e fattorizza (LU sparsa) la matrice A per Eulero implicito 2D:
-    # A * T_next = rhs
-    #
-    # Le incognite sono i centri cella T[i, j] con i lungo x e j lungo z,
-    # appiattiti in ordine C (indice = i * n_z + j). Tutti e quattro i lati
-    # della sezione includono la convezione verso l'ambiente; su ogni cella
-    # agisce inoltre la conduzione circonferenziale (y) verso il resto della
-    # fascetta, assunto a temperatura ambiente. Il termine è ricavato
-    # considerando l'intero volume della fascia cilindrica (equazione
-    # dell'aletta): il calore conduce lungo y attraverso l'intero spessore
-    # mentre le superfici esterna e interna dell'intero cilindro scambiano
-    # per convezione, dando un sink distribuito uniformemente su ogni cella
-    # q_y = -(h_esterno + h_interno) / spessore * (T - T_amb).
-    k = materiale["conducibilita_termica_W_mK"]
-    rho = materiale["densita_kg_m3"]
-    cp = materiale["calore_specifico_J_kgK"]
-    alpha = k / (rho * cp)
-
-    r_x = alpha * dt_s / (dx_m * dx_m)
-    r_z = alpha * dt_s / (dz_m * dz_m)
-    b_esterno = h_esterno_W_m2K * dt_s / (rho * cp * dz_m)
-    b_interno = h_interno_W_m2K * dt_s / (rho * cp * dz_m)
-    b_bordo = h_bordi_W_m2K * dt_s / (rho * cp * dx_m)
-    c_y = (h_esterno_W_m2K + h_interno_W_m2K) * dt_s / (rho * cp * spessore_m)
-
-    n = n_x * n_z
-    scambio = np.full(n, c_y)
-    scambio[0::n_z] += b_esterno
-    scambio[n_z - 1::n_z] += b_interno
-    scambio[:n_z] += b_bordo
-    scambio[n - n_z:] += b_bordo
-
-    A = (
-        sp.identity(n)
-        + r_x * sp.kron(_laplaciano_1d(n_x), sp.identity(n_z))
-        + r_z * sp.kron(sp.identity(n_x), _laplaciano_1d(n_z))
-        + sp.diags(scambio)
+    # confronta il tempo di diffusione attraverso lo spessore con il tempo di
+    # transito della sorgente. Per Fo >> 1 lo spessore si equilibra molto
+    # prima che la sorgente sia passata e una shell a una temperatura per
+    # nodo è adeguata; scendendo verso Fo ~ 1 servirebbe una formulazione
+    # multistrato.
+    alpha = materiale["conducibilita_termica_W_mK"] / (
+        materiale["densita_kg_m3"] * materiale["calore_specifico_J_kgK"]
     )
-    return splu(sp.csc_matrix(A))
+    velocita = abs(sorgente["velocita_m_s"])
+    if velocita <= 0.0:
+        return float("inf")
+    tempo_transito = sorgente["sigma_punto_m"] / velocita
+    return alpha * tempo_transito / (spessore_m * spessore_m)
 
 
-def prepara_stato_termico(fascetta: dict, aria: dict, sorgente: dict) -> dict:
-    # Prepara griglia, coefficienti e solutore fattorizzato per un run:
+def prepara_stato_termico(
+    fascetta: dict, aria: dict, sorgente: dict, sensore: dict
+) -> dict:
+    # Costruisce mesh, matrici globali e solutore fattorizzato per un run:
     # tutto ciò che resta costante durante l'integrazione temporale.
     materiale = MATERIALI[fascetta["materiale"]]
 
-    lunghezza = fascetta["lunghezza_mm"] / 1000.0
-    spessore = fascetta["spessore_mm"] / 1000.0
-    n_x = fascetta["n_nodi_x"]
-    n_z = fascetta["n_nodi_z"]
-    dx = lunghezza / n_x
-    dz = spessore / n_z
+    lunghezza_m = fascetta["lunghezza_mm"] / 1000.0
+    spessore_m = fascetta["spessore_mm"] / 1000.0
+    raggio_medio_m = fascetta["diametro_medio_mm"] / 2000.0
 
-    dt = SIMULAZIONE["dt_interno_s"]
+    mesh = costruisci_mesh_cilindrica(
+        lunghezza_m=lunghezza_m,
+        raggio_medio_m=raggio_medio_m,
+        n_elementi_x=fascetta["n_elementi_x"],
+        n_elementi_theta=fascetta["n_elementi_theta"],
+    )
+
+    C = assemblaggio.assembla_capacita(mesh, materiale, spessore_m)
+    K_cond = assemblaggio.assembla_conduzione(mesh, materiale, spessore_m)
+    K_conv, f_ambiente = assemblaggio.assembla_convezione(mesh, aria, spessore_m)
+
+    dt_s = SIMULAZIONE["dt_interno_s"]
+    solutore = sol.costruisci_solutore(C, K_cond + K_conv, dt_s)
+
+    assemblatore = assemblaggio.prepara_assemblatore_sorgente(mesh)
+
+    # Il fattore circonferenziale del flusso non dipende dal tempo: le
+    # sorgenti traslano solo lungo x.
+    fattore_arco = src.fattore_circonferenziale(
+        sorgente, raggio_medio_m, assemblatore["arco_gauss_m"]
+    )
+    attenuazione_al_sensore = float(
+        src.fattore_circonferenziale(sorgente, raggio_medio_m, np.zeros(1))[0]
+    )
 
     if sorgente["skin_depth_fissa_m"] is None:
-        skin_depth = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
+        skin_depth_m = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
     else:
-        skin_depth = float(sorgente["skin_depth_fissa_m"])
+        skin_depth_m = float(sorgente["skin_depth_fissa_m"])
 
-    rho = materiale["densita_kg_m3"]
-    cp = materiale["calore_specifico_J_kgK"]
-    alpha = materiale["conducibilita_termica_W_mK"] / (rho * cp)
-
-    z_centri = (np.arange(n_z) + 0.5) * dz
+    x_sensore_m = sensore["x_mm"] / 1000.0
 
     return {
-        "n_x": n_x,
-        "n_z": n_z,
-        "dx_m": dx,
-        "dz_m": dz,
-        "dt_s": dt,
-        "x_centri_m": (np.arange(n_x) + 0.5) * dx,
-        "z_centri_m": z_centri,
-        "skin_depth_m": skin_depth,
-        "rho": rho,
-        "cp": cp,
-        "profilo_z": profilo_deposizione_z_1_m(z_centri, spessore, skin_depth),
-        "b_esterno": aria["h_esterno_W_m2K"] * dt / (rho * cp * dz),
-        "b_interno": aria["h_interno_W_m2K"] * dt / (rho * cp * dz),
-        "b_bordo": aria["h_bordi_W_m2K"] * dt / (rho * cp * dx),
-        "c_y": (aria["h_esterno_W_m2K"] + aria["h_interno_W_m2K"]) * dt / (rho * cp * spessore),
+        "mesh": mesh,
+        "materiale": materiale,
+        "spessore_m": spessore_m,
+        "raggio_medio_m": raggio_medio_m,
+        "dt_s": dt_s,
+        "solutore": solutore,
+        "assemblatore": assemblatore,
+        "f_ambiente": f_ambiente,
+        "fattore_arco": fattore_arco,
+        "attenuazione_al_sensore": attenuazione_al_sensore,
+        "x_sensore_m": x_sensore_m,
+        "interpolatore_sensore": sens.prepara_interpolatore(mesh, x_sensore_m, 0.0),
+        "sorgente": sorgente,
         "T_ambiente_C": aria["temperatura_ambiente_C"],
-        "solutore": costruisci_solutore_implicito_2d(
-            n_x=n_x,
-            n_z=n_z,
-            dt_s=dt,
-            dx_m=dx,
-            dz_m=dz,
-            spessore_m=spessore,
-            materiale=materiale,
-            h_esterno_W_m2K=aria["h_esterno_W_m2K"],
-            h_interno_W_m2K=aria["h_interno_W_m2K"],
-            h_bordi_W_m2K=aria["h_bordi_W_m2K"],
+        "skin_depth_m": skin_depth_m,
+        "numero_fourier_spessore": numero_fourier_spessore(
+            materiale, spessore_m, sorgente
         ),
     }
 
 
-def passo_implicito(stato: dict, T: np.ndarray, q_x: np.ndarray) -> np.ndarray:
-    # Avanza il campo di temperatura di un passo dt: assembla il termine noto
-    # (sorgente volumetrica, convezione sui quattro lati, conduzione
-    # circonferenziale verso l'ambiente) e risolve il sistema implicito.
-    T_amb = stato["T_ambiente_C"]
-    rhs = T + (stato["dt_s"] / (stato["rho"] * stato["cp"])) * (
-        q_x[:, None] * stato["profilo_z"][None, :]
+def campo_iniziale(stato: dict) -> np.ndarray:
+    return np.full(stato["mesh"]["n_nodi"], stato["T_ambiente_C"], dtype=float)
+
+
+def passo_termico(stato: dict, T: np.ndarray, t_s: float) -> tuple[np.ndarray, dict]:
+    # Avanza il campo di temperatura di un passo dt e restituisce anche lo
+    # stato della sorgente in quell'istante.
+    sorgente = stato["sorgente"]
+    x_riferimento_m, flusso_x = src.profilo_flusso_x_W_m2(
+        sorgente, stato["x_sensore_m"], t_s, stato["assemblatore"]["x_gauss_m"]
+    )
+    f_sorgente = assemblaggio.assembla_sorgente(
+        stato["assemblatore"], flusso_x, stato["fattore_arco"]
+    )
+    T_next = sol.passo_implicito(
+        stato["solutore"], T, stato["f_ambiente"], f_sorgente
     )
-    rhs += stato["c_y"] * T_amb
-    rhs[:, 0] += stato["b_esterno"] * T_amb
-    rhs[:, -1] += stato["b_interno"] * T_amb
-    rhs[0, :] += stato["b_bordo"] * T_amb
-    rhs[-1, :] += stato["b_bordo"] * T_amb
-    return stato["solutore"].solve(rhs.ravel()).reshape(T.shape)
 
-
-def quantizza(valore: float, passo: float) -> float:
-    if passo <= 0.0:
-        return valore
-    return round(valore / passo) * passo
+    _, flusso_al_sensore = src.profilo_flusso_x_W_m2(
+        sorgente, stato["x_sensore_m"], t_s, np.array([stato["x_sensore_m"]])
+    )
+    return T_next, {
+        "x_riferimento_m": x_riferimento_m,
+        "flusso_al_sensore_W_m2": float(
+            flusso_al_sensore[0] * stato["attenuazione_al_sensore"]
+        ),
+    }
 
 
 def configurazione_randomizzata(indice_run: int, rng: random.Random) -> dict:
@@ -323,25 +213,22 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
 
     nome_materiale = fascetta["materiale"]
 
-    stato = prepara_stato_termico(fascetta, aria, sorgente)
-    n_x = stato["n_x"]
-    n_z = stato["n_z"]
-    x_centri = stato["x_centri_m"]
+    stato = prepara_stato_termico(fascetta, aria, sorgente, sensore)
+    mesh = stato["mesh"]
+    interpolatore = stato["interpolatore_sensore"]
     skin_depth = stato["skin_depth_m"]
 
-    x_sensore = sensore["x_mm"] / 1000.0
-    i_sensore = min(n_x - 1, max(0, int(x_sensore / stato["dx_m"])))
-
     dt = stato["dt_s"]
     durata = SIMULAZIONE["durata_s"]
     periodo_campionamento = 1.0 / SIMULAZIONE["frequenza_campionamento_hz"]
 
-    T = np.full((n_x, n_z), aria["temperatura_ambiente_C"], dtype=float)
-    T_sensore = T[i_sensore, -1]
+    T = campo_iniziale(stato)
+    T_letta = sens.temperatura_osservata_C(interpolatore, T)
 
     prossimo_campione_t = 0.0
-    T_vera_max = T[i_sensore, -1]
-    T_misurata_max = T_sensore
+    T_vera_max = T_letta
+    T_misurata_max = T_letta
+    T_fascetta_max = float(T.max())
 
     output_csv.parent.mkdir(parents=True, exist_ok=True)
 
@@ -356,7 +243,8 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
             "skin_depth_m",
             "T_vera_lato_sensore_C",
             "T_misurata_sensore_C",
-            "T_lato_caldo_C",
+            "T_max_fascetta_C",
+            "T_media_fascetta_C",
             "T_ambiente_C",
             "velocita_m_s",
             "sigma_punto_m",
@@ -366,37 +254,41 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
 
         t = 0.0
         while t <= durata + 1e-12:
-            x_sorgente, q_x = profilo_flusso_incidente_W_m2(
-                sorgente, x_sensore, t, x_centri
-            )
-            T = passo_implicito(stato, T, q_x)
+            T, info = passo_termico(stato, T, t)
 
             # Temperatura vera della superficie interna nel punto osservato
             # dal sensore infrarosso.
-            T_vera_lato_sensore = T[i_sensore, -1]
+            T_vera = sens.temperatura_osservata_C(interpolatore, T)
+            T_letta = sens.aggiorna_inerzia_C(
+                T_letta, T_vera, dt, sensore["costante_tempo_s"]
+            )
 
-            # Inerzia del sensore del primo ordine.
-            tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
-            T_sensore += (T_vera_lato_sensore - T_sensore) * dt / tau_sensore
+            # I massimi delle temperature vere sono aggiornati a ogni passo
+            # interno, non solo agli istanti di campionamento: servono come
+            # ground truth e il picco è stretto rispetto al periodo di
+            # campionamento CSV, che da solo lo taglierebbe di oltre un grado.
+            T_max = float(T.max())
+            T_vera_max = max(T_vera_max, T_vera)
+            T_fascetta_max = max(T_fascetta_max, T_max)
 
             # Campionamento CSV.
             if t + 1e-12 >= prossimo_campione_t:
-                misurata = T_sensore + rng.gauss(0.0, sensore["rumore_std_C"])
-                misurata = quantizza(misurata, sensore["quantizzazione_C"])
+                misurata = T_letta + rng.gauss(0.0, sensore["rumore_std_C"])
+                misurata = sens.quantizza(misurata, sensore["quantizzazione_C"])
 
-                T_vera_max = max(T_vera_max, T_vera_lato_sensore)
                 T_misurata_max = max(T_misurata_max, misurata)
 
                 writer.writerow([
                     cfg_run["id_run"],
                     f"{t:.6f}",
-                    f"{x_sorgente:.9f}",
+                    f"{info['x_riferimento_m']:.9f}",
                     f"{sorgente['offset_y_percorso_m']:.9f}",
-                    f"{q_x[i_sensore]:.6f}",
+                    f"{info['flusso_al_sensore_W_m2']:.6f}",
                     f"{skin_depth:.9e}",
-                    f"{T_vera_lato_sensore:.6f}",
+                    f"{T_vera:.6f}",
                     f"{misurata:.6f}",
-                    f"{T[i_sensore, 0]:.6f}",
+                    f"{T_max:.6f}",
+                    f"{float(T.mean()):.6f}",
                     f"{aria['temperatura_ambiente_C']:.6f}",
                     f"{sorgente['velocita_m_s']:.9f}",
                     f"{sorgente['sigma_punto_m']:.9f}",
@@ -411,11 +303,14 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
         "id_run": cfg_run["id_run"],
         "file_csv": str(output_csv.name),
         "materiale": nome_materiale,
-        "diametro_m": fascetta["diametro_mm"] / 1000.0,
+        "diametro_medio_m": fascetta["diametro_medio_mm"] / 1000.0,
         "lunghezza_m": fascetta["lunghezza_mm"] / 1000.0,
         "spessore_m": fascetta["spessore_mm"] / 1000.0,
-        "n_nodi_x": n_x,
-        "n_nodi_z": n_z,
+        "n_elementi_x": mesh["n_elementi_x"],
+        "n_elementi_theta": mesh["n_elementi_theta"],
+        "n_nodi": mesh["n_nodi"],
+        "lato_elemento_x_m": mesh["lato_x_m"],
+        "lato_elemento_arco_m": mesh["lato_arco_m"],
         "durata_s": durata,
         "frequenza_campionamento_hz": SIMULAZIONE["frequenza_campionamento_hz"],
         "dt_interno_s": dt,
@@ -425,9 +320,12 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
         "h_bordi_W_m2K": aria["h_bordi_W_m2K"],
         "x_inizio_m": sorgente["x_inizio_m"],
         "x_fine_m": sorgente["x_fine_m"],
-        "x_sensore_m": x_sensore,
+        "x_sensore_m": stato["x_sensore_m"],
         "distanza_sensore_parete_m": sensore["distanza_parete_mm"] / 1000.0,
         "offset_y_percorso_m": sorgente["offset_y_percorso_m"],
+        "theta_percorso_rad": src.theta_percorso_rad(
+            sorgente, stato["raggio_medio_m"]
+        ),
         "velocita_m_s": sorgente["velocita_m_s"],
         "numero_sorgenti": sorgente.get("numero_sorgenti", 1),
         "distanza_sorgenti_m": sorgente.get("distanza_sorgenti_m", 0.0),
@@ -436,11 +334,13 @@ def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
         "efficienza_riscaldamento": sorgente["efficienza_riscaldamento"],
         "frequenza_hz": sorgente["frequenza_hz"],
         "skin_depth_m": skin_depth,
+        "numero_fourier_spessore": stato["numero_fourier_spessore"],
         "costante_tempo_sensore_s": sensore["costante_tempo_s"],
         "rumore_std_sensore_C": sensore["rumore_std_C"],
         "quantizzazione_sensore_C": sensore["quantizzazione_C"],
         "T_vera_max_lato_sensore_C": T_vera_max,
         "T_misurata_max_sensore_C": T_misurata_max,
+        "T_max_fascetta_C": T_fascetta_max,
     }
 
 

solutore.py

@@ -0,0 +1,39 @@
+# Integrazione temporale implicita del sistema semidiscreto della shell.
+#
+# Dopo l'assemblaggio il problema è
+#
+#     C * dT/dt + K * T = f(t)        con K = K_cond + K_conv
+#
+# Con Eulero implicito:
+#
+#     (C/dt + K) * T^(n+1) = (C/dt) * T^n + f^(n+1)
+#
+# La matrice A = C/dt + K dipende solo da mesh, materiale, passo temporale e
+# condizioni convettive: nessuno di questi cambia durante un run, quindi si
+# fattorizza LU una volta sola e a ogni passo si risolvono solo i sistemi
+# triangolari. Se in futuro si introducessero proprietà dipendenti dalla
+# temperatura, A andrebbe ricostruita e rifattorizzata a ogni passo.
+
+import scipy.sparse as sp
+from scipy.sparse.linalg import splu
+
+
+def costruisci_solutore(C, K, dt_s: float) -> dict:
+    # A è simmetrica e definita positiva: capacità, conduzione e convezione
+    # lo sono tutte. L'ordinamento MMD sul grafo di A + A^T sfrutta la
+    # simmetria e riduce il riempimento della fattorizzazione rispetto al
+    # COLAMD di default, pensato per matrici non simmetriche.
+    A = (C / dt_s + K).tocsc()
+    return {
+        "lu": splu(A, permc_spec="MMD_AT_PLUS_A"),
+        "C_su_dt": (C / dt_s).tocsr(),
+        "dt_s": dt_s,
+    }
+
+
+def passo_implicito(solutore: dict, T, f_costante, f_sorgente):
+    # Avanza il campo di un passo dt. f_costante raccoglie i termini noti
+    # indipendenti dal tempo (convezione verso l'ambiente su facce e bordi),
+    # f_sorgente il carico della sorgente all'istante di arrivo.
+    rhs = solutore["C_su_dt"] @ T + f_costante + f_sorgente
+    return solutore["lu"].solve(rhs)

sorgente.py

@@ -0,0 +1,125 @@
+# Sorgenti a induzione in movimento sulla superficie esterna della fascetta.
+#
+# Il gruppo è formato da più sorgenti equidistanti che traslano insieme lungo
+# x con la stessa velocità. Tutte percorrono la stessa generatrice, cioè la
+# stessa posizione circonferenziale, definita dall'offset rispetto al punto
+# osservato dal sensore.
+#
+# Il flusso termico incidente su un punto della superficie è
+#
+#     q''(x, arco, t) = q_picco * somma_i exp(-((x - x_i(t))^2) / (2 sigma^2))
+#                                * exp(-(delta_arco^2) / (2 sigma^2))
+#
+# dove delta_arco è la distanza circonferenziale minima tra il punto e il
+# percorso delle sorgenti, misurata lungo la circonferenza e riavvolta su
+# 2*pi*R: passare da un lato o dall'altro della fascetta è equivalente.
+#
+# Il flusso è quindi il prodotto di un fattore che dipende solo da x e varia
+# nel tempo e di un fattore che dipende solo dall'arco ed è costante: questa
+# separazione è sfruttata dall'assemblaggio del vettore sorgente.
+
+import math
+
+import numpy as np
+
+
+def _spread_sorgenti_m(sorgente: dict) -> float:
+    # Distanza lungo x tra la prima e l'ultima sorgente del gruppo.
+    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
+    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
+    return (numero_sorgenti - 1) * distanza
+
+
+def x_riferimento_iniziale_m(sorgente: dict, x_sensore_m: float) -> float:
+    # Posizione a t=0 della sorgente di indice 0. Le sorgenti i sono a
+    # x_i = riferimento + i * distanza, quindi per v >= 0 l'indice 0 è la
+    # più arretrata nel verso di marcia, per v < 0 è la più avanzata.
+    # x_inizio_m è la distanza dal sensore della sorgente più avanzata
+    # (quella che lo raggiunge per prima).
+    spread = _spread_sorgenti_m(sorgente)
+    x_inizio = sorgente["x_inizio_m"]
+    if sorgente["velocita_m_s"] >= 0:
+        return (x_sensore_m - x_inizio) - spread
+    return x_sensore_m + x_inizio
+
+
+def x_riferimento_finale_m(sorgente: dict, x_sensore_m: float) -> float:
+    # Posizione di fine corsa della sorgente di indice 0. x_fine_m è la
+    # distanza dal sensore della sorgente più arretrata nel verso di marcia
+    # (quella che lo supera per ultima).
+    spread = _spread_sorgenti_m(sorgente)
+    x_fine = sorgente["x_fine_m"]
+    if sorgente["velocita_m_s"] >= 0:
+        return x_sensore_m + x_fine
+    return (x_sensore_m - x_fine) - spread
+
+
+def _intervallo_attivo(inizio: float, fine: float, v: float, x_m: float) -> bool:
+    if v >= 0:
+        return inizio <= x_m <= fine
+    return fine <= x_m <= inizio
+
+
+def theta_percorso_rad(sorgente: dict, raggio_medio_m: float) -> float:
+    # Posizione circonferenziale del percorso delle sorgenti, ricavata
+    # dall'offset in arco rispetto al sensore, che sta per convenzione a
+    # theta = 0.
+    return sorgente["offset_y_percorso_m"] / raggio_medio_m
+
+
+def fattore_circonferenziale(
+    sorgente: dict, raggio_medio_m: float, arco_m: np.ndarray
+) -> np.ndarray:
+    # Attenuazione gaussiana del flusso in funzione della distanza
+    # circonferenziale dal percorso delle sorgenti. La distanza è riavvolta
+    # sulla circonferenza, così la sorgente scalda anche i punti raggiunti
+    # "dall'altra parte" della fascetta.
+    circonferenza = 2.0 * math.pi * raggio_medio_m
+    delta = arco_m - sorgente["offset_y_percorso_m"]
+    delta = (delta + circonferenza / 2.0) % circonferenza - circonferenza / 2.0
+    sigma = sorgente["sigma_punto_m"]
+    return np.exp(-0.5 * (delta * delta) / (sigma * sigma))
+
+
+def profilo_flusso_x_W_m2(
+    sorgente: dict, x_sensore_m: float, t_s: float, x_m: np.ndarray
+) -> tuple[float, np.ndarray]:
+    # Restituisce la posizione della sorgente di riferimento e il fattore
+    # assiale del flusso efficace [W/m²], somma dei contributi di tutte le
+    # sorgenti attive valutata nei punti x richiesti.
+    x_rif_iniziale = x_riferimento_iniziale_m(sorgente, x_sensore_m)
+    x_rif_finale = x_riferimento_finale_m(sorgente, x_sensore_m)
+    x_riferimento = x_rif_iniziale + sorgente["velocita_m_s"] * t_s
+
+    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
+    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
+    v = sorgente["velocita_m_s"]
+    zero_dopo_fine = sorgente.get("zero_dopo_fine", True)
+    sigma = sorgente["sigma_punto_m"]
+
+    q_picco = (
+        sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"]
+    )
+
+    q_x = np.zeros_like(x_m)
+    for i in range(numero_sorgenti):
+        x_i = x_riferimento + i * distanza
+
+        if zero_dopo_fine:
+            inizio_i = x_rif_iniziale + i * distanza
+            fine_i = x_rif_finale + i * distanza
+            if not _intervallo_attivo(inizio_i, fine_i, v, x_i):
+                continue
+
+        dx = x_m - x_i
+        q_x += q_picco * np.exp(-0.5 * (dx * dx) / (sigma * sigma))
+
+    return x_riferimento, q_x
+
+
+def posizioni_sorgenti_m(sorgente: dict, x_riferimento_m: float) -> np.ndarray:
+    # Posizioni lungo x di tutte le sorgenti del gruppo a partire da quella
+    # di riferimento.
+    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
+    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
+    return x_riferimento_m + np.arange(numero_sorgenti) * distanza

valida_shell.py

@@ -0,0 +1,252 @@
+# Validazione del solutore shell su casi con soluzione nota.
+#
+# Non è una suite di test automatica: stampa i risultati e i confronti perché
+# molti sono verifiche di convergenza, dove conta l'andamento dell'errore più
+# di una soglia di pass/fail. Va eseguita a mano dopo ogni modifica alla
+# fisica o all'assemblaggio:
+#
+#     python valida_shell.py
+
+import math
+
+import numpy as np
+
+import assemblaggio as asm
+import solutore as sol
+import sorgente as src
+from materials import MATERIALI
+from mesh import campo_su_griglia, costruisci_mesh_cilindrica
+
+
+MATERIALE = MATERIALI["banda_stagnata"]
+K_TERM = MATERIALE["conducibilita_termica_W_mK"]
+RHO = MATERIALE["densita_kg_m3"]
+CP = MATERIALE["calore_specifico_J_kgK"]
+ALPHA = K_TERM / (RHO * CP)
+
+LUNGHEZZA_M = 0.100
+RAGGIO_M = 0.035
+SPESSORE_M = 0.18e-3
+DT_S = 0.01
+
+
+def _esito(condizione: bool) -> str:
+    return "OK" if condizione else "FALLITO"
+
+
+def _campo_nullo(mesh: dict) -> np.ndarray:
+    return np.zeros(mesh["n_nodi"])
+
+
+def verifica_matrici(mesh: dict) -> None:
+    print("=== 1. proprietà delle matrici assemblate ===")
+    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
+    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
+
+    attesa = RHO * CP * SPESSORE_M * mesh["area_totale_m2"]
+    errore = abs(C.sum() / attesa - 1.0)
+    print(f"  capacità totale {C.sum():.6f} vs rho*cp*s*A {attesa:.6f}"
+          f"  err={errore:.1e} {_esito(errore < 1e-12)}")
+
+    # A temperatura costante il flusso deve essere nullo, quindi ogni riga
+    # della matrice di conduzione somma a zero.
+    somma_righe = np.abs(np.asarray(K.sum(axis=1))).max()
+    print(f"  max |somma righe K_cond| = {somma_righe:.1e}"
+          f" {_esito(somma_righe < 1e-9 * K_TERM)}")
+
+    asimmetria = abs(K - K.T).max()
+    print(f"  simmetria di K_cond: scarto {asimmetria:.1e}"
+          f" {_esito(asimmetria < 1e-12 * abs(K).max())}")
+
+    autovalori = np.linalg.eigvalsh(K.toarray()[:400, :400])
+    print(f"  blocco campione semidefinito positivo:"
+          f" autovalore minimo {autovalori.min():.1e}"
+          f" {_esito(autovalori.min() > -1e-9)}")
+
+
+def verifica_raffreddamento(mesh: dict) -> None:
+    # Campo uniforme senza sorgente: deve restare uniforme e decadere verso la
+    # temperatura ambiente con la costante di tempo rho*cp*s/(h_est+h_int).
+    # I bordi assiali sono esclusi, altrimenti il decadimento non è uniforme.
+    print("\n=== 2. raffreddamento convettivo uniforme ===")
+    aria = {
+        "temperatura_ambiente_C": 20.0,
+        "h_esterno_W_m2K": 12.0,
+        "h_interno_W_m2K": 8.0,
+        "h_bordi_W_m2K": 0.0,
+    }
+    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
+    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
+    K_conv, f_ambiente = asm.assembla_convezione(mesh, aria, SPESSORE_M)
+    solutore = sol.costruisci_solutore(C, K + K_conv, DT_S)
+
+    sovratemperatura = 200.0
+    T = np.full(mesh["n_nodi"], aria["temperatura_ambiente_C"] + sovratemperatura)
+    n_passi = 3000
+    for _ in range(n_passi):
+        T = sol.passo_implicito(solutore, T, f_ambiente, _campo_nullo(mesh))
+
+    t = n_passi * DT_S
+    tau = RHO * CP * SPESSORE_M / (aria["h_esterno_W_m2K"] + aria["h_interno_W_m2K"])
+    attesa = aria["temperatura_ambiente_C"] + sovratemperatura * math.exp(-t / tau)
+    errore = abs(T.mean() - attesa)
+    # Eulero implicito è del primo ordine: l'errore atteso è dell'ordine di
+    # dt/(2*tau) sulla sovratemperatura residua.
+    tolleranza = 3.0 * (attesa - aria["temperatura_ambiente_C"]) * DT_S / (2 * tau)
+
+    print(f"  tau = {tau:.3f} s, t = {t:.1f} s")
+    print(f"  numerica {T.mean():.4f} °C vs analitica {attesa:.4f} °C"
+          f"  err={errore:.4f} °C {_esito(errore < tolleranza)}")
+    dispersione = T.max() - T.min()
+    print(f"  campo ancora uniforme: max-min = {dispersione:.1e} °C"
+          f" {_esito(dispersione < 1e-9)}")
+
+
+def _sorgente_prova() -> dict:
+    return {
+        "x_inizio_m": 0.5,
+        "x_fine_m": 0.5,
+        "offset_y_percorso_m": 0.0,
+        "velocita_m_s": -1.0,
+        "numero_sorgenti": 1,
+        "distanza_sorgenti_m": 0.0,
+        "sigma_punto_m": 0.012,
+        "flusso_termico_picco_W_m2": 1e6,
+        "efficienza_riscaldamento": 1.0,
+        "zero_dopo_fine": True,
+    }
+
+
+def verifica_sorgente(mesh: dict) -> None:
+    print("\n=== 3. vettore della sorgente e conservazione dell'energia ===")
+    sorgente = _sorgente_prova()
+    assemblatore = asm.prepara_assemblatore_sorgente(mesh)
+    fattore_arco = src.fattore_circonferenziale(
+        sorgente, RAGGIO_M, assemblatore["arco_gauss_m"]
+    )
+
+    # Con la gaussiana interamente dentro la fascetta la potenza assemblata
+    # deve valere l'integrale analitico 2*pi*sigma^2*q_max; con la sorgente
+    # centrata su un bordo assiale deve valerne esattamente la metà.
+    potenza_analitica = (
+        sorgente["flusso_termico_picco_W_m2"]
+        * 2.0 * math.pi * sorgente["sigma_punto_m"] ** 2
+    )
+    # x_riferimento(t) = x_sensore + x_inizio + v*t, con x_sensore = 0.05.
+    for t, descrizione, atteso in [
+        (0.50, "centrata a metà fascetta", 1.0),
+        (0.55, "centrata sul bordo x = 0", 0.5),
+    ]:
+        _, flusso_x = src.profilo_flusso_x_W_m2(
+            sorgente, 0.05, t, assemblatore["x_gauss_m"]
+        )
+        f = asm.assembla_sorgente(assemblatore, flusso_x, fattore_arco)
+        rapporto = f.sum() / potenza_analitica
+        print(f"  {descrizione:26s}: {f.sum():9.4f} W / {potenza_analitica:.4f} W"
+              f" = {rapporto:.4f} (atteso {atteso:.1f})"
+              f" {_esito(abs(rapporto - atteso) < 1e-3)}")
+
+    # Senza convezione tutta l'energia iniettata deve finire nel campo.
+    aria = {
+        "temperatura_ambiente_C": 20.0,
+        "h_esterno_W_m2K": 0.0,
+        "h_interno_W_m2K": 0.0,
+        "h_bordi_W_m2K": 0.0,
+    }
+    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
+    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
+    _, f_ambiente = asm.assembla_convezione(mesh, aria, SPESSORE_M)
+    solutore = sol.costruisci_solutore(C, K, DT_S)
+
+    T = np.full(mesh["n_nodi"], 20.0)
+    energia_iniziale = float((C @ T).sum())
+    energia_iniettata = 0.0
+    t = 0.0
+    for _ in range(150):
+        _, flusso_x = src.profilo_flusso_x_W_m2(
+            sorgente, 0.05, t, assemblatore["x_gauss_m"]
+        )
+        f = asm.assembla_sorgente(assemblatore, flusso_x, fattore_arco)
+        energia_iniettata += f.sum() * DT_S
+        T = sol.passo_implicito(solutore, T, f_ambiente, f)
+        t += DT_S
+
+    incremento = float((C @ T).sum()) - energia_iniziale
+    errore = abs(incremento / energia_iniettata - 1.0)
+    print(f"  iniettata {energia_iniettata:.4f} J, accumulata {incremento:.4f} J"
+          f"  err={errore:.1e} {_esito(errore < 1e-10)}")
+
+
+def verifica_modi_circonferenziali(mesh: dict) -> None:
+    # Un modo sin(n*theta) su un cilindro adiabatico decade come
+    # exp(-alpha * n^2 * t / R^2): confronto diretto con l'analitico.
+    print("\n=== 4. modi sinusoidali circonferenziali ===")
+    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
+    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
+    solutore = sol.costruisci_solutore(C, K, DT_S)
+    theta = np.tile(mesh["theta_nodi_rad"], mesh["n_nodi_x"])
+
+    for n in (1, 2, 4):
+        T = 100.0 * np.sin(n * theta)
+        for _ in range(200):
+            T = sol.passo_implicito(
+                solutore, T, _campo_nullo(mesh), _campo_nullo(mesh)
+            )
+        t = 200 * DT_S
+        attesa = 100.0 * math.exp(-ALPHA * n**2 / RAGGIO_M**2 * t)
+        print(f"  n={n}: numerica {np.abs(T).max():8.4f} vs analitica {attesa:8.4f}"
+              f"  err={abs(np.abs(T).max() / attesa - 1):.3%}")
+
+
+def verifica_convergenza() -> None:
+    print("\n=== 5. convergenza al raffinamento della mesh (modo n=4) ===")
+    for n_theta in (30, 60, 110, 220):
+        mesh = costruisci_mesh_cilindrica(LUNGHEZZA_M, RAGGIO_M, 20, n_theta)
+        C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
+        K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
+        solutore = sol.costruisci_solutore(C, K, DT_S)
+        theta = np.tile(mesh["theta_nodi_rad"], mesh["n_nodi_x"])
+        T = 100.0 * np.sin(4 * theta)
+        for _ in range(200):
+            T = sol.passo_implicito(
+                solutore, T, _campo_nullo(mesh), _campo_nullo(mesh)
+            )
+        attesa = 100.0 * math.exp(-ALPHA * 16 / RAGGIO_M**2 * 200 * DT_S)
+        print(f"  n_elementi_theta={n_theta:4d} (arco {mesh['lato_arco_m']*1e3:5.2f} mm):"
+              f" err={abs(np.abs(T).max() / attesa - 1):.4%}")
+
+
+def verifica_periodicita(mesh: dict) -> None:
+    # Un impulso su theta = 0 deve diffondere in modo identico nei due versi:
+    # è la prova che la connettività riavvolge davvero la circonferenza.
+    print("\n=== 6. periodicità circonferenziale ===")
+    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
+    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
+    solutore = sol.costruisci_solutore(C, K, DT_S)
+
+    T = np.full(mesh["n_nodi"], 20.0)
+    T[50 * mesh["n_elementi_theta"]] = 500.0
+    for _ in range(50):
+        T = sol.passo_implicito(solutore, T, _campo_nullo(mesh), _campo_nullo(mesh))
+
+    griglia = campo_su_griglia(mesh, T)
+    verso_positivo = griglia[50, 1:6]
+    verso_negativo = griglia[50, -1:-6:-1]
+    scarto = np.abs(verso_positivo - verso_negativo).max()
+    print(f"  verso theta+ : {np.array2string(verso_positivo, precision=4)}")
+    print(f"  verso theta- : {np.array2string(verso_negativo, precision=4)}")
+    print(f"  scarto massimo = {scarto:.1e} {_esito(scarto < 1e-9)}")
+
+
+def main() -> None:
+    mesh = costruisci_mesh_cilindrica(LUNGHEZZA_M, RAGGIO_M, 100, 110)
+    verifica_matrici(mesh)
+    verifica_raffreddamento(mesh)
+    verifica_sorgente(mesh)
+    verifica_modi_circonferenziali(mesh)
+    verifica_convergenza()
+    verifica_periodicita(mesh)
+
+
+if __name__ == "__main__":
+    main()