Esporta in CSV la temperatura simulata e la lettura del sensore

Nuova azione `python main.py csv`: esporta_csv.py integra il campo FEM e scrive
output/csv/run_XXXX.csv con la temperatura vera della parete nel punto
osservato, la lettura del sensore reale e il massimo nodale istante per
istante, più un metadata.csv con una riga per analisi. La cartella è ricreata a
ogni esecuzione.

Le righe sono campionate a FEM["frequenza_campionamento_hz"],
indipendente dal passo di integrazione, e FEM["num_run"] analisi vengono
esportate di seguito: con la randomizzazione abilitata ognuna ha parametri
diversi ma riproducibili dallo stesso seed.

simula_campo_fem accetta salva_campi=False, perché all'export servono solo le
serie scalari e non i campi nodali di ogni campione, e restituisce anche la
temperatura nodale massima per campione.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

.gitignore

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

CLAUDE.md

@@ -18,45 +18,39 @@ source .venv/bin/activate
 # Installare le dipendenze
 pip install -r requirements.txt
 
-# Generare il dataset (scrive dataset/run_XXXX.csv + dataset/metadata.csv)
-python simulate.py
-
-# Visualizzare il primo run
-python plot_csv.py
+# 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
 ```
 
+Ogni modulo resta eseguibile anche direttamente (`python plot_mesh.py`, `python plot_animazione_fem.py`, `python esporta_csv.py`).
+
 Attivare sempre il venv (`source .venv/bin/activate`) prima di eseguire qualsiasi comando Python.
 
 Non sono configurati test o linter.
 
 ## Architettura
 
-Generatore di dataset per misurazioni termiche pseudo-realistiche di una fascetta (anello cilindrico sottile) riscaldata da sorgenti a induzione in movimento.
+Analisi termica transitoria a elementi finiti di una fascetta (anello cilindrico sottile) riscaldata da sorgenti a induzione in movimento.
 
-**Geometria:** la fascetta ha diametro, spessore e lunghezza configurabili. Il dominio simulato è la sezione rettangolare lunghezza × spessore, con origine (0, 0) nel vertice in alto a sinistra: x = lunghezza (le sorgenti si muovono in direzione -x sul lato esterno), z = spessore (0 = lato esterno, spessore = lato interno). La coordinata circonferenziale y non è risolta: l'offset y delle sorgenti è collassato in un'attenuazione gaussiana del flusso. Il sensore è un pirometro a infrarossi dentro la fascetta che misura la superficie interna in un punto x fisso.
+**Geometria:** la fascetta ha diametro, 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.
 
 **Flusso dei dati:**
 
-1. `config.py` — tutti i parametri configurabili (dizionari SIMULAZIONE, FASCETTA, ARIA, SORGENTE, SENSORE, RANDOMIZZAZIONE)
-2. `materials.py` — dizionario MATERIALI con proprietà termofisiche ed elettriche per materiale
-3. `simulate.py` — motore principale: genera N run randomizzati, scrive i CSV, scrive `metadata.csv`
-4. `plot_csv.py` — visualizzazione autonoma per un singolo run
+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`
 
-**Pipeline fisica dentro `simula_singolo()` in [simulate.py](simulate.py):**
+**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.
 
-- La skin depth è calcolata dalla resistività elettrica del materiale e dalla frequenza di induzione (`calcola_skin_depth_m`)
-- Le sorgenti gaussiane in movimento producono un profilo di flusso termico superficiale q(x) sul lato esterno, variabile nel tempo (`profilo_flusso_incidente_W_m2`)
-- Quel flusso è ridistribuito volumetricamente attraverso lo spessore con decadimento esponenziale in z (`profilo_deposizione_z_1_m`): q_vol(x, z) = q(x) · p(z)
-- Uno schema 2D a volumi finiti con Eulero implicito integra l'equazione del calore su `n_nodi_x × n_nodi_z` celle: `prepara_stato_termico` costruisce griglia, coefficienti e matrice sparsa fattorizzata LU una volta per run (`costruisci_solutore_implicito_2d`, che restituisce l'oggetto `splu`), poi `passo_implicito` avanza il campo risolvendo solo il sistema triangolare
-- Le condizioni al contorno sono incorporate nella matrice: convezione su tutti e quattro i lati della sezione, più un termine di conduzione circonferenziale (y) verso il resto della fascetta assunto a temperatura ambiente. Il termine è un'equazione di aletta ricavata sull'intero volume del cilindro: il calore conduce lungo y attraverso l'intero spessore mentre le superfici esterna e interna dell'intero cilindro scambiano per convezione, dando q_y = -(h_esterno + h_interno)/spessore · (T - T_amb), senza parametri di conduzione y configurabili a parte
-- La temperatura iniziale del campo è la temperatura ambiente (randomizzata per run)
-- L'output del sensore aggiunge inerzia del primo ordine, rumore gaussiano e quantizzazione
+**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` e `passo_implicito` sono condivisi con `plot_animazione.py`, che riproduce la fisica di run_0001 per animare la sezione: ogni modifica alla fisica va fatta lì, non duplicata.
+**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`.
 
