temp

.gitignore

@@ -1,3 +1,3 @@
 .venv/
 __pycache__/
-output/
+dataset/

CLAUDE.md

@@ -18,39 +18,64 @@ source .venv/bin/activate
 # Installare le dipendenze
 pip install -r requirements.txt
 
-# Punto di ingresso unico: `python main.py` elenca le azioni disponibili
-python main.py mesh      # disegna la sola mesh a elementi shell
-python main.py fem       # animazione 3D del campo FEM sulla mesh shell
-python main.py csv       # esporta output/csv/run_XXXX.csv + metadata.csv
-```
+# Generare il dataset (scrive dataset/run_XXXX.csv + dataset/metadata.csv)
+python simulate.py
 
-Ogni modulo resta eseguibile anche direttamente (`python plot_mesh.py`, `python plot_animazione_fem.py`, `python esporta_csv.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
 
-Analisi termica transitoria a elementi finiti di una fascetta (anello cilindrico sottile) riscaldata da sorgenti a induzione in movimento.
+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 discretizzato è la superficie media del cilindro: x = asse (le sorgenti si muovono lungo x sulla superficie esterna), theta = coordinata circonferenziale con theta = 0 sul piano del sensore, s = R·theta la lunghezza d'arco. Lo spessore non è discretizzato. Il sensore è un pirometro a infrarossi dentro la fascetta che misura la parete in un punto x fisso sul piano theta = 0.
+**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. `main.py` — dispatcher da riga di comando che seleziona l'azione (mesh, fem, csv)
-2. `config.py` — tutti i parametri configurabili (dizionari FASCETTA, MESH, FEM, ARIA, SORGENTE, SENSORE, RANDOMIZZAZIONE, USCITA)
-3. `materials.py` — dizionario MATERIALI con proprietà termofisiche ed elettriche per materiale
-4. `mesh.py` — mesh a elementi shell quadrangolari della superficie media del cilindro, visualizzata da `plot_mesh.py`
-5. `fem.py` — solutore termico transitorio a elementi finiti sulla mesh shell, animato da `plot_animazione_fem.py` ed esportato in CSV da `esporta_csv.py`
+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. `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
 
-**Mesh shell (`mesh.py`):** griglia strutturata `n_elementi_x × n_elementi_circonferenza` di quadrilateri a 4 nodi sulla superficie media (raggio = (diametro − spessore)/2); lo spessore non è discretizzato, è un attributo della shell. La mesh è chiusa lungo la circonferenza (nessun nodo duplicato sulla cucitura) e usa coordinate globali x = asse, y = R·sin(theta), z = R·cos(theta), con theta = 0 sul piano del sensore.
+**Pipeline fisica in [simulate.py](simulate.py):**
 
-**Pipeline FEM dentro `fem.py`:** equazione risolta per unità di superficie media, `rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) − (h_est + h_int)(T − T_amb)`, con convezione `h_bordi` sui due bordi anulari (area = spessore × perimetro). Il collasso in spessore è lecito perché Biot ≈ 1e-7 e il tempo di diffusione nella parete è di pochi ms: la profondità di penetrazione del riscaldamento non entra nel modello, conta solo il flusso assorbito. La skin depth è comunque calcolata dalla frequenza dell'induttore (`calcola_skin_depth_m`) e riportata in `stato["skin_depth_m"]` come sola diagnostica. L'impronta della sorgente è una gaussiana isotropa in (x, s) con distanza circonferenziale valutata sull'immagine più vicina (superficie chiusa). Gli elementi sono bilineari a 4 nodi e tutti identici (rettangoli dx × ds), quindi le matrici di elemento (`MASSA_RIF`, `RIGIDEZZA_X`, `RIGIDEZZA_S`) sono calcolate una volta e assemblate in forma vettorizzata; la matrice di massa di superficie è riusata per capacità, convezione e carico della sorgente. Eulero implicito con matrice costante fattorizzata LU una volta per analisi (`prepara_stato_fem`), poi `passo_implicito_fem` assembla solo il termine noto. In `fem.py` stanno anche la cinematica del gruppo di sorgenti (`_x_riferimento_iniziale_m`, `_x_riferimento_finale_m`, `_intervallo_attivo`) e la lettura della configurazione (`configurazione_randomizzata`).
+- `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
 
-**Sensore e randomizzazione:** il campo nodale è la temperatura vera della parete; `simula_campo_fem` restituisce anche `T_lette`, la lettura del sensore degradata da inerzia del primo ordine (a ogni passo dt), rumore gaussiano e quantizzazione (a ogni campionamento di frame). `configurazione_randomizzata` perturba velocità, flusso di picco, sigma dello spot, offset y, temperatura ambiente e rumore del sensore con un RNG con seed `FEM["seed"]`, così due analisi restano riproducibili ma diverse; con `RANDOMIZZAZIONE["abilitata"] = False` restituisce i valori nominali di `config.py`.
+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.
 
-**Export CSV (`esporta_csv.py`):** `FEM["num_run"]` analisi (una per file `output/csv/run_XXXX.csv`), righe campionate a `FEM["frequenza_campionamento_hz"]` indipendentemente da `FEM["dt_s"]`. Colonne: `id_run, tempo_s, x_sorgente_m, T_vera_sensore_C, T_misurata_sensore_C, T_massima_fascetta_C, T_ambiente_C, offset_y_sorgente_m, velocita_m_s, sigma_punto_m, flusso_picco_W_m2, skin_depth_m, materiale`. `metadata.csv` ha una riga per analisi con tutti i parametri e i picchi. L'export chiama `simula_campo_fem(salva_campi=False)`: usa solo le serie scalari, i campi nodali non vengono accumulati. La cartella `output/csv` è cancellata e ricreata a ogni esecuzione (`shutil.rmtree`).
+`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.
+
+**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`
 
@@ -58,8 +83,11 @@ Ogni parametro in [config.py](config.py) ha un commento che spiega solo cos'è (
 
 ## Vincoli progettuali chiave
 
-- Il modello è una shell: la circonferenza è risolta, lo spessore è collassato (temperatura uniforme nella parete). Non introdurre effetti attraverso lo spessore (skin depth, gradiente radiale) senza prima discretizzarlo.
+- 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 sola (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.
-- Le visualizzazioni salvano in `USCITA["cartella"]` solo quando il backend matplotlib non è interattivo; altrimenti aprono una finestra.
+- 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,135 +1,199 @@
-# Analisi Termica FEM — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso
+# Simulatore Termico a Elementi Shell — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso
 
-Questo progetto calcola il transitorio termico di una fascetta (anello cilindrico
-sottile) riscaldata da un gruppo di sorgenti a induzione in movimento, osservata da un
-sensore a infrarossi fisso. La fascetta è discretizzata con una mesh a elementi shell
-quadrangolari e il campo di temperatura è risolto a elementi finiti, poi visualizzato in
-un'animazione 3D isometrica o esportato in CSV come serie temporale (temperatura vera
-della parete e lettura del sensore).
+Questo progetto genera misurazioni CSV pseudo-realistiche della temperatura di una
+fascetta (anello cilindrico sottile) riscaldata da un gruppo di sorgenti a induzione
+in movimento, osservata da un sensore a infrarossi fisso. Lo scopo è produrre dataset
+per l'addestramento e la validazione di modelli di stima/regressione termica.
 
 ## Geometria
 
 La fascetta è un anello cilindrico definito da tre dimensioni:
 