-**Randomizzazione per run** (`configurazione_randomizzata`): ogni run perturba velocità, flusso di picco, sigma del punto, offset y, temperatura ambiente e rumore del sensore con estrazioni gaussiane/uniformi da un RNG con seed fisso, garantendo riproducibilità.
-
-**Schema di output** (`dataset/run_XXXX.csv`): serie temporale con colonne `id_run, tempo_s, x_sorgente_m, offset_y_sorgente_m, flusso_termico_sorgente_W_m2, skin_depth_m, T_vera_lato_sensore_C, T_misurata_sensore_C, T_lato_caldo_C, T_ambiente_C, velocita_m_s, sigma_punto_m, flusso_picco_W_m2, materiale`. `metadata.csv` ha una riga per run con tutti i parametri e le temperature di picco.
+**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`).
 
 ## Convenzioni su `config.py`
 
@@ -64,8 +58,8 @@ Ogni parametro in [config.py](config.py) ha un commento che spiega solo cos'è (
 
 ## Vincoli progettuali chiave
 
-- Il modello è 2D nella sezione (x = lunghezza, z = spessore). La coordinata circonferenziale y non è risolta spazialmente — l'offset y del percorso delle sorgenti è collassato in un'attenuazione gaussiana del flusso, e la conduzione lungo y è un termine di scambio lineare verso la temperatura ambiente; il diametro è registrato solo come geometria del setup.
+- Il modello è una shell: la circonferenza è risolta, lo spessore è collassato (temperatura uniforme nella parete). Non introdurre effetti attraverso lo spessore (skin depth, gradiente radiale) senza prima discretizzarlo.
 - Le posizioni di inizio/fine corsa delle sorgenti (`x_inizio_m`, `x_fine_m`) sono distanze dal punto x del sensore lungo il verso di marcia; il segno di `velocita_m_s` determina il verso (negativo = -x).
-- La matrice implicita è costruita e fattorizzata una volta per run (proprietà del materiale costanti, nessun coefficiente dipendente dalla temperatura). Se si aggiungono proprietà dipendenti dalla temperatura, la matrice deve essere ricostruita e rifattorizzata ad ogni passo temporale.
-- `simulate.py` cancella e ricrea l'intera cartella di output ad ogni esecuzione (`shutil.rmtree`).
+- 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.
 - 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,113 +1,135 @@
-# Simulatore Termico 2D — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso
+# Analisi Termica FEM — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso
 
-Questo progetto genera misurazioni CSV pseudo-realistiche della temperatura di una
-fascetta (anello cilindrico sottile) riscaldata da un gruppo di sorgenti a induzione
-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.
+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).
 
 ## Geometria
 
 La fascetta è un anello cilindrico definito da tre dimensioni:
 
-- **diametro** (default 70 mm) — il diametro del cilindro;
-- **spessore** (default 0.12 mm) — lo spessore della parete;
+- **diametro** (default 70 mm) — il diametro esterno del cilindro;
+- **spessore** (default 0.18 mm) — lo spessore della parete;
 - **lunghezza** (default 100 mm) — l'estensione assiale lungo `x`.
 
-Il dominio simulato è la **sezione rettangolare lunghezza × spessore**. Il sistema di
-coordinate ha l'origine `(0, 0)` nel vertice in alto a sinistra della sezione:
+Il dominio discretizzato è la **superficie media** del cilindro (raggio =
+(diametro − spessore)/2). Lo spessore non è discretizzato: è un attributo degli
+elementi shell.
 
 ```text
-        sorgenti (induttori), in moto verso -x
+        sorgenti (induttori), in moto lungo x sulla superficie esterna
               ▼           ▼           ▼
-  (0,0) ─────────────────────────────────────► x
-    │   ┌───────────────────────────────────┐  z = 0        lato ESTERNO
-    │   │  sezione della fascetta           │               (flusso termico)
-    │   └───────────────────────────────────┘  z = spessore lato INTERNO
-    ▼                     ┆
-    z                     ┆  linea di vista
-                          ▲
-                     sensore IR (fisso, x = 50 mm,
-                     a 10 mm dalla parete interna)
+              ┌───────────────────────────────────┐
+        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)
 ```
 
-- **x** = direzione della lunghezza, da `0` a `lunghezza`. Le sorgenti viaggiano in
-  direzione `-x` sul lato esterno.
-- **z** = direzione dello spessore, da `0` (lato esterno, dove arriva il flusso
-  termico) a `spessore` (lato interno, osservato dal sensore).
-- **y** = coordinata circonferenziale (lungo la circonferenza π·diametro). Non è
-  risolta spazialmente: vedi sotto come viene trattata.
+- **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)`.
 
 Il **sensore** è un pirometro a infrarossi posto all'interno della fascetta, a una
-distanza configurabile dalla parete interna (default 10 mm). Essendo senza contatto,
-la distanza non influenza la misura: il sensore legge la temperatura della superficie
-interna nel punto `x` configurato (default 50 mm, al centro della lunghezza).
+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`.
 
 ## Modello fisico
 
-Non è una simulazione FEM elettromagnetica + termica completa: è un generatore
-pratico di dataset. La catena di approssimazioni è la seguente.
+L'equazione risolta è quella del calore per unità di superficie media:
 
-1. **Sorgenti gaussiane in moto** — ogni sorgente ha un'impronta gaussiana di raggio
-   `sigma_punto_m`. Un gruppo di `numero_sorgenti` sorgenti equidistanti
-   (`distanza_sorgenti_m`) si muove rigidamente a velocità costante. Il profilo di
-   flusso sul lato esterno è la somma dei contributi:
-   `q(x, t) = Σᵢ q_picco · efficienza · exp(-((x - xᵢ(t))² + Δy²) / (2σ²))`.
-   L'offset circonferenziale `Δy` tra il percorso delle sorgenti e il punto osservato
-   dal sensore non è risolto spazialmente: entra come attenuazione gaussiana del flusso.
+```text
+rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) − (h_est + h_int)·(T − T_amb)
+```
 
-2. **Skin depth** — il riscaldamento a induzione è approssimato come riscaldamento
-   volumetrico che decade esponenzialmente con la profondità `z`:
-   `q_vol(x, z) = q(x) · exp(-z/δ) / (δ·(1 - exp(-spessore/δ)))`, normalizzato in modo
-   da conservare il flusso superficiale. La skin depth `δ = √(2ρₑ/(ωμ))` è calcolata
-   dalla resistività elettrica e dalla permeabilità del materiale alla frequenza di
-   induzione, oppure può essere imposta con `skin_depth_fissa_m`. Per la banda
-   stagnata a 20 kHz risulta ≈ 0.1 mm, confrontabile con lo spessore: la parete è
-   quasi isoterma attraverso lo spessore.
+con `t` = spessore della parete e `∇²` il laplaciano nel piano `(x, s)`.
 
-3. **Diffusione 2D del calore** — l'equazione del calore è integrata nella sezione
-   `(x, z)` con volumi finiti ed Eulero implicito (incondizionatamente stabile).
+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.
 
-4. **Scambi con l'esterno** — la sezione scambia calore con l'ambiente su tutto il
-   contorno:
-   - convezione sul lato esterno (`h_esterno_W_m2K`), sul lato interno
-     (`h_interno_W_m2K`) e sui due bordi in x (`h_bordi_W_m2K`);
-   - **conduzione circonferenziale**: la sezione cede calore per conduzione lungo `y`
-     al resto della fascetta, assunto a temperatura ambiente. Il termine è
-     un'equazione di aletta ricavata sull'intero volume del cilindro: il calore
-     conduce lungo `y` attraverso l'intero spessore mentre le superfici esterna e
-     interna dell'intero cilindro perdono calore per convezione, dando
-     `q_y = -(h_esterno + h_interno)/spessore · (T - T_amb)` — nessun parametro di
-     conduzione `y` aggiuntivo da configurare.
+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`.
 