-- **diametro** (default 70 mm) — il diametro esterno del cilindro;
+- **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 discretizzato è la **superficie media** del cilindro (raggio =
-(diametro − spessore)/2). Lo spessore non è discretizzato: è un attributo degli
-elementi shell.
+Il dominio simulato è l'**intera superficie cilindrica**, discretizzata con elementi
+shell quadrilateri. La superficie è parametrizzata da `(x, θ)`:
 
 ```text
-        sorgenti (induttori), in moto lungo x sulla superficie esterna
-              ▼           ▼           ▼
-              ┌───────────────────────────────────┐
-        x = 0 │   superficie media del cilindro   │ x = lunghezza
-              └───────────────────────────────────┘
-                         ┆  linea di vista
-                         ▲
-                    sensore IR (fisso, x = 50 mm, theta = 0,
-                    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** = asse del cilindro, da `0` a `lunghezza`.
-- **theta** = coordinata circonferenziale, con `theta = 0` sul piano del sensore;
-  nelle formule si usa la lunghezza d'arco `s = R·theta`.
-- Coordinate globali per la visualizzazione: `x`, `y = R·sin(theta)`,
-  `z = R·cos(theta)`.
+- **x** = direzione della lunghezza, da `0` a `lunghezza`. Le sorgenti viaggiano in
+  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 parete nel nodo
-più vicino a `x = 50 mm` sul piano `theta = 0`.
+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) sulla
+generatrice `θ = 0`.
 
 ## Modello fisico
 
-L'equazione risolta è quella del calore per unità di superficie media:
+Il campo elettromagnetico non è simulato: la parte termica è però una FEM vera. La
+catena di modellazione è la seguente.
+
+1. **Shell termica** — per una parete sottile con temperatura uniforme nello spessore
+   l'equazione del calore sulla superficie cilindrica è
+
+   ```text
+   ρ cp s ∂T/∂t = k s ( ∂²T/∂x² + (1/R²) ∂²T/∂θ² ) + q'' − q_conv
+   ```
+
+   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
+
+   ```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 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
-rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) − (h_est + h_int)·(T − T_amb)
+Fo = α · (σ/|v|) / s²
 ```
 
-con `t` = spessore della parete e `∇²` il laplaciano nel piano `(x, 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.
 
-1. **Collasso dello spessore** — la temperatura è uniforme attraverso la parete. Con
-   spessore 0.18 mm il numero di Biot `h·t/k` vale ~1e-7 e il tempo di diffusione
-   attraverso la parete `t²/alpha` è di pochi millisecondi, molto più rapido del
-   transito delle sorgenti. Di conseguenza la profondità di penetrazione del
-   riscaldamento a induzione non entra nel modello: conta solo il flusso totale
-   assorbito per unità di superficie. La skin depth
-   `δ = √(2ρₑ/(ωμ))` è comunque calcolata dalla frequenza dell'induttore
-   (`frequenza_hz`, oppure imposta con `skin_depth_fissa_m`) e riportata come
-   grandezza diagnostica: per la banda stagnata a 20 kHz risulta ≈ 0.14 mm,
-   confrontabile con lo spessore, il che conferma che la parete è quasi isoterma
-   attraverso lo spessore.
+Il criterio per sapere quando servirebbe una formulazione multistrato è `Fo ≈ 1`, cioè
+uno spessore critico `s* = √(α σ/|v|)`:
 
-2. **Sorgenti gaussiane in moto** — ogni sorgente ha un'impronta gaussiana isotropa nel
-   piano `(x, s)` di raggio `sigma_punto_m`. Un gruppo di `numero_sorgenti` sorgenti
-   equidistanti (`distanza_sorgenti_m`) si muove rigidamente a velocità costante:
-   `q(x, s, t) = Σᵢ q_picco · efficienza · exp(-((x - xᵢ(t))² + (s - s₀)²) / (2σ²))`.
-   La distanza circonferenziale è valutata sull'immagine più vicina, perché la
-   superficie è chiusa. L'offset `s₀` del percorso rispetto al piano del sensore è
-   `offset_y_percorso_m`.
-
-3. **Scambi con l'ambiente** — convezione sulla faccia esterna (`h_esterno_W_m2K`) e su
-   quella interna (`h_interno_W_m2K`) su tutta la superficie, più convezione sui due
-   bordi anulari `x = 0` e `x = lunghezza` (`h_bordi_W_m2K`, su un'area pari a
-   spessore × perimetro).
-
-4. **Temperatura iniziale** — il campo parte uniformemente alla temperatura ambiente
-   dell'analisi (randomizzabile).
-
-5. **Sensore reale** — la lettura aggiunge alla temperatura vera della parete nel punto
-   osservato: inerzia del primo ordine (`costante_tempo_s`), rumore gaussiano
-   (`rumore_std_C`) e quantizzazione (`quantizzazione_C`), gli ultimi due applicati a
-   ogni campionamento.
-
-6. **Randomizzazione** — con `RANDOMIZZAZIONE["abilitata"]` ogni analisi perturba
-   velocità, flusso di picco, sigma dello spot, offset y del percorso, temperatura
-   ambiente e rumore del sensore con estrazioni gaussiane/uniformi da un RNG con seed
-   fisso (`FEM["seed"]`): analisi diverse dello stesso setup differiscono come
-   differirebbero due passaggi reali, restando riproducibili.
-
-Irraggiamento non modellato e proprietà dei materiali costanti con la temperatura: a
-~210 °C le perdite radiative non sono del tutto trascurabili rispetto alla convezione, e
-per gli acciai ferromagnetici la skin depth reale varia fortemente con temperatura e
-campo (punto di Curie non modellato).
+| 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
 
-- **Mesh** (`mesh.py`): griglia strutturata `n_elementi_x × n_elementi_circonferenza`
-  di quadrilateri a 4 nodi (default 40 × 48), chiusa lungo la circonferenza senza nodi
-  duplicati sulla cucitura.
-- **Elementi**: shell bilineari a 4 nodi. Tutti gli elementi sono rettangoli identici
-  `dx × ds`, quindi le matrici di elemento (massa e rigidezza) sono calcolate una volta
-  sola e assemblate in forma vettorizzata. La matrice di massa di superficie è riusata
-  per capacità termica, convezione sulle facce e carico della sorgente; la convezione
-  sui bordi anulari usa una matrice di massa 1D sull'anello di nodi.
-- **Integrazione temporale**: Eulero implicito con passo `FEM["dt_s"]` (default 1 ms).
-  Tutti i termini sono lineari e costanti, quindi la matrice di sistema è **assemblata e
-  fattorizzata LU una sola volta** (`scipy.sparse.linalg.splu`) in `prepara_stato_fem`;
-  ogni passo assembla solo il termine noto e risolve i sistemi triangolari.
-- **Flusso nodale**: l'impronta gaussiana isotropa è separabile,
-  `q(x, s) = q_x(x) · q_s(s)`, e il fattore circonferenziale è costante nel tempo
-  (precalcolato). A ogni passo si valuta quindi solo il fattore assiale sui `n_x + 1`
-  nodi distinti in x, espanso con un prodotto esterno.
+- 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
-main.py                 punto di ingresso da riga di comando
-config.py               tutti i parametri (geometria, mesh, FEM, aria, sorgente, sensore)
+config.py               tutti i parametri di simulazione
 materials.py            proprietà termofisiche ed elettriche dei materiali
-mesh.py                 generazione della mesh a elementi shell
-fem.py                  solutore termico transitorio a elementi finiti
-plot_mesh.py            disegno della sola mesh
-plot_animazione_fem.py  animazione 3D del campo di temperatura
-esporta_csv.py          export CSV della serie temporale
-output/                 immagini e GIF salvate quando il backend non è interattivo
-output/csv/             CSV esportati (ricreata a ogni esecuzione)
+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
@@ -143,67 +207,69 @@ pip install -r requirements.txt
 ## Uso
 
 ```bash
-python main.py          # elenco delle azioni disponibili
-python main.py mesh     # disegna la sola mesh a elementi shell
-python main.py fem      # integra il campo FEM e lo anima in 3D
-python main.py csv      # esporta la serie temporale in output/csv/
+# genera il dataset (ATTENZIONE: cancella e ricrea la cartella dataset/)
+python simulate.py
+
+# grafici statici del primo run (temperatura e flusso nel tempo)
+python plot_csv.py
+
+# animazione della superficie sviluppata durante il passaggio delle sorgenti
+python plot_animazione.py
+
+# animazione 3D isometrica del cilindro
+python plot_animazione_3d.py
 ```
 
-Ogni modulo resta eseguibile anche direttamente (`python plot_mesh.py`,
-`python plot_animazione_fem.py`, `python esporta_csv.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/`.
 
-Gli script aprono finestre interattive (backend Qt); se il backend non è interattivo
-(es. sessione senza display) salvano automaticamente PNG/GIF in `output/`.
+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`.
 
-L'animazione FEM mostra due pannelli: la fascetta in vista isometrica con il campo di
-temperatura dipinto sugli elementi shell e i marker delle sorgenti in transito, e la
-storia della temperatura nel punto osservato dal sensore (valore nodale vero e lettura
-del sensore reale). Finestra temporale iniziale e cadenza dei fotogrammi si
-regolano con le costanti in testa a `plot_animazione_fem.py`; la durata simulata è
-`FEM["durata_s"]`.
+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.
 
-## Output CSV
+### Vedere gli elementi shell
 
-`python main.py csv` cancella e ricrea `output/csv/`, poi esporta `FEM["num_run"]` analisi
-(una per file `run_XXXX.csv`) più un `metadata.csv` con una riga per analisi. Le righe
-sono campionate a `FEM["frequenza_campionamento_hz"]`, indipendente dal passo di
-integrazione.
+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:
 
-### `output/csv/run_XXXX.csv` — serie temporale
+```python
+# in plot_animazione.py e/o plot_animazione_3d.py
+MOSTRA_MESH = True
 
-| Colonna                 | Significato                                                     |
-|-------------------------|-----------------------------------------------------------------|
-| `id_run`                | identificativo dell'analisi                                     |
-| `tempo_s`               | tempo simulato                                                  |
-| `x_sorgente_m`          | posizione della sorgente di riferimento del gruppo              |
-| `T_vera_sensore_C`      | temperatura vera della parete nel punto osservato (valore nodale FEM) |
-| `T_misurata_sensore_C`  | lettura del sensore (inerzia + rumore + quantizzazione)         |
-| `T_massima_fascetta_C`  | temperatura nodale massima su tutta la fascetta in quell'istante |
-| `T_ambiente_C`          | temperatura ambiente dell'analisi                               |
-| `offset_y_sorgente_m`, `velocita_m_s`, `sigma_punto_m`, `flusso_picco_W_m2` | parametri effettivi dell'analisi (variano se la randomizzazione è attiva) |
-| `skin_depth_m`          | skin depth diagnostica                                          |
-| `materiale`             | chiave del materiale                                            |
+# in config.py, dentro FASCETTA
+"n_elementi_x": 25,
+"n_elementi_theta": 28,
+```
 
-### `output/csv/metadata.csv` — una riga per analisi
-
-Tutti i parametri effettivi (geometria, mesh, coefficienti di scambio, sorgenti, sensore)
-e i valori di picco: `T_vera_max_sensore_C`, `T_misurata_max_sensore_C`,
-`T_massima_fascetta_C`.
+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
 
 Tutto si modifica in `config.py`. I dizionari principali:
 
-| Dizionario  | Contenuto                                                                 |
-|-------------|---------------------------------------------------------------------------|
-| `FASCETTA`        | diametro, lunghezza, spessore, materiale                            |
-| `MESH`            | numero di elementi shell lungo x e lungo la circonferenza           |
-| `FEM`             | passo temporale, durata, seed, campionamento CSV, numero di analisi |
-| `ARIA`            | temperatura ambiente e coefficienti di convezione (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 analisi                              |
-| `USCITA`          | cartella per immagini e animazioni salvate                          |
+| Dizionario       | Contenuto                                                                 |
+|------------------|---------------------------------------------------------------------------|
+| `SIMULAZIONE`    | numero di run, durata, campionamento CSV, passo interno, seed, cartella   |
+| `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                                        |
 
 Punti da conoscere:
 
@@ -214,14 +280,58 @@ Punti da conoscere:
   alla fine della propria corsa.
 - **Materiale**: `FASCETTA["materiale"]` deve essere una chiave di `MATERIALI` in
   `materials.py`. Per aggiungere un materiale basta una nuova voce nel dizionario
-  (conducibilità termica, densità, calore specifico, resistività elettrica,
-  permeabilità relativa).
+  (conducibilità, densità, calore specifico, resistività elettrica, permeabilità).
 - **Proprietà di `banda_stagnata`**: la banda stagnata è un nastro di acciaio a basso
   tenore di carbonio (0,15–0,25% C) rivestito su entrambe le facce da un sottile
   strato di stagno elettrolitico, dello spessore di pochi micrometri — trascurabile
   rispetto allo spessore tipico della parete (es. 0,18 mm) e quindi ininfluente sulle
-  proprietà termiche, elettriche e magnetiche in massa. I valori in `materials.py` sono
-  quindi quelli dell'acciaio dolce sottostante:
+  proprietà termiche, elettriche e magnetiche in massa. I valori in `materials.py`
+  sono quindi quelli dell'acciaio dolce sottostante:
   [Banda stagnata: tutto quello che c'è da sapere (MUNDOLATAS)](https://mundolatas.com/it/banda-stagnata-tutto-quello-che-ce-da-sapere/),
   [Bande stagnate elettrolitiche (EUROPERF)](https://www.europerf.it/it/banda-stagnata-elettrolitica.php),
   [Differenza tra banda stagnata e acciaio inossidabile (Wuxi Bright Packing)](https://it.brightmetalcan.com/info/difference-between-tinplate-and-stainless-stee-48700260.html).
+- **Randomizzazione**: ogni run perturba velocità, flusso di picco, sigma, offset y,
+  temperatura ambiente e rumore del sensore con estrazioni da un RNG a seed fisso
+  (`SIMULAZIONE["seed"]`): il dataset è riproducibile.
+
+## Output
+
+### `dataset/run_XXXX.csv` — serie temporale del run
+
+| Colonna                        | Significato                                                        |
+|--------------------------------|--------------------------------------------------------------------|
+| `id_run`                       | identificativo del run                                             |
+| `tempo_s`                      | tempo simulato                                                     |
+| `x_sorgente_m`                 | posizione della sorgente di riferimento del gruppo                 |
+| `offset_y_sorgente_m`          | offset circonferenziale del percorso (costante nel run)            |
+| `flusso_termico_sorgente_W_m2` | flusso efficace nel punto x del sensore                            |
+| `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_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, 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
+   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 ~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,21 +1,29 @@
-# Configurazione dell'analisi termica a elementi finiti della 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 discretizzato è la superficie media del cilindro (mesh a
-#   elementi shell quadrangolari, vedi mesh.py): x = asse della fascetta,
-#   s = R·theta = coordinata circonferenziale. Lo spessore non è
-#   discretizzato, la temperatura è uniforme attraverso la parete.
-# - Le sorgenti a induzione si muovono lungo x sulla superficie esterna, con
-#   un'impronta gaussiana isotropa nel piano (x, s).
-# - La superficie scambia per convezione con l'aria sulla faccia esterna e su
-#   quella interna, più i due bordi anulari x = 0 e x = lunghezza.
+# - 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
-#   contatto la temperatura della superficie interna nel punto x = "x_mm",
-#   sul piano circonferenziale theta = 0.
+#   contatto la temperatura della superficie interna nel punto x = "x_mm".
 #
 # Unità di misura:
 # - lunghezza: m (mm dove indicato dal suffisso)
@@ -24,9 +32,36 @@
 # - flusso termico: W/m²
 # - coefficiente di convezione: W/(m² K)
 
+SIMULAZIONE = {
+    # Numero di file CSV da generare.
+    "num_run": 1,
+
+    # Tempo simulato totale.
+    "durata_s": 30.0,
+
+    # Frequenza di campionamento CSV.
+    # Esempio: 2 Hz significa una riga ogni 0.5 s.
+    "frequenza_campionamento_hz": 10.0,
+
+    # Passo di integrazione numerica interna.
+    # Può essere inferiore al periodo di campionamento CSV.
+    "dt_interno_s": 0.001,
+
+    # Seed per la riproducibilità.
+    "seed": 42,
+
+    # Cartella di output.
+    "cartella_output": "dataset",
+
+    # Numero di processi paralleli per la generazione dei run.
+    # None = usa tutti i core disponibili.
+    "num_processi": None,
+}
+
 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,
@@ -34,36 +69,15 @@ FASCETTA = {
     # Spessore della parete [mm].
     "spessore_mm": 0.18,
 
-    # Deve corrispondere a una chiave in materials.py.
-    "materiale": "banda_stagnata",
-}
-
-MESH = {
-    # Numero di elementi shell lungo l'asse x (lunghezza della fascetta).
-    "n_elementi_x": 40,
+    # Numero di elementi shell lungo x (lunghezza).
+    "n_elementi_x": 100,
 
     # Numero di elementi shell lungo la circonferenza.
-    # La mesh è chiusa su se stessa: non c'è una riga di nodi duplicata.
-    "n_elementi_circonferenza": 48,
-}
+    # Più elementi = maggiore risoluzione spaziale, simulazione più lenta.
+    "n_elementi_theta": 110,
 
-FEM = {
-    # Passo di integrazione temporale del solutore a elementi finiti.
-    "dt_s": 0.001,
-
-    # Tempo simulato totale dell'analisi transitoria.
-    "durata_s": 30.0,
-
-    # Seed per la riproducibilità di randomizzazione e rumore del sensore.
-    "seed": 42,
-
-    # Frequenza di campionamento delle righe del CSV esportato.
-    # Esempio: 10 Hz significa una riga ogni 0.1 s.
-    "frequenza_campionamento_hz": 10.0,
-
-    # Numero di analisi da esportare in CSV, una per file.
-    # Con la randomizzazione abilitata ogni analisi ha parametri diversi.
-    "num_run": 1,
+    # Deve corrispondere a una chiave in materials.py.
+    "materiale": "banda_stagnata",
 }
 
 ARIA = {
@@ -76,7 +90,8 @@ ARIA = {
     # Coefficiente di convezione sulla faccia interna (lato sensore).
     "h_interno_W_m2K": 8.0,
 
-    # Coefficiente di convezione sui bordi anulari (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,
 }
 
@@ -89,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.
@@ -113,8 +129,10 @@ SORGENTE = {
     # Frazione del flusso incidente che diventa effettivamente calore nella fascetta.
     "efficienza_riscaldamento": 0.35,
 
-    # Frequenza di induzione dell'induttore, usata per stimare la skin depth
-    # se skin_depth_fissa_m è None.
+    # 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.
@@ -127,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 lo spessore [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,
@@ -148,7 +168,7 @@ SENSORE = {
 }
 
 RANDOMIZZAZIONE = {
-    # Se abilitata, ogni analisi varia leggermente alcuni parametri.
+    # Se abilitata, ogni run varia leggermente alcuni parametri.
     "abilitata": False,
 
     # Deviazioni standard relative.
@@ -165,9 +185,3 @@ RANDOMIZZAZIONE = {
     # dalla linea ideale allineata con il sensore.
     "offset_y_max_assoluto_m": 0.001,
 }
-
-USCITA = {
-    # Cartella in cui salvare immagini e animazioni quando il backend
-    # matplotlib non è interattivo.
-    "cartella": "output",
-}

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

esporta_csv.py

@@ -1,153 +0,0 @@
-# Export CSV della serie temporale calcolata dal solutore FEM (fem.py).
-#
-# Per ogni analisi viene scritto un file run_XXXX.csv con la temperatura vera
-# della parete nel punto osservato dal sensore e la lettura del sensore reale
-# (inerzia, rumore, quantizzazione), più i parametri effettivi dell'analisi.
-# `metadata.csv` riassume un'analisi per riga.
-#
-# Le righe sono campionate a FEM["frequenza_campionamento_hz"], indipendente
-# dal passo di integrazione FEM["dt_s"]. I campi nodali non vengono accumulati:
-# servono solo le serie scalari.
-
-import csv
-import random
-import shutil
-from pathlib import Path
-
-from config import FEM, USCITA
-from fem import configurazione_randomizzata, simula_campo_fem
-from mesh import genera_mesh, riepilogo_mesh
-
-INTESTAZIONE = [
-    "id_run",
-    "tempo_s",
-    "x_sorgente_m",
-    "T_vera_sensore_C",
-    "T_misurata_sensore_C",
-    "T_massima_fascetta_C",
-    "T_ambiente_C",
-    "offset_y_sorgente_m",
-    "velocita_m_s",
-    "sigma_punto_m",
-    "flusso_picco_W_m2",
-    "skin_depth_m",
-    "materiale",
-]
-
-
-def esporta_run(
-    id_run: str,
-    percorso_csv: Path,
-    mesh_dati: dict,
-    rng: random.Random,
-) -> dict:
-    """Integra un'analisi FEM, ne scrive il CSV e restituisce i suoi metadati."""
-    cfg = configurazione_randomizzata(rng)
-    dati = simula_campo_fem(
-        dt_frame_s=1.0 / FEM["frequenza_campionamento_hz"],
-        cfg=cfg,
-        mesh_dati=mesh_dati,
-        rng=rng,
-        salva_campi=False,
-    )
-
-    fascetta = cfg["fascetta"]
-    aria = cfg["aria"]
-    sorgente = cfg["sorgente"]
-    sensore = cfg["sensore"]
-    skin_depth = dati["skin_depth_m"]
-
-    percorso_csv.parent.mkdir(parents=True, exist_ok=True)
-    with percorso_csv.open("w", newline="") as f:
-        writer = csv.writer(f)
-        writer.writerow(INTESTAZIONE)
-        for k, t in enumerate(dati["tempi"]):
-            writer.writerow([
-                id_run,
-                f"{t:.6f}",
-                f"{dati['x_riferimenti'][k]:.9f}",
-                f"{dati['T_vere'][k]:.6f}",
-                f"{dati['T_lette'][k]:.6f}",
-                f"{dati['T_massime'][k]:.6f}",
-                f"{aria['temperatura_ambiente_C']:.6f}",
-                f"{sorgente['offset_y_percorso_m']:.9f}",
-                f"{sorgente['velocita_m_s']:.9f}",
-                f"{sorgente['sigma_punto_m']:.9f}",
-                f"{sorgente['flusso_termico_picco_W_m2']:.6f}",
-                f"{skin_depth:.9e}",
-                fascetta["materiale"],
-            ])
-
-    return {
-        "id_run": id_run,
-        "file_csv": percorso_csv.name,
-        "materiale": fascetta["materiale"],
-        "diametro_m": fascetta["diametro_mm"] / 1000.0,
-        "lunghezza_m": fascetta["lunghezza_mm"] / 1000.0,
-        "spessore_m": fascetta["spessore_mm"] / 1000.0,
-        "n_elementi_x": mesh_dati["n_elementi_x"],
-        "n_elementi_circonferenza": mesh_dati["n_elementi_circonferenza"],
-        "durata_s": FEM["durata_s"],
-        "dt_s": FEM["dt_s"],
-        "frequenza_campionamento_hz": FEM["frequenza_campionamento_hz"],
-        "temperatura_ambiente_C": aria["temperatura_ambiente_C"],
-        "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"],
-        "x_inizio_m": sorgente["x_inizio_m"],
-        "x_fine_m": sorgente["x_fine_m"],
-        "x_sensore_m": dati["stato"]["x_sensore_m"],
-        "distanza_sensore_parete_m": sensore["distanza_parete_mm"] / 1000.0,
-        "offset_y_percorso_m": sorgente["offset_y_percorso_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),
-        "sigma_punto_m": sorgente["sigma_punto_m"],
-        "flusso_termico_picco_W_m2": sorgente["flusso_termico_picco_W_m2"],
-        "efficienza_riscaldamento": sorgente["efficienza_riscaldamento"],
-        "frequenza_hz": sorgente["frequenza_hz"],
-        "skin_depth_m": skin_depth,
-        "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_sensore_C": float(dati["T_vere"].max()),
-        "T_misurata_max_sensore_C": float(dati["T_lette"].max()),
-        "T_massima_fascetta_C": float(dati["T_massime"].max()),
-    }
-
-
-def main() -> None:
-    mesh_dati = genera_mesh()
-    print(riepilogo_mesh(mesh_dati))
-
-    # La cartella dei CSV è ricreata da zero a ogni esecuzione, così non
-    # restano run di esecuzioni precedenti con parametri diversi.
-    cartella = Path(USCITA["cartella"]) / "csv"
-    if cartella.exists():
-        shutil.rmtree(cartella)
-    cartella.mkdir(parents=True, exist_ok=True)
-
-    rng = random.Random(FEM["seed"])
-    righe_metadata = []
-    for indice in range(1, FEM["num_run"] + 1):
-        id_run = f"run_{indice:04d}"
-        print(f"Integrazione FEM di {id_run} in corso...")
-        riga = esporta_run(id_run, cartella / f"{id_run}.csv", mesh_dati, rng)
-        print(
-            f"  T vera max al sensore: {riga['T_vera_max_sensore_C']:.1f} °C — "
-            f"T misurata max: {riga['T_misurata_max_sensore_C']:.1f} °C"
-        )
-        righe_metadata.append(riga)
-
-    percorso_metadata = cartella / "metadata.csv"
-    with percorso_metadata.open("w", newline="") as f:
-        writer = csv.DictWriter(f, fieldnames=list(righe_metadata[0].keys()))
-        writer.writeheader()
-        writer.writerows(righe_metadata)
-
-    print(f"Esportate {len(righe_metadata)} analisi in: {cartella.resolve()}")
-    print(f"Metadata: {percorso_metadata.resolve()}")
-
-
-if __name__ == "__main__":
-    main()

fem.py

@@ -1,460 +0,0 @@
-# Solutore termico transitorio a elementi finiti sulla mesh shell della
-# fascetta (vedi mesh.py per la geometria).
-#
-# Il dominio è la superficie media del cilindro, con x e la coordinata
-# circonferenziale s = R·theta entrambe risolte, mentre lo spessore è
-# collassato (temperatura uniforme attraverso la parete).
-#
-# Il collasso in spessore è lecito: con spessore 0.18 mm il numero di Biot
-# h·t/k vale ~1e-7 e il tempo di diffusione attraverso la parete t²/alpha è
-# di pochi millisecondi, molto più rapido del transito delle sorgenti. Di
-# conseguenza la profondità di penetrazione del riscaldamento a induzione non
-# entra nel modello: conta solo il flusso totale assorbito per unità di
-# superficie. La skin depth è comunque calcolata dalla frequenza dell'induttore
-# come grandezza diagnostica (`calcola_skin_depth_m`).
-#
-# Equazione risolta (per unità di superficie media):
-#
-#   rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) - (h_est + h_int)·(T - T_amb)
-#
-# con convezione aggiuntiva sui due bordi anulari x = 0 e x = L (coefficiente
-# h_bordi su un'area pari a spessore × perimetro). L'impronta della sorgente
-# è una gaussiana isotropa nel piano (x, s), con la distanza circonferenziale
-# valutata sull'immagine più vicina perché la superficie è chiusa.
-#
-# Discretizzazione: elementi shell quadrangolari bilineari a 4 nodi. Gli
-# elementi della mesh sono tutti rettangoli identici dx × ds, quindi le
-# matrici di elemento sono calcolate una volta sola e assemblate in forma
-# vettorizzata. Integrazione temporale con Eulero implicito: la matrice di
-# sistema è costante e viene fattorizzata LU una volta per run.
-
-import math
-import random
-from copy import deepcopy
-
-import numpy as np
-import scipy.sparse as sp
-from scipy.sparse.linalg import splu
-
-from config import ARIA, FASCETTA, FEM, RANDOMIZZAZIONE, SENSORE, SORGENTE
-from materials import MATERIALI
-from mesh import genera_mesh
-
-MU0 = 4.0 * math.pi * 1e-7
-
-# Matrici di riferimento dell'elemento rettangolare bilineare a 4 nodi, con
-# nodi locali in ordine antiorario: (0,0), (a,0), (a,b), (0,b).
-#
-# MASSA_RIF va moltiplicata per l'area a·b e dà l'integrale di N_i·N_j;
-# RIGIDEZZA_X per b/a e RIGIDEZZA_S per a/b, e la loro somma dà l'integrale
-# di grad(N_i)·grad(N_j).
-MASSA_RIF = np.array(
-    [[4.0, 2.0, 1.0, 2.0],
-     [2.0, 4.0, 2.0, 1.0],
-     [1.0, 2.0, 4.0, 2.0],
-     [2.0, 1.0, 2.0, 4.0]]
-) / 36.0
-
-RIGIDEZZA_X = np.array(
-    [[2.0, -2.0, -1.0, 1.0],
-     [-2.0, 2.0, 1.0, -1.0],
-     [-1.0, 1.0, 2.0, -2.0],
-     [1.0, -1.0, -2.0, 2.0]]
-) / 6.0
-
-RIGIDEZZA_S = np.array(
-    [[2.0, 1.0, -1.0, -2.0],
-     [1.0, 2.0, -2.0, -1.0],
-     [-1.0, -2.0, 2.0, 1.0],
-     [-2.0, -1.0, 1.0, 2.0]]
-) / 6.0
-
-
-def calcola_skin_depth_m(materiale: dict, frequenza_hz: float) -> float:
-    # Skin depth elettromagnetica approssimata:
-    # delta = sqrt(2 * rho_e / (omega * mu))
-    #
-    # Semplificata. Per acciai ferromagnetici il comportamento reale
-    # è fortemente non lineare con temperatura e campo magnetico.
-    #
-    # Nel modello shell è una grandezza diagnostica: serve a verificare che la
-    # deposizione del calore resti confinata entro uno spessore confrontabile
-    # con la parete, non entra nell'equazione risolta.
-    rho_e = materiale["resistivita_elettrica_ohm_m"]
-    mu_r = materiale["permeabilita_relativa"]
-    omega = 2.0 * math.pi * frequenza_hz
-    mu = MU0 * mu_r
-    return math.sqrt(2.0 * rho_e / (omega * mu))
-
-
-def quantizza(valore: float, passo: float) -> float:
-    if passo <= 0.0:
-        return valore
-    return round(valore / passo) * passo
-
-
-def configurazione_randomizzata(rng: random.Random | None = None) -> dict:
-    """Copia dei dizionari di configurazione usati da un'analisi FEM.
-
-    Se `RANDOMIZZAZIONE["abilitata"]`, velocità, flusso di picco, sigma dello
-    spot, offset y del percorso, temperatura ambiente e rumore del sensore sono
-    perturbati con estrazioni dall'RNG passato: analisi diverse dello stesso
-    setup differiscono come differirebbero due passaggi reali.
-    """
-    fascetta = deepcopy(FASCETTA)
-    aria = deepcopy(ARIA)
-    sorgente = deepcopy(SORGENTE)
-    sensore = deepcopy(SENSORE)
-
-    if RANDOMIZZAZIONE.get("abilitata", False):
-        if rng is None:
-            rng = random.Random(FEM["seed"])
-
-        def perturba_rel(valore: float, std_rel: float, fattore_min: float = 0.1) -> float:
-            fattore = rng.gauss(1.0, std_rel)
-            fattore = max(fattore_min, fattore)
-            return valore * fattore
-
-        sorgente["velocita_m_s"] = perturba_rel(
-            sorgente["velocita_m_s"],
-            RANDOMIZZAZIONE["velocita_std_rel"],
-        )
-        sorgente["flusso_termico_picco_W_m2"] = perturba_rel(
-            sorgente["flusso_termico_picco_W_m2"],
-            RANDOMIZZAZIONE["flusso_picco_std_rel"],
-        )
-        sorgente["sigma_punto_m"] = perturba_rel(
-            sorgente["sigma_punto_m"],
-            RANDOMIZZAZIONE["sigma_punto_std_rel"],
-        )
-        sorgente["offset_y_percorso_m"] = rng.uniform(
-            -RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
-            RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
-        )
-        aria["temperatura_ambiente_C"] += rng.gauss(
-            0.0,
-            RANDOMIZZAZIONE["temperatura_ambiente_std_C"],
-        )
-        sensore["rumore_std_C"] = perturba_rel(
-            sensore["rumore_std_C"],
-            RANDOMIZZAZIONE["rumore_sensore_std_rel"],
-            fattore_min=0.0,
-        )
-
-    return {
-        "fascetta": fascetta,
-        "aria": aria,
-        "sorgente": sorgente,
-        "sensore": sensore,
-    }
-
-
-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 _assembla(elementi: np.ndarray, matrice_elemento: np.ndarray, n_nodi: int):
-    # Assembla una matrice globale sparsa a partire da un'unica matrice di
-    # elemento 4×4, uguale per tutti gli elementi. I contributi ripetuti sullo
-    # stesso (riga, colonna) vengono sommati dal formato COO.
-    n_elementi = len(elementi)
-    righe = np.broadcast_to(elementi[:, :, None], (n_elementi, 4, 4)).ravel()
-    colonne = np.broadcast_to(elementi[:, None, :], (n_elementi, 4, 4)).ravel()
-    valori = np.broadcast_to(matrice_elemento, (n_elementi, 4, 4)).ravel()
-    return sp.coo_matrix(
-        (valori, (righe, colonne)), shape=(n_nodi, n_nodi)
-    ).tocsr()
-
-
-def _massa_anello(indici: np.ndarray, lunghezza_segmento_m: float, n_nodi: int):
-    # Matrice di massa 1D su un anello chiuso di nodi (integrale di N_i·N_j
-    # lungo la circonferenza), usata per la convezione sui bordi x = 0 e x = L.
-    n = len(indici)
-    successivi = indici[(np.arange(n) + 1) % n]
-    coefficiente = lunghezza_segmento_m / 6.0
-    righe = np.concatenate([indici, indici, successivi, successivi])
-    colonne = np.concatenate([indici, successivi, indici, successivi])
-    valori = coefficiente * np.concatenate(
-        [np.full(n, 2.0), np.full(n, 1.0), np.full(n, 1.0), np.full(n, 2.0)]
-    )
-    return sp.coo_matrix(
-        (valori, (righe, colonne)), shape=(n_nodi, n_nodi)
-    ).tocsr()
-
-
-def prepara_stato_fem(
-    fascetta: dict,
-    aria: dict,
-    sorgente: dict,
-    sensore: dict,
-    mesh_dati: dict | None = None,
-    dt_s: float | None = None,
-) -> dict:
-    """Assembla le matrici FEM e fattorizza il sistema implicito.
-
-    Restituisce lo stato costante per l'integrazione temporale: matrici,
-    solutore LU, coordinate nodali, peso circonferenziale della sorgente e
-    indice del nodo osservato dal sensore.
-    """
-    if mesh_dati is None:
-        mesh_dati = genera_mesh(fascetta)
-    if dt_s is None:
-        dt_s = FEM["dt_s"]
-
-    materiale = MATERIALI[fascetta["materiale"]]
-    if sorgente["skin_depth_fissa_m"] is None:
-        skin_depth_m = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
-    else:
-        skin_depth_m = float(sorgente["skin_depth_fissa_m"])
-    k = materiale["conducibilita_termica_W_mK"]
-    rho = materiale["densita_kg_m3"]
-    cp = materiale["calore_specifico_J_kgK"]
-
-    spessore_m = mesh_dati["spessore_m"]
-    raggio_m = mesh_dati["raggio_m"]
-    circonferenza_m = 2.0 * math.pi * raggio_m
-    n_x = mesh_dati["n_elementi_x"]
-    n_theta = mesh_dati["n_elementi_circonferenza"]
-    elementi = mesh_dati["elementi"]
-    indice_nodo = mesh_dati["indice_nodo"]
-    n_nodi = len(mesh_dati["nodi"])
-
-    dx_m = mesh_dati["lunghezza_m"] / n_x
-    ds_m = circonferenza_m / n_theta
-
-    h_esterno = aria["h_esterno_W_m2K"]
-    h_interno = aria["h_interno_W_m2K"]
-    h_bordi = aria["h_bordi_W_m2K"]
-    T_ambiente = aria["temperatura_ambiente_C"]
-
-    # Matrice di massa di superficie: riusata per la capacità termica, per la
-    # convezione sulle due facce e per il carico della sorgente.
-    M_area = _assembla(elementi, MASSA_RIF * (dx_m * ds_m), n_nodi)
-    K_cond = _assembla(
-        elementi,
-        k * spessore_m * (RIGIDEZZA_X * (ds_m / dx_m) + RIGIDEZZA_S * (dx_m / ds_m)),
-        n_nodi,
-    )
-
-    # Convezione sui bordi anulari: l'area di scambio di ogni segmento è
-    # spessore × lunghezza del segmento circonferenziale.
-    H_bordi = h_bordi * spessore_m * (
-        _massa_anello(indice_nodo[0, :], ds_m, n_nodi)
-        + _massa_anello(indice_nodo[-1, :], ds_m, n_nodi)
-    )
-
-    C_su_dt = (rho * cp * spessore_m / dt_s) * M_area
-    A = C_su_dt + K_cond + (h_esterno + h_interno) * M_area + H_bordi
-
-    # Termine noto costante della convezione: (somma delle righe) × h × T_amb,
-    # perché la somma delle righe della matrice di massa è l'integrale di N_i.
-    area_nodale = np.asarray(M_area.sum(axis=1)).ravel()
-    bordo_nodale = np.asarray(H_bordi.sum(axis=1)).ravel()
-    carico_convezione = (h_esterno + h_interno) * area_nodale * T_ambiente
-    carico_convezione += bordo_nodale * T_ambiente
-
-    # Peso circonferenziale dell'impronta della sorgente: dipende solo
-    # dall'offset y del percorso e dalla sigma dello spot, entrambi costanti
-    # durante il run, quindi è valutato una volta sola qui. La distanza è
-    # quella dell'immagine più vicina, perché la superficie è chiusa.
-    sigma_m = sorgente["sigma_punto_m"]
-    ds_sorgente = raggio_m * mesh_dati["theta_nodi_rad"] - sorgente["offset_y_percorso_m"]
-    ds_sorgente -= circonferenza_m * np.round(ds_sorgente / circonferenza_m)
-    peso_circonferenziale = np.exp(-0.5 * (ds_sorgente ** 2) / (sigma_m * sigma_m))
-
-    # Nodo osservato dal sensore: x più vicino a x_mm sul piano theta = 0.
-    x_sensore_m = sensore["x_mm"] / 1000.0
-    i_sensore = int(np.argmin(np.abs(mesh_dati["x_nodi_m"] - x_sensore_m)))
-    indice_sensore = int(indice_nodo[i_sensore, 0])
-
-    return {
-        "mesh": mesh_dati,
-        "dt_s": dt_s,
-        "M_area": M_area,
-        "C_su_dt": C_su_dt,
-        "carico_convezione": carico_convezione,
-        "solutore": splu(sp.csc_matrix(A)),
-        "x_nodi_m": mesh_dati["x_nodi_m"],
-        "peso_circonferenziale": peso_circonferenziale,
-        "forma_griglia": (n_x + 1, n_theta),
-        "circonferenza_m": circonferenza_m,
-        "area_nodale_m2": area_nodale,
-        "T_ambiente_C": T_ambiente,
-        "x_sensore_m": x_sensore_m,
-        "indice_sensore": indice_sensore,
-        "n_nodi": n_nodi,
-        # Diagnostica: non entra nell'equazione risolta.
-        "skin_depth_m": skin_depth_m,
-    }
-
-
-def flusso_nodale_W_m2(
-    sorgente: dict,
-    stato: dict,
-    t_s: float,
-) -> tuple[float, np.ndarray]:
-    """Flusso termico assorbito nei nodi all'istante t_s.
-
-    Restituisce (x_riferimento_m, q_nodale) dove x_riferimento_m è la
-    posizione della sorgente di indice 0. Ogni sorgente attiva contribuisce
-    con un'impronta gaussiana isotropa in (x, s).
-
-    L'impronta gaussiana isotropa è separabile, q(x, s) = q_x(x) · q_s(s), e
-    il fattore circonferenziale q_s è costante nel tempo (precalcolato in
-    `prepara_stato_fem`). Qui si valuta quindi solo il fattore assiale sui
-    n_x + 1 nodi distinti in x e si espande con un prodotto esterno, invece di
-    valutare l'esponenziale su tutti i nodi della mesh.
-    """
-    x_sensore_m = stato["x_sensore_m"]
-    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"]
-
-    peso_assiale = np.zeros(len(stato["x_nodi_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 = stato["x_nodi_m"] - x_i
-        peso_assiale += np.exp(-0.5 * (dx * dx) / (sigma * sigma))
-
-    q_nodale = q_picco * np.outer(peso_assiale, stato["peso_circonferenziale"])
-    return x_riferimento, q_nodale.ravel()
-
-
-def passo_implicito_fem(stato: dict, T: np.ndarray, q_nodale: np.ndarray) -> np.ndarray:
-    # Avanza il campo nodale di un passo dt: termine noto = capacità sul campo
-    # precedente + carico della sorgente + carico di convezione, poi risoluzione
-    # del sistema già fattorizzato.
-    rhs = stato["C_su_dt"] @ T
-    rhs += stato["M_area"] @ q_nodale
-    rhs += stato["carico_convezione"]
-    return stato["solutore"].solve(rhs)
-
-
-def campo_iniziale(stato: dict) -> np.ndarray:
-    return np.full(stato["n_nodi"], stato["T_ambiente_C"], dtype=float)
-
-
-def simula_campo_fem(
-    dt_frame_s: float,
-    cfg: dict | None = None,
-    durata_s: float | None = None,
-    mesh_dati: dict | None = None,
-    rng: random.Random | None = None,
-    salva_campi: bool = True,
-) -> dict:
-    """Integra il campo FEM fino a durata_s, campionando ogni dt_frame_s.
-
-    Il campo nodale è la temperatura vera della parete; la lettura del sensore
-    ne è la versione degradata: inerzia del primo ordine, rumore gaussiano e
-    quantizzazione, questi ultimi due applicati a ogni campionamento.
-
-    Con `salva_campi = False` i campi nodali non vengono accumulati (utile per
-    l'export CSV, che usa solo le serie scalari): `campi` resta una lista vuota.
-    """
-    if rng is None:
-        rng = random.Random(FEM["seed"])
-    if cfg is None:
-        cfg = configurazione_randomizzata(rng)
-    if durata_s is None:
-        durata_s = FEM["durata_s"]
-
-    fascetta = cfg["fascetta"]
-    aria = cfg["aria"]
-    sorgente = cfg["sorgente"]
-    sensore = cfg["sensore"]
-
-    stato = prepara_stato_fem(fascetta, aria, sorgente, sensore, mesh_dati=mesh_dati)
-    dt = stato["dt_s"]
-    indice_sensore = stato["indice_sensore"]
-
-    T = campo_iniziale(stato)
-    T_sensore = float(T[indice_sensore])
-    tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
-
-    tempi, campi, x_riferimenti = [], [], []
-    T_vere, T_lette, T_massime = [], [], []
-
-    prossimo_frame_t = 0.0
-    t = 0.0
-    while t <= durata_s + 1e-12:
-        x_rif, q_nodale = flusso_nodale_W_m2(sorgente, stato, t)
-        T = passo_implicito_fem(stato, T, q_nodale)
-
-        T_sensore += (T[indice_sensore] - T_sensore) * dt / tau_sensore
-
-        if t + 1e-12 >= prossimo_frame_t:
-            letta = T_sensore + rng.gauss(0.0, sensore["rumore_std_C"])
-            letta = quantizza(letta, sensore["quantizzazione_C"])
-
-            tempi.append(t)
-            if salva_campi:
-                campi.append(T.copy())
-            x_riferimenti.append(x_rif)
-            T_vere.append(float(T[indice_sensore]))
-            T_lette.append(letta)
-            T_massime.append(float(T.max()))
-            prossimo_frame_t += dt_frame_s
-
-        t += dt
-
-    return {
-        "stato": stato,
-        "mesh": stato["mesh"],
-        "tempi": np.array(tempi),
-        "campi": campi,
-        "x_riferimenti": np.array(x_riferimenti),
-        "T_vere": np.array(T_vere),
-        "T_lette": np.array(T_lette),
-        "T_massime": np.array(T_massime),
-        "cfg": cfg,
-        "sorgente": sorgente,
-        "T_ambiente_C": stato["T_ambiente_C"],
-        "skin_depth_m": stato["skin_depth_m"],
-    }

main.py

@@ -1,61 +0,0 @@
-# Punto di ingresso unico del progetto: seleziona da riga di comando quale
-# azione eseguire.
-#
-#   python main.py mesh          # disegna la sola mesh a elementi shell
-#   python main.py fem           # animazione 3D del campo FEM sulla mesh shell
-#   python main.py csv           # esporta in CSV la serie temporale del FEM
-#
-# Senza argomenti stampa l'elenco delle azioni disponibili.
-
-import argparse
-import sys
-
-
-def _azione_mesh() -> None:
-    import plot_mesh
-
-    plot_mesh.main()
-
-
-def _azione_fem() -> None:
-    import plot_animazione_fem
-
-    plot_animazione_fem.main()
-
-
-def _azione_csv() -> None:
-    import esporta_csv
-
-    esporta_csv.main()
-
-
-# Chiave da riga di comando -> (funzione, descrizione mostrata nell'help).
-AZIONI = {
-    "mesh": (_azione_mesh, "Disegna la mesh a elementi shell della fascetta"),
-    "fem": (_azione_fem, "Animazione 3D del campo FEM sulla mesh a elementi shell"),
-    "csv": (_azione_csv, "Esporta in CSV la temperatura simulata e la lettura del sensore"),
-}
-
-
-def main() -> None:
-    parser = argparse.ArgumentParser(
-        description="Analisi termica FEM della fascetta riscaldata a induzione.",
-        formatter_class=argparse.RawTextHelpFormatter,
-    )
-    parser.add_argument(
-        "azione",
-        nargs="?",
-        choices=list(AZIONI),
-        help="\n".join(f"{nome}: {descrizione}" for nome, (_, descrizione) in AZIONI.items()),
-    )
-    argomenti = parser.parse_args()
-
-    if argomenti.azione is None:
-        parser.print_help()
-        sys.exit(1)
-
-    AZIONI[argomenti.azione][0]()
-
-
-if __name__ == "__main__":
-    main()

materials.py

@@ -1,4 +1,4 @@
-# Database dei materiali per l'analisi termica.
+# Database dei materiali per il simulatore termico.
 #
 # Tutte le unità sono SI:
 # - conducibilita_termica_W_mK

mesh.py

@@ -1,114 +1,119 @@
-# Generazione della mesh a elementi shell quadrangolari della fascetta.
+# Generazione della mesh a elementi shell quadrilateri sulla superficie
+# cilindrica della fascetta.
 #
-# La fascetta è modellata come la superficie media di un cilindro: gli
-# elementi sono quadrilateri a 4 nodi disposti su una griglia strutturata
-# n_elementi_x × n_elementi_circonferenza. Lo spessore della parete non è
-# discretizzato, è un attributo degli elementi shell.
+# 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.
 #
-# Sistema di coordinate globale (coerente con plot_animazione_3d.py):
-# - x = asse del cilindro, da 0 a lunghezza
-# - y, z = piano della sezione circolare, con y = R·sin(theta), z = R·cos(theta)
-# - theta = 0 corrisponde al piano del sensore (y = 0), cresce in verso
-#   antiorario nel piano y-z
+# 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.
 #
-# La mesh è chiusa lungo la circonferenza: l'ultima colonna di elementi
-# richiude sulla prima riga di nodi, senza nodi duplicati sulla cucitura.
+# Ordine dei nodi locali di ogni elemento (antiorario nel piano x-arco):
+#
+#     4 ---- 3      arco
+#     |      |       ^
+#     |      |       |
+#     1 ---- 2       +---> x
 
 import numpy as np
 
-from config import FASCETTA, MESH
 
+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
 
-def raggio_medio_m(fascetta: dict = FASCETTA) -> float:
-    # Raggio della superficie media della parete: il diametro configurato è
-    # quello esterno, la shell sta a metà dello spessore.
-    diametro_m = fascetta["diametro_mm"] / 1000.0
-    spessore_m = fascetta["spessore_mm"] / 1000.0
-    return (diametro_m - spessore_m) / 2.0
+    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
 
-def genera_mesh(fascetta: dict = FASCETTA, mesh: dict = MESH) -> dict:
-    """Costruisce la mesh shell quadrangolare della fascetta.
+    # 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
 
-    Restituisce un dizionario con:
-    - nodi: array (n_nodi, 3) con le coordinate globali [m]
-    - elementi: array (n_elementi, 4) con gli indici dei nodi, in verso
-      antiorario visto dall'esterno del cilindro
-    - indice_nodo: array (n_x + 1, n_theta) che mappa (i, j) sull'indice
-      globale del nodo
-    - x_nodi_m, theta_nodi_rad: coordinate parametriche della griglia
-    - spessore_m, raggio_m, lunghezza_m: geometria della shell
-    """
-    n_x = int(mesh["n_elementi_x"])
-    n_theta = int(mesh["n_elementi_circonferenza"])
-    if n_x < 1 or n_theta < 3:
-        raise ValueError(
-            "Servono almeno 1 elemento lungo x e 3 lungo la circonferenza."
-        )
+    def indice(i, j):
+        return i * n_elementi_theta + j
 
-    lunghezza_m = fascetta["lunghezza_mm"] / 1000.0
-    spessore_m = fascetta["spessore_mm"] / 1000.0
-    raggio_m = raggio_medio_m(fascetta)
-
-    # Lungo x la griglia è aperta (n_x + 1 file di nodi), lungo theta è
-    # chiusa (n_theta file, l'ultima si ricongiunge alla prima).
-    x_nodi_m = np.linspace(0.0, lunghezza_m, n_x + 1)
-    theta_nodi_rad = np.linspace(0.0, 2.0 * np.pi, n_theta, endpoint=False)
-
-    X, Theta = np.meshgrid(x_nodi_m, theta_nodi_rad, indexing="ij")
-    nodi = np.column_stack(
-        (
-            X.ravel(),
-            (raggio_m * np.sin(Theta)).ravel(),
-            (raggio_m * np.cos(Theta)).ravel(),
-        )
+    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,
     )
 
-    indice_nodo = np.arange((n_x + 1) * n_theta).reshape(n_x + 1, n_theta)
-
-    i = np.arange(n_x)[:, None]
-    j = np.arange(n_theta)[None, :]
-    j_succ = (j + 1) % n_theta
-    elementi = np.stack(
-        (
-            np.broadcast_to(indice_nodo[i, j], (n_x, n_theta)),
-            np.broadcast_to(indice_nodo[i + 1, j], (n_x, n_theta)),
-            np.broadcast_to(indice_nodo[i + 1, j_succ], (n_x, n_theta)),
-            np.broadcast_to(indice_nodo[i, j_succ], (n_x, n_theta)),
-        ),
-        axis=-1,
-    ).reshape(-1, 4)
+    # 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 {
-        "nodi": nodi,
-        "elementi": elementi,
-        "indice_nodo": indice_nodo,
+        "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,
-        "n_elementi_x": n_x,
-        "n_elementi_circonferenza": n_theta,
-        "spessore_m": spessore_m,
-        "raggio_m": raggio_m,
-        "lunghezza_m": lunghezza_m,
+        "arco_nodi_m": arco_nodi_m,
+        "connettivita": connettivita,
+        "segmenti_bordo": segmenti_bordo,
+        "area_totale_m2": lunghezza_m * circonferenza_m,
     }
 
 
-def riepilogo_mesh(mesh_dati: dict) -> str:
-    n_nodi = len(mesh_dati["nodi"])
-    n_elementi = len(mesh_dati["elementi"])
-    passo_x_mm = 1000.0 * mesh_dati["lunghezza_m"] / mesh_dati["n_elementi_x"]
-    passo_circ_mm = (
-        1000.0
-        * 2.0
-        * np.pi
-        * mesh_dati["raggio_m"]
-        / mesh_dati["n_elementi_circonferenza"]
-    )
-    return (
-        f"Mesh shell: {n_elementi} elementi quadrangolari, {n_nodi} nodi\n"
-        f"  griglia: {mesh_dati['n_elementi_x']} (x) × "
-        f"{mesh_dati['n_elementi_circonferenza']} (circonferenza)\n"
-        f"  passo elemento: {passo_x_mm:.2f} mm (x), {passo_circ_mm:.2f} mm (circonferenza)\n"
-        f"  raggio medio: {1000.0 * mesh_dati['raggio_m']:.2f} mm, "
-        f"spessore shell: {1000.0 * mesh_dati['spessore_m']:.3f} mm"
-    )
+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

@@ -0,0 +1,281 @@
+# 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) 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
+
+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_termico,
+    prepara_stato_termico,
+)
+
+# Istante di inizio dei fotogrammi mostrati (la simulazione parte comunque da 0).
+T_INIZIO_ANIMAZIONE_S = 0.40
+
+# Istante di fine dell'animazione.
+T_FINE_ANIMAZIONE_S = 30
+
+# Tempo simulato tra un fotogramma e il successivo.
+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, sensore)
+    mesh = stato["mesh"]
+    dt = stato["dt_s"]
+    interpolatore = stato["interpolatore_sensore"]
+
+    # 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 = campo_iniziale(stato)
+    T_letta = sens.temperatura_osservata_C(interpolatore, T)
+
+    tempi, campi, flussi, x_riferimenti = [], [], [], []
+    T_vere, T_lette = [], []
+
+    prossimo_frame_t = 0.0
+    t = 0.0
+    while t <= T_FINE_ANIMAZIONE_S + 1e-12:
+        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(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
+
+    return {
+        "tempi": np.array(tempi),
+        "campi": campi,
+        "flussi": flussi,
+        "x_riferimenti": np.array(x_riferimenti),
+        "T_vere": np.array(T_vere),
+        "T_lette": np.array(T_lette),
+        "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"],
+        "offset_sorgente_mm": sorgente["offset_y_percorso_m"] * 1000.0,
+        "sorgente": sorgente,
+        "T_ambiente_C": aria["temperatura_ambiente_C"],
+    }
+
+
+def main() -> None:
+    rng = random.Random(SIMULAZIONE["seed"])
+    cfg_run = configurazione_randomizzata(1, rng)
+    dati = simula_campi(cfg_run)
+
+    tempi = dati["tempi"]
+    indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S))
+    n_frame = len(tempi) - indice_inizio
+
+    lunghezza_mm = dati["lunghezza_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]
+
+    # 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_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("Superficie sviluppata della fascetta: sorgenti in transito e sensore")
+
+    # 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
+    )
+    ax_flusso.annotate(
+        "verso di marcia",
+        xy=(0.28, 0.85), xytext=(0.55, 0.85), xycoords="axes fraction",
+        arrowprops={"arrowstyle": "->", "color": "gray"},
+        color="gray", va="center",
+    )
+    ax_flusso.set_xlim(*x_vista_mm)
+    ax_flusso.set_ylim(0.0, q_max_MW * 1.25)
+    ax_flusso.set_ylabel("q(x) [MW/m²]")
+    ax_flusso.grid(True, alpha=0.3)
+    ax_flusso.set_xticklabels([])
+
+    # 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=dati["T_ambiente_C"],
+        vmax=T_max,
+    )
+    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_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 superficie, mantenendo allineati gli assi x.
+    barra = fig.colorbar(
+        mappa, ax=(ax_flusso, ax_superficie, ax_storia), pad=0.02, aspect=35
+    )
+    barra.set_label("T [°C]")
+
+    # Pannello 3: temperatura nel punto osservato dal sensore.
+    linea_vera, = ax_storia.plot([], [], label="T vera lato interno")
+    linea_letta, = ax_storia.plot([], [], label="T sensore (con inerzia)")
+    cursore = ax_storia.axvline(tempi[indice_inizio], color="gray", linewidth=0.8)
+    ax_storia.set_xlim(0.0, T_FINE_ANIMAZIONE_S)
+    ax_storia.set_ylim(15.0, max(dati["T_vere"].max(), dati["T_lette"].max()) * 1.08)
+    ax_storia.set_xlabel("Tempo [s]")
+    ax_storia.set_ylabel("T [°C]")
+    ax_storia.legend(loc="upper left")
+    ax_storia.grid(True, alpha=0.3)
+
+    testo_tempo = ax_flusso.set_title(f"t = {tempi[indice_inizio]:.3f} s", loc="right")
+
+    def aggiorna(frame: int):
+        k = indice_inizio + frame
+        t = tempi[k]
+
+        linea_flusso.set_data(dati["x_nodi_mm"], dati["flussi"][k] / 1e6)
+
+        x_sorgenti_mm = (
+            dati["x_riferimenti"][k] * 1000.0
+            + np.arange(numero_sorgenti) * distanza_mm
+        )
+        visibili = (x_sorgenti_mm >= x_vista_mm[0]) & (x_sorgenti_mm <= x_vista_mm[1])
+        marker_sorgenti.set_data(
+            x_sorgenti_mm[visibili],
+            np.full(int(visibili.sum()), q_max_MW * 1.12),
+        )
+
+        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])
+        cursore.set_xdata([t, t])
+        testo_tempo.set_text(f"t = {t:.3f} s")
+
+        return (
+            linea_flusso, marker_sorgenti, mappa,
+            linea_vera, linea_letta, cursore, testo_tempo,
+        )
+
+    animazione = FuncAnimation(
+        fig, aggiorna, frames=n_frame, interval=INTERVALLO_RIPRODUZIONE_MS, blit=False
+    )
+
+    # Se il backend non è interattivo si salva una GIF invece di mostrare la finestra.
+    if matplotlib.get_backend().lower() == "agg":
+        percorso = Path("dataset") / "animazione_superficie.gif"
+        animazione.save(
+            percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS)
+        )
+        print(f"Backend non interattivo: animazione salvata in {percorso}")
+        return
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

plot_animazione_3d.py

@@ -0,0 +1,140 @@
+# Animazione 3D del barattolo in vista isometrica: colormap della temperatura
+# sulla superficie cilindrica.
+#
+# 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
+
+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_INIZIO_ANIMAZIONE_S, simula_campi
+from simulate import configurazione_randomizzata
+
+# Tempo simulato tra un fotogramma e il successivo dell'animazione 3D. Più
+# largo di DT_FRAME_S_BASE perché ricreare una superficie 3D a ogni
+# fotogramma è più costoso della semplice imshow 2D.
+DT_FRAME_S = 0.15
+
+# Millisecondi tra i fotogrammi in riproduzione.
+INTERVALLO_RIPRODUZIONE_MS = 60
+
+# 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:
+    rng = random.Random(SIMULAZIONE["seed"])
+    cfg_run = configurazione_randomizzata(1, rng)
+    dati = simula_campi(cfg_run)
+
+    tempi = dati["tempi"]
+    indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S))
+    passo = max(1, round(DT_FRAME_S / DT_FRAME_S_BASE))
+    indici_frame = list(range(indice_inizio, len(tempi), passo))
+
+    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"]
+
+    T_ambiente = dati["T_ambiente_C"]
+    sorgente = dati["sorgente"]
+    raggio_m = mesh["raggio_medio_m"]
+    lunghezza_m = mesh["lunghezza_m"]
+
+    Xm, Ym, Zm = coordinate_3d(mesh)
+
+    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.
+    norm = matplotlib.colors.Normalize(
+        vmin=T_ambiente - 0.4 * (T_max - T_ambiente), vmax=T_max
+    )
+    cmap = matplotlib.colormaps["inferno"]
+    lightsource = matplotlib.colors.LightSource(azdeg=315, altdeg=45)
+
+    fig = plt.figure(figsize=(9, 7))
+    ax = fig.add_subplot(projection="3d")
+    ax.view_init(elev=35.264, azim=45)
+    ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m))
+    ax.set_xlabel("x [m]")
+    ax.set_axis_off()
+
+    # La colorbar mostra il range reale delle temperature: il norm esteso
+    # verso il basso serve solo a schiarire il colore di base della
+    # superficie fredda, non deve comparire nella scala mostrata all'utente.
+    norm_colorbar = matplotlib.colors.Normalize(vmin=T_ambiente, vmax=T_max)
+    mappabile = cm.ScalarMappable(cmap=cmap, norm=norm_colorbar)
+    mappabile.set_array([])
+    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()
+
+        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,
+            **(bordi if MOSTRA_MESH else {"linewidth": 0}),
+        )
+
+        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():
+            ax.scatter(
+                x_sorgenti_m[visibili],
+                np.full(int(visibili.sum()), raggio_m * np.sin(theta_sorgente) * 1.05),
+                np.full(int(visibili.sum()), raggio_m * np.cos(theta_sorgente) * 1.05),
+                color="cyan", s=25, depthshade=False,
+            )
+
+        ax.set_title(f"t = {tempi[k]:.3f} s")
+        return ()
+
+    animazione = FuncAnimation(
+        fig, disegna_frame, frames=indici_frame,
+        interval=INTERVALLO_RIPRODUZIONE_MS, blit=False,
+    )
+
+    if matplotlib.get_backend().lower() == "agg":
+        percorso = Path("dataset") / "animazione_3d.gif"
+        animazione.save(
+            percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS)
+        )
+        print(f"Backend non interattivo: animazione salvata in {percorso}")
+        return
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

plot_animazione_fem.py

@@ -1,198 +0,0 @@
-# Animazione 3D isometrica della fascetta con il campo di temperatura calcolato
-# a elementi finiti sulla mesh shell (fem.py).
-#
-# La coordinata circonferenziale è una direzione risolta del modello FEM,
-# quindi la temperatura dipinta su ogni elemento è quella effettivamente
-# calcolata, senza ricostruzioni.
-#
-# Ogni elemento shell è disegnato come faccia piana con colore pari alla media
-# dei suoi quattro valori nodali. L'ombreggiatura è calcolata a mano dalla
-# normale radiale dell'elemento: Poly3DCollection ombreggia solo alla
-# creazione, mentre qui i colori delle facce cambiano a ogni fotogramma.
-
-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 mpl_toolkits.mplot3d.art3d import Poly3DCollection
-
-from config import USCITA
-from fem import simula_campo_fem
-from mesh import genera_mesh, riepilogo_mesh
-
-# Istante del primo fotogramma mostrato: i primi istanti sono ancora uniformi
-# alla temperatura ambiente.
-T_INIZIO_ANIMAZIONE_S = 0.40
-
-# Tempo simulato tra un fotogramma e il successivo.
-DT_FRAME_S = 0.15
-
-# Millisecondi tra i fotogrammi in riproduzione.
-INTERVALLO_RIPRODUZIONE_MS = 60
-
-# Direzione della luce (x = asse del cilindro, y-z = piano della sezione) e
-# quota di luce ambiente, per dare volume alla superficie.
-DIREZIONE_LUCE = np.array([-0.25, 0.45, 0.85])
-LUCE_AMBIENTE = 0.45
-
-# Se True disegna il reticolo degli elementi sopra il campo. Con la mesh fitta
-# i bordi degli elementi sul lato nascosto rendono la superficie confusa,
-# quindi per default i bordi prendono il colore della faccia: così le facce
-# adiacenti combaciano senza lasciare fessure di antialiasing.
-MOSTRA_BORDI_ELEMENTI = False
-
-
-def _illuminazione(mesh_dati: dict) -> np.ndarray:
-    # Fattore moltiplicativo di luminosità per elemento, dalla normale
-    # radiale uscente valutata nel baricentro dell'elemento.
-    facce = mesh_dati["nodi"][mesh_dati["elementi"]]
-    baricentri = facce.mean(axis=1)
-    normali = baricentri.copy()
-    # La normale della superficie cilindrica è radiale: nessuna componente
-    # lungo l'asse x.
-    normali[:, 0] = 0.0
-    normali /= np.linalg.norm(normali, axis=1, keepdims=True)
-
-    luce = DIREZIONE_LUCE / np.linalg.norm(DIREZIONE_LUCE)
-    diffusa = np.clip(normali @ luce, 0.0, None)
-    return LUCE_AMBIENTE + (1.0 - LUCE_AMBIENTE) * diffusa
-
-
-def main() -> None:
-    mesh_dati = genera_mesh()
-    print(riepilogo_mesh(mesh_dati))
-
-    print("Integrazione FEM in corso...")
-    dati = simula_campo_fem(dt_frame_s=DT_FRAME_S, mesh_dati=mesh_dati)
-
-    tempi = dati["tempi"]
-    indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S))
-    indici_frame = list(range(indice_inizio, len(tempi)))
-
-    T_ambiente = dati["T_ambiente_C"]
-    T_max = max(campo.max() for campo in dati["campi"])
-    print(
-        f"T massima nodale: {T_max:.1f} °C — "
-        f"T massima nel punto del sensore: {dati['T_vere'].max():.1f} °C\n"
-        f"Skin depth (diagnostica): {1000.0 * dati['skin_depth_m']:.3f} mm"
-    )
-
-    elementi = mesh_dati["elementi"]
-    facce = mesh_dati["nodi"][elementi]
-    illuminazione = _illuminazione(mesh_dati)[:, None]
-
-    raggio_m = mesh_dati["raggio_m"]
-    lunghezza_m = mesh_dati["lunghezza_m"]
-
-    sorgente = dati["sorgente"]
-    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
-    distanza_m = sorgente.get("distanza_sorgenti_m", 0.0)
-    theta_sorgente = sorgente["offset_y_percorso_m"] / raggio_m
-
-    norm = matplotlib.colors.Normalize(vmin=T_ambiente, vmax=T_max)
-    cmap = matplotlib.colormaps["inferno"]
-
-    fig = plt.figure(figsize=(9, 8.5))
-    griglia = fig.add_gridspec(2, 1, height_ratios=[2.6, 1.0], hspace=0.05)
-    ax = fig.add_subplot(griglia[0], projection="3d")
-    ax_storia = fig.add_subplot(griglia[1])
-    ax.view_init(elev=35.264, azim=45)
-    # zoom > 1 riempie il riquadro: con set_axis_off gli assi 3D lascerebbero
-    # altrimenti molto margine vuoto sopra il pannello della storia.
-    ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m), zoom=1.12)
-    ax.set_axis_off()
-    ax.set_xlim(0.0, lunghezza_m)
-    ax.set_ylim(-raggio_m, raggio_m)
-    ax.set_zlim(-raggio_m, raggio_m)
-
-    collezione = Poly3DCollection(facce, linewidths=0.25, shade=False)
-    ax.add_collection3d(collezione)
-
-    marker_sorgenti = ax.plot(
-        [], [], [], "o", color="cyan", markersize=6, zorder=10
-    )[0]
-
-    mappabile = cm.ScalarMappable(cmap=cmap, norm=norm)
-    mappabile.set_array([])
-    # La colorbar è agganciata a entrambi i pannelli per non restringere solo
-    # quello 3D: così i due riquadri restano della stessa larghezza.
-    barra = fig.colorbar(mappabile, ax=(ax, ax_storia), shrink=0.6, pad=0.02, aspect=30)
-    barra.set_label("T [°C]")
-
-    titolo = ax.set_title("")
-
-    # Pannello inferiore: storia della temperatura nel punto osservato dal
-    # sensore. T_vere è il valore nodale calcolato dal FEM, T_lette lo stesso
-    # segnale degradato dal sensore reale (inerzia, rumore, quantizzazione).
-    # A parte il rumore le due curve sono quasi sovrapposte: la prima è
-    # tracciata spessa e trasparente perché la seconda resti leggibile sopra
-    # di essa.
-    linea_vera, = ax_storia.plot(
-        [], [], color="tab:blue", linewidth=3.0, alpha=0.4,
-        label="T nodo del sensore (FEM)",
-    )
-    linea_letta, = ax_storia.plot(
-        [], [], color="tab:red", linewidth=1.0,
-        label="T misurata dal sensore (inerzia + rumore + quantizzazione)",
-    )
-    cursore = ax_storia.axvline(tempi[indice_inizio], color="gray", linewidth=0.8)
-    ax_storia.set_xlim(0.0, tempi[-1])
-    ax_storia.set_ylim(
-        T_ambiente - 0.05 * (T_max - T_ambiente),
-        max(dati["T_vere"].max(), dati["T_lette"].max()) * 1.08,
-    )
-    ax_storia.set_xlabel("Tempo [s]")
-    ax_storia.set_ylabel("T [°C]")
-    ax_storia.legend(loc="upper left")
-    ax_storia.grid(True, alpha=0.3)
-
-    def disegna_frame(k: int):
-        T_elementi = dati["campi"][k][elementi].mean(axis=1)
-        colori = cmap(norm(T_elementi))
-        # Ombreggiatura moltiplicativa sui soli canali RGB.
-        colori[:, :3] *= illuminazione
-        collezione.set_facecolor(colori)
-        collezione.set_edgecolor("#40404060" if MOSTRA_BORDI_ELEMENTI else colori)
-
-        x_sorgenti_m = dati["x_riferimenti"][k] + np.arange(numero_sorgenti) * distanza_m
-        visibili = (x_sorgenti_m >= 0.0) & (x_sorgenti_m <= lunghezza_m)
-        n_visibili = int(visibili.sum())
-        marker_sorgenti.set_data_3d(
-            x_sorgenti_m[visibili],
-            np.full(n_visibili, raggio_m * np.sin(theta_sorgente) * 1.06),
-            np.full(n_visibili, raggio_m * np.cos(theta_sorgente) * 1.06),
-        )
-
-        linea_vera.set_data(tempi[: k + 1], dati["T_vere"][: k + 1])
-        linea_letta.set_data(tempi[: k + 1], dati["T_lette"][: k + 1])
-        cursore.set_xdata([tempi[k], tempi[k]])
-
-        titolo.set_text(
-            f"FEM shell — t = {tempi[k]:.2f} s — "
-            f"T sensore = {dati['T_lette'][k]:.1f} °C"
-        )
-        return collezione, marker_sorgenti, linea_vera, linea_letta, cursore, titolo
-
-    animazione = FuncAnimation(
-        fig, disegna_frame, frames=indici_frame,
-        interval=INTERVALLO_RIPRODUZIONE_MS, blit=False,
-    )
-
-    if matplotlib.get_backend().lower() == "agg":
-        cartella = Path(USCITA["cartella"])
-        cartella.mkdir(parents=True, exist_ok=True)
-        percorso = cartella / "animazione_fem.gif"
-        animazione.save(
-            percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS)
-        )
-        print(f"Backend non interattivo: animazione salvata in {percorso}")
-        return
-
-    plt.show()
-
-
-if __name__ == "__main__":
-    main()

plot_csv.py

@@ -0,0 +1,48 @@
+from pathlib import Path
+
+import matplotlib
+import matplotlib.pyplot as plt
+import pandas as pd
+
+
+def main() -> None:
+    percorso_csv = Path("dataset/run_0001.csv")
+    if not percorso_csv.exists():
+        raise FileNotFoundError(
+            "File CSV non trovato. Esegui prima `python simulate.py`."
+        )
+
+    df = pd.read_csv(percorso_csv)
+
+    fig_temperature = plt.figure()
+    plt.plot(df["tempo_s"], df["T_vera_lato_sensore_C"], label="Temperatura vera lato interno")
+    plt.plot(df["tempo_s"], df["T_misurata_sensore_C"], label="Temperatura misurata dal sensore")
+    plt.xlabel("Tempo [s]")
+    plt.ylabel("Temperatura [°C]")
+    plt.title("Sensore fisso, sorgenti in moto")
+    plt.legend()
+    plt.grid(True)
+    plt.tight_layout()
+
+    fig_flusso = plt.figure()
+    plt.plot(df["tempo_s"], df["flusso_termico_sorgente_W_m2"])
+    plt.xlabel("Tempo [s]")
+    plt.ylabel("Flusso termico efficace [W/m²]")
+    plt.title("Flusso termico nel punto x osservato dal sensore")
+    plt.grid(True)
+    plt.tight_layout()
+
+    # Se il backend non è interattivo (nessun display o backend GUI
+    # disponibile), plt.show() non aprirebbe nulla: si salvano i PNG.
+    if matplotlib.get_backend().lower() == "agg":
+        for fig, nome in [(fig_temperature, "temperature"), (fig_flusso, "flusso")]:
+            percorso = Path("dataset") / f"grafico_{nome}.png"
+            fig.savefig(percorso, dpi=150)
+            print(f"Backend non interattivo: grafico salvato in {percorso}")
+        return
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

plot_mesh.py

@@ -1,60 +0,0 @@
-# Visualizzazione della sola mesh a elementi shell della fascetta, in vista
-# isometrica: facce chiare con i bordi degli elementi in evidenza.
-
-from pathlib import Path
-
-import matplotlib
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d.art3d import Poly3DCollection
-
-from config import USCITA
-from mesh import genera_mesh, riepilogo_mesh
-
-
-def main() -> None:
-    dati = genera_mesh()
-    print(riepilogo_mesh(dati))
-
-    nodi = dati["nodi"]
-    facce = nodi[dati["elementi"]]
-
-    raggio_m = dati["raggio_m"]
-    lunghezza_m = dati["lunghezza_m"]
-
-    fig = plt.figure(figsize=(9, 7))
-    ax = fig.add_subplot(projection="3d")
-    ax.view_init(elev=35.264, azim=45)
-    ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m))
-    ax.set_axis_off()
-
-    collezione = Poly3DCollection(
-        facce,
-        facecolors="#d9e3f0",
-        edgecolors="#2b3a55",
-        linewidths=0.4,
-        shade=True,
-        lightsource=matplotlib.colors.LightSource(azdeg=315, altdeg=45),
-    )
-    ax.add_collection3d(collezione)
-
-    ax.set_xlim(0.0, lunghezza_m)
-    ax.set_ylim(-raggio_m, raggio_m)
-    ax.set_zlim(-raggio_m, raggio_m)
-    ax.set_title(
-        f"Mesh shell: {dati['n_elementi_x']} × "
-        f"{dati['n_elementi_circonferenza']} elementi quadrangolari"
-    )
-
-    if matplotlib.get_backend().lower() == "agg":
-        cartella = Path(USCITA["cartella"])
-        cartella.mkdir(parents=True, exist_ok=True)
-        percorso = cartella / "mesh.png"
-        fig.savefig(percorso, dpi=150)
-        print(f"Backend non interattivo: immagine salvata in {percorso}")
-        return
-
-    plt.show()
-
-
-if __name__ == "__main__":
-    main()

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

@@ -0,0 +1,386 @@
+# 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
+import random
+import shutil
+from concurrent.futures import ProcessPoolExecutor
+from copy import deepcopy
+from pathlib import Path
+
+import numpy as np
+
+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
+
+
+def calcola_skin_depth_m(materiale: dict, frequenza_hz: float) -> float:
+    # Skin depth elettromagnetica approssimata:
+    # delta = sqrt(2 * rho_e / (omega * mu))
+    #
+    # Semplificata. Per acciai ferromagnetici il comportamento reale
+    # è fortemente non lineare con temperatura e campo magnetico.
+    rho_e = materiale["resistivita_elettrica_ohm_m"]
+    mu_r = materiale["permeabilita_relativa"]
+    omega = 2.0 * math.pi * frequenza_hz
+    mu = MU0 * mu_r
+    return math.sqrt(2.0 * rho_e / (omega * mu))
+
+
+def numero_fourier_spessore(
+    materiale: dict, spessore_m: float, sorgente: dict
+) -> float:
+    # Verifica quantitativa dell'ipotesi di parete termicamente sottile:
+    #
+    #     Fo = alpha * (sigma / |v|) / spessore^2
+    #
+    # 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"]
+    )
+    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, 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_m = fascetta["lunghezza_mm"] / 1000.0
+    spessore_m = fascetta["spessore_mm"] / 1000.0
+    raggio_medio_m = fascetta["diametro_medio_mm"] / 2000.0
+
+    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_m = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
+    else:
+        skin_depth_m = float(sorgente["skin_depth_fissa_m"])
+
+    x_sensore_m = sensore["x_mm"] / 1000.0
+
+    return {
+        "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"],
+        "skin_depth_m": skin_depth_m,
+        "numero_fourier_spessore": numero_fourier_spessore(
+            materiale, spessore_m, sorgente
+        ),
+    }
+
+
+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
+    )
+
+    _, 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:
+    fascetta = deepcopy(FASCETTA)
+    aria = deepcopy(ARIA)
+    sorgente = deepcopy(SORGENTE)
+    sensore = deepcopy(SENSORE)
+
+    if RANDOMIZZAZIONE.get("abilitata", False):
+        def perturba_rel(valore: float, std_rel: float, fattore_min: float = 0.1) -> float:
+            fattore = rng.gauss(1.0, std_rel)
+            fattore = max(fattore_min, fattore)
+            return valore * fattore
+
+        sorgente["velocita_m_s"] = perturba_rel(
+            sorgente["velocita_m_s"],
+            RANDOMIZZAZIONE["velocita_std_rel"],
+        )
+        sorgente["flusso_termico_picco_W_m2"] = perturba_rel(
+            sorgente["flusso_termico_picco_W_m2"],
+            RANDOMIZZAZIONE["flusso_picco_std_rel"],
+        )
+        sorgente["sigma_punto_m"] = perturba_rel(
+            sorgente["sigma_punto_m"],
+            RANDOMIZZAZIONE["sigma_punto_std_rel"],
+        )
+        sorgente["offset_y_percorso_m"] = rng.uniform(
+            -RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
+            RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
+        )
+        aria["temperatura_ambiente_C"] += rng.gauss(
+            0.0,
+            RANDOMIZZAZIONE["temperatura_ambiente_std_C"],
+        )
+        sensore["rumore_std_C"] = perturba_rel(
+            sensore["rumore_std_C"],
+            RANDOMIZZAZIONE["rumore_sensore_std_rel"],
+            fattore_min=0.0,
+        )
+
+    return {
+        "id_run": f"run_{indice_run:04d}",
+        "fascetta": fascetta,
+        "aria": aria,
+        "sorgente": sorgente,
+        "sensore": sensore,
+    }
+
+
+def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
+    fascetta = cfg_run["fascetta"]
+    aria = cfg_run["aria"]
+    sorgente = cfg_run["sorgente"]
+    sensore = cfg_run["sensore"]
+
+    nome_materiale = fascetta["materiale"]
+
+    stato = prepara_stato_termico(fascetta, aria, sorgente, sensore)
+    mesh = stato["mesh"]
+    interpolatore = stato["interpolatore_sensore"]
+    skin_depth = stato["skin_depth_m"]
+
+    dt = stato["dt_s"]
+    durata = SIMULAZIONE["durata_s"]
+    periodo_campionamento = 1.0 / SIMULAZIONE["frequenza_campionamento_hz"]
+
+    T = campo_iniziale(stato)
+    T_letta = sens.temperatura_osservata_C(interpolatore, T)
+
+    prossimo_campione_t = 0.0
+    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)
+
+    with output_csv.open("w", newline="") as f:
+        writer = csv.writer(f)
+        writer.writerow([
+            "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",
+        ])
+
+        t = 0.0
+        while t <= durata + 1e-12:
+            T, info = passo_termico(stato, T, t)
+
+            # Temperatura vera della superficie interna nel punto osservato
+            # dal sensore infrarosso.
+            T_vera = sens.temperatura_osservata_C(interpolatore, T)
+            T_letta = sens.aggiorna_inerzia_C(
+                T_letta, T_vera, dt, sensore["costante_tempo_s"]
+            )
+
+            # 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_letta + rng.gauss(0.0, sensore["rumore_std_C"])
+                misurata = sens.quantizza(misurata, sensore["quantizzazione_C"])
+
+                T_misurata_max = max(T_misurata_max, misurata)
+
+                writer.writerow([
+                    cfg_run["id_run"],
+                    f"{t:.6f}",
+                    f"{info['x_riferimento_m']:.9f}",
+                    f"{sorgente['offset_y_percorso_m']:.9f}",
+                    f"{info['flusso_al_sensore_W_m2']:.6f}",
+                    f"{skin_depth:.9e}",
+                    f"{T_vera:.6f}",
+                    f"{misurata:.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}",
+                    f"{sorgente['flusso_termico_picco_W_m2']:.6f}",
+                    nome_materiale,
+                ])
+                prossimo_campione_t += periodo_campionamento
+
+            t += dt
+
+    return {
+        "id_run": cfg_run["id_run"],
+        "file_csv": str(output_csv.name),
+        "materiale": nome_materiale,
+        "diametro_medio_m": fascetta["diametro_medio_mm"] / 1000.0,
+        "lunghezza_m": fascetta["lunghezza_mm"] / 1000.0,
+        "spessore_m": fascetta["spessore_mm"] / 1000.0,
+        "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,
+        "temperatura_ambiente_C": aria["temperatura_ambiente_C"],
+        "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"],
+        "x_inizio_m": sorgente["x_inizio_m"],
+        "x_fine_m": sorgente["x_fine_m"],
+        "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),
+        "sigma_punto_m": sorgente["sigma_punto_m"],
+        "flusso_termico_picco_W_m2": sorgente["flusso_termico_picco_W_m2"],
+        "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,
+    }
+
+
+def _esegui_run(indice_e_seme: tuple[int, int]) -> dict:
+    # Ogni run riceve un seme indipendente derivato dal seed globale, così
+    # l'esecuzione in parallelo resta riproducibile indipendentemente
+    # dall'ordine in cui i processi la completano.
+    indice, seme = indice_e_seme
+    rng = random.Random(seme)
+    cfg_run = configurazione_randomizzata(indice, rng)
+    cartella_output = Path(SIMULAZIONE["cartella_output"])
+    percorso_csv = cartella_output / f"{cfg_run['id_run']}.csv"
+    return simula_singolo(cfg_run, percorso_csv, rng)
+
+
+def main() -> None:
+    cartella_output = Path(SIMULAZIONE["cartella_output"])
+    if cartella_output.exists():
+        shutil.rmtree(cartella_output)
+    cartella_output.mkdir(parents=True, exist_ok=True)
+
+    rng_semi = random.Random(SIMULAZIONE["seed"])
+    num_run = SIMULAZIONE["num_run"]
+    semi = [rng_semi.randrange(2**63) for _ in range(num_run)]
+
+    num_processi = SIMULAZIONE["num_processi"] or os.cpu_count() or 1
+    with ProcessPoolExecutor(max_workers=num_processi) as executor:
+        righe_metadata = list(
+            executor.map(_esegui_run, enumerate(semi, start=1))
+        )
+
+    percorso_metadata = cartella_output / "metadata.csv"
+    with percorso_metadata.open("w", newline="") as f:
+        writer = csv.DictWriter(f, fieldnames=list(righe_metadata[0].keys()))
+        writer.writeheader()
+        writer.writerows(righe_metadata)
+
+    print(f"Generati {len(righe_metadata)} run in: {cartella_output.resolve()}")
+    print(f"Metadata: {percorso_metadata.resolve()}")
+
+
+if __name__ == "__main__":
+    main()

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()