-5. **Temperatura iniziale** — il campo parte uniformemente alla temperatura ambiente
-   del run (che è randomizzata, quindi varia run per run).
+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).
 
-6. **Sensore realistico** — la lettura aggiunge alla temperatura vera della superficie
-   interna: inerzia del primo ordine (`costante_tempo_s`), rumore gaussiano
-   (`rumore_std_C`) e quantizzazione (`quantizzazione_C`).
+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).
 
 ## Metodo numerico
 
-- Griglia a volumi finiti `n_nodi_x × n_nodi_z` (default 100 × 15); le incognite sono
-  i centri cella.
-- Eulero implicito con passo `dt_interno_s` (default 0.2 ms), più fine del periodo di
-  campionamento CSV.
-- Tutti i termini (diffusione, convezione, conduzione circonferenziale) sono lineari e
-  costanti nel run: la matrice sparsa viene costruita e **fattorizzata LU una sola
-  volta per run** (`scipy.sparse.linalg.splu`); ogni passo temporale risolve solo i
-  sistemi triangolari. Un run da 30 s simulati richiede ~20 s di calcolo.
+- **Mesh** (`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.
 - Se in futuro si introducessero proprietà dipendenti dalla temperatura, la matrice
   andrebbe ricostruita e rifattorizzata a ogni passo.
 
 ## File
 
 ```text
-config.py           tutti i parametri di simulazione
-materials.py        proprietà termofisiche ed elettriche dei materiali
-simulate.py         motore fisico + generazione dei CSV
-plot_csv.py         grafici rapidi (temperature e flusso) del primo run
-plot_animazione.py  animazione della sezione: campo T(x,z), sorgenti, sensore
-dataset/            output generato da simulate.py (ricreato a ogni esecuzione)
+main.py                 punto di ingresso da riga di comando
+config.py               tutti i parametri (geometria, mesh, FEM, aria, sorgente, sensore)
+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)
 ```
 
 ## Installazione
@@ -121,38 +143,67 @@ pip install -r requirements.txt
 ## Uso
 
 ```bash
-# 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 sezione durante il passaggio delle sorgenti
-python plot_animazione.py
+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/
 ```
 
-Gli script di visualizzazione aprono finestre interattive (backend Qt); se il backend
-non è interattivo (es. sessione senza display) salvano automaticamente PNG/GIF in
-`dataset/`.
+Ogni modulo resta eseguibile anche direttamente (`python plot_mesh.py`,
+`python plot_animazione_fem.py`, `python esporta_csv.py`).
 
-L'animazione riproduce esattamente la fisica di `run_0001` (stesso seed) e mostra tre
-pannelli allineati: il profilo di flusso `q(x)` con le sorgenti in transito, il campo
-di temperatura nella sezione con il sensore IR, e la temperatura nel punto osservato
-(vera e con inerzia del sensore). Finestra temporale e cadenza dei fotogrammi si
-regolano con le costanti in testa a `plot_animazione.py`.
+Gli script aprono finestre interattive (backend Qt); se il backend non è interattivo
+(es. sessione senza display) salvano automaticamente PNG/GIF in `output/`.
+
+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"]`.
+
+## Output CSV
+
+`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.
+
+### `output/csv/run_XXXX.csv` — serie temporale
+
+| 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                                            |
+
+### `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`.
 
 ## Configurazione
 
 Tutto si modifica in `config.py`. I dizionari principali:
 
-| Dizionario       | Contenuto                                                                 |
-|------------------|---------------------------------------------------------------------------|
-| `SIMULAZIONE`    | numero di run, durata, campionamento CSV, passo interno, seed, cartella   |
-| `FASCETTA`       | diametro, lunghezza, spessore, griglia, conduzione circonferenziale, materiale |
-| `ARIA`           | temperatura ambiente e coefficienti di convezione dei quattro lati        |
-| `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                                        |
+| 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                          |
 
 Punti da conoscere:
 
@@ -163,54 +214,14 @@ 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à, densità, calore specifico, resistività elettrica, permeabilità).
+  (conducibilità termica, densità, calore specifico, resistività elettrica,
+  permeabilità relativa).
 - **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_lato_caldo_C`               | temperatura della superficie esterna nello stesso punto x          |
-| `T_ambiente_C`                 | temperatura ambiente del run                                       |
-| `velocita_m_s`, `sigma_punto_m`, `flusso_picco_W_m2` | parametri randomizzati del run              |
-| `materiale`                    | chiave del materiale                                               |
-
-### `dataset/metadata.csv` — una riga per run
-
-Contiene tutti i parametri effettivi del run (geometria, griglia, coefficienti di
-scambio, parametri delle sorgenti e del sensore, valori randomizzati) e le temperature
-di picco vera e misurata: utile come ground truth e per filtrare i run.
-
-## Limitazioni
-
-1. Il campo elettromagnetico non è simulato: l'accoppiamento induttivo è ridotto a
-   impronta gaussiana × efficienza × decadimento esponenziale in z.
-2. La coordinata circonferenziale y non è risolta: offset del percorso e conduzione
-   verso il resto della fascetta sono modelli collassati (attenuazione gaussiana e
-   scambio lineare verso T ambiente).
-3. Le proprietà dei materiali sono costanti con la temperatura; per gli acciai
-   ferromagnetici la skin depth reale varia fortemente con temperatura e campo
-   (punto di Curie non modellato).
-4. Irraggiamento non modellato: a ~220 °C le perdite radiative non sono del tutto
-   trascurabili rispetto alla convezione.

config.py

@@ -1,30 +1,21 @@
-# Configurazione per il simulatore termico 2D della sezione di una fascetta.
+# Configurazione dell'analisi termica a elementi finiti della fascetta.
 #
 # Geometria e modello fisico:
 # - La fascetta è un anello cilindrico con diametro "diametro_mm", spessore
 #   "spessore_mm" e lunghezza "lunghezza_mm".
-# - Il dominio simulato è la sezione rettangolare lunghezza × spessore.
-# - Sistema di coordinate: origine (0, 0) nel vertice in alto a sinistra
-#   della sezione. x = direzione della lunghezza (da 0 a lunghezza),
-#   z = direzione dello spessore (0 = lato esterno, dove agiscono le
-#   sorgenti; spessore = lato interno, osservato dal sensore).
-# - y è la coordinata circonferenziale: non è risolta spazialmente, l'offset
-#   y del percorso delle sorgenti è collassato in un'attenuazione gaussiana
-#   del flusso.
-# - Le sorgenti a induzione si muovono in direzione -x sul lato esterno.
-# - Il riscaldamento a induzione è approssimato come riscaldamento volumetrico
-#   che decade esponenzialmente con la profondità z secondo la skin depth.
-# - La sezione scambia per convezione con l'aria su tutti e quattro i lati
-#   (esterno, interno e i due bordi in x). Scambia inoltre per conduzione
-#   lungo y con il resto della fascetta, assunto a temperatura ambiente:
-#   il calore conduce attraverso l'intero volume dello spessore mentre le
-#   superfici esterna e interna dell'intero cilindro perdono calore per
-#   convezione (equazione dell'aletta), derivato da h_esterno, h_interno e
-#   spessore_mm senza parametri di conduzione y aggiuntivi.
+# - 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 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".
+#   contatto la temperatura della superficie interna nel punto x = "x_mm",
+#   sul piano circonferenziale theta = 0.
 #
 # Unità di misura:
 # - lunghezza: m (mm dove indicato dal suffisso)
@@ -33,32 +24,6 @@
 # - 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.0002,
-
-    # 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,
@@ -69,28 +34,49 @@ FASCETTA = {
     # Spessore della parete [mm].
     "spessore_mm": 0.18,
 
-    # Numero di celle del volume finito lungo x (lunghezza).
-    "n_nodi_x": 100,
-
-    # Numero di celle del volume finito lungo z (spessore).
-    # Più nodi = maggiore risoluzione spaziale, simulazione più lenta.
-    "n_nodi_z": 15,
-
     # 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 la circonferenza.
+    # La mesh è chiusa su se stessa: non c'è una riga di nodi duplicata.
+    "n_elementi_circonferenza": 48,
+}
+
+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,
+}
+
 ARIA = {
     # Temperatura dell'aria ambiente.
     "temperatura_ambiente_C": 25.0,
 
-    # Coefficiente di convezione sul lato esterno (z = 0, lato sorgenti).
+    # Coefficiente di convezione sulla faccia esterna (lato sorgenti).
     "h_esterno_W_m2K": 12.0,
 
-    # Coefficiente di convezione sul lato interno (z = spessore, lato sensore).
+    # Coefficiente di convezione sulla faccia interna (lato sensore).
     "h_interno_W_m2K": 8.0,
 
-    # Coefficiente di convezione sui bordi laterali (x = 0 e x = lunghezza).
+    # Coefficiente di convezione sui bordi anulari (x = 0 e x = lunghezza).
     "h_bordi_W_m2K": 10.0,
 }
 
@@ -127,7 +113,8 @@ SORGENTE = {
     # Frazione del flusso incidente che diventa effettivamente calore nella fascetta.
     "efficienza_riscaldamento": 0.35,
 
-    # Frequenza di induzione usata per stimare la skin depth se skin_depth_fissa_m è None.
+    # Frequenza di induzione dell'induttore, usata per stimare la skin depth
+    # se skin_depth_fissa_m è None.
     "frequenza_hz": 20000.0,
 
     # Override della skin depth. Usare None per calcolarla dalle proprietà elettriche del materiale.
@@ -142,7 +129,7 @@ SENSORE = {
     # Coordinata x del punto della superficie interna osservato dal sensore [mm].
     "x_mm": 50.0,
 
-    # Distanza del sensore dalla parete interna lungo z [mm].
+    # Distanza del sensore dalla parete interna lungo lo spessore [mm].
     # Il sensore è a infrarossi: la distanza non influenza la misura,
     # è registrata solo come geometria del setup.
     "distanza_parete_mm": 10.0,
@@ -161,7 +148,7 @@ SENSORE = {
 }
 
 RANDOMIZZAZIONE = {
-    # Se abilitata, ogni run varia leggermente alcuni parametri.
+    # Se abilitata, ogni analisi varia leggermente alcuni parametri.
     "abilitata": False,
 
     # Deviazioni standard relative.
@@ -178,3 +165,9 @@ 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",
+}

esporta_csv.py

@@ -0,0 +1,153 @@
+# 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

@@ -0,0 +1,460 @@
+# 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

@@ -0,0 +1,61 @@
+# 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 il simulatore termico.
+# Database dei materiali per l'analisi termica.
 #
 # Tutte le unità sono SI:
 # - conducibilita_termica_W_mK

mesh.py

@@ -0,0 +1,114 @@
+# Generazione della mesh a elementi shell quadrangolari 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.
+#
+# 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
+#
+# La mesh è chiusa lungo la circonferenza: l'ultima colonna di elementi
+# richiude sulla prima riga di nodi, senza nodi duplicati sulla cucitura.
+
+import numpy as np
+
+from config import FASCETTA, MESH
+
+
+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
+
+
+def genera_mesh(fascetta: dict = FASCETTA, mesh: dict = MESH) -> dict:
+    """Costruisce la mesh shell quadrangolare della fascetta.
+
+    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."
+        )
+
+    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(),
+        )
+    )
+
+    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)
+
+    return {
+        "nodi": nodi,
+        "elementi": elementi,
+        "indice_nodo": indice_nodo,
+        "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,
+    }
+
+
+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"
+    )

plot_animazione.py

@@ -1,230 +0,0 @@
-# Animazione della sezione della fascetta durante il passaggio delle sorgenti.
-#
-# Riproduce la fisica di run_0001 (stesso seed di simulate.py) e mostra:
-# - il profilo di flusso termico q(x) sul lato esterno e le sorgenti in moto;
-# - il campo di temperatura T(x, z) nella sezione lunghezza × spessore;
-# - il sensore infrarosso e la temperatura nel punto osservato.
-
-import random
-from pathlib import Path
-
-import matplotlib
-import matplotlib.pyplot as plt
-import numpy as np
-from matplotlib.animation import FuncAnimation, PillowWriter
-
-from config import SIMULAZIONE
-from simulate import (
-    configurazione_randomizzata,
-    passo_implicito,
-    prepara_stato_termico,
-    profilo_flusso_incidente_W_m2,
-)
-
-# 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
-
-
-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.
-    fascetta = cfg_run["fascetta"]
-    aria = cfg_run["aria"]
-    sorgente = cfg_run["sorgente"]
-    sensore = cfg_run["sensore"]
-
-    stato = prepara_stato_termico(fascetta, aria, sorgente)
-    n_x = stato["n_x"]
-    n_z = stato["n_z"]
-    x_centri = stato["x_centri_m"]
-    dt = stato["dt_s"]
-
-    x_sensore = sensore["x_mm"] / 1000.0
-    i_sensore = min(n_x - 1, max(0, int(x_sensore / stato["dx_m"])))
-
-    T = np.full((n_x, n_z), aria["temperatura_ambiente_C"], dtype=float)
-    T_sensore = T[i_sensore, -1]
-    tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
-
-    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:
-        x_rif, q_x = profilo_flusso_incidente_W_m2(sorgente, x_sensore, t, x_centri)
-        T = passo_implicito(stato, T, q_x)
-
-        T_sensore += (T[i_sensore, -1] - T_sensore) * dt / tau_sensore
-
-        if t + 1e-12 >= prossimo_frame_t:
-            tempi.append(t)
-            campi.append(T.copy())
-            flussi.append(q_x.copy())
-            x_riferimenti.append(x_rif)
-            T_vere.append(T[i_sensore, -1])
-            T_lette.append(T_sensore)
-            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),
-        "x_centri_mm": x_centri * 1000.0,
-        "spessore_mm": fascetta["spessore_mm"],
-        "lunghezza_mm": fascetta["lunghezza_mm"],
-        "x_sensore_mm": sensore["x_mm"],
-        "sorgente": cfg_run["sorgente"],
-    }
-
-
-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"]
-    spessore_mm = dati["spessore_mm"]
-    x_vista_mm = (-10.0, lunghezza_mm + 10.0)
-
-    q_max_MW = max(q.max() for q in dati["flussi"]) / 1e6
-    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
-
-    fig, (ax_flusso, ax_sezione, ax_storia) = plt.subplots(
-        3, 1, figsize=(10, 8), height_ratios=[1.0, 1.6, 1.2],
-        gridspec_kw={"hspace": 0.45},
-    )
-    fig.suptitle("Sezione della fascetta: sorgenti in transito e sensore")
-
-    # Pannello 1: profilo di flusso sul lato esterno 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: campo di temperatura nella sezione (z verso il basso,
-    # origine nel vertice in alto a sinistra come nel modello).
-    immagine = ax_sezione.imshow(
-        dati["campi"][indice_inizio].T,
-        extent=(0.0, lunghezza_mm, spessore_mm, 0.0),
-        aspect="auto",
-        cmap="inferno",
-        vmin=cfg_run["aria"]["temperatura_ambiente_C"],
-        vmax=T_max,
-        interpolation="bilinear",
-    )
-    ax_sezione.set_xlim(*x_vista_mm)
-    ax_sezione.set_ylim(3.2 * spessore_mm, -0.6 * spessore_mm)
-    ax_sezione.set_ylabel("z [mm]")
-    ax_sezione.set_xlabel("x [mm]")
-    # Sensore infrarosso sotto la parete interna (posizione schematica,
-    # non in scala) con linea di vista tratteggiata.
-    x_sens = dati["x_sensore_mm"]
-    ax_sezione.plot([x_sens], [2.4 * spessore_mm], "^", color="tab:blue", markersize=12)
-    ax_sezione.plot(
-        [x_sens, x_sens], [1.1 * spessore_mm, 2.1 * spessore_mm],
-        linestyle="--", color="tab:blue", linewidth=1,
-    )
-    ax_sezione.text(
-        x_sens + 3, 2.4 * spessore_mm, "sensore IR",
-        color="tab:blue", va="center",
-    )
-    # La colorbar è agganciata a tutti i pannelli per non restringere solo
-    # quello della sezione, mantenendo allineati gli assi x.
-    barra = fig.colorbar(
-        immagine, ax=(ax_flusso, ax_sezione, 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_centri_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),
-        )
-
-        immagine.set_data(dati["campi"][k].T)
-
-        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, immagine,
-            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_sezione.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

@@ -1,150 +0,0 @@
-# Animazione 3D del barattolo in vista isometrica: colormap della temperatura
-# sulla superficie esterna.
-#
-# Il modello risolve solo la sezione x-z (vedi CLAUDE.md): la coordinata
-# circonferenziale y non è simulata, è collassata in un'attenuazione
-# gaussiana del flusso. Per estrudere il campo attorno alla circonferenza si
-# riusa la stessa gaussiana: la sovratemperatura rispetto al piano y=0 (dove
-# si trova il sensore) viene scalata in funzione della distanza
-# circonferenziale dal percorso delle sorgenti. È quindi una ricostruzione
-# visiva, non un calcolo di diffusione in y.
-
-import random
-
-import matplotlib
-import matplotlib.pyplot as plt
-import numpy as np
-from matplotlib import cm
-from matplotlib.animation import FuncAnimation, PillowWriter
-from pathlib import Path
-
-from config import FASCETTA, SIMULAZIONE
-from plot_animazione import DT_FRAME_S as DT_FRAME_S_BASE
-from plot_animazione import T_FINE_ANIMAZIONE_S, T_INIZIO_ANIMAZIONE_S, simula_campi
-from 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
-
-# Numero di punti lungo la circonferenza per l'estrusione della superficie.
-N_THETA = 72
-
-
-def attenuazione_circonferenziale(
-    y_m: np.ndarray, offset_y_m: float, sigma_m: float
-) -> np.ndarray:
-    # Fattore che scala la sovratemperatura (T - T_ambiente) del piano y=0
-    # in funzione della distanza circonferenziale y dal piano stesso,
-    # normalizzato in modo da valere 1 in y=0.
-    sigma = max(sigma_m, 1e-9)
-    esponente = -0.5 * ((y_m - offset_y_m) ** 2 - offset_y_m ** 2) / (sigma * sigma)
-    return np.exp(esponente)
-
-
-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))
-
-    T_ambiente = cfg_run["aria"]["temperatura_ambiente_C"]
-    sorgente = cfg_run["sorgente"]
-    sigma_m = sorgente["sigma_punto_m"]
-    offset_y_m = sorgente["offset_y_percorso_m"]
-    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
-    distanza_m = sorgente.get("distanza_sorgenti_m", 0.0)
-
-    raggio_m = (FASCETTA["diametro_mm"] / 1000.0) / 2.0
-    x_centri_m = dati["x_centri_mm"] / 1000.0
-    lunghezza_m = dati["lunghezza_mm"] / 1000.0
-
-    theta = np.linspace(-np.pi, np.pi, N_THETA)
-    y_circ_m = theta * raggio_m
-    attenuazione = attenuazione_circonferenziale(y_circ_m, offset_y_m, sigma_m)
-
-    Xm, Thetam = np.meshgrid(x_centri_m, theta)
-    Ym = raggio_m * np.sin(Thetam)
-    Zm = raggio_m * np.cos(Thetam)
-
-    T_max = max(c[:, 0].max() for c in dati["campi"])
-    # 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]")
-
-    def disegna_frame(k: int):
-        ax.cla()
-        ax.view_init(elev=35.264, azim=45)
-        ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m))
-        ax.set_axis_off()
-
-        T_lato_esterno = dati["campi"][k][:, 0]
-        T_superficie = T_ambiente + (T_lato_esterno[None, :] - T_ambiente) * attenuazione[:, None]
-        colori = cmap(norm(T_superficie))
-
-        ax.plot_surface(
-            Xm, Ym, Zm, facecolors=colori, rstride=1, cstride=1,
-            antialiased=False, shade=True, lightsource=lightsource, linewidth=0,
-        )
-
-        x_sorgenti_m = dati["x_riferimenti"][k] + np.arange(numero_sorgenti) * distanza_m
-        visibili = (x_sorgenti_m >= 0.0) & (x_sorgenti_m <= lunghezza_m)
-        if visibili.any():
-            theta_sorgente = offset_y_m / raggio_m
-            ax.scatter(
-                x_sorgenti_m[visibili],
-                np.full(int(visibili.sum()), raggio_m * np.sin(theta_sorgente) * 1.05),
-                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

@@ -0,0 +1,198 @@
+# 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

@@ -1,48 +0,0 @@
-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

@@ -0,0 +1,60 @@
+# 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()

simulate.py

@@ -1,486 +0,0 @@
-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 scipy.sparse as sp
-from scipy.sparse.linalg import splu
-
-from config import ARIA, FASCETTA, RANDOMIZZAZIONE, SENSORE, SIMULAZIONE, SORGENTE
-from materials import MATERIALI
-
-
-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 _spread_sorgenti_m(sorgente: dict) -> float:
-    # Distanza lungo x tra la prima e l'ultima sorgente del gruppo.
-    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
-    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
-    return (numero_sorgenti - 1) * distanza
-
-
-def _x_riferimento_iniziale_m(sorgente: dict, x_sensore_m: float) -> float:
-    # Posizione a t=0 della sorgente di indice 0. Le sorgenti i sono a
-    # x_i = riferimento + i * distanza, quindi per v >= 0 l'indice 0 è la
-    # più arretrata nel verso di marcia, per v < 0 è la più avanzata.
-    # x_inizio_m è la distanza dal sensore della sorgente più avanzata
-    # (quella che lo raggiunge per prima).
-    spread = _spread_sorgenti_m(sorgente)
-    x_inizio = sorgente["x_inizio_m"]
-    if sorgente["velocita_m_s"] >= 0:
-        return (x_sensore_m - x_inizio) - spread
-    return x_sensore_m + x_inizio
-
-
-def _x_riferimento_finale_m(sorgente: dict, x_sensore_m: float) -> float:
-    # Posizione di fine corsa della sorgente di indice 0. x_fine_m è la
-    # distanza dal sensore della sorgente più arretrata nel verso di marcia
-    # (quella che lo supera per ultima).
-    spread = _spread_sorgenti_m(sorgente)
-    x_fine = sorgente["x_fine_m"]
-    if sorgente["velocita_m_s"] >= 0:
-        return x_sensore_m + x_fine
-    return (x_sensore_m - x_fine) - spread
-
-
-def _intervallo_attivo(inizio: float, fine: float, v: float, x_m: float) -> bool:
-    if v >= 0:
-        return inizio <= x_m <= fine
-    return fine <= x_m <= inizio
-
-
-def profilo_flusso_incidente_W_m2(
-    sorgente: dict,
-    x_sensore_m: float,
-    t_s: float,
-    x_centri_m: np.ndarray,
-) -> tuple[float, np.ndarray]:
-    # Restituisce x_sorgente_m (posizione della sorgente di riferimento) e
-    # il profilo di flusso termico efficace q(x) [W/m²] sul lato esterno,
-    # somma dei contributi di tutte le sorgenti attive.
-    #
-    # Ogni sorgente in moto ha un'impronta gaussiana lungo x, valutata sui
-    # centri cella della sezione. L'offset circonferenziale y non è risolto
-    # spazialmente: entra come attenuazione gaussiana del flusso.
-    x_rif_iniziale = _x_riferimento_iniziale_m(sorgente, x_sensore_m)
-    x_rif_finale = _x_riferimento_finale_m(sorgente, x_sensore_m)
-    x_riferimento = x_rif_iniziale + sorgente["velocita_m_s"] * t_s
-
-    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
-    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
-    v = sorgente["velocita_m_s"]
-    zero_dopo_fine = sorgente.get("zero_dopo_fine", True)
-    dy = sorgente["offset_y_percorso_m"]
-    sigma = sorgente["sigma_punto_m"]
-
-    q_picco = sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"]
-    attenuazione_y = math.exp(-0.5 * (dy * dy) / (sigma * sigma))
-
-    q_x = np.zeros_like(x_centri_m)
-    for i in range(numero_sorgenti):
-        x_i = x_riferimento + i * distanza
-
-        if zero_dopo_fine:
-            inizio_i = x_rif_iniziale + i * distanza
-            fine_i = x_rif_finale + i * distanza
-            if not _intervallo_attivo(inizio_i, fine_i, v, x_i):
-                continue
-
-        dx = x_centri_m - x_i
-        q_x += q_picco * attenuazione_y * np.exp(-0.5 * (dx * dx) / (sigma * sigma))
-
-    return x_riferimento, q_x
-
-
-def profilo_deposizione_z_1_m(
-    z_centri_m: np.ndarray,
-    spessore_m: float,
-    skin_depth_m: float,
-) -> np.ndarray:
-    # Profilo di deposizione volumetrica del flusso superficiale [1/m]:
-    #
-    # p(z) = exp(-z / delta) / (delta * (1 - exp(-spessore / delta)))
-    #
-    # Normalizzato in modo che integrale_0^spessore p(z) dz = 1, così che
-    # q_vol(x, z) = q(x) * p(z) conservi il flusso superficiale.
-    delta = max(skin_depth_m, 1e-9)
-    normalizzazione = delta * (1.0 - math.exp(-spessore_m / delta))
-    return np.exp(-z_centri_m / delta) / normalizzazione
-
-
-def _laplaciano_1d(n: int) -> sp.spmatrix:
-    # Operatore alle differenze -T'' su n celle con bordi adiabatici (Neumann).
-    diagonale = np.full(n, 2.0)
-    diagonale[0] = 1.0
-    diagonale[-1] = 1.0
-    fuori = -np.ones(n - 1)
-    return sp.diags([fuori, diagonale, fuori], [-1, 0, 1])
-
-
-def costruisci_solutore_implicito_2d(
-    n_x: int,
-    n_z: int,
-    dt_s: float,
-    dx_m: float,
-    dz_m: float,
-    spessore_m: float,
-    materiale: dict,
-    h_esterno_W_m2K: float,
-    h_interno_W_m2K: float,
-    h_bordi_W_m2K: float,
-):
-    # Costruisce e fattorizza (LU sparsa) la matrice A per Eulero implicito 2D:
-    # A * T_next = rhs
-    #
-    # Le incognite sono i centri cella T[i, j] con i lungo x e j lungo z,
-    # appiattiti in ordine C (indice = i * n_z + j). Tutti e quattro i lati
-    # della sezione includono la convezione verso l'ambiente; su ogni cella
-    # agisce inoltre la conduzione circonferenziale (y) verso il resto della
-    # fascetta, assunto a temperatura ambiente. Il termine è ricavato
-    # considerando l'intero volume della fascia cilindrica (equazione
-    # dell'aletta): il calore conduce lungo y attraverso l'intero spessore
-    # mentre le superfici esterna e interna dell'intero cilindro scambiano
-    # per convezione, dando un sink distribuito uniformemente su ogni cella
-    # q_y = -(h_esterno + h_interno) / spessore * (T - T_amb).
-    k = materiale["conducibilita_termica_W_mK"]
-    rho = materiale["densita_kg_m3"]
-    cp = materiale["calore_specifico_J_kgK"]
-    alpha = k / (rho * cp)
-
-    r_x = alpha * dt_s / (dx_m * dx_m)
-    r_z = alpha * dt_s / (dz_m * dz_m)
-    b_esterno = h_esterno_W_m2K * dt_s / (rho * cp * dz_m)
-    b_interno = h_interno_W_m2K * dt_s / (rho * cp * dz_m)
-    b_bordo = h_bordi_W_m2K * dt_s / (rho * cp * dx_m)
-    c_y = (h_esterno_W_m2K + h_interno_W_m2K) * dt_s / (rho * cp * spessore_m)
-
-    n = n_x * n_z
-    scambio = np.full(n, c_y)
-    scambio[0::n_z] += b_esterno
-    scambio[n_z - 1::n_z] += b_interno
-    scambio[:n_z] += b_bordo
-    scambio[n - n_z:] += b_bordo
-
-    A = (
-        sp.identity(n)
-        + r_x * sp.kron(_laplaciano_1d(n_x), sp.identity(n_z))
-        + r_z * sp.kron(sp.identity(n_x), _laplaciano_1d(n_z))
-        + sp.diags(scambio)
-    )
-    return splu(sp.csc_matrix(A))
-
-
-def prepara_stato_termico(fascetta: dict, aria: dict, sorgente: dict) -> dict:
-    # Prepara griglia, coefficienti e solutore fattorizzato per un run:
-    # tutto ciò che resta costante durante l'integrazione temporale.
-    materiale = MATERIALI[fascetta["materiale"]]
-
-    lunghezza = fascetta["lunghezza_mm"] / 1000.0
-    spessore = fascetta["spessore_mm"] / 1000.0
-    n_x = fascetta["n_nodi_x"]
-    n_z = fascetta["n_nodi_z"]
-    dx = lunghezza / n_x
-    dz = spessore / n_z
-
-    dt = SIMULAZIONE["dt_interno_s"]
-
-    if sorgente["skin_depth_fissa_m"] is None:
-        skin_depth = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
-    else:
-        skin_depth = float(sorgente["skin_depth_fissa_m"])
-
-    rho = materiale["densita_kg_m3"]
-    cp = materiale["calore_specifico_J_kgK"]
-    alpha = materiale["conducibilita_termica_W_mK"] / (rho * cp)
-
-    z_centri = (np.arange(n_z) + 0.5) * dz
-
-    return {
-        "n_x": n_x,
-        "n_z": n_z,
-        "dx_m": dx,
-        "dz_m": dz,
-        "dt_s": dt,
-        "x_centri_m": (np.arange(n_x) + 0.5) * dx,
-        "z_centri_m": z_centri,
-        "skin_depth_m": skin_depth,
-        "rho": rho,
-        "cp": cp,
-        "profilo_z": profilo_deposizione_z_1_m(z_centri, spessore, skin_depth),
-        "b_esterno": aria["h_esterno_W_m2K"] * dt / (rho * cp * dz),
-        "b_interno": aria["h_interno_W_m2K"] * dt / (rho * cp * dz),
-        "b_bordo": aria["h_bordi_W_m2K"] * dt / (rho * cp * dx),
-        "c_y": (aria["h_esterno_W_m2K"] + aria["h_interno_W_m2K"]) * dt / (rho * cp * spessore),
-        "T_ambiente_C": aria["temperatura_ambiente_C"],
-        "solutore": costruisci_solutore_implicito_2d(
-            n_x=n_x,
-            n_z=n_z,
-            dt_s=dt,
-            dx_m=dx,
-            dz_m=dz,
-            spessore_m=spessore,
-            materiale=materiale,
-            h_esterno_W_m2K=aria["h_esterno_W_m2K"],
-            h_interno_W_m2K=aria["h_interno_W_m2K"],
-            h_bordi_W_m2K=aria["h_bordi_W_m2K"],
-        ),
-    }
-
-
-def passo_implicito(stato: dict, T: np.ndarray, q_x: np.ndarray) -> np.ndarray:
-    # Avanza il campo di temperatura di un passo dt: assembla il termine noto
-    # (sorgente volumetrica, convezione sui quattro lati, conduzione
-    # circonferenziale verso l'ambiente) e risolve il sistema implicito.
-    T_amb = stato["T_ambiente_C"]
-    rhs = T + (stato["dt_s"] / (stato["rho"] * stato["cp"])) * (
-        q_x[:, None] * stato["profilo_z"][None, :]
-    )
-    rhs += stato["c_y"] * T_amb
-    rhs[:, 0] += stato["b_esterno"] * T_amb
-    rhs[:, -1] += stato["b_interno"] * T_amb
-    rhs[0, :] += stato["b_bordo"] * T_amb
-    rhs[-1, :] += stato["b_bordo"] * T_amb
-    return stato["solutore"].solve(rhs.ravel()).reshape(T.shape)
-
-
-def quantizza(valore: float, passo: float) -> float:
-    if passo <= 0.0:
-        return valore
-    return round(valore / passo) * passo
-
-
-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)
-    n_x = stato["n_x"]
-    n_z = stato["n_z"]
-    x_centri = stato["x_centri_m"]
-    skin_depth = stato["skin_depth_m"]
-
-    x_sensore = sensore["x_mm"] / 1000.0
-    i_sensore = min(n_x - 1, max(0, int(x_sensore / stato["dx_m"])))
-
-    dt = stato["dt_s"]
-    durata = SIMULAZIONE["durata_s"]
-    periodo_campionamento = 1.0 / SIMULAZIONE["frequenza_campionamento_hz"]
-
-    T = np.full((n_x, n_z), aria["temperatura_ambiente_C"], dtype=float)
-    T_sensore = T[i_sensore, -1]
-
-    prossimo_campione_t = 0.0
-    T_vera_max = T[i_sensore, -1]
-    T_misurata_max = T_sensore
-
-    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_lato_caldo_C",
-            "T_ambiente_C",
-            "velocita_m_s",
-            "sigma_punto_m",
-            "flusso_picco_W_m2",
-            "materiale",
-        ])
-
-        t = 0.0
-        while t <= durata + 1e-12:
-            x_sorgente, q_x = profilo_flusso_incidente_W_m2(
-                sorgente, x_sensore, t, x_centri
-            )
-            T = passo_implicito(stato, T, q_x)
-
-            # Temperatura vera della superficie interna nel punto osservato
-            # dal sensore infrarosso.
-            T_vera_lato_sensore = T[i_sensore, -1]
-
-            # Inerzia del sensore del primo ordine.
-            tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
-            T_sensore += (T_vera_lato_sensore - T_sensore) * dt / tau_sensore
-
-            # Campionamento CSV.
-            if t + 1e-12 >= prossimo_campione_t:
-                misurata = T_sensore + rng.gauss(0.0, sensore["rumore_std_C"])
-                misurata = quantizza(misurata, sensore["quantizzazione_C"])
-
-                T_vera_max = max(T_vera_max, T_vera_lato_sensore)
-                T_misurata_max = max(T_misurata_max, misurata)
-
-                writer.writerow([
-                    cfg_run["id_run"],
-                    f"{t:.6f}",
-                    f"{x_sorgente:.9f}",
-                    f"{sorgente['offset_y_percorso_m']:.9f}",
-                    f"{q_x[i_sensore]:.6f}",
-                    f"{skin_depth:.9e}",
-                    f"{T_vera_lato_sensore:.6f}",
-                    f"{misurata:.6f}",
-                    f"{T[i_sensore, 0]:.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_m": fascetta["diametro_mm"] / 1000.0,
-        "lunghezza_m": fascetta["lunghezza_mm"] / 1000.0,
-        "spessore_m": fascetta["spessore_mm"] / 1000.0,
-        "n_nodi_x": n_x,
-        "n_nodi_z": n_z,
-        "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": x_sensore,
-        "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_lato_sensore_C": T_vera_max,
-        "T_misurata_max_sensore_C": T_misurata_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()