temp

2026-08-03 00:18:00 +02:00
21 changed files with 2095 additions and 1253 deletions
@@ -1,3 +1,3 @@
 .venv/
 __pycache__/
-output/
+dataset/
@@ -18,39 +18,64 @@ source .venv/bin/activate
 # Installare le dipendenze
 pip install -r requirements.txt
-# Punto di ingresso unico: `python main.py` elenca le azioni disponibili
+# Generare il dataset (scrive dataset/run_XXXX.csv + dataset/metadata.csv)
-python main.py mesh      # disegna la sola mesh a elementi shell
+python simulate.py
 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`).
+# Visualizzare il primo run
 python plot_csv.py
 # Animazione della superficie sviluppata e vista 3D isometrica
 python plot_animazione.py
 python plot_animazione_3d.py
 ```
 Attivare sempre il venv (`source .venv/bin/activate`) prima di eseguire qualsiasi comando Python.
-Non sono configurati test o linter.
+Non sono configurati linter. La validazione del solutore si esegue a mano:
 ```bash
 python valida_shell.py
 ```
 Va rieseguita dopo ogni modifica alla fisica o all'assemblaggio: confronta le matrici e il solutore con soluzioni analitiche (raffreddamento convettivo, modi sinusoidali circonferenziali, conservazione dell'energia, potenza della sorgente, periodicità) e stampa gli errori.
 ## Architettura
-Analisi termica transitoria a elementi finiti di una fascetta (anello cilindrico sottile) riscaldata da sorgenti a induzione in movimento.
+Generatore di dataset per misurazioni termiche pseudo-realistiche di una fascetta (anello cilindrico sottile) riscaldata da sorgenti a induzione in movimento.
-**Geometria:** la fascetta ha diametro, spessore e lunghezza configurabili. Il dominio discretizzato è la superficie media del cilindro: x = asse (le sorgenti si muovono lungo x sulla superficie esterna), theta = coordinata circonferenziale con theta = 0 sul piano del sensore, s = R·theta la lunghezza d'arco. Lo spessore non è discretizzato. Il sensore è un pirometro a infrarossi dentro la fascetta che misura la parete in un punto x fisso sul piano theta = 0.
+**Geometria:** la fascetta ha diametro medio, spessore e lunghezza configurabili. Il dominio simulato è l'intera superficie cilindrica, parametrizzata da x (lunghezza, le sorgenti si muovono in direzione -x) e theta (circonferenza, periodica). Lo spessore non è discretizzato: la parete è termicamente sottile, quindi faccia esterna e faccia interna hanno la stessa temperatura. Il sensore è un pirometro a infrarossi dentro la fascetta che misura la superficie interna in un punto x fisso, sulla generatrice theta = 0; l'offset circonferenziale delle sorgenti è misurato rispetto a quella generatrice.
 **Flusso dei dati:**
-1. `main.py` — dispatcher da riga di comando che seleziona l'azione (mesh, fem, csv)
+1. `config.py` — tutti i parametri configurabili (dizionari SIMULAZIONE, FASCETTA, ARIA, SORGENTE, SENSORE, RANDOMIZZAZIONE)
-2. `config.py` — tutti i parametri configurabili (dizionari FASCETTA, MESH, FEM, ARIA, SORGENTE, SENSORE, RANDOMIZZAZIONE, USCITA)
+2. `materials.py` — dizionario MATERIALI con proprietà termofisiche ed elettriche per materiale
-3. `materials.py` — dizionario MATERIALI con proprietà termofisiche ed elettriche per materiale
+3. `mesh.py` — nodi, connettività Q4 e periodicità circonferenziale della shell cilindrica
-4. `mesh.py` — mesh a elementi shell quadrangolari della superficie media del cilindro, visualizzata da `plot_mesh.py`
+4. `elementi_shell.py` — funzioni di forma bilineari, quadratura di Gauss 2×2, matrici del singolo elemento
-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`
+5. `assemblaggio.py` — assemblaggio delle matrici globali e del vettore della sorgente
 6. `sorgente.py` — moto del gruppo di sorgenti e distribuzione gaussiana del flusso
 7. `solutore.py` — Eulero implicito e fattorizzazione LU sparsa
 8. `sensore.py` — interpolazione della temperatura, inerzia, rumore, quantizzazione
 9. `simulate.py` — motore principale: genera N run randomizzati, scrive i CSV, scrive `metadata.csv`
 10. `plot_csv.py`, `plot_animazione.py`, `plot_animazione_3d.py` — visualizzazioni autonome
-**Mesh shell (`mesh.py`):** griglia strutturata `n_elementi_x × n_elementi_circonferenza` di quadrilateri a 4 nodi sulla superficie media (raggio = (diametro − spessore)/2); lo spessore non è discretizzato, è un attributo della shell. La mesh è chiusa lungo la circonferenza (nessun nodo duplicato sulla cucitura) e usa coordinate globali x = asse, y = R·sin(theta), z = R·cos(theta), con theta = 0 sul piano del sensore.
+**Pipeline fisica in [simulate.py](simulate.py):**
-**Pipeline FEM dentro `fem.py`:** equazione risolta per unità di superficie media, `rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) − (h_est + h_int)(T − T_amb)`, con convezione `h_bordi` sui due bordi anulari (area = spessore × perimetro). Il collasso in spessore è lecito perché Biot ≈ 1e-7 e il tempo di diffusione nella parete è di pochi ms: la profondità di penetrazione del riscaldamento non entra nel modello, conta solo il flusso assorbito. La skin depth è comunque calcolata dalla frequenza dell'induttore (`calcola_skin_depth_m`) e riportata in `stato["skin_depth_m"]` come sola diagnostica. L'impronta della sorgente è una gaussiana isotropa in (x, s) con distanza circonferenziale valutata sull'immagine più vicina (superficie chiusa). Gli elementi sono bilineari a 4 nodi e tutti identici (rettangoli dx × ds), quindi le matrici di elemento (`MASSA_RIF`, `RIGIDEZZA_X`, `RIGIDEZZA_S`) sono calcolate una volta e assemblate in forma vettorizzata; la matrice di massa di superficie è riusata per capacità, convezione e carico della sorgente. Eulero implicito con matrice costante fattorizzata LU una volta per analisi (`prepara_stato_fem`), poi `passo_implicito_fem` assembla solo il termine noto. In `fem.py` stanno anche la cinematica del gruppo di sorgenti (`_x_riferimento_iniziale_m`, `_x_riferimento_finale_m`, `_intervallo_attivo`) e la lettura della configurazione (`configurazione_randomizzata`).
+- `prepara_stato_termico` costruisce la mesh, assembla C, K_cond e K_conv, fattorizza LU una volta per run e prepara l'assemblatore della sorgente e l'interpolatore del sensore
 - La conduzione assiale e quella circonferenziale nascono dall'assemblaggio: elementi adiacenti condividono i nodi del bordo comune, quindi non serve alcun accoppiamento esplicito tra elementi
 - Convezione su entrambe le facce (h_esterno + h_interno su tutta l'area) e sui due bordi assiali (h_bordi · spessore lungo il perimetro), entrambe incorporate in K_conv e nel termine noto costante `f_ambiente`
 - Le sorgenti gaussiane in movimento producono q''(x, theta, t); la distanza circonferenziale è riavvolta su 2πR, quindi la sorgente scalda anche i punti raggiunti dall'altro verso. Il vettore di carico è ∫ N^T q'' dA, integrato con la stessa quadratura di Gauss delle matrici
 - Il flusso è separabile in un fattore assiale variabile nel tempo e uno circonferenziale costante: `assembla_sorgente` sfrutta la separazione, non è un'approssimazione
 - `passo_termico` avanza il campo risolvendo (C/dt + K) T = C/dt · T_prec + f_ambiente + f_sorgente
 - La temperatura iniziale del campo è la temperatura ambiente (randomizzata per run)
 - L'output del sensore aggiunge inerzia del primo ordine, rumore gaussiano e quantizzazione
-**Sensore e randomizzazione:** il campo nodale è la temperatura vera della parete; `simula_campo_fem` restituisce anche `T_lette`, la lettura del sensore degradata da inerzia del primo ordine (a ogni passo dt), rumore gaussiano e quantizzazione (a ogni campionamento di frame). `configurazione_randomizzata` perturba velocità, flusso di picco, sigma dello spot, offset y, temperatura ambiente e rumore del sensore con un RNG con seed `FEM["seed"]`, così due analisi restano riproducibili ma diverse; con `RANDOMIZZAZIONE["abilitata"] = False` restituisce i valori nominali di `config.py`.
+La skin depth (`calcola_skin_depth_m`) e il numero di Fourier sullo spessore (`numero_fourier_spessore`) sono calcolati e registrati in `metadata.csv` come diagnostiche: non entrano nella soluzione, servono a verificare che l'ipotesi di parete termicamente sottile regga.
-**Export CSV (`esporta_csv.py`):** `FEM["num_run"]` analisi (una per file `output/csv/run_XXXX.csv`), righe campionate a `FEM["frequenza_campionamento_hz"]` indipendentemente da `FEM["dt_s"]`. Colonne: `id_run, tempo_s, x_sorgente_m, T_vera_sensore_C, T_misurata_sensore_C, T_massima_fascetta_C, T_ambiente_C, offset_y_sorgente_m, velocita_m_s, sigma_punto_m, flusso_picco_W_m2, skin_depth_m, materiale`. `metadata.csv` ha una riga per analisi con tutti i parametri e i picchi. L'export chiama `simula_campo_fem(salva_campi=False)`: usa solo le serie scalari, i campi nodali non vengono accumulati. La cartella `output/csv` è cancellata e ricreata a ogni esecuzione (`shutil.rmtree`).
+`prepara_stato_termico` e `passo_termico` sono condivisi con `plot_animazione.py`, che riproduce la fisica di run_0001 per animare la superficie, e da lì con `plot_animazione_3d.py`: ogni modifica alla fisica va fatta lì, non duplicata.
 **Randomizzazione per run** (`configurazione_randomizzata`): ogni run perturba velocità, flusso di picco, sigma del punto, offset circonferenziale, temperatura ambiente e rumore del sensore con estrazioni gaussiane/uniformi da un RNG con seed fisso, garantendo riproducibilità.
 **Schema di output** (`dataset/run_XXXX.csv`): serie temporale con colonne `id_run, tempo_s, x_sorgente_m, offset_y_sorgente_m, flusso_termico_sorgente_W_m2, skin_depth_m, T_vera_lato_sensore_C, T_misurata_sensore_C, T_max_fascetta_C, T_media_fascetta_C, T_ambiente_C, velocita_m_s, sigma_punto_m, flusso_picco_W_m2, materiale`. `metadata.csv` ha una riga per run con tutti i parametri, le dimensioni della mesh, le diagnostiche e le temperature di picco.
 ## Convenzioni su `config.py`
@@ -58,8 +83,11 @@ Ogni parametro in [config.py](config.py) ha un commento che spiega solo cos'è (
 ## Vincoli progettuali chiave
- Il modello è una shell: la circonferenza è risolta, lo spessore è collassato (temperatura uniforme nella parete). Non introdurre effetti attraverso lo spessore (skin depth, gradiente radiale) senza prima discretizzarlo.
+- Il modello è una shell con una temperatura per nodo: lo spessore non è discretizzato e la temperatura della faccia esterna coincide con quella della faccia interna. L'ipotesi regge finché `numero_fourier_spessore` (α · σ/|v| / spessore²) resta molto maggiore di 1; sotto ~1 servirebbe una formulazione multistrato. Con i valori di default vale circa 5.4, e il gradiente misurato nello spessore col vecchio modello 2D era di 0.002 °C al picco.
 - La mesh cilindrica è strutturata e uniforme, quindi tutti gli elementi sono lo stesso rettangolo: le matrici elementari si calcolano una volta e si replicano. Se si introducessero mesh non uniformi o spessore variabile, `assemblaggio.py` dovrebbe ciclare sugli elementi.
 - La periodicità circonferenziale è strutturale, non un vincolo imposto: non esiste un nodo separato per theta = 2π, la connettività riavvolge l'ultimo elemento sul primo.
 - Le posizioni di inizio/fine corsa delle sorgenti (`x_inizio_m`, `x_fine_m`) sono distanze dal punto x del sensore lungo il verso di marcia; il segno di `velocita_m_s` determina il verso (negativo = -x).
- La matrice implicita è costruita e fattorizzata una volta sola (proprietà del materiale costanti, nessun coefficiente dipendente dalla temperatura). Se si aggiungono proprietà dipendenti dalla temperatura, la matrice deve essere ricostruita e rifattorizzata ad ogni passo temporale.
+- 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.
- Le visualizzazioni salvano in `USCITA["cartella"]` solo quando il backend matplotlib non è interattivo; altrimenti aprono una finestra.
+- La matrice è simmetrica definita positiva: `splu` usa `permc_spec="MMD_AT_PLUS_A"`, che dimezza circa il riempimento rispetto al COLAMD di default.
 - `simulate.py` cancella e ricrea l'intera cartella di output ad ogni esecuzione (`shutil.rmtree`).
 - Aggiungere un nuovo materiale richiede solo una nuova voce nel dizionario `MATERIALI` in [materials.py](materials.py); la chiave del materiale va poi impostata in `FASCETTA["materiale"]` in [config.py](config.py).
@@ -1,135 +1,199 @@
-# Analisi Termica FEM — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso
+# Simulatore Termico a Elementi Shell — Fascetta, Sorgenti a Induzione Mobili, Sensore IR Fisso
-Questo progetto calcola il transitorio termico di una fascetta (anello cilindrico
+Questo progetto genera misurazioni CSV pseudo-realistiche della temperatura di una
-sottile) riscaldata da un gruppo di sorgenti a induzione in movimento, osservata da un
+fascetta (anello cilindrico sottile) riscaldata da un gruppo di sorgenti a induzione
-sensore a infrarossi fisso. La fascetta è discretizzata con una mesh a elementi shell
+in movimento, osservata da un sensore a infrarossi fisso. Lo scopo è produrre dataset
-quadrangolari e il campo di temperatura è risolto a elementi finiti, poi visualizzato in
+per l'addestramento e la validazione di modelli di stima/regressione termica.
 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 esterno del cilindro;
+- **diametro medio** (default 70 mm) — il diametro della superficie media, su cui
  giace la mesh;
 - **spessore** (default 0.18 mm) — lo spessore della parete;
 - **lunghezza** (default 100 mm) — l'estensione assiale lungo `x`.
-Il dominio discretizzato è la **superficie media** del cilindro (raggio =
+Il dominio simulato è l'**intera superficie cilindrica**, discretizzata con elementi
-(diametro − spessore)/2). Lo spessore non è discretizzato: è un attributo degli
+shell quadrilateri. La superficie è parametrizzata da `(x, θ)`:
 elementi shell.
 ```text
-        sorgenti (induttori), in moto lungo x sulla superficie esterna
+                 sorgenti (induttori), in moto verso -x
-              ▼           ▼           ▼
+                       ▼           ▼           ▼
-              ┌───────────────────────────────────┐
+                  ╭────────────────────────────────╮
-        x = 0 │   superficie media del cilindro   │ x = lunghezza
+     θ = 0  ─ ─ ─ ┤   ●  ← punto osservato dal     ├ ─ ─ ─  generatrice del sensore
-              └───────────────────────────────────┘
+                  │      sensore IR (x = 50 mm)    │
-                         ┆  linea di vista
+                  ╰────────────────────────────────╯
-                         ▲
+                  x = 0                     x = lunghezza
-                    sensore IR (fisso, x = 50 mm, theta = 0,
+
-                    a 10 mm dalla parete interna)
+     superficie sviluppata in piano; il bordo superiore e quello inferiore
     sono lo stesso punto della fascetta (θ = 0 ≡ θ = 2π)
 ```
- **x** = asse del cilindro, da `0` a `lunghezza`.
+- **x** = direzione della lunghezza, da `0` a `lunghezza`. Le sorgenti viaggiano in
- **theta** = coordinata circonferenziale, con `theta = 0` sul piano del sensore;
+  direzione `-x` sulla superficie esterna.
-  nelle formule si usa la lunghezza d'arco `s = R·theta`.
+- **θ** = coordinata circonferenziale, **periodica** e risolta spazialmente. La
- Coordinate globali per la visualizzazione: `x`, `y = R·sin(theta)`,
+  generatrice `θ = 0` è quella del sensore; l'offset circonferenziale delle sorgenti
-  `z = R·cos(theta)`.
+  (`offset_y_percorso_m`, misurato lungo l'arco) è riferito ad essa.
 - Lo **spessore non è discretizzato**: la parete è termicamente sottile, quindi la
  faccia esterna e quella interna hanno la stessa temperatura.
 Il **sensore** è un pirometro a infrarossi posto all'interno della fascetta, a una
-distanza configurabile dalla parete interna (default 10 mm). Essendo senza contatto, la
+distanza configurabile dalla parete interna (default 10 mm). Essendo senza contatto,
-distanza non influenza la misura: il sensore legge la temperatura della parete nel nodo
+la distanza non influenza la misura: il sensore legge la temperatura della superficie
-più vicino a `x = 50 mm` sul piano `theta = 0`.
+interna nel punto `x` configurato (default 50 mm, al centro della lunghezza) sulla
 generatrice `θ = 0`.
 ## Modello fisico
-L'equazione risolta è quella del calore per unità di superficie media:
+Il campo elettromagnetico non è simulato: la parte termica è però una FEM vera. La
 catena di modellazione è la seguente.
 1. **Shell termica** — per una parete sottile con temperatura uniforme nello spessore
   l'equazione del calore sulla superficie cilindrica è
   ```text
   ρ cp s ∂T/∂t = k s ( ∂²T/∂x² + (1/R²) ∂²T/∂θ² ) + q'' − q_conv
   ```
   La conduzione **assiale** e quella **circonferenziale** sono entrambe risolte.
   Discretizzando con elementi quadrilateri a 4 nodi e una temperatura per nodo si
   ottiene `C Ṫ + K T = f(t)`, con
   ```text
   Cₑ      = ∫ ρ cp s Nᵀ N dA          (matrice consistente)
   Kₑ_cond = ∫ Bᵀ D B s dA             D = k·I  (materiale isotropo)
   ```
   La conduzione tra elementi adiacenti nasce dall'assemblaggio: due elementi
   confinanti condividono i nodi del bordo comune, quindi condividono gli stessi
   gradi di libertà termici. Non serve alcuna legge di scambio tra elementi dello
   stesso materiale.
 2. **Convezione** — nelle shell la faccia esterna e quella interna sono **entrambe**
   superfici fisiche esposte, quindi ogni elemento scambia su tutta la sua area:
   ```text
   Kₑ_conv = ∫ (h_est + h_int) Nᵀ N dA        fₑ_conv = ∫ (h_est + h_int) T∞ Nᵀ dA
   ```
   I due bordi assiali espongono solo lo spessore della lamiera:
   `Kₑ_bordo = ∫ h_b s Nᵀ N dl`. Con `s = 0.18 mm` su `100 mm` di lunghezza il
   contributo dei bordi vale ~0.2% dell'area ed è quindi marginale.
 3. **Sorgenti gaussiane in moto** — un gruppo di `numero_sorgenti` sorgenti
   equidistanti (`distanza_sorgenti_m`) si muove rigidamente a velocità costante lungo
   la generatrice individuata da `offset_y_percorso_m`. Il flusso sulla superficie
   esterna è
   ```text
   q''(x, θ, t) = q_picco · efficienza · Σᵢ exp( −((x − xᵢ(t))² + Δy²) / (2σ²) )
   Δy = R · wrap(θ − θ_s)
   ```
   La distanza circonferenziale è **riavvolta** su `2πR`: la sorgente scalda anche i
   punti raggiunti dall'altro verso della circonferenza. Il carico nodale è
   `fₑ_src = ∫ Nᵀ q'' dA`, integrato con la stessa quadratura di Gauss delle matrici:
   la posizione circonferenziale dell'induttore è quindi simulata davvero, non
   rappresentata da un coefficiente di attenuazione.
 4. **Skin depth** — `δ = √(2ρₑ/(ωμ))`, calcolata dalla resistività elettrica e dalla
   permeabilità del materiale, oppure imposta con `skin_depth_fissa_m`. Poiché la
   shell non discretizza lo spessore, la skin depth **non entra nella soluzione**: è
   registrata in `metadata.csv` come diagnostica. Per la banda stagnata a 20 kHz vale
   0.142 mm, confrontabile con lo spessore di 0.18 mm — cioè la deposizione è già
   quasi uniforme attraverso la parete.
 5. **Temperatura iniziale** — il campo parte uniformemente alla temperatura ambiente
   del run (che è randomizzata, quindi varia run per run).
 6. **Sensore realistico** — la temperatura nel punto osservato è interpolata con le
   funzioni di forma dell'elemento che lo contiene; la lettura vi aggiunge inerzia del
   primo ordine (`costante_tempo_s`), rumore gaussiano (`rumore_std_C`) e
   quantizzazione (`quantizzazione_C`).
 ### Perché una sola temperatura nello spessore
 La parete è termicamente sottile con ampio margine. Il numero di Fourier che confronta
 il tempo di diffusione attraverso lo spessore con il tempo di transito della sorgente,
 ```text
-rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) − (h_est + h_int)·(T − T_amb)
+Fo = α · (σ/|v|) / s²
 ```
-con `t` = spessore della parete e `∇²` il laplaciano nel piano `(x, s)`.
+vale **5.4** con i parametri di default (ed è riportato per ogni run in
 `metadata.csv`). Il modello 2D `(x, z)` precedente, che risolveva lo spessore con 15
 nodi, dava un salto esterno–interno di **0.002 °C al picco** su una sovratemperatura
 di ~200 °C, e 0.53 °C come massimo su tutto il run. Risolvere lo spessore non aggiunge
 informazione a questi parametri.
-1. **Collasso dello spessore** — la temperatura è uniforme attraverso la parete. Con
+Il criterio per sapere quando servirebbe una formulazione multistrato è `Fo ≈ 1`, cioè
-   spessore 0.18 mm il numero di Biot `h·t/k` vale ~1e-7 e il tempo di diffusione
+uno spessore critico `s* = √(α σ/|v|)`:
   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.
-2. **Sorgenti gaussiane in moto** — ogni sorgente ha un'impronta gaussiana isotropa nel
+| materiale                | Fo a 0.18 mm | spessore critico |
-   piano `(x, s)` di raggio `sigma_punto_m`. Un gruppo di `numero_sorgenti` sorgenti
+|--------------------------|--------------|------------------|
-   equidistanti (`distanza_sorgenti_m`) si muove rigidamente a velocità costante:
+| alluminio                | 31.3         | 1.01 mm          |
-   `q(x, s, t) = Σᵢ q_picco · efficienza · exp(-((x - xᵢ(t))² + (s - s₀)²) / (2σ²))`.
+| ferro                    | 8.4          | 0.52 mm          |
-   La distanza circonferenziale è valutata sull'immagine più vicina, perché la
+| banda stagnata / acciaio dolce | 5.4    | 0.42 mm          |
-   superficie è chiusa. L'offset `s₀` del percorso rispetto al piano del sensore è
+| inox 304                 | 1.5          | 0.22 mm          |
   `offset_y_percorso_m`.
 3. **Scambi con l'ambiente** — convezione sulla faccia esterna (`h_esterno_W_m2K`) e su
   quella interna (`h_interno_W_m2K`) su tutta la superficie, più convezione sui due
   bordi anulari `x = 0` e `x = lunghezza` (`h_bordi_W_m2K`, su un'area pari a
   spessore × perimetro).
 4. **Temperatura iniziale** — il campo parte uniformemente alla temperatura ambiente
   dell'analisi (randomizzabile).
 5. **Sensore reale** — la lettura aggiunge alla temperatura vera della parete nel punto
   osservato: inerzia del primo ordine (`costante_tempo_s`), rumore gaussiano
   (`rumore_std_C`) e quantizzazione (`quantizzazione_C`), gli ultimi due applicati a
   ogni campionamento.
 6. **Randomizzazione** — con `RANDOMIZZAZIONE["abilitata"]` ogni analisi perturba
   velocità, flusso di picco, sigma dello spot, offset y del percorso, temperatura
   ambiente e rumore del sensore con estrazioni gaussiane/uniformi da un RNG con seed
   fisso (`FEM["seed"]`): analisi diverse dello stesso setup differiscono come
   differirebbero due passaggi reali, restando riproducibili.
 Irraggiamento non modellato e proprietà dei materiali costanti con la temperatura: a
 ~210 °C le perdite radiative non sono del tutto trascurabili rispetto alla convezione, e
 per gli acciai ferromagnetici la skin depth reale varia fortemente con temperatura e
 campo (punto di Curie non modellato).
 ## Metodo numerico
- **Mesh** (`mesh.py`): griglia strutturata `n_elementi_x × n_elementi_circonferenza`
+- Mesh strutturata `n_elementi_x × n_elementi_theta` (default 100 × 110, cioè elementi
-  di quadrilateri a 4 nodi (default 40 × 48), chiusa lungo la circonferenza senza nodi
+  di 1.0 × 2.0 mm) per **11110 nodi**; la periodicità circonferenziale è strutturale,
-  duplicati sulla cucitura.
+  non un vincolo imposto: non esiste un nodo separato per `θ = 2π`.
- **Elementi**: shell bilineari a 4 nodi. Tutti gli elementi sono rettangoli identici
+- Quadratura di Gauss 2×2 per tutte le matrici elementari e per il vettore della
-  `dx × ds`, quindi le matrici di elemento (massa e rigidezza) sono calcolate una volta
+  sorgente. La mesh è uniforme, quindi tutti gli elementi sono lo stesso rettangolo:
-  sola e assemblate in forma vettorizzata. La matrice di massa di superficie è riusata
+  le matrici elementari si calcolano una volta e si replicano.
-  per capacità termica, convezione sulle facce e carico della sorgente; la convezione
+- Eulero implicito (incondizionatamente stabile) con passo `dt_interno_s`
-  sui bordi anulari usa una matrice di massa 1D sull'anello di nodi.
+  (default 1 ms), più fine del periodo di campionamento CSV.
- **Integrazione temporale**: Eulero implicito con passo `FEM["dt_s"]` (default 1 ms).
+- Tutti i termini sono lineari e costanti nel run: `A = C/dt + K_cond + K_conv` viene
-  Tutti i termini sono lineari e costanti, quindi la matrice di sistema è **assemblata e
+  **fattorizzata LU una sola volta per run** (`scipy.sparse.linalg.splu`); ogni passo
-  fattorizzata LU una sola volta** (`scipy.sparse.linalg.splu`) in `prepara_stato_fem`;
+  risolve solo i sistemi triangolari. `A` è simmetrica definita positiva, quindi si usa
-  ogni passo assembla solo il termine noto e risolve i sistemi triangolari.
+  l'ordinamento `MMD_AT_PLUS_A`, che dimezza circa il riempimento rispetto al COLAMD di
- **Flusso nodale**: l'impronta gaussiana isotropa è separabile,
+  default. Un run da 30 s simulati richiede ~65 s di calcolo; i run sono generati in
-  `q(x, s) = q_x(x) · q_s(s)`, e il fattore circonferenziale è costante nel tempo
+  parallelo su più processi.
  (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.
 ### Validazione
 `python valida_shell.py` verifica il solutore su casi con soluzione nota. Va rieseguito
 dopo ogni modifica alla fisica o all'assemblaggio. Risultati con la mesh di default:
 | Verifica                                        | Esito                                    |
 |-------------------------------------------------|------------------------------------------|
 | capacità totale = `ρ cp s A`                    | esatta a precisione macchina             |
 | somma righe di `K_cond` = 0 (T costante → flusso nullo) | 4·10⁻¹⁸                          |
 | simmetria e semidefinita positività di `K_cond` | verificate                               |
 | raffreddamento convettivo uniforme vs `exp(−t/τ)` | errore 0.011 °C, pari a `dt/(2τ)`      |
 | conservazione dell'energia (adiabatico)         | errore relativo 4·10⁻¹⁴                  |
 | `Σ f_src` vs potenza analitica `2πσ² q_max`     | rapporto 1.0000; **0.5000** con sorgente centrata sul bordo |
 | modi sinusoidali circonferenziali `n = 1, 2, 4` vs `exp(−α n²t/R²)` | errore 0.04% – 0.17%   |
 | convergenza al raffinamento (modo `n = 4`)      | 2.75% → 1.07% → 0.17% → 0.05% da 30 a 220 elementi |
 | simmetria della diffusione attorno a `θ = 0`    | 2·10⁻¹⁴                                  |
 L'ultima riga è la prova che la periodicità funziona davvero: un impulso su `θ = 0`
 diffonde in modo identico nei due versi solo se la connettività riavvolge la
 circonferenza.
 ## File
 ```text
-main.py                 punto di ingresso da riga di comando
+config.py               tutti i parametri di simulazione
 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
+mesh.py                 nodi, connettività Q4, periodicità circonferenziale
-fem.py                  solutore termico transitorio a elementi finiti
+elementi_shell.py       funzioni di forma, quadratura, matrici del singolo elemento
-plot_mesh.py            disegno della sola mesh
+assemblaggio.py         assemblaggio delle matrici globali e del vettore sorgente
-plot_animazione_fem.py  animazione 3D del campo di temperatura
+sorgente.py             moto degli induttori e distribuzione gaussiana del flusso
-esporta_csv.py          export CSV della serie temporale
+solutore.py             integrazione implicita e fattorizzazione sparsa
-output/                 immagini e GIF salvate quando il backend non è interattivo
+sensore.py              interpolazione, inerzia, rumore, quantizzazione
-output/csv/             CSV esportati (ricreata a ogni esecuzione)
+simulate.py             orchestrazione dei run + generazione dei CSV
 valida_shell.py         verifiche del solutore su casi con soluzione nota
 plot_csv.py             grafici rapidi (temperature e flusso) del primo run
 plot_animazione.py      animazione della superficie sviluppata T(x, θ)
 plot_animazione_3d.py   animazione 3D isometrica del cilindro
 dataset/                output generato da simulate.py (ricreato a ogni esecuzione)
 ```
 ## Installazione
@@ -143,67 +207,69 @@ pip install -r requirements.txt
 ## Uso
 ```bash
-python main.py          # elenco delle azioni disponibili
+# genera il dataset (ATTENZIONE: cancella e ricrea la cartella dataset/)
-python main.py mesh     # disegna la sola mesh a elementi shell
+python simulate.py
-python main.py fem      # integra il campo FEM e lo anima in 3D
+
-python main.py csv      # esporta la serie temporale in output/csv/
+# grafici statici del primo run (temperatura e flusso nel tempo)
 python plot_csv.py
 # animazione della superficie sviluppata durante il passaggio delle sorgenti
 python plot_animazione.py
 # animazione 3D isometrica del cilindro
 python plot_animazione_3d.py
 ```
-Ogni modulo resta eseguibile anche direttamente (`python plot_mesh.py`,
+Gli script di visualizzazione aprono finestre interattive (backend Qt); se il backend
-`python plot_animazione_fem.py`, `python esporta_csv.py`).
+non è interattivo (es. sessione senza display) salvano automaticamente PNG/GIF in
 `dataset/`.
-Gli script aprono finestre interattive (backend Qt); se il backend non è interattivo
+Entrambe le animazioni riproducono esattamente la fisica di `run_0001` (stesso seed).
-(es. sessione senza display) salvano automaticamente PNG/GIF in `output/`.
+Quella 2D mostra tre pannelli allineati: il profilo di flusso `q(x)` con le sorgenti in
 transito, il campo di temperatura sulla superficie sviluppata in piano (asse verticale
 = arco circonferenziale, con il sensore a quota 0), e la temperatura nel punto
 osservato (vera e con inerzia del sensore). Quella 3D disegna il campo calcolato sulla
 superficie del cilindro — non è più una ricostruzione, il modello risolve davvero la
 distribuzione circonferenziale. Finestra temporale e cadenza dei fotogrammi si regolano
 con le costanti in testa a `plot_animazione.py`.
-L'animazione FEM mostra due pannelli: la fascetta in vista isometrica con il campo di
+La mappa sviluppata usa `pcolormesh` con shading Gouraud, che interpola linearmente tra
-temperatura dipinto sugli elementi shell e i marker delle sorgenti in transito, e la
+i valori nodali: è esattamente ciò che fanno le funzioni di forma bilineari del Q4,
-storia della temperatura nel punto osservato dal sensore (valore nodale vero e lettura
+quindi si guarda la soluzione dell'elemento finito e non un'interpolazione grafica.
 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
+### Vedere gli elementi shell
-`python main.py csv` cancella e ricrea `output/csv/`, poi esporta `FEM["num_run"]` analisi
+Entrambe le animazioni hanno una costante `MOSTRA_MESH` che disegna i bordi degli
-(una per file `run_XXXX.csv`) più un `metadata.csv` con una riga per analisi. Le righe
+elementi. È `False` di default perché con la mesh di produzione gli elementi sono
-sono campionate a `FEM["frequenza_campionamento_hz"]`, indipendente dal passo di
+1 × 2 mm, più fitti dei pixel disponibili, e il reticolo coprirebbe il campo. Per
-integrazione.
+ispezionare la mesh servono due modifiche insieme:
-### `output/csv/run_XXXX.csv` — serie temporale
+```python
 # in plot_animazione.py e/o plot_animazione_3d.py
 MOSTRA_MESH = True
-| Colonna                 | Significato                                                     |
+# in config.py, dentro FASCETTA
-|-------------------------|-----------------------------------------------------------------|
+"n_elementi_x": 25,
-| `id_run`                | identificativo dell'analisi                                     |
+"n_elementi_theta": 28,
-| `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
+Così si vedono i singoli quadrilateri, in piano e avvolti sul cilindro. Attenzione: con
-
+28 elementi sulla circonferenza l'arco vale 7.9 mm contro un `sigma_punto_m` di 12 mm,
-Tutti i parametri effettivi (geometria, mesh, coefficienti di scambio, sorgenti, sensore)
+cioè meno di due elementi per sigma — la mesh grossolana serve a **guardare la
-e i valori di picco: `T_vera_max_sensore_C`, `T_misurata_max_sensore_C`,
+discretizzazione**, non a produrre risultati (vedi la tabella di convergenza sopra).
 `T_massima_fascetta_C`.
 ## Configurazione
 Tutto si modifica in `config.py`. I dizionari principali:
-| Dizionario  | Contenuto                                                                 |
+| Dizionario       | Contenuto                                                                 |
-|-------------|---------------------------------------------------------------------------|
+|------------------|---------------------------------------------------------------------------|
-| `FASCETTA`        | diametro, lunghezza, spessore, materiale                            |
+| `SIMULAZIONE`    | numero di run, durata, campionamento CSV, passo interno, seed, cartella   |
-| `MESH`            | numero di elementi shell lungo x e lungo la circonferenza           |
+| `FASCETTA`       | diametro medio, lunghezza, spessore, densità della mesh, materiale        |
-| `FEM`             | passo temporale, durata, seed, campionamento CSV, numero di analisi |
+| `ARIA`           | temperatura ambiente e coefficienti di convezione di facce e bordi        |
-| `ARIA`            | temperatura ambiente e coefficienti di convezione (facce e bordi)   |
+| `SORGENTE`       | corsa, velocità (il segno dà il verso), gruppo di sorgenti, gaussiana, flusso, frequenza |
-| `SORGENTE`        | corsa, velocità (il segno dà il verso), gruppo di sorgenti, gaussiana, flusso, frequenza |
+| `SENSORE`        | posizione (x e distanza dalla parete), inerzia, rumore, quantizzazione    |
-| `SENSORE`         | posizione (x e distanza dalla parete), inerzia, rumore, quantizzazione |
+| `RANDOMIZZAZIONE`| entità delle perturbazioni per run                                        |
 | `RANDOMIZZAZIONE` | entità delle perturbazioni per analisi                              |
 | `USCITA`          | cartella per immagini e animazioni salvate                          |
 Punti da conoscere:
@@ -214,14 +280,58 @@ Punti da conoscere:
  alla fine della propria corsa.
 - **Materiale**: `FASCETTA["materiale"]` deve essere una chiave di `MATERIALI` in
  `materials.py`. Per aggiungere un materiale basta una nuova voce nel dizionario
-  (conducibilità termica, densità, calore specifico, resistività elettrica,
+  (conducibilità, densità, calore specifico, resistività elettrica, permeabilità).
  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
+  proprietà termiche, elettriche e magnetiche in massa. I valori in `materials.py`
-  quindi quelli dell'acciaio dolce sottostante:
+  sono quindi quelli dell'acciaio dolce sottostante:
  [Banda stagnata: tutto quello che c'è da sapere (MUNDOLATAS)](https://mundolatas.com/it/banda-stagnata-tutto-quello-che-ce-da-sapere/),
  [Bande stagnate elettrolitiche (EUROPERF)](https://www.europerf.it/it/banda-stagnata-elettrolitica.php),
  [Differenza tra banda stagnata e acciaio inossidabile (Wuxi Bright Packing)](https://it.brightmetalcan.com/info/difference-between-tinplate-and-stainless-stee-48700260.html).
 - **Randomizzazione**: ogni run perturba velocità, flusso di picco, sigma, offset y,
  temperatura ambiente e rumore del sensore con estrazioni da un RNG a seed fisso
  (`SIMULAZIONE["seed"]`): il dataset è riproducibile.
 ## Output
 ### `dataset/run_XXXX.csv` — serie temporale del run
 | Colonna                        | Significato                                                        |
 |--------------------------------|--------------------------------------------------------------------|
 | `id_run`                       | identificativo del run                                             |
 | `tempo_s`                      | tempo simulato                                                     |
 | `x_sorgente_m`                 | posizione della sorgente di riferimento del gruppo                 |
 | `offset_y_sorgente_m`          | offset circonferenziale del percorso (costante nel run)            |
 | `flusso_termico_sorgente_W_m2` | flusso efficace nel punto x del sensore                            |
 | `skin_depth_m`                 | skin depth usata (costante nel run)                                |
 | `T_vera_lato_sensore_C`        | temperatura vera della superficie interna nel punto del sensore    |
 | `T_misurata_sensore_C`         | lettura del sensore (inerzia + rumore + quantizzazione)            |
 | `T_max_fascetta_C`             | temperatura massima su tutta la superficie in quell'istante        |
 | `T_media_fascetta_C`           | temperatura media su tutta la superficie in quell'istante          |
 | `T_ambiente_C`                 | temperatura ambiente del run                                       |
 | `velocita_m_s`, `sigma_punto_m`, `flusso_picco_W_m2` | parametri randomizzati del run              |
 | `materiale`                    | chiave del materiale                                               |
 ### `dataset/metadata.csv` — una riga per run
 Contiene tutti i parametri effettivi del run (geometria, dimensioni della mesh,
 coefficienti di scambio, parametri delle sorgenti e del sensore, valori randomizzati),
 le diagnostiche `skin_depth_m` e `numero_fourier_spessore`, e le temperature di picco:
 utile come ground truth e per filtrare i run.
 ## Limitazioni
 1. Il campo elettromagnetico non è simulato: l'accoppiamento induttivo è ridotto a
   un'impronta gaussiana moltiplicata per un'efficienza costante.
 2. Lo spessore non è discretizzato: il modello non può rappresentare un gradiente tra
   faccia esterna e faccia interna. Con i parametri di default il gradiente reale è di
   pochi millesimi di grado (vedi `numero_fourier_spessore`), ma con pareti più spesse
   di ~0.4 mm, materiali a bassa diffusività o sorgenti molto più veloci servirebbe una
   shell multistrato.
 3. Le proprietà dei materiali sono costanti con la temperatura; per gli acciai
   ferromagnetici la skin depth reale varia fortemente con temperatura e campo
   (punto di Curie non modellato).
 4. Irraggiamento non modellato: a ~210 °C le perdite radiative non sono del tutto
   trascurabili rispetto alla convezione.
@@ -0,0 +1,151 @@
 # Assemblaggio delle matrici globali della shell termica.
 #
 # Tutti gli elementi della mesh cilindrica sono lo stesso rettangolo, quindi
 # ogni matrice elementare si calcola una volta e si replica su tutti gli
 # elementi. L'assemblaggio è la somma dei contributi sui gradi di libertà
 # condivisi: due elementi confinanti condividono i due nodi del bordo comune,
 # ed è questo che genera la conduzione tra elementi adiacenti. Non serve
 # nessun accoppiamento aggiuntivo tra elementi dello stesso materiale.
 import numpy as np
 import scipy.sparse as sp
 import elementi_shell as el
 def _assembla_matrice(matrice_elemento: np.ndarray, connettivita: np.ndarray, n_nodi: int):
    # Replica una matrice elementare costante su tutti gli elementi e somma i
    # contributi sui gradi di libertà condivisi.
    n_elementi, n_locali = connettivita.shape
    righe = np.repeat(connettivita, n_locali, axis=1).ravel()
    colonne = np.tile(connettivita, (1, n_locali)).ravel()
    dati = np.tile(matrice_elemento.ravel(), n_elementi)
    return sp.coo_matrix(
        (dati, (righe, colonne)), shape=(n_nodi, n_nodi)
    ).tocsr()
 def _assembla_vettore(vettore_elemento: np.ndarray, connettivita: np.ndarray, n_nodi: int):
    dati = np.tile(vettore_elemento, connettivita.shape[0])
    return np.bincount(connettivita.ravel(), weights=dati, minlength=n_nodi)
 def assembla_capacita(mesh: dict, materiale: dict, spessore_m: float):
    # C = integrale di rho * cp * spessore * N^T N dA.
    #
    # Matrice consistente: con integrazione temporale implicita conserva
    # meglio l'energia rispetto alla versione concentrata sulla diagonale.
    capacita_superficiale = (
        materiale["densita_kg_m3"] * materiale["calore_specifico_J_kgK"] * spessore_m
    )
    M_e = el.matrice_massa_elemento(mesh["lato_x_m"], mesh["lato_arco_m"])
    return _assembla_matrice(
        capacita_superficiale * M_e, mesh["connettivita"], mesh["n_nodi"]
    )
 def assembla_conduzione(mesh: dict, materiale: dict, spessore_m: float):
    # K_cond = integrale di B^T D B * spessore dA, con D isotropa.
    K_e = el.matrice_conduzione_elemento(
        mesh["lato_x_m"],
        mesh["lato_arco_m"],
        materiale["conducibilita_termica_W_mK"],
        spessore_m,
    )
    return _assembla_matrice(K_e, mesh["connettivita"], mesh["n_nodi"])
 def assembla_convezione(mesh: dict, aria: dict, spessore_m: float):
    # Convezione sulle facce e sui bordi assiali.
    #
    # Nelle shell la faccia esterna e quella interna sono entrambe superfici
    # fisiche esposte all'aria, quindi ogni elemento scambia su tutta la sua
    # area con il coefficiente combinato h_esterno + h_interno. I bordi
    # assiali (x = 0 e x = lunghezza) espongono invece solo lo spessore della
    # lamiera: il loro contributo è proporzionale a h_bordi * spessore ed è
    # marginale rispetto a quello delle facce.
    T_ambiente = aria["temperatura_ambiente_C"]
    h_facce = aria["h_esterno_W_m2K"] + aria["h_interno_W_m2K"]
    M_e = el.matrice_massa_elemento(mesh["lato_x_m"], mesh["lato_arco_m"])
    f_e = el.vettore_carico_uniforme_elemento(mesh["lato_x_m"], mesh["lato_arco_m"])
    K_conv = _assembla_matrice(h_facce * M_e, mesh["connettivita"], mesh["n_nodi"])
    f_ambiente = _assembla_vettore(
        h_facce * T_ambiente * f_e, mesh["connettivita"], mesh["n_nodi"]
    )
    h_bordi_efficace = aria["h_bordi_W_m2K"] * spessore_m
    M_bordo = el.matrice_massa_segmento(mesh["lato_arco_m"])
    f_bordo = el.vettore_carico_segmento(mesh["lato_arco_m"])
    K_conv = K_conv + _assembla_matrice(
        h_bordi_efficace * M_bordo, mesh["segmenti_bordo"], mesh["n_nodi"]
    )
    f_ambiente = f_ambiente + _assembla_vettore(
        h_bordi_efficace * T_ambiente * f_bordo,
        mesh["segmenti_bordo"],
        mesh["n_nodi"],
    )
    return K_conv, f_ambiente
 def prepara_assemblatore_sorgente(mesh: dict) -> dict:
    # Precalcola quanto serve per assemblare il vettore della sorgente a ogni
    # passo temporale: le funzioni di forma nei punti di Gauss pesate per il
    # peso di quadratura e lo jacobiano, e le coordinate dei punti di Gauss
    # separate per direzione.
    det_j = mesh["lato_x_m"] * mesh["lato_arco_m"] / 4.0
    pesi = el.pesi_ai_punti_gauss()
    N_gauss = el.matrice_forme_ai_punti_gauss()
    offset_x, offset_arco = el.coordinate_locali_punti_gauss(
        mesh["lato_x_m"], mesh["lato_arco_m"]
    )
    # Coordinate assolute dei punti di Gauss, per direzione: (n_elementi_x, 2)
    # e (n_elementi_theta, 2).
    x_gauss = mesh["x_nodi_m"][: mesh["n_elementi_x"], None] + offset_x[None, :]
    arco_gauss = mesh["arco_nodi_m"][:, None] + offset_arco[None, :]
    return {
        "peso_forme": (pesi * det_j)[:, None] * N_gauss,
        "connettivita": mesh["connettivita"],
        "n_nodi": mesh["n_nodi"],
        "n_elementi_x": mesh["n_elementi_x"],
        "n_elementi_theta": mesh["n_elementi_theta"],
        "x_gauss_m": x_gauss,
        "arco_gauss_m": arco_gauss,
    }
 def assembla_sorgente(
    assemblatore: dict, flusso_x_W_m2: np.ndarray, fattore_arco: np.ndarray
 ) -> np.ndarray:
    # f_src = integrale di N^T q'' dA, con q'' fattorizzato come
    # q''(x, arco) = flusso_x(x) * fattore_arco(arco).
    #
    # La fattorizzazione non è un'approssimazione: la gaussiana della
    # sorgente è separabile in x e in arco, e tutte le sorgenti del gruppo
    # condividono la stessa posizione circonferenziale, quindi il fattore
    # circonferenziale è comune e costante nel tempo.
    #
    # flusso_x_W_m2 ha forma (n_elementi_x, 2) e fattore_arco
    # (n_elementi_theta, 2): sono i valori nei punti di Gauss di ciascuna
    # direzione.
    #
    # Per gran parte di un run le sorgenti sono fuori dalla fascetta e il
    # fattore assiale è identicamente nullo: in quel caso il vettore è nullo
    # e non serve percorrere la mesh.
    if not flusso_x_W_m2.any():
        return np.zeros(assemblatore["n_nodi"])
    q_gauss = (
        flusso_x_W_m2[:, None, :, None] * fattore_arco[None, :, None, :]
    ).reshape(-1, 4)
    contributi = q_gauss @ assemblatore["peso_forme"]
    return np.bincount(
        assemblatore["connettivita"].ravel(),
        weights=contributi.ravel(),
        minlength=assemblatore["n_nodi"],
    )
@@ -1,21 +1,29 @@
-# Configurazione dell'analisi termica a elementi finiti della fascetta.
+# Configurazione per il simulatore termico a elementi shell di una fascetta.
 #
 # Geometria e modello fisico:
-# - La fascetta è un anello cilindrico con diametro "diametro_mm", spessore
+# - La fascetta è un anello cilindrico con diametro medio "diametro_medio_mm",
-#   "spessore_mm" e lunghezza "lunghezza_mm".
+#   spessore "spessore_mm" e lunghezza "lunghezza_mm".
-# - Il dominio discretizzato è la superficie media del cilindro (mesh a
+# - Il dominio simulato è l'intera superficie cilindrica, discretizzata con
-#   elementi shell quadrangolari, vedi mesh.py): x = asse della fascetta,
+#   elementi shell quadrilateri a 4 nodi e una temperatura per nodo.
-#   s = R·theta = coordinata circonferenziale. Lo spessore non è
+# - Sistema di coordinate: x = direzione della lunghezza (da 0 a lunghezza),
-#   discretizzato, la temperatura è uniforme attraverso la parete.
+#   theta = coordinata circonferenziale, periodica (theta = 0 coincide con
-# - Le sorgenti a induzione si muovono lungo x sulla superficie esterna, con
+#   theta = 2*pi). Il sensore sta per convenzione su theta = 0.
-#   un'impronta gaussiana isotropa nel piano (x, s).
+# - Lo spessore non è discretizzato: la parete è termicamente sottile, quindi
-# - La superficie scambia per convezione con l'aria sulla faccia esterna e su
+#   la faccia esterna e quella interna hanno la stessa temperatura. La
-#   quella interna, più i due bordi anulari x = 0 e x = lunghezza.
+#   validità dell'ipotesi è riportata a ogni run come numero di Fourier sullo
 #   spessore in metadata.csv.
 # - Le sorgenti a induzione si muovono in direzione -x sulla superficie
 #   esterna, lungo la generatrice individuata da "offset_y_percorso_m".
 # - La conduzione assiale e quella circonferenziale sono entrambe risolte e
 #   nascono dall'assemblaggio: elementi adiacenti condividono i nodi del
 #   bordo comune.
 # - Faccia esterna e faccia interna scambiano entrambe per convezione con
 #   l'aria su tutta l'area; i due bordi assiali scambiano attraverso lo
 #   spessore della lamiera.
 # - La temperatura iniziale della fascetta è pari alla temperatura ambiente.
 # - Il sensore è un pirometro a infrarossi posto all'interno della fascetta,
 #   a distanza "distanza_parete_mm" dalla parete interna: misura senza
-#   contatto la temperatura della superficie interna nel punto x = "x_mm",
+#   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)
@@ -24,9 +32,36 @@
 # - flusso termico: W/m²
 # - coefficiente di convezione: W/(m² K)
 SIMULAZIONE = {
    # Numero di file CSV da generare.
    "num_run": 1,
    # Tempo simulato totale.
    "durata_s": 30.0,
    # Frequenza di campionamento CSV.
    # Esempio: 2 Hz significa una riga ogni 0.5 s.
    "frequenza_campionamento_hz": 10.0,
    # Passo di integrazione numerica interna.
    # Può essere inferiore al periodo di campionamento CSV.
    "dt_interno_s": 0.001,
    # Seed per la riproducibilità.
    "seed": 42,
    # Cartella di output.
    "cartella_output": "dataset",
    # Numero di processi paralleli per la generazione dei run.
    # None = usa tutti i core disponibili.
    "num_processi": None,
 }
 FASCETTA = {
-    # Diametro della fascetta [mm].
+    # Diametro medio della fascetta [mm], cioè il diametro della superficie
-    "diametro_mm": 70.0,
+    # su cui giace la mesh shell.
    "diametro_medio_mm": 70.0,
    # Lunghezza della fascetta lungo x [mm].
    "lunghezza_mm": 100.0,
@@ -34,36 +69,15 @@ FASCETTA = {
    # Spessore della parete [mm].
    "spessore_mm": 0.18,
-    # Deve corrispondere a una chiave in materials.py.
+    # Numero di elementi shell lungo x (lunghezza).
-    "materiale": "banda_stagnata",
+    "n_elementi_x": 100,
 }
 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.
+    # Più elementi = maggiore risoluzione spaziale, simulazione più lenta.
-    "n_elementi_circonferenza": 48,
+    "n_elementi_theta": 110,
 }
-FEM = {
+    # Deve corrispondere a una chiave in materials.py.
-    # Passo di integrazione temporale del solutore a elementi finiti.
+    "materiale": "banda_stagnata",
    "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 = {
@@ -76,7 +90,8 @@ ARIA = {
    # Coefficiente di convezione sulla faccia interna (lato sensore).
    "h_interno_W_m2K": 8.0,
-    # Coefficiente di convezione sui bordi anulari (x = 0 e x = lunghezza).
+    # Coefficiente di convezione sui bordi assiali (x = 0 e x = lunghezza),
    # che espongono all'aria solo lo spessore della lamiera.
    "h_bordi_W_m2K": 10.0,
 }
@@ -89,9 +104,10 @@ SORGENTE = {
    # gruppo più lontana dal sensore (quella che lo supera per ultima).
    "x_fine_m": 5.0,
-    # Offset circonferenziale (y) tra il percorso delle sorgenti e il punto
+    # Offset circonferenziale, misurato lungo l'arco, tra il percorso delle
-    # osservato dal sensore. Se diverso da zero, la sorgente passa a lato,
+    # sorgenti e la generatrice osservata dal sensore. Se diverso da zero, la
-    # riducendo il picco.
+    # sorgente passa a lato e il calore deve diffondere lungo la
    # circonferenza per raggiungere il sensore.
    "offset_y_percorso_m": 0.0,
    # Velocità delle sorgenti lungo x. Il segno indica il verso di marcia.
@@ -113,8 +129,10 @@ SORGENTE = {
    # Frazione del flusso incidente che diventa effettivamente calore nella fascetta.
    "efficienza_riscaldamento": 0.35,
-    # Frequenza di induzione dell'induttore, usata per stimare la skin depth
+    # Frequenza di induzione usata per stimare la skin depth se skin_depth_fissa_m è None.
-    # se skin_depth_fissa_m è None.
+    # La shell non discretizza lo spessore, quindi la skin depth non entra
    # nella soluzione: è calcolata e registrata come diagnostica, per
    # verificare che resti confrontabile con lo spessore della lamiera.
    "frequenza_hz": 20000.0,
    # Override della skin depth. Usare None per calcolarla dalle proprietà elettriche del materiale.
@@ -127,9 +145,11 @@ SORGENTE = {
 SENSORE = {
    # Coordinata x del punto della superficie interna osservato dal sensore [mm].
    # La posizione circonferenziale è theta = 0 per convenzione: è l'offset
    # delle sorgenti a essere misurato rispetto al sensore.
    "x_mm": 50.0,
-    # Distanza del sensore dalla parete interna lungo lo spessore [mm].
+    # Distanza del sensore dalla parete interna [mm].
    # Il sensore è a infrarossi: la distanza non influenza la misura,
    # è registrata solo come geometria del setup.
    "distanza_parete_mm": 10.0,
@@ -148,7 +168,7 @@ SENSORE = {
 }
 RANDOMIZZAZIONE = {
-    # Se abilitata, ogni analisi varia leggermente alcuni parametri.
+    # Se abilitata, ogni run varia leggermente alcuni parametri.
    "abilitata": False,
    # Deviazioni standard relative.
@@ -165,9 +185,3 @@ RANDOMIZZAZIONE = {
    # dalla linea ideale allineata con il sensore.
    "offset_y_max_assoluto_m": 0.001,
 }
 USCITA = {
    # Cartella in cui salvare immagini e animazioni quando il backend
    # matplotlib non è interattivo.
    "cartella": "output",
 }
@@ -0,0 +1,135 @@
 # Elemento shell termico quadrilatero a 4 nodi (Q4) con una temperatura per
 # nodo, cioè temperatura uniforme nello spessore.
 #
 # L'elemento vive nel piano (x, arco) della superficie sviluppata. Poiché la
 # mesh cilindrica è strutturata e uniforme, tutti gli elementi sono lo stesso
 # rettangolo lato_x × lato_arco: le matrici elementari si calcolano una volta
 # sola e si replicano in fase di assemblaggio.
 #
 # Le funzioni di forma bilineari sono
 #
 #     N_a(xi, eta) = (1 + xi_a * xi) * (1 + eta_a * eta) / 4
 #
 # con (xi_a, eta_a) i vertici dell'elemento di riferimento. La mappa verso
 # l'elemento fisico è affine, quindi lo jacobiano è costante:
 # det(J) = lato_x * lato_arco / 4.
 #
 # Integrazione con quadratura di Gauss 2x2, esatta per i polinomi che
 # compaiono nelle matrici di massa e conduzione dell'elemento Q4.
 import numpy as np
 # Vertici dell'elemento di riferimento, nello stesso ordine dei nodi locali
 # usati dalla connettività in mesh.py.
 NODI_RIFERIMENTO = np.array([(-1.0, -1.0), (1.0, -1.0), (1.0, 1.0), (-1.0, 1.0)])
 # Quadratura di Gauss a 2 punti per direzione.
 PUNTI_GAUSS_1D = np.array([-1.0 / np.sqrt(3.0), 1.0 / np.sqrt(3.0)])
 PESI_GAUSS_1D = np.array([1.0, 1.0])
 def funzioni_forma(xi: float, eta: float) -> np.ndarray:
    xi_a, eta_a = NODI_RIFERIMENTO[:, 0], NODI_RIFERIMENTO[:, 1]
    return 0.25 * (1.0 + xi_a * xi) * (1.0 + eta_a * eta)
 def derivate_funzioni_forma(xi: float, eta: float) -> tuple[np.ndarray, np.ndarray]:
    xi_a, eta_a = NODI_RIFERIMENTO[:, 0], NODI_RIFERIMENTO[:, 1]
    dN_dxi = 0.25 * xi_a * (1.0 + eta_a * eta)
    dN_deta = 0.25 * eta_a * (1.0 + xi_a * xi)
    return dN_dxi, dN_deta
 def _punti_quadratura_2d():
    # Punti (xi, eta) e pesi della quadratura 2x2, nell'ordine
    # g = p * 2 + q con xi = punto p (direzione x) ed eta = punto q
    # (direzione circonferenziale). Lo stesso ordine è atteso dai valori del
    # flusso passati a vettore_sorgente_elementi().
    for p, xi in enumerate(PUNTI_GAUSS_1D):
        for q, eta in enumerate(PUNTI_GAUSS_1D):
            yield xi, eta, PESI_GAUSS_1D[p] * PESI_GAUSS_1D[q]
 def matrice_forme_ai_punti_gauss() -> np.ndarray:
    # Matrice (4 punti di Gauss, 4 nodi) con le funzioni di forma valutate
    # nei punti di quadratura.
    return np.array([funzioni_forma(xi, eta) for xi, eta, _ in _punti_quadratura_2d()])
 def pesi_ai_punti_gauss() -> np.ndarray:
    # Pesi di quadratura nei 4 punti, nello stesso ordine.
    return np.array([peso for _, _, peso in _punti_quadratura_2d()])
 def coordinate_locali_punti_gauss(lato_x_m: float, lato_arco_m: float):
    # Offset dei punti di Gauss rispetto al vertice di riferimento
    # dell'elemento, separati per direzione: 2 valori lungo x e 2 lungo
    # l'arco. La separazione è possibile perché l'elemento è un rettangolo.
    offset_x = lato_x_m * (1.0 + PUNTI_GAUSS_1D) / 2.0
    offset_arco = lato_arco_m * (1.0 + PUNTI_GAUSS_1D) / 2.0
    return offset_x, offset_arco
 def matrice_massa_elemento(lato_x_m: float, lato_arco_m: float) -> np.ndarray:
    # Integrale di N^T N sull'area dell'elemento. Matrice consistente (non
    # concentrata): è la base sia della capacità termica sia della convezione
    # sulle facce.
    det_j = lato_x_m * lato_arco_m / 4.0
    M = np.zeros((4, 4))
    for xi, eta, peso in _punti_quadratura_2d():
        N = funzioni_forma(xi, eta)
        M += peso * det_j * np.outer(N, N)
    return M
 def matrice_conduzione_elemento(
    lato_x_m: float, lato_arco_m: float, conducibilita_W_mK: float, spessore_m: float
 ) -> np.ndarray:
    # Integrale di B^T D B * spessore sull'area, con D = k * identità
    # (materiale isotropo). Rappresenta la conduzione tangenziale, cioè sia
    # quella assiale sia quella circonferenziale: sul cilindro sviluppato le
    # due direzioni sono ortogonali e la derivata circonferenziale rispetto
    # all'arco equivale a (1/R) d/dtheta.
    det_j = lato_x_m * lato_arco_m / 4.0
    K = np.zeros((4, 4))
    for xi, eta, peso in _punti_quadratura_2d():
        dN_dxi, dN_deta = derivate_funzioni_forma(xi, eta)
        # Jacobiano diagonale e costante: la derivata fisica è quella
        # naturale riscalata dal semilato dell'elemento.
        B = np.stack([dN_dxi * (2.0 / lato_x_m), dN_deta * (2.0 / lato_arco_m)])
        K += peso * det_j * conducibilita_W_mK * spessore_m * (B.T @ B)
    return K
 def vettore_carico_uniforme_elemento(lato_x_m: float, lato_arco_m: float) -> np.ndarray:
    # Integrale di N^T sull'area: distribuisce ai nodi un carico superficiale
    # uniforme unitario.
    det_j = lato_x_m * lato_arco_m / 4.0
    f = np.zeros(4)
    for xi, eta, peso in _punti_quadratura_2d():
        f += peso * det_j * funzioni_forma(xi, eta)
    return f
 def _funzioni_forma_segmento(xi: float) -> np.ndarray:
    return np.array([0.5 * (1.0 - xi), 0.5 * (1.0 + xi)])
 def matrice_massa_segmento(lunghezza_m: float) -> np.ndarray:
    # Integrale di N^T N su un segmento a 2 nodi, per i bordi assiali.
    det_j = lunghezza_m / 2.0
    M = np.zeros((2, 2))
    for xi, peso in zip(PUNTI_GAUSS_1D, PESI_GAUSS_1D):
        N = _funzioni_forma_segmento(xi)
        M += peso * det_j * np.outer(N, N)
    return M
 def vettore_carico_segmento(lunghezza_m: float) -> np.ndarray:
    # Integrale di N^T su un segmento a 2 nodi.
    det_j = lunghezza_m / 2.0
    f = np.zeros(2)
    for xi, peso in zip(PUNTI_GAUSS_1D, PESI_GAUSS_1D):
        f += peso * det_j * _funzioni_forma_segmento(xi)
    return f
@@ -1,153 +0,0 @@
 # Export CSV della serie temporale calcolata dal solutore FEM (fem.py).
 #
 # Per ogni analisi viene scritto un file run_XXXX.csv con la temperatura vera
 # della parete nel punto osservato dal sensore e la lettura del sensore reale
 # (inerzia, rumore, quantizzazione), più i parametri effettivi dell'analisi.
 # `metadata.csv` riassume un'analisi per riga.
 #
 # Le righe sono campionate a FEM["frequenza_campionamento_hz"], indipendente
 # dal passo di integrazione FEM["dt_s"]. I campi nodali non vengono accumulati:
 # servono solo le serie scalari.
 import csv
 import random
 import shutil
 from pathlib import Path
 from config import FEM, USCITA
 from fem import configurazione_randomizzata, simula_campo_fem
 from mesh import genera_mesh, riepilogo_mesh
 INTESTAZIONE = [
    "id_run",
    "tempo_s",
    "x_sorgente_m",
    "T_vera_sensore_C",
    "T_misurata_sensore_C",
    "T_massima_fascetta_C",
    "T_ambiente_C",
    "offset_y_sorgente_m",
    "velocita_m_s",
    "sigma_punto_m",
    "flusso_picco_W_m2",
    "skin_depth_m",
    "materiale",
 ]
 def esporta_run(
    id_run: str,
    percorso_csv: Path,
    mesh_dati: dict,
    rng: random.Random,
 ) -> dict:
    """Integra un'analisi FEM, ne scrive il CSV e restituisce i suoi metadati."""
    cfg = configurazione_randomizzata(rng)
    dati = simula_campo_fem(
        dt_frame_s=1.0 / FEM["frequenza_campionamento_hz"],
        cfg=cfg,
        mesh_dati=mesh_dati,
        rng=rng,
        salva_campi=False,
    )
    fascetta = cfg["fascetta"]
    aria = cfg["aria"]
    sorgente = cfg["sorgente"]
    sensore = cfg["sensore"]
    skin_depth = dati["skin_depth_m"]
    percorso_csv.parent.mkdir(parents=True, exist_ok=True)
    with percorso_csv.open("w", newline="") as f:
        writer = csv.writer(f)
        writer.writerow(INTESTAZIONE)
        for k, t in enumerate(dati["tempi"]):
            writer.writerow([
                id_run,
                f"{t:.6f}",
                f"{dati['x_riferimenti'][k]:.9f}",
                f"{dati['T_vere'][k]:.6f}",
                f"{dati['T_lette'][k]:.6f}",
                f"{dati['T_massime'][k]:.6f}",
                f"{aria['temperatura_ambiente_C']:.6f}",
                f"{sorgente['offset_y_percorso_m']:.9f}",
                f"{sorgente['velocita_m_s']:.9f}",
                f"{sorgente['sigma_punto_m']:.9f}",
                f"{sorgente['flusso_termico_picco_W_m2']:.6f}",
                f"{skin_depth:.9e}",
                fascetta["materiale"],
            ])
    return {
        "id_run": id_run,
        "file_csv": percorso_csv.name,
        "materiale": fascetta["materiale"],
        "diametro_m": fascetta["diametro_mm"] / 1000.0,
        "lunghezza_m": fascetta["lunghezza_mm"] / 1000.0,
        "spessore_m": fascetta["spessore_mm"] / 1000.0,
        "n_elementi_x": mesh_dati["n_elementi_x"],
        "n_elementi_circonferenza": mesh_dati["n_elementi_circonferenza"],
        "durata_s": FEM["durata_s"],
        "dt_s": FEM["dt_s"],
        "frequenza_campionamento_hz": FEM["frequenza_campionamento_hz"],
        "temperatura_ambiente_C": aria["temperatura_ambiente_C"],
        "h_esterno_W_m2K": aria["h_esterno_W_m2K"],
        "h_interno_W_m2K": aria["h_interno_W_m2K"],
        "h_bordi_W_m2K": aria["h_bordi_W_m2K"],
        "x_inizio_m": sorgente["x_inizio_m"],
        "x_fine_m": sorgente["x_fine_m"],
        "x_sensore_m": dati["stato"]["x_sensore_m"],
        "distanza_sensore_parete_m": sensore["distanza_parete_mm"] / 1000.0,
        "offset_y_percorso_m": sorgente["offset_y_percorso_m"],
        "velocita_m_s": sorgente["velocita_m_s"],
        "numero_sorgenti": sorgente.get("numero_sorgenti", 1),
        "distanza_sorgenti_m": sorgente.get("distanza_sorgenti_m", 0.0),
        "sigma_punto_m": sorgente["sigma_punto_m"],
        "flusso_termico_picco_W_m2": sorgente["flusso_termico_picco_W_m2"],
        "efficienza_riscaldamento": sorgente["efficienza_riscaldamento"],
        "frequenza_hz": sorgente["frequenza_hz"],
        "skin_depth_m": skin_depth,
        "costante_tempo_sensore_s": sensore["costante_tempo_s"],
        "rumore_std_sensore_C": sensore["rumore_std_C"],
        "quantizzazione_sensore_C": sensore["quantizzazione_C"],
        "T_vera_max_sensore_C": float(dati["T_vere"].max()),
        "T_misurata_max_sensore_C": float(dati["T_lette"].max()),
        "T_massima_fascetta_C": float(dati["T_massime"].max()),
    }
 def main() -> None:
    mesh_dati = genera_mesh()
    print(riepilogo_mesh(mesh_dati))
    # La cartella dei CSV è ricreata da zero a ogni esecuzione, così non
    # restano run di esecuzioni precedenti con parametri diversi.
    cartella = Path(USCITA["cartella"]) / "csv"
    if cartella.exists():
        shutil.rmtree(cartella)
    cartella.mkdir(parents=True, exist_ok=True)
    rng = random.Random(FEM["seed"])
    righe_metadata = []
    for indice in range(1, FEM["num_run"] + 1):
        id_run = f"run_{indice:04d}"
        print(f"Integrazione FEM di {id_run} in corso...")
        riga = esporta_run(id_run, cartella / f"{id_run}.csv", mesh_dati, rng)
        print(
            f"  T vera max al sensore: {riga['T_vera_max_sensore_C']:.1f} °C — "
            f"T misurata max: {riga['T_misurata_max_sensore_C']:.1f} °C"
        )
        righe_metadata.append(riga)
    percorso_metadata = cartella / "metadata.csv"
    with percorso_metadata.open("w", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=list(righe_metadata[0].keys()))
        writer.writeheader()
        writer.writerows(righe_metadata)
    print(f"Esportate {len(righe_metadata)} analisi in: {cartella.resolve()}")
    print(f"Metadata: {percorso_metadata.resolve()}")
 if __name__ == "__main__":
    main()
@@ -1,460 +0,0 @@
 # Solutore termico transitorio a elementi finiti sulla mesh shell della
 # fascetta (vedi mesh.py per la geometria).
 #
 # Il dominio è la superficie media del cilindro, con x e la coordinata
 # circonferenziale s = R·theta entrambe risolte, mentre lo spessore è
 # collassato (temperatura uniforme attraverso la parete).
 #
 # Il collasso in spessore è lecito: con spessore 0.18 mm il numero di Biot
 # h·t/k vale ~1e-7 e il tempo di diffusione attraverso la parete t²/alpha è
 # di pochi millisecondi, molto più rapido del transito delle sorgenti. Di
 # conseguenza la profondità di penetrazione del riscaldamento a induzione non
 # entra nel modello: conta solo il flusso totale assorbito per unità di
 # superficie. La skin depth è comunque calcolata dalla frequenza dell'induttore
 # come grandezza diagnostica (`calcola_skin_depth_m`).
 #
 # Equazione risolta (per unità di superficie media):
 #
 #   rho·cp·t ∂T/∂t = k·t ∇²T + q(x, s, t) - (h_est + h_int)·(T - T_amb)
 #
 # con convezione aggiuntiva sui due bordi anulari x = 0 e x = L (coefficiente
 # h_bordi su un'area pari a spessore × perimetro). L'impronta della sorgente
 # è una gaussiana isotropa nel piano (x, s), con la distanza circonferenziale
 # valutata sull'immagine più vicina perché la superficie è chiusa.
 #
 # Discretizzazione: elementi shell quadrangolari bilineari a 4 nodi. Gli
 # elementi della mesh sono tutti rettangoli identici dx × ds, quindi le
 # matrici di elemento sono calcolate una volta sola e assemblate in forma
 # vettorizzata. Integrazione temporale con Eulero implicito: la matrice di
 # sistema è costante e viene fattorizzata LU una volta per run.
 import math
 import random
 from copy import deepcopy
 import numpy as np
 import scipy.sparse as sp
 from scipy.sparse.linalg import splu
 from config import ARIA, FASCETTA, FEM, RANDOMIZZAZIONE, SENSORE, SORGENTE
 from materials import MATERIALI
 from mesh import genera_mesh
 MU0 = 4.0 * math.pi * 1e-7
 # Matrici di riferimento dell'elemento rettangolare bilineare a 4 nodi, con
 # nodi locali in ordine antiorario: (0,0), (a,0), (a,b), (0,b).
 #
 # MASSA_RIF va moltiplicata per l'area a·b e dà l'integrale di N_i·N_j;
 # RIGIDEZZA_X per b/a e RIGIDEZZA_S per a/b, e la loro somma dà l'integrale
 # di grad(N_i)·grad(N_j).
 MASSA_RIF = np.array(
    [[4.0, 2.0, 1.0, 2.0],
     [2.0, 4.0, 2.0, 1.0],
     [1.0, 2.0, 4.0, 2.0],
     [2.0, 1.0, 2.0, 4.0]]
 ) / 36.0
 RIGIDEZZA_X = np.array(
    [[2.0, -2.0, -1.0, 1.0],
     [-2.0, 2.0, 1.0, -1.0],
     [-1.0, 1.0, 2.0, -2.0],
     [1.0, -1.0, -2.0, 2.0]]
 ) / 6.0
 RIGIDEZZA_S = np.array(
    [[2.0, 1.0, -1.0, -2.0],
     [1.0, 2.0, -2.0, -1.0],
     [-1.0, -2.0, 2.0, 1.0],
     [-2.0, -1.0, 1.0, 2.0]]
 ) / 6.0
 def calcola_skin_depth_m(materiale: dict, frequenza_hz: float) -> float:
    # Skin depth elettromagnetica approssimata:
    # delta = sqrt(2 * rho_e / (omega * mu))
    #
    # Semplificata. Per acciai ferromagnetici il comportamento reale
    # è fortemente non lineare con temperatura e campo magnetico.
    #
    # Nel modello shell è una grandezza diagnostica: serve a verificare che la
    # deposizione del calore resti confinata entro uno spessore confrontabile
    # con la parete, non entra nell'equazione risolta.
    rho_e = materiale["resistivita_elettrica_ohm_m"]
    mu_r = materiale["permeabilita_relativa"]
    omega = 2.0 * math.pi * frequenza_hz
    mu = MU0 * mu_r
    return math.sqrt(2.0 * rho_e / (omega * mu))
 def quantizza(valore: float, passo: float) -> float:
    if passo <= 0.0:
        return valore
    return round(valore / passo) * passo
 def configurazione_randomizzata(rng: random.Random | None = None) -> dict:
    """Copia dei dizionari di configurazione usati da un'analisi FEM.
    Se `RANDOMIZZAZIONE["abilitata"]`, velocità, flusso di picco, sigma dello
    spot, offset y del percorso, temperatura ambiente e rumore del sensore sono
    perturbati con estrazioni dall'RNG passato: analisi diverse dello stesso
    setup differiscono come differirebbero due passaggi reali.
    """
    fascetta = deepcopy(FASCETTA)
    aria = deepcopy(ARIA)
    sorgente = deepcopy(SORGENTE)
    sensore = deepcopy(SENSORE)
    if RANDOMIZZAZIONE.get("abilitata", False):
        if rng is None:
            rng = random.Random(FEM["seed"])
        def perturba_rel(valore: float, std_rel: float, fattore_min: float = 0.1) -> float:
            fattore = rng.gauss(1.0, std_rel)
            fattore = max(fattore_min, fattore)
            return valore * fattore
        sorgente["velocita_m_s"] = perturba_rel(
            sorgente["velocita_m_s"],
            RANDOMIZZAZIONE["velocita_std_rel"],
        )
        sorgente["flusso_termico_picco_W_m2"] = perturba_rel(
            sorgente["flusso_termico_picco_W_m2"],
            RANDOMIZZAZIONE["flusso_picco_std_rel"],
        )
        sorgente["sigma_punto_m"] = perturba_rel(
            sorgente["sigma_punto_m"],
            RANDOMIZZAZIONE["sigma_punto_std_rel"],
        )
        sorgente["offset_y_percorso_m"] = rng.uniform(
            -RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
            RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
        )
        aria["temperatura_ambiente_C"] += rng.gauss(
            0.0,
            RANDOMIZZAZIONE["temperatura_ambiente_std_C"],
        )
        sensore["rumore_std_C"] = perturba_rel(
            sensore["rumore_std_C"],
            RANDOMIZZAZIONE["rumore_sensore_std_rel"],
            fattore_min=0.0,
        )
    return {
        "fascetta": fascetta,
        "aria": aria,
        "sorgente": sorgente,
        "sensore": sensore,
    }
 def _spread_sorgenti_m(sorgente: dict) -> float:
    # Distanza lungo x tra la prima e l'ultima sorgente del gruppo.
    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
    return (numero_sorgenti - 1) * distanza
 def _x_riferimento_iniziale_m(sorgente: dict, x_sensore_m: float) -> float:
    # Posizione a t=0 della sorgente di indice 0. Le sorgenti i sono a
    # x_i = riferimento + i * distanza, quindi per v >= 0 l'indice 0 è la
    # più arretrata nel verso di marcia, per v < 0 è la più avanzata.
    # x_inizio_m è la distanza dal sensore della sorgente più avanzata
    # (quella che lo raggiunge per prima).
    spread = _spread_sorgenti_m(sorgente)
    x_inizio = sorgente["x_inizio_m"]
    if sorgente["velocita_m_s"] >= 0:
        return (x_sensore_m - x_inizio) - spread
    return x_sensore_m + x_inizio
 def _x_riferimento_finale_m(sorgente: dict, x_sensore_m: float) -> float:
    # Posizione di fine corsa della sorgente di indice 0. x_fine_m è la
    # distanza dal sensore della sorgente più arretrata nel verso di marcia
    # (quella che lo supera per ultima).
    spread = _spread_sorgenti_m(sorgente)
    x_fine = sorgente["x_fine_m"]
    if sorgente["velocita_m_s"] >= 0:
        return x_sensore_m + x_fine
    return (x_sensore_m - x_fine) - spread
 def _intervallo_attivo(inizio: float, fine: float, v: float, x_m: float) -> bool:
    if v >= 0:
        return inizio <= x_m <= fine
    return fine <= x_m <= inizio
 def _assembla(elementi: np.ndarray, matrice_elemento: np.ndarray, n_nodi: int):
    # Assembla una matrice globale sparsa a partire da un'unica matrice di
    # elemento 4×4, uguale per tutti gli elementi. I contributi ripetuti sullo
    # stesso (riga, colonna) vengono sommati dal formato COO.
    n_elementi = len(elementi)
    righe = np.broadcast_to(elementi[:, :, None], (n_elementi, 4, 4)).ravel()
    colonne = np.broadcast_to(elementi[:, None, :], (n_elementi, 4, 4)).ravel()
    valori = np.broadcast_to(matrice_elemento, (n_elementi, 4, 4)).ravel()
    return sp.coo_matrix(
        (valori, (righe, colonne)), shape=(n_nodi, n_nodi)
    ).tocsr()
 def _massa_anello(indici: np.ndarray, lunghezza_segmento_m: float, n_nodi: int):
    # Matrice di massa 1D su un anello chiuso di nodi (integrale di N_i·N_j
    # lungo la circonferenza), usata per la convezione sui bordi x = 0 e x = L.
    n = len(indici)
    successivi = indici[(np.arange(n) + 1) % n]
    coefficiente = lunghezza_segmento_m / 6.0
    righe = np.concatenate([indici, indici, successivi, successivi])
    colonne = np.concatenate([indici, successivi, indici, successivi])
    valori = coefficiente * np.concatenate(
        [np.full(n, 2.0), np.full(n, 1.0), np.full(n, 1.0), np.full(n, 2.0)]
    )
    return sp.coo_matrix(
        (valori, (righe, colonne)), shape=(n_nodi, n_nodi)
    ).tocsr()
 def prepara_stato_fem(
    fascetta: dict,
    aria: dict,
    sorgente: dict,
    sensore: dict,
    mesh_dati: dict | None = None,
    dt_s: float | None = None,
 ) -> dict:
    """Assembla le matrici FEM e fattorizza il sistema implicito.
    Restituisce lo stato costante per l'integrazione temporale: matrici,
    solutore LU, coordinate nodali, peso circonferenziale della sorgente e
    indice del nodo osservato dal sensore.
    """
    if mesh_dati is None:
        mesh_dati = genera_mesh(fascetta)
    if dt_s is None:
        dt_s = FEM["dt_s"]
    materiale = MATERIALI[fascetta["materiale"]]
    if sorgente["skin_depth_fissa_m"] is None:
        skin_depth_m = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
    else:
        skin_depth_m = float(sorgente["skin_depth_fissa_m"])
    k = materiale["conducibilita_termica_W_mK"]
    rho = materiale["densita_kg_m3"]
    cp = materiale["calore_specifico_J_kgK"]
    spessore_m = mesh_dati["spessore_m"]
    raggio_m = mesh_dati["raggio_m"]
    circonferenza_m = 2.0 * math.pi * raggio_m
    n_x = mesh_dati["n_elementi_x"]
    n_theta = mesh_dati["n_elementi_circonferenza"]
    elementi = mesh_dati["elementi"]
    indice_nodo = mesh_dati["indice_nodo"]
    n_nodi = len(mesh_dati["nodi"])
    dx_m = mesh_dati["lunghezza_m"] / n_x
    ds_m = circonferenza_m / n_theta
    h_esterno = aria["h_esterno_W_m2K"]
    h_interno = aria["h_interno_W_m2K"]
    h_bordi = aria["h_bordi_W_m2K"]
    T_ambiente = aria["temperatura_ambiente_C"]
    # Matrice di massa di superficie: riusata per la capacità termica, per la
    # convezione sulle due facce e per il carico della sorgente.
    M_area = _assembla(elementi, MASSA_RIF * (dx_m * ds_m), n_nodi)
    K_cond = _assembla(
        elementi,
        k * spessore_m * (RIGIDEZZA_X * (ds_m / dx_m) + RIGIDEZZA_S * (dx_m / ds_m)),
        n_nodi,
    )
    # Convezione sui bordi anulari: l'area di scambio di ogni segmento è
    # spessore × lunghezza del segmento circonferenziale.
    H_bordi = h_bordi * spessore_m * (
        _massa_anello(indice_nodo[0, :], ds_m, n_nodi)
        + _massa_anello(indice_nodo[-1, :], ds_m, n_nodi)
    )
    C_su_dt = (rho * cp * spessore_m / dt_s) * M_area
    A = C_su_dt + K_cond + (h_esterno + h_interno) * M_area + H_bordi
    # Termine noto costante della convezione: (somma delle righe) × h × T_amb,
    # perché la somma delle righe della matrice di massa è l'integrale di N_i.
    area_nodale = np.asarray(M_area.sum(axis=1)).ravel()
    bordo_nodale = np.asarray(H_bordi.sum(axis=1)).ravel()
    carico_convezione = (h_esterno + h_interno) * area_nodale * T_ambiente
    carico_convezione += bordo_nodale * T_ambiente
    # Peso circonferenziale dell'impronta della sorgente: dipende solo
    # dall'offset y del percorso e dalla sigma dello spot, entrambi costanti
    # durante il run, quindi è valutato una volta sola qui. La distanza è
    # quella dell'immagine più vicina, perché la superficie è chiusa.
    sigma_m = sorgente["sigma_punto_m"]
    ds_sorgente = raggio_m * mesh_dati["theta_nodi_rad"] - sorgente["offset_y_percorso_m"]
    ds_sorgente -= circonferenza_m * np.round(ds_sorgente / circonferenza_m)
    peso_circonferenziale = np.exp(-0.5 * (ds_sorgente ** 2) / (sigma_m * sigma_m))
    # Nodo osservato dal sensore: x più vicino a x_mm sul piano theta = 0.
    x_sensore_m = sensore["x_mm"] / 1000.0
    i_sensore = int(np.argmin(np.abs(mesh_dati["x_nodi_m"] - x_sensore_m)))
    indice_sensore = int(indice_nodo[i_sensore, 0])
    return {
        "mesh": mesh_dati,
        "dt_s": dt_s,
        "M_area": M_area,
        "C_su_dt": C_su_dt,
        "carico_convezione": carico_convezione,
        "solutore": splu(sp.csc_matrix(A)),
        "x_nodi_m": mesh_dati["x_nodi_m"],
        "peso_circonferenziale": peso_circonferenziale,
        "forma_griglia": (n_x + 1, n_theta),
        "circonferenza_m": circonferenza_m,
        "area_nodale_m2": area_nodale,
        "T_ambiente_C": T_ambiente,
        "x_sensore_m": x_sensore_m,
        "indice_sensore": indice_sensore,
        "n_nodi": n_nodi,
        # Diagnostica: non entra nell'equazione risolta.
        "skin_depth_m": skin_depth_m,
    }
 def flusso_nodale_W_m2(
    sorgente: dict,
    stato: dict,
    t_s: float,
 ) -> tuple[float, np.ndarray]:
    """Flusso termico assorbito nei nodi all'istante t_s.
    Restituisce (x_riferimento_m, q_nodale) dove x_riferimento_m è la
    posizione della sorgente di indice 0. Ogni sorgente attiva contribuisce
    con un'impronta gaussiana isotropa in (x, s).
    L'impronta gaussiana isotropa è separabile, q(x, s) = q_x(x) · q_s(s), e
    il fattore circonferenziale q_s è costante nel tempo (precalcolato in
    `prepara_stato_fem`). Qui si valuta quindi solo il fattore assiale sui
    n_x + 1 nodi distinti in x e si espande con un prodotto esterno, invece di
    valutare l'esponenziale su tutti i nodi della mesh.
    """
    x_sensore_m = stato["x_sensore_m"]
    x_rif_iniziale = _x_riferimento_iniziale_m(sorgente, x_sensore_m)
    x_rif_finale = _x_riferimento_finale_m(sorgente, x_sensore_m)
    x_riferimento = x_rif_iniziale + sorgente["velocita_m_s"] * t_s
    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
    v = sorgente["velocita_m_s"]
    zero_dopo_fine = sorgente.get("zero_dopo_fine", True)
    sigma = sorgente["sigma_punto_m"]
    q_picco = sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"]
    peso_assiale = np.zeros(len(stato["x_nodi_m"]))
    for i in range(numero_sorgenti):
        x_i = x_riferimento + i * distanza
        if zero_dopo_fine:
            inizio_i = x_rif_iniziale + i * distanza
            fine_i = x_rif_finale + i * distanza
            if not _intervallo_attivo(inizio_i, fine_i, v, x_i):
                continue
        dx = stato["x_nodi_m"] - x_i
        peso_assiale += np.exp(-0.5 * (dx * dx) / (sigma * sigma))
    q_nodale = q_picco * np.outer(peso_assiale, stato["peso_circonferenziale"])
    return x_riferimento, q_nodale.ravel()
 def passo_implicito_fem(stato: dict, T: np.ndarray, q_nodale: np.ndarray) -> np.ndarray:
    # Avanza il campo nodale di un passo dt: termine noto = capacità sul campo
    # precedente + carico della sorgente + carico di convezione, poi risoluzione
    # del sistema già fattorizzato.
    rhs = stato["C_su_dt"] @ T
    rhs += stato["M_area"] @ q_nodale
    rhs += stato["carico_convezione"]
    return stato["solutore"].solve(rhs)
 def campo_iniziale(stato: dict) -> np.ndarray:
    return np.full(stato["n_nodi"], stato["T_ambiente_C"], dtype=float)
 def simula_campo_fem(
    dt_frame_s: float,
    cfg: dict | None = None,
    durata_s: float | None = None,
    mesh_dati: dict | None = None,
    rng: random.Random | None = None,
    salva_campi: bool = True,
 ) -> dict:
    """Integra il campo FEM fino a durata_s, campionando ogni dt_frame_s.
    Il campo nodale è la temperatura vera della parete; la lettura del sensore
    ne è la versione degradata: inerzia del primo ordine, rumore gaussiano e
    quantizzazione, questi ultimi due applicati a ogni campionamento.
    Con `salva_campi = False` i campi nodali non vengono accumulati (utile per
    l'export CSV, che usa solo le serie scalari): `campi` resta una lista vuota.
    """
    if rng is None:
        rng = random.Random(FEM["seed"])
    if cfg is None:
        cfg = configurazione_randomizzata(rng)
    if durata_s is None:
        durata_s = FEM["durata_s"]
    fascetta = cfg["fascetta"]
    aria = cfg["aria"]
    sorgente = cfg["sorgente"]
    sensore = cfg["sensore"]
    stato = prepara_stato_fem(fascetta, aria, sorgente, sensore, mesh_dati=mesh_dati)
    dt = stato["dt_s"]
    indice_sensore = stato["indice_sensore"]
    T = campo_iniziale(stato)
    T_sensore = float(T[indice_sensore])
    tau_sensore = max(sensore["costante_tempo_s"], 1e-9)
    tempi, campi, x_riferimenti = [], [], []
    T_vere, T_lette, T_massime = [], [], []
    prossimo_frame_t = 0.0
    t = 0.0
    while t <= durata_s + 1e-12:
        x_rif, q_nodale = flusso_nodale_W_m2(sorgente, stato, t)
        T = passo_implicito_fem(stato, T, q_nodale)
        T_sensore += (T[indice_sensore] - T_sensore) * dt / tau_sensore
        if t + 1e-12 >= prossimo_frame_t:
            letta = T_sensore + rng.gauss(0.0, sensore["rumore_std_C"])
            letta = quantizza(letta, sensore["quantizzazione_C"])
            tempi.append(t)
            if salva_campi:
                campi.append(T.copy())
            x_riferimenti.append(x_rif)
            T_vere.append(float(T[indice_sensore]))
            T_lette.append(letta)
            T_massime.append(float(T.max()))
            prossimo_frame_t += dt_frame_s
        t += dt
    return {
        "stato": stato,
        "mesh": stato["mesh"],
        "tempi": np.array(tempi),
        "campi": campi,
        "x_riferimenti": np.array(x_riferimenti),
        "T_vere": np.array(T_vere),
        "T_lette": np.array(T_lette),
        "T_massime": np.array(T_massime),
        "cfg": cfg,
        "sorgente": sorgente,
        "T_ambiente_C": stato["T_ambiente_C"],
        "skin_depth_m": stato["skin_depth_m"],
    }
@@ -1,61 +0,0 @@
 # Punto di ingresso unico del progetto: seleziona da riga di comando quale
 # azione eseguire.
 #
 #   python main.py mesh          # disegna la sola mesh a elementi shell
 #   python main.py fem           # animazione 3D del campo FEM sulla mesh shell
 #   python main.py csv           # esporta in CSV la serie temporale del FEM
 #
 # Senza argomenti stampa l'elenco delle azioni disponibili.
 import argparse
 import sys
 def _azione_mesh() -> None:
    import plot_mesh
    plot_mesh.main()
 def _azione_fem() -> None:
    import plot_animazione_fem
    plot_animazione_fem.main()
 def _azione_csv() -> None:
    import esporta_csv
    esporta_csv.main()
 # Chiave da riga di comando -> (funzione, descrizione mostrata nell'help).
 AZIONI = {
    "mesh": (_azione_mesh, "Disegna la mesh a elementi shell della fascetta"),
    "fem": (_azione_fem, "Animazione 3D del campo FEM sulla mesh a elementi shell"),
    "csv": (_azione_csv, "Esporta in CSV la temperatura simulata e la lettura del sensore"),
 }
 def main() -> None:
    parser = argparse.ArgumentParser(
        description="Analisi termica FEM della fascetta riscaldata a induzione.",
        formatter_class=argparse.RawTextHelpFormatter,
    )
    parser.add_argument(
        "azione",
        nargs="?",
        choices=list(AZIONI),
        help="\n".join(f"{nome}: {descrizione}" for nome, (_, descrizione) in AZIONI.items()),
    )
    argomenti = parser.parse_args()
    if argomenti.azione is None:
        parser.print_help()
        sys.exit(1)
    AZIONI[argomenti.azione][0]()
 if __name__ == "__main__":
    main()
@@ -1,4 +1,4 @@
-# Database dei materiali per l'analisi termica.
+# Database dei materiali per il simulatore termico.
 #
 # Tutte le unità sono SI:
 # - conducibilita_termica_W_mK
@@ -1,114 +1,119 @@
-# Generazione della mesh a elementi shell quadrangolari della fascetta.
+# Generazione della mesh a elementi shell quadrilateri sulla superficie
 # cilindrica della fascetta.
 #
-# La fascetta è modellata come la superficie media di un cilindro: gli
+# La superficie è parametrizzata da (x, theta): x è la coordinata assiale,
-# elementi sono quadrilateri a 4 nodi disposti su una griglia strutturata
+# theta quella circonferenziale. Al posto di theta si usa quasi ovunque
-# n_elementi_x × n_elementi_circonferenza. Lo spessore della parete non è
+# l'ascissa curvilinea arco = raggio_medio * theta, così la shell diventa un
-# discretizzato, è un attributo degli elementi shell.
+# rettangolo lunghezza × circonferenza e gli elementi sono tutti identici.
 #
-# Sistema di coordinate globale (coerente con plot_animazione_3d.py):
+# I nodi formano una griglia (n_elementi_x + 1) × n_elementi_theta: in
-# - x = asse del cilindro, da 0 a lunghezza
+# direzione circonferenziale non c'è un nodo finale distinto perché theta = 0
-# - y, z = piano della sezione circolare, con y = R·sin(theta), z = R·cos(theta)
+# e theta = 2*pi sono lo stesso nodo. La periodicità è quindi strutturale,
-# - theta = 0 corrisponde al piano del sensore (y = 0), cresce in verso
+# nasce dalla connettività e non da vincoli imposti a posteriori.
 #   antiorario nel piano y-z
 #
-# La mesh è chiusa lungo la circonferenza: l'ultima colonna di elementi
+# Ordine dei nodi locali di ogni elemento (antiorario nel piano x-arco):
-# richiude sulla prima riga di nodi, senza nodi duplicati sulla cucitura.
+#
 #     4 ---- 3      arco
 #     |      |       ^
 #     |      |       |
 #     1 ---- 2       +---> x
 import numpy as np
 from config import FASCETTA, MESH
 def costruisci_mesh_cilindrica(
    lunghezza_m: float,
    raggio_medio_m: float,
    n_elementi_x: int,
    n_elementi_theta: int,
 ) -> dict:
    # Costruisce nodi e connettività della shell cilindrica.
    n_nodi_x = n_elementi_x + 1
    n_nodi = n_nodi_x * n_elementi_theta
    n_elementi = n_elementi_x * n_elementi_theta
-def raggio_medio_m(fascetta: dict = FASCETTA) -> float:
+    lato_x_m = lunghezza_m / n_elementi_x
-    # Raggio della superficie media della parete: il diametro configurato è
+    circonferenza_m = 2.0 * np.pi * raggio_medio_m
-    # quello esterno, la shell sta a metà dello spessore.
+    lato_arco_m = circonferenza_m / n_elementi_theta
    diametro_m = fascetta["diametro_mm"] / 1000.0
    spessore_m = fascetta["spessore_mm"] / 1000.0
    return (diametro_m - spessore_m) / 2.0
    x_nodi_m = np.arange(n_nodi_x) * lato_x_m
    theta_nodi_rad = np.arange(n_elementi_theta) * (2.0 * np.pi / n_elementi_theta)
    arco_nodi_m = theta_nodi_rad * raggio_medio_m
-def genera_mesh(fascetta: dict = FASCETTA, mesh: dict = MESH) -> dict:
+    # Connettività: elemento (i, j) collega i nodi (i, j), (i+1, j),
-    """Costruisce la mesh shell quadrangolare della fascetta.
+    # (i+1, j+1), (i, j+1), con j+1 riavvolto sulla circonferenza.
    i_elem = np.repeat(np.arange(n_elementi_x), n_elementi_theta)
    j_elem = np.tile(np.arange(n_elementi_theta), n_elementi_x)
    j_succ = (j_elem + 1) % n_elementi_theta
-    Restituisce un dizionario con:
+    def indice(i, j):
-    - nodi: array (n_nodi, 3) con le coordinate globali [m]
+        return i * n_elementi_theta + j
    - 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
+    connettivita = np.stack(
-    spessore_m = fascetta["spessore_mm"] / 1000.0
+        [
-    raggio_m = raggio_medio_m(fascetta)
+            indice(i_elem, j_elem),
-
+            indice(i_elem + 1, j_elem),
-    # Lungo x la griglia è aperta (n_x + 1 file di nodi), lungo theta è
+            indice(i_elem + 1, j_succ),
-    # chiusa (n_theta file, l'ultima si ricongiunge alla prima).
+            indice(i_elem, j_succ),
-    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)
+        axis=1,
    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)
+    # Segmenti dei due bordi assiali (x = 0 e x = lunghezza): anelli chiusi di
-
+    # n_elementi_theta segmenti, usati per la convezione sullo spessore.
-    i = np.arange(n_x)[:, None]
+    j = np.arange(n_elementi_theta)
-    j = np.arange(n_theta)[None, :]
+    j_dopo = (j + 1) % n_elementi_theta
-    j_succ = (j + 1) % n_theta
+    segmenti_bordo = np.concatenate(
-    elementi = np.stack(
+        [
-        (
+            np.stack([indice(0, j), indice(0, j_dopo)], axis=1),
-            np.broadcast_to(indice_nodo[i, j], (n_x, n_theta)),
+            np.stack(
-            np.broadcast_to(indice_nodo[i + 1, j], (n_x, n_theta)),
+                [indice(n_elementi_x, j), indice(n_elementi_x, j_dopo)], axis=1
-            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,
+        "n_elementi_x": n_elementi_x,
-        "elementi": elementi,
+        "n_elementi_theta": n_elementi_theta,
-        "indice_nodo": indice_nodo,
+        "n_elementi": n_elementi,
        "n_nodi_x": n_nodi_x,
        "n_nodi": n_nodi,
        "lunghezza_m": lunghezza_m,
        "raggio_medio_m": raggio_medio_m,
        "circonferenza_m": circonferenza_m,
        "lato_x_m": lato_x_m,
        "lato_arco_m": lato_arco_m,
        "x_nodi_m": x_nodi_m,
        "theta_nodi_rad": theta_nodi_rad,
-        "n_elementi_x": n_x,
+        "arco_nodi_m": arco_nodi_m,
-        "n_elementi_circonferenza": n_theta,
+        "connettivita": connettivita,
-        "spessore_m": spessore_m,
+        "segmenti_bordo": segmenti_bordo,
-        "raggio_m": raggio_m,
+        "area_totale_m2": lunghezza_m * circonferenza_m,
        "lunghezza_m": lunghezza_m,
    }
-def riepilogo_mesh(mesh_dati: dict) -> str:
+def campo_su_griglia(mesh: dict, T: np.ndarray) -> np.ndarray:
-    n_nodi = len(mesh_dati["nodi"])
+    # Rimappa il vettore nodale sulla griglia (n_nodi_x, n_elementi_theta),
-    n_elementi = len(mesh_dati["elementi"])
+    # comoda per le mappe sviluppate e per le superfici 3D.
-    passo_x_mm = 1000.0 * mesh_dati["lunghezza_m"] / mesh_dati["n_elementi_x"]
+    return T.reshape(mesh["n_nodi_x"], mesh["n_elementi_theta"])
-    passo_circ_mm = (
+
-        1000.0
+
-        * 2.0
+def coordinate_3d(mesh: dict, chiudi_circonferenza: bool = True):
-        * np.pi
+    # Coordinate cartesiane dei nodi per la vista 3D: x lungo l'asse del
-        * mesh_dati["raggio_m"]
+    # cilindro, y e z sulla sezione circolare.
-        / mesh_dati["n_elementi_circonferenza"]
+    #
-    )
+    # Con chiudi_circonferenza si ripete la prima colonna in coda, così la
-    return (
+    # superficie disegnata non mostra una fessura in theta = 0. È solo una
-        f"Mesh shell: {n_elementi} elementi quadrangolari, {n_nodi} nodi\n"
+    # necessità di disegno: il nodo ripetuto non è un grado di libertà.
-        f"  griglia: {mesh_dati['n_elementi_x']} (x) × "
+    theta = mesh["theta_nodi_rad"]
-        f"{mesh_dati['n_elementi_circonferenza']} (circonferenza)\n"
+    if chiudi_circonferenza:
-        f"  passo elemento: {passo_x_mm:.2f} mm (x), {passo_circ_mm:.2f} mm (circonferenza)\n"
+        theta = np.append(theta, 2.0 * np.pi)
-        f"  raggio medio: {1000.0 * mesh_dati['raggio_m']:.2f} mm, "
+    X, Theta = np.meshgrid(mesh["x_nodi_m"], theta, indexing="ij")
-        f"spessore shell: {1000.0 * mesh_dati['spessore_m']:.3f} mm"
+    R = mesh["raggio_medio_m"]
-    )
+    return X, R * np.sin(Theta), R * np.cos(Theta)
 def chiudi_campo(campo: np.ndarray) -> np.ndarray:
    # Ripete la prima colonna circonferenziale in coda, in accordo con
    # coordinate_3d(chiudi_circonferenza=True).
    return np.concatenate([campo, campo[:, :1]], axis=1)
@@ -0,0 +1,281 @@
 # Animazione della superficie della fascetta durante il passaggio delle
 # sorgenti.
 #
 # Riproduce la fisica di run_0001 (stesso seed di simulate.py) e mostra:
 # - il profilo di flusso termico q(x) lungo la generatrice delle sorgenti;
 # - il campo di temperatura sulla superficie cilindrica sviluppata in piano,
 #   cioè T(x, arco) con l'arco centrato sulla generatrice del sensore;
 # - il sensore infrarosso e la temperatura nel punto osservato.
 #
 # La mappa sviluppata è periodica in verticale: il bordo superiore e quello
 # inferiore sono lo stesso punto della fascetta.
 import random
 from pathlib import Path
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib.animation import FuncAnimation, PillowWriter
 from matplotlib.collections import LineCollection
 import sensore as sens
 import sorgente as src
 from config import SIMULAZIONE
 from mesh import campo_su_griglia, chiudi_campo
 from simulate import (
    campo_iniziale,
    configurazione_randomizzata,
    passo_termico,
    prepara_stato_termico,
 )
 # Istante di inizio dei fotogrammi mostrati (la simulazione parte comunque da 0).
 T_INIZIO_ANIMAZIONE_S = 0.40
 # Istante di fine dell'animazione.
 T_FINE_ANIMAZIONE_S = 30
 # Tempo simulato tra un fotogramma e il successivo.
 DT_FRAME_S = 0.05
 # Millisecondi tra i fotogrammi in riproduzione.
 INTERVALLO_RIPRODUZIONE_MS = 30
 # Disegna i bordi degli elementi shell sopra il campo di temperatura.
 # Con la mesh di default (100 x 110, elementi di 1 x 2 mm) il reticolo è più
 # fitto dei pixel disponibili e copre il campo: per ispezionare la mesh
 # conviene attivarlo e insieme abbassare n_elementi_x e n_elementi_theta in
 # config.py, per esempio a 25 x 28.
 MOSTRA_MESH = False
 def simula_campi(cfg_run: dict) -> dict:
    # Esegue la simulazione fino a T_FINE_ANIMAZIONE_S salvando, a ogni
    # fotogramma, campo di temperatura, profilo di flusso e stato del sensore.
    #
    # I campi sono restituiti come griglie (n_nodi_x, n_elementi_theta) con la
    # colonna dell'arco ruotata in modo da portare la generatrice del sensore
    # (theta = 0) al centro: è solo una scelta di visualizzazione.
    fascetta = cfg_run["fascetta"]
    aria = cfg_run["aria"]
    sorgente = cfg_run["sorgente"]
    sensore = cfg_run["sensore"]
    stato = prepara_stato_termico(fascetta, aria, sorgente, sensore)
    mesh = stato["mesh"]
    dt = stato["dt_s"]
    interpolatore = stato["interpolatore_sensore"]
    # Dopo np.roll(..., rotazione) la colonna theta = 0 finisce all'indice
    # "rotazione", quindi l'arco con segno della colonna m vale
    # (m - rotazione) * lato_arco.
    rotazione = mesh["n_elementi_theta"] // 2
    arco_centrato_m = (
        np.arange(mesh["n_elementi_theta"]) - rotazione
    ) * mesh["lato_arco_m"]
    T = campo_iniziale(stato)
    T_letta = sens.temperatura_osservata_C(interpolatore, T)
    tempi, campi, flussi, x_riferimenti = [], [], [], []
    T_vere, T_lette = [], []
    prossimo_frame_t = 0.0
    t = 0.0
    while t <= T_FINE_ANIMAZIONE_S + 1e-12:
        T, info = passo_termico(stato, T, t)
        T_vera = sens.temperatura_osservata_C(interpolatore, T)
        T_letta = sens.aggiorna_inerzia_C(
            T_letta, T_vera, dt, sensore["costante_tempo_s"]
        )
        if t + 1e-12 >= prossimo_frame_t:
            _, q_nodi = src.profilo_flusso_x_W_m2(
                sorgente, stato["x_sensore_m"], t, mesh["x_nodi_m"]
            )
            tempi.append(t)
            campi.append(np.roll(campo_su_griglia(mesh, T), rotazione, axis=1))
            flussi.append(q_nodi)
            x_riferimenti.append(info["x_riferimento_m"])
            T_vere.append(T_vera)
            T_lette.append(T_letta)
            prossimo_frame_t += DT_FRAME_S
        t += dt
    return {
        "tempi": np.array(tempi),
        "campi": campi,
        "flussi": flussi,
        "x_riferimenti": np.array(x_riferimenti),
        "T_vere": np.array(T_vere),
        "T_lette": np.array(T_lette),
        "mesh": mesh,
        "rotazione": rotazione,
        "x_nodi_mm": mesh["x_nodi_m"] * 1000.0,
        "arco_centrato_mm": arco_centrato_m * 1000.0,
        "arco_chiuso_mm": np.append(
            arco_centrato_m, arco_centrato_m[-1] + mesh["lato_arco_m"]
        ) * 1000.0,
        "lunghezza_mm": fascetta["lunghezza_mm"],
        "x_sensore_mm": sensore["x_mm"],
        "offset_sorgente_mm": sorgente["offset_y_percorso_m"] * 1000.0,
        "sorgente": sorgente,
        "T_ambiente_C": aria["temperatura_ambiente_C"],
    }
 def main() -> None:
    rng = random.Random(SIMULAZIONE["seed"])
    cfg_run = configurazione_randomizzata(1, rng)
    dati = simula_campi(cfg_run)
    tempi = dati["tempi"]
    indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S))
    n_frame = len(tempi) - indice_inizio
    lunghezza_mm = dati["lunghezza_mm"]
    x_vista_mm = (-10.0, lunghezza_mm + 10.0)
    arco_mm = dati["arco_chiuso_mm"]
    arco_min_mm, arco_max_mm = arco_mm[0], arco_mm[-1]
    # Il minimo evita un asse degenere se la finestra dell'animazione si
    # chiude prima che le sorgenti entrino nella fascetta.
    q_max_MW = max(max(q.max() for q in dati["flussi"]) / 1e6, 1e-3)
    T_max = max(c.max() for c in dati["campi"])
    sorgente = dati["sorgente"]
    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza_mm = sorgente.get("distanza_sorgenti_m", 0.0) * 1000.0
    offset_mm = dati["offset_sorgente_mm"]
    fig, (ax_flusso, ax_superficie, ax_storia) = plt.subplots(
        3, 1, figsize=(10, 8), height_ratios=[1.0, 1.6, 1.2],
        gridspec_kw={"hspace": 0.45},
    )
    fig.suptitle("Superficie sviluppata della fascetta: sorgenti in transito e sensore")
    # Pannello 1: profilo di flusso lungo x e posizioni delle sorgenti.
    linea_flusso, = ax_flusso.plot([], [], color="tab:red")
    marker_sorgenti, = ax_flusso.plot(
        [], [], "v", color="tab:red", markersize=10, clip_on=False
    )
    ax_flusso.annotate(
        "verso di marcia",
        xy=(0.28, 0.85), xytext=(0.55, 0.85), xycoords="axes fraction",
        arrowprops={"arrowstyle": "->", "color": "gray"},
        color="gray", va="center",
    )
    ax_flusso.set_xlim(*x_vista_mm)
    ax_flusso.set_ylim(0.0, q_max_MW * 1.25)
    ax_flusso.set_ylabel("q(x) [MW/m²]")
    ax_flusso.grid(True, alpha=0.3)
    ax_flusso.set_xticklabels([])
    # Pannello 2: temperatura sulla superficie sviluppata. L'asse verticale è
    # l'arco circonferenziale, con 0 sulla generatrice del sensore.
    #
    # pcolormesh con shading "gouraud" interpola linearmente tra i valori
    # nodali: è esattamente ciò che fanno le funzioni di forma bilineari
    # dell'elemento Q4, quindi il disegno è la soluzione dell'elemento finito
    # e non un'interpolazione grafica arbitraria.
    Xg, Yg = np.meshgrid(dati["x_nodi_mm"], dati["arco_chiuso_mm"], indexing="ij")
    mappa = ax_superficie.pcolormesh(
        Xg, Yg,
        chiudi_campo(dati["campi"][indice_inizio]),
        shading="gouraud",
        cmap="inferno",
        vmin=dati["T_ambiente_C"],
        vmax=T_max,
    )
    if MOSTRA_MESH:
        # Bordi degli elementi shell: linee lungo x a ogni nodo
        # circonferenziale e viceversa.
        segmenti = [
            [(x, arco_min_mm), (x, arco_max_mm)] for x in dati["x_nodi_mm"]
        ] + [
            [(0.0, a), (lunghezza_mm, a)] for a in dati["arco_chiuso_mm"]
        ]
        ax_superficie.add_collection(
            LineCollection(segmenti, colors="white", linewidths=0.5, alpha=0.55)
        )
    ax_superficie.set_xlim(*x_vista_mm)
    ax_superficie.set_ylim(arco_min_mm, arco_max_mm)
    ax_superficie.set_ylabel("arco circonferenziale [mm]")
    ax_superficie.set_xlabel("x [mm]")
    ax_superficie.plot(
        [dati["x_sensore_mm"]], [0.0], "^", color="tab:blue", markersize=11
    )
    ax_superficie.text(
        dati["x_sensore_mm"] + 3, 0.0, "sensore IR", color="tab:blue", va="center"
    )
    ax_superficie.axhline(
        offset_mm, color="tab:cyan", linewidth=0.8, linestyle="--", alpha=0.7
    )
    # La colorbar è agganciata a tutti i pannelli per non restringere solo
    # quello della superficie, mantenendo allineati gli assi x.
    barra = fig.colorbar(
        mappa, ax=(ax_flusso, ax_superficie, ax_storia), pad=0.02, aspect=35
    )
    barra.set_label("T [°C]")
    # Pannello 3: temperatura nel punto osservato dal sensore.
    linea_vera, = ax_storia.plot([], [], label="T vera lato interno")
    linea_letta, = ax_storia.plot([], [], label="T sensore (con inerzia)")
    cursore = ax_storia.axvline(tempi[indice_inizio], color="gray", linewidth=0.8)
    ax_storia.set_xlim(0.0, T_FINE_ANIMAZIONE_S)
    ax_storia.set_ylim(15.0, max(dati["T_vere"].max(), dati["T_lette"].max()) * 1.08)
    ax_storia.set_xlabel("Tempo [s]")
    ax_storia.set_ylabel("T [°C]")
    ax_storia.legend(loc="upper left")
    ax_storia.grid(True, alpha=0.3)
    testo_tempo = ax_flusso.set_title(f"t = {tempi[indice_inizio]:.3f} s", loc="right")
    def aggiorna(frame: int):
        k = indice_inizio + frame
        t = tempi[k]
        linea_flusso.set_data(dati["x_nodi_mm"], dati["flussi"][k] / 1e6)
        x_sorgenti_mm = (
            dati["x_riferimenti"][k] * 1000.0
            + np.arange(numero_sorgenti) * distanza_mm
        )
        visibili = (x_sorgenti_mm >= x_vista_mm[0]) & (x_sorgenti_mm <= x_vista_mm[1])
        marker_sorgenti.set_data(
            x_sorgenti_mm[visibili],
            np.full(int(visibili.sum()), q_max_MW * 1.12),
        )
        mappa.set_array(chiudi_campo(dati["campi"][k]))
        linea_vera.set_data(tempi[: k + 1], dati["T_vere"][: k + 1])
        linea_letta.set_data(tempi[: k + 1], dati["T_lette"][: k + 1])
        cursore.set_xdata([t, t])
        testo_tempo.set_text(f"t = {t:.3f} s")
        return (
            linea_flusso, marker_sorgenti, mappa,
            linea_vera, linea_letta, cursore, testo_tempo,
        )
    animazione = FuncAnimation(
        fig, aggiorna, frames=n_frame, interval=INTERVALLO_RIPRODUZIONE_MS, blit=False
    )
    # Se il backend non è interattivo si salva una GIF invece di mostrare la finestra.
    if matplotlib.get_backend().lower() == "agg":
        percorso = Path("dataset") / "animazione_superficie.gif"
        animazione.save(
            percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS)
        )
        print(f"Backend non interattivo: animazione salvata in {percorso}")
        return
    plt.show()
 if __name__ == "__main__":
    main()
@@ -0,0 +1,140 @@
 # Animazione 3D del barattolo in vista isometrica: colormap della temperatura
 # sulla superficie cilindrica.
 #
 # Il modello shell risolve la temperatura sull'intera superficie, quindi il
 # campo disegnato qui è quello calcolato: non c'è più nessuna ricostruzione
 # della distribuzione circonferenziale. La faccia esterna e quella interna
 # hanno la stessa temperatura (una temperatura per nodo), quindi la superficie
 # mostrata vale per entrambe.
 import random
 from pathlib import Path
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib import cm
 from matplotlib.animation import FuncAnimation, PillowWriter
 import sorgente as src
 from config import SIMULAZIONE
 from mesh import chiudi_campo, coordinate_3d
 from plot_animazione import DT_FRAME_S as DT_FRAME_S_BASE
 from plot_animazione import T_INIZIO_ANIMAZIONE_S, simula_campi
 from simulate import configurazione_randomizzata
 # Tempo simulato tra un fotogramma e il successivo dell'animazione 3D. Più
 # largo di DT_FRAME_S_BASE perché ricreare una superficie 3D a ogni
 # fotogramma è più costoso della semplice imshow 2D.
 DT_FRAME_S = 0.15
 # Millisecondi tra i fotogrammi in riproduzione.
 INTERVALLO_RIPRODUZIONE_MS = 60
 # Disegna i bordi degli elementi shell sulla superficie del cilindro.
 # Con la mesh di default (100 x 110, elementi di 1 x 2 mm) il reticolo satura
 # il disegno e sbiadisce i colori: per ispezionare la mesh conviene attivarlo
 # e insieme abbassare n_elementi_x e n_elementi_theta in config.py, per
 # esempio a 25 x 28.
 MOSTRA_MESH = False
 def main() -> None:
    rng = random.Random(SIMULAZIONE["seed"])
    cfg_run = configurazione_randomizzata(1, rng)
    dati = simula_campi(cfg_run)
    tempi = dati["tempi"]
    indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S))
    passo = max(1, round(DT_FRAME_S / DT_FRAME_S_BASE))
    indici_frame = list(range(indice_inizio, len(tempi), passo))
    mesh = dati["mesh"]
    # simula_campi ruota le colonne per centrare il sensore nella mappa
    # sviluppata: qui serve l'ordine originale della mesh.
    rotazione = dati["rotazione"]
    T_ambiente = dati["T_ambiente_C"]
    sorgente = dati["sorgente"]
    raggio_m = mesh["raggio_medio_m"]
    lunghezza_m = mesh["lunghezza_m"]
    Xm, Ym, Zm = coordinate_3d(mesh)
    T_max = max(c.max() for c in dati["campi"])
    # vmin più basso della temperatura ambiente reale: altrimenti le zone
    # fredde cadrebbero sul nero puro di "inferno" e, essendo lo shading
    # moltiplicativo, nessuna illuminazione basterebbe a renderle visibili.
    norm = matplotlib.colors.Normalize(
        vmin=T_ambiente - 0.4 * (T_max - T_ambiente), vmax=T_max
    )
    cmap = matplotlib.colormaps["inferno"]
    lightsource = matplotlib.colors.LightSource(azdeg=315, altdeg=45)
    fig = plt.figure(figsize=(9, 7))
    ax = fig.add_subplot(projection="3d")
    ax.view_init(elev=35.264, azim=45)
    ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m))
    ax.set_xlabel("x [m]")
    ax.set_axis_off()
    # La colorbar mostra il range reale delle temperature: il norm esteso
    # verso il basso serve solo a schiarire il colore di base della
    # superficie fredda, non deve comparire nella scala mostrata all'utente.
    norm_colorbar = matplotlib.colors.Normalize(vmin=T_ambiente, vmax=T_max)
    mappabile = cm.ScalarMappable(cmap=cmap, norm=norm_colorbar)
    mappabile.set_array([])
    barra = fig.colorbar(mappabile, ax=ax, shrink=0.6, pad=0.05)
    barra.set_label("T [°C]")
    theta_sorgente = src.theta_percorso_rad(sorgente, raggio_m)
    def disegna_frame(k: int):
        ax.cla()
        ax.view_init(elev=35.264, azim=45)
        ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m))
        ax.set_axis_off()
        campo = np.roll(dati["campi"][k], -rotazione, axis=1)
        colori = cmap(norm(chiudi_campo(campo)))
        # Ogni faccia disegnata è un elemento shell della mesh: con
        # MOSTRA_MESH se ne vedono i bordi.
        bordi = {"edgecolor": (1.0, 1.0, 1.0, 0.25), "linewidth": 0.15}
        ax.plot_surface(
            Xm, Ym, Zm, facecolors=colori, rstride=1, cstride=1,
            antialiased=False, shade=True, lightsource=lightsource,
            **(bordi if MOSTRA_MESH else {"linewidth": 0}),
        )
        x_sorgenti_m = src.posizioni_sorgenti_m(sorgente, dati["x_riferimenti"][k])
        visibili = (x_sorgenti_m >= 0.0) & (x_sorgenti_m <= lunghezza_m)
        if visibili.any():
            ax.scatter(
                x_sorgenti_m[visibili],
                np.full(int(visibili.sum()), raggio_m * np.sin(theta_sorgente) * 1.05),
                np.full(int(visibili.sum()), raggio_m * np.cos(theta_sorgente) * 1.05),
                color="cyan", s=25, depthshade=False,
            )
        ax.set_title(f"t = {tempi[k]:.3f} s")
        return ()
    animazione = FuncAnimation(
        fig, disegna_frame, frames=indici_frame,
        interval=INTERVALLO_RIPRODUZIONE_MS, blit=False,
    )
    if matplotlib.get_backend().lower() == "agg":
        percorso = Path("dataset") / "animazione_3d.gif"
        animazione.save(
            percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS)
        )
        print(f"Backend non interattivo: animazione salvata in {percorso}")
        return
    plt.show()
 if __name__ == "__main__":
    main()
@@ -1,198 +0,0 @@
 # Animazione 3D isometrica della fascetta con il campo di temperatura calcolato
 # a elementi finiti sulla mesh shell (fem.py).
 #
 # La coordinata circonferenziale è una direzione risolta del modello FEM,
 # quindi la temperatura dipinta su ogni elemento è quella effettivamente
 # calcolata, senza ricostruzioni.
 #
 # Ogni elemento shell è disegnato come faccia piana con colore pari alla media
 # dei suoi quattro valori nodali. L'ombreggiatura è calcolata a mano dalla
 # normale radiale dell'elemento: Poly3DCollection ombreggia solo alla
 # creazione, mentre qui i colori delle facce cambiano a ogni fotogramma.
 from pathlib import Path
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib import cm
 from matplotlib.animation import FuncAnimation, PillowWriter
 from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 from config import USCITA
 from fem import simula_campo_fem
 from mesh import genera_mesh, riepilogo_mesh
 # Istante del primo fotogramma mostrato: i primi istanti sono ancora uniformi
 # alla temperatura ambiente.
 T_INIZIO_ANIMAZIONE_S = 0.40
 # Tempo simulato tra un fotogramma e il successivo.
 DT_FRAME_S = 0.15
 # Millisecondi tra i fotogrammi in riproduzione.
 INTERVALLO_RIPRODUZIONE_MS = 60
 # Direzione della luce (x = asse del cilindro, y-z = piano della sezione) e
 # quota di luce ambiente, per dare volume alla superficie.
 DIREZIONE_LUCE = np.array([-0.25, 0.45, 0.85])
 LUCE_AMBIENTE = 0.45
 # Se True disegna il reticolo degli elementi sopra il campo. Con la mesh fitta
 # i bordi degli elementi sul lato nascosto rendono la superficie confusa,
 # quindi per default i bordi prendono il colore della faccia: così le facce
 # adiacenti combaciano senza lasciare fessure di antialiasing.
 MOSTRA_BORDI_ELEMENTI = False
 def _illuminazione(mesh_dati: dict) -> np.ndarray:
    # Fattore moltiplicativo di luminosità per elemento, dalla normale
    # radiale uscente valutata nel baricentro dell'elemento.
    facce = mesh_dati["nodi"][mesh_dati["elementi"]]
    baricentri = facce.mean(axis=1)
    normali = baricentri.copy()
    # La normale della superficie cilindrica è radiale: nessuna componente
    # lungo l'asse x.
    normali[:, 0] = 0.0
    normali /= np.linalg.norm(normali, axis=1, keepdims=True)
    luce = DIREZIONE_LUCE / np.linalg.norm(DIREZIONE_LUCE)
    diffusa = np.clip(normali @ luce, 0.0, None)
    return LUCE_AMBIENTE + (1.0 - LUCE_AMBIENTE) * diffusa
 def main() -> None:
    mesh_dati = genera_mesh()
    print(riepilogo_mesh(mesh_dati))
    print("Integrazione FEM in corso...")
    dati = simula_campo_fem(dt_frame_s=DT_FRAME_S, mesh_dati=mesh_dati)
    tempi = dati["tempi"]
    indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S))
    indici_frame = list(range(indice_inizio, len(tempi)))
    T_ambiente = dati["T_ambiente_C"]
    T_max = max(campo.max() for campo in dati["campi"])
    print(
        f"T massima nodale: {T_max:.1f} °C — "
        f"T massima nel punto del sensore: {dati['T_vere'].max():.1f} °C\n"
        f"Skin depth (diagnostica): {1000.0 * dati['skin_depth_m']:.3f} mm"
    )
    elementi = mesh_dati["elementi"]
    facce = mesh_dati["nodi"][elementi]
    illuminazione = _illuminazione(mesh_dati)[:, None]
    raggio_m = mesh_dati["raggio_m"]
    lunghezza_m = mesh_dati["lunghezza_m"]
    sorgente = dati["sorgente"]
    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza_m = sorgente.get("distanza_sorgenti_m", 0.0)
    theta_sorgente = sorgente["offset_y_percorso_m"] / raggio_m
    norm = matplotlib.colors.Normalize(vmin=T_ambiente, vmax=T_max)
    cmap = matplotlib.colormaps["inferno"]
    fig = plt.figure(figsize=(9, 8.5))
    griglia = fig.add_gridspec(2, 1, height_ratios=[2.6, 1.0], hspace=0.05)
    ax = fig.add_subplot(griglia[0], projection="3d")
    ax_storia = fig.add_subplot(griglia[1])
    ax.view_init(elev=35.264, azim=45)
    # zoom > 1 riempie il riquadro: con set_axis_off gli assi 3D lascerebbero
    # altrimenti molto margine vuoto sopra il pannello della storia.
    ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m), zoom=1.12)
    ax.set_axis_off()
    ax.set_xlim(0.0, lunghezza_m)
    ax.set_ylim(-raggio_m, raggio_m)
    ax.set_zlim(-raggio_m, raggio_m)
    collezione = Poly3DCollection(facce, linewidths=0.25, shade=False)
    ax.add_collection3d(collezione)
    marker_sorgenti = ax.plot(
        [], [], [], "o", color="cyan", markersize=6, zorder=10
    )[0]
    mappabile = cm.ScalarMappable(cmap=cmap, norm=norm)
    mappabile.set_array([])
    # La colorbar è agganciata a entrambi i pannelli per non restringere solo
    # quello 3D: così i due riquadri restano della stessa larghezza.
    barra = fig.colorbar(mappabile, ax=(ax, ax_storia), shrink=0.6, pad=0.02, aspect=30)
    barra.set_label("T [°C]")
    titolo = ax.set_title("")
    # Pannello inferiore: storia della temperatura nel punto osservato dal
    # sensore. T_vere è il valore nodale calcolato dal FEM, T_lette lo stesso
    # segnale degradato dal sensore reale (inerzia, rumore, quantizzazione).
    # A parte il rumore le due curve sono quasi sovrapposte: la prima è
    # tracciata spessa e trasparente perché la seconda resti leggibile sopra
    # di essa.
    linea_vera, = ax_storia.plot(
        [], [], color="tab:blue", linewidth=3.0, alpha=0.4,
        label="T nodo del sensore (FEM)",
    )
    linea_letta, = ax_storia.plot(
        [], [], color="tab:red", linewidth=1.0,
        label="T misurata dal sensore (inerzia + rumore + quantizzazione)",
    )
    cursore = ax_storia.axvline(tempi[indice_inizio], color="gray", linewidth=0.8)
    ax_storia.set_xlim(0.0, tempi[-1])
    ax_storia.set_ylim(
        T_ambiente - 0.05 * (T_max - T_ambiente),
        max(dati["T_vere"].max(), dati["T_lette"].max()) * 1.08,
    )
    ax_storia.set_xlabel("Tempo [s]")
    ax_storia.set_ylabel("T [°C]")
    ax_storia.legend(loc="upper left")
    ax_storia.grid(True, alpha=0.3)
    def disegna_frame(k: int):
        T_elementi = dati["campi"][k][elementi].mean(axis=1)
        colori = cmap(norm(T_elementi))
        # Ombreggiatura moltiplicativa sui soli canali RGB.
        colori[:, :3] *= illuminazione
        collezione.set_facecolor(colori)
        collezione.set_edgecolor("#40404060" if MOSTRA_BORDI_ELEMENTI else colori)
        x_sorgenti_m = dati["x_riferimenti"][k] + np.arange(numero_sorgenti) * distanza_m
        visibili = (x_sorgenti_m >= 0.0) & (x_sorgenti_m <= lunghezza_m)
        n_visibili = int(visibili.sum())
        marker_sorgenti.set_data_3d(
            x_sorgenti_m[visibili],
            np.full(n_visibili, raggio_m * np.sin(theta_sorgente) * 1.06),
            np.full(n_visibili, raggio_m * np.cos(theta_sorgente) * 1.06),
        )
        linea_vera.set_data(tempi[: k + 1], dati["T_vere"][: k + 1])
        linea_letta.set_data(tempi[: k + 1], dati["T_lette"][: k + 1])
        cursore.set_xdata([tempi[k], tempi[k]])
        titolo.set_text(
            f"FEM shell — t = {tempi[k]:.2f} s — "
            f"T sensore = {dati['T_lette'][k]:.1f} °C"
        )
        return collezione, marker_sorgenti, linea_vera, linea_letta, cursore, titolo
    animazione = FuncAnimation(
        fig, disegna_frame, frames=indici_frame,
        interval=INTERVALLO_RIPRODUZIONE_MS, blit=False,
    )
    if matplotlib.get_backend().lower() == "agg":
        cartella = Path(USCITA["cartella"])
        cartella.mkdir(parents=True, exist_ok=True)
        percorso = cartella / "animazione_fem.gif"
        animazione.save(
            percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS)
        )
        print(f"Backend non interattivo: animazione salvata in {percorso}")
        return
    plt.show()
 if __name__ == "__main__":
    main()
@@ -0,0 +1,48 @@
 from pathlib import Path
 import matplotlib
 import matplotlib.pyplot as plt
 import pandas as pd
 def main() -> None:
    percorso_csv = Path("dataset/run_0001.csv")
    if not percorso_csv.exists():
        raise FileNotFoundError(
            "File CSV non trovato. Esegui prima `python simulate.py`."
        )
    df = pd.read_csv(percorso_csv)
    fig_temperature = plt.figure()
    plt.plot(df["tempo_s"], df["T_vera_lato_sensore_C"], label="Temperatura vera lato interno")
    plt.plot(df["tempo_s"], df["T_misurata_sensore_C"], label="Temperatura misurata dal sensore")
    plt.xlabel("Tempo [s]")
    plt.ylabel("Temperatura [°C]")
    plt.title("Sensore fisso, sorgenti in moto")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    fig_flusso = plt.figure()
    plt.plot(df["tempo_s"], df["flusso_termico_sorgente_W_m2"])
    plt.xlabel("Tempo [s]")
    plt.ylabel("Flusso termico efficace [W/m²]")
    plt.title("Flusso termico nel punto x osservato dal sensore")
    plt.grid(True)
    plt.tight_layout()
    # Se il backend non è interattivo (nessun display o backend GUI
    # disponibile), plt.show() non aprirebbe nulla: si salvano i PNG.
    if matplotlib.get_backend().lower() == "agg":
        for fig, nome in [(fig_temperature, "temperature"), (fig_flusso, "flusso")]:
            percorso = Path("dataset") / f"grafico_{nome}.png"
            fig.savefig(percorso, dpi=150)
            print(f"Backend non interattivo: grafico salvato in {percorso}")
        return
    plt.show()
 if __name__ == "__main__":
    main()
@@ -1,60 +0,0 @@
 # Visualizzazione della sola mesh a elementi shell della fascetta, in vista
 # isometrica: facce chiare con i bordi degli elementi in evidenza.
 from pathlib import Path
 import matplotlib
 import matplotlib.pyplot as plt
 from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 from config import USCITA
 from mesh import genera_mesh, riepilogo_mesh
 def main() -> None:
    dati = genera_mesh()
    print(riepilogo_mesh(dati))
    nodi = dati["nodi"]
    facce = nodi[dati["elementi"]]
    raggio_m = dati["raggio_m"]
    lunghezza_m = dati["lunghezza_m"]
    fig = plt.figure(figsize=(9, 7))
    ax = fig.add_subplot(projection="3d")
    ax.view_init(elev=35.264, azim=45)
    ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m))
    ax.set_axis_off()
    collezione = Poly3DCollection(
        facce,
        facecolors="#d9e3f0",
        edgecolors="#2b3a55",
        linewidths=0.4,
        shade=True,
        lightsource=matplotlib.colors.LightSource(azdeg=315, altdeg=45),
    )
    ax.add_collection3d(collezione)
    ax.set_xlim(0.0, lunghezza_m)
    ax.set_ylim(-raggio_m, raggio_m)
    ax.set_zlim(-raggio_m, raggio_m)
    ax.set_title(
        f"Mesh shell: {dati['n_elementi_x']} × "
        f"{dati['n_elementi_circonferenza']} elementi quadrangolari"
    )
    if matplotlib.get_backend().lower() == "agg":
        cartella = Path(USCITA["cartella"])
        cartella.mkdir(parents=True, exist_ok=True)
        percorso = cartella / "mesh.png"
        fig.savefig(percorso, dpi=150)
        print(f"Backend non interattivo: immagine salvata in {percorso}")
        return
    plt.show()
 if __name__ == "__main__":
    main()
@@ -0,0 +1,60 @@
 # Pirometro a infrarossi montato all'interno della fascetta.
 #
 # Il sensore osserva senza contatto un punto fisso della superficie interna,
 # individuato dalla coordinata assiale x e dalla posizione circonferenziale
 # theta = 0, che per convenzione è la generatrice di riferimento del modello.
 #
 # Con una shell a una temperatura per nodo la faccia interna e quella esterna
 # hanno la stessa temperatura: il valore osservato è la temperatura della
 # shell nel punto, interpolata con le funzioni di forma dell'elemento che lo
 # contiene. L'approssimazione regge finché la parete è termicamente sottile
 # rispetto al transito della sorgente (vedi numero_fourier_spessore in
 # simulate.py).
 #
 # Alla temperatura vera si sovrappone la catena di misura: inerzia del primo
 # ordine, rumore gaussiano e quantizzazione della lettura.
 import numpy as np
 import elementi_shell as el
 def prepara_interpolatore(mesh: dict, x_m: float, theta_rad: float) -> dict:
    # Individua l'elemento che contiene il punto osservato e i pesi nodali
    # con cui ricostruirvi la temperatura.
    n_x = mesh["n_elementi_x"]
    n_theta = mesh["n_elementi_theta"]
    x_limitato = min(max(x_m, 0.0), mesh["lunghezza_m"])
    i = min(n_x - 1, max(0, int(x_limitato / mesh["lato_x_m"])))
    xi = 2.0 * (x_limitato - mesh["x_nodi_m"][i]) / mesh["lato_x_m"] - 1.0
    arco = (theta_rad * mesh["raggio_medio_m"]) % mesh["circonferenza_m"]
    j = int(arco / mesh["lato_arco_m"]) % n_theta
    eta = 2.0 * (arco - mesh["arco_nodi_m"][j]) / mesh["lato_arco_m"] - 1.0
    elemento = i * n_theta + j
    return {
        "nodi": mesh["connettivita"][elemento],
        "pesi": el.funzioni_forma(xi, eta),
        "x_m": x_limitato,
        "theta_rad": theta_rad,
    }
 def temperatura_osservata_C(interpolatore: dict, T: np.ndarray) -> float:
    return float(T[interpolatore["nodi"]] @ interpolatore["pesi"])
 def aggiorna_inerzia_C(
    T_letta_C: float, T_vera_C: float, dt_s: float, costante_tempo_s: float
 ) -> float:
    # Risposta del primo ordine del sensore.
    tau = max(costante_tempo_s, 1e-9)
    return T_letta_C + (T_vera_C - T_letta_C) * dt_s / tau
 def quantizza(valore: float, passo: float) -> float:
    if passo <= 0.0:
        return valore
    return round(valore / passo) * passo
@@ -0,0 +1,386 @@
 # Motore principale: genera N run randomizzati della fascetta riscaldata da
 # sorgenti a induzione in movimento e scrive i CSV del dataset.
 #
 # Il modello termico è una shell a elementi finiti quadrilateri sull'intera
 # superficie cilindrica, con una temperatura per nodo. La pipeline di un run è
 # in prepara_stato_termico() e passo_termico(), condivise con i moduli di
 # animazione: ogni modifica alla fisica va fatta lì, non duplicata.
 import csv
 import math
 import os
 import random
 import shutil
 from concurrent.futures import ProcessPoolExecutor
 from copy import deepcopy
 from pathlib import Path
 import numpy as np
 import assemblaggio
 import sensore as sens
 import solutore as sol
 import sorgente as src
 from config import ARIA, FASCETTA, RANDOMIZZAZIONE, SENSORE, SIMULAZIONE, SORGENTE
 from materials import MATERIALI
 from mesh import costruisci_mesh_cilindrica
 MU0 = 4.0 * math.pi * 1e-7
 def calcola_skin_depth_m(materiale: dict, frequenza_hz: float) -> float:
    # Skin depth elettromagnetica approssimata:
    # delta = sqrt(2 * rho_e / (omega * mu))
    #
    # Semplificata. Per acciai ferromagnetici il comportamento reale
    # è fortemente non lineare con temperatura e campo magnetico.
    rho_e = materiale["resistivita_elettrica_ohm_m"]
    mu_r = materiale["permeabilita_relativa"]
    omega = 2.0 * math.pi * frequenza_hz
    mu = MU0 * mu_r
    return math.sqrt(2.0 * rho_e / (omega * mu))
 def numero_fourier_spessore(
    materiale: dict, spessore_m: float, sorgente: dict
 ) -> float:
    # Verifica quantitativa dell'ipotesi di parete termicamente sottile:
    #
    #     Fo = alpha * (sigma / |v|) / spessore^2
    #
    # confronta il tempo di diffusione attraverso lo spessore con il tempo di
    # transito della sorgente. Per Fo >> 1 lo spessore si equilibra molto
    # prima che la sorgente sia passata e una shell a una temperatura per
    # nodo è adeguata; scendendo verso Fo ~ 1 servirebbe una formulazione
    # multistrato.
    alpha = materiale["conducibilita_termica_W_mK"] / (
        materiale["densita_kg_m3"] * materiale["calore_specifico_J_kgK"]
    )
    velocita = abs(sorgente["velocita_m_s"])
    if velocita <= 0.0:
        return float("inf")
    tempo_transito = sorgente["sigma_punto_m"] / velocita
    return alpha * tempo_transito / (spessore_m * spessore_m)
 def prepara_stato_termico(
    fascetta: dict, aria: dict, sorgente: dict, sensore: dict
 ) -> dict:
    # Costruisce mesh, matrici globali e solutore fattorizzato per un run:
    # tutto ciò che resta costante durante l'integrazione temporale.
    materiale = MATERIALI[fascetta["materiale"]]
    lunghezza_m = fascetta["lunghezza_mm"] / 1000.0
    spessore_m = fascetta["spessore_mm"] / 1000.0
    raggio_medio_m = fascetta["diametro_medio_mm"] / 2000.0
    mesh = costruisci_mesh_cilindrica(
        lunghezza_m=lunghezza_m,
        raggio_medio_m=raggio_medio_m,
        n_elementi_x=fascetta["n_elementi_x"],
        n_elementi_theta=fascetta["n_elementi_theta"],
    )
    C = assemblaggio.assembla_capacita(mesh, materiale, spessore_m)
    K_cond = assemblaggio.assembla_conduzione(mesh, materiale, spessore_m)
    K_conv, f_ambiente = assemblaggio.assembla_convezione(mesh, aria, spessore_m)
    dt_s = SIMULAZIONE["dt_interno_s"]
    solutore = sol.costruisci_solutore(C, K_cond + K_conv, dt_s)
    assemblatore = assemblaggio.prepara_assemblatore_sorgente(mesh)
    # Il fattore circonferenziale del flusso non dipende dal tempo: le
    # sorgenti traslano solo lungo x.
    fattore_arco = src.fattore_circonferenziale(
        sorgente, raggio_medio_m, assemblatore["arco_gauss_m"]
    )
    attenuazione_al_sensore = float(
        src.fattore_circonferenziale(sorgente, raggio_medio_m, np.zeros(1))[0]
    )
    if sorgente["skin_depth_fissa_m"] is None:
        skin_depth_m = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
    else:
        skin_depth_m = float(sorgente["skin_depth_fissa_m"])
    x_sensore_m = sensore["x_mm"] / 1000.0
    return {
        "mesh": mesh,
        "materiale": materiale,
        "spessore_m": spessore_m,
        "raggio_medio_m": raggio_medio_m,
        "dt_s": dt_s,
        "solutore": solutore,
        "assemblatore": assemblatore,
        "f_ambiente": f_ambiente,
        "fattore_arco": fattore_arco,
        "attenuazione_al_sensore": attenuazione_al_sensore,
        "x_sensore_m": x_sensore_m,
        "interpolatore_sensore": sens.prepara_interpolatore(mesh, x_sensore_m, 0.0),
        "sorgente": sorgente,
        "T_ambiente_C": aria["temperatura_ambiente_C"],
        "skin_depth_m": skin_depth_m,
        "numero_fourier_spessore": numero_fourier_spessore(
            materiale, spessore_m, sorgente
        ),
    }
 def campo_iniziale(stato: dict) -> np.ndarray:
    return np.full(stato["mesh"]["n_nodi"], stato["T_ambiente_C"], dtype=float)
 def passo_termico(stato: dict, T: np.ndarray, t_s: float) -> tuple[np.ndarray, dict]:
    # Avanza il campo di temperatura di un passo dt e restituisce anche lo
    # stato della sorgente in quell'istante.
    sorgente = stato["sorgente"]
    x_riferimento_m, flusso_x = src.profilo_flusso_x_W_m2(
        sorgente, stato["x_sensore_m"], t_s, stato["assemblatore"]["x_gauss_m"]
    )
    f_sorgente = assemblaggio.assembla_sorgente(
        stato["assemblatore"], flusso_x, stato["fattore_arco"]
    )
    T_next = sol.passo_implicito(
        stato["solutore"], T, stato["f_ambiente"], f_sorgente
    )
    _, flusso_al_sensore = src.profilo_flusso_x_W_m2(
        sorgente, stato["x_sensore_m"], t_s, np.array([stato["x_sensore_m"]])
    )
    return T_next, {
        "x_riferimento_m": x_riferimento_m,
        "flusso_al_sensore_W_m2": float(
            flusso_al_sensore[0] * stato["attenuazione_al_sensore"]
        ),
    }
 def configurazione_randomizzata(indice_run: int, rng: random.Random) -> dict:
    fascetta = deepcopy(FASCETTA)
    aria = deepcopy(ARIA)
    sorgente = deepcopy(SORGENTE)
    sensore = deepcopy(SENSORE)
    if RANDOMIZZAZIONE.get("abilitata", False):
        def perturba_rel(valore: float, std_rel: float, fattore_min: float = 0.1) -> float:
            fattore = rng.gauss(1.0, std_rel)
            fattore = max(fattore_min, fattore)
            return valore * fattore
        sorgente["velocita_m_s"] = perturba_rel(
            sorgente["velocita_m_s"],
            RANDOMIZZAZIONE["velocita_std_rel"],
        )
        sorgente["flusso_termico_picco_W_m2"] = perturba_rel(
            sorgente["flusso_termico_picco_W_m2"],
            RANDOMIZZAZIONE["flusso_picco_std_rel"],
        )
        sorgente["sigma_punto_m"] = perturba_rel(
            sorgente["sigma_punto_m"],
            RANDOMIZZAZIONE["sigma_punto_std_rel"],
        )
        sorgente["offset_y_percorso_m"] = rng.uniform(
            -RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
            RANDOMIZZAZIONE["offset_y_max_assoluto_m"],
        )
        aria["temperatura_ambiente_C"] += rng.gauss(
            0.0,
            RANDOMIZZAZIONE["temperatura_ambiente_std_C"],
        )
        sensore["rumore_std_C"] = perturba_rel(
            sensore["rumore_std_C"],
            RANDOMIZZAZIONE["rumore_sensore_std_rel"],
            fattore_min=0.0,
        )
    return {
        "id_run": f"run_{indice_run:04d}",
        "fascetta": fascetta,
        "aria": aria,
        "sorgente": sorgente,
        "sensore": sensore,
    }
 def simula_singolo(cfg_run: dict, output_csv: Path, rng: random.Random) -> dict:
    fascetta = cfg_run["fascetta"]
    aria = cfg_run["aria"]
    sorgente = cfg_run["sorgente"]
    sensore = cfg_run["sensore"]
    nome_materiale = fascetta["materiale"]
    stato = prepara_stato_termico(fascetta, aria, sorgente, sensore)
    mesh = stato["mesh"]
    interpolatore = stato["interpolatore_sensore"]
    skin_depth = stato["skin_depth_m"]
    dt = stato["dt_s"]
    durata = SIMULAZIONE["durata_s"]
    periodo_campionamento = 1.0 / SIMULAZIONE["frequenza_campionamento_hz"]
    T = campo_iniziale(stato)
    T_letta = sens.temperatura_osservata_C(interpolatore, T)
    prossimo_campione_t = 0.0
    T_vera_max = T_letta
    T_misurata_max = T_letta
    T_fascetta_max = float(T.max())
    output_csv.parent.mkdir(parents=True, exist_ok=True)
    with output_csv.open("w", newline="") as f:
        writer = csv.writer(f)
        writer.writerow([
            "id_run",
            "tempo_s",
            "x_sorgente_m",
            "offset_y_sorgente_m",
            "flusso_termico_sorgente_W_m2",
            "skin_depth_m",
            "T_vera_lato_sensore_C",
            "T_misurata_sensore_C",
            "T_max_fascetta_C",
            "T_media_fascetta_C",
            "T_ambiente_C",
            "velocita_m_s",
            "sigma_punto_m",
            "flusso_picco_W_m2",
            "materiale",
        ])
        t = 0.0
        while t <= durata + 1e-12:
            T, info = passo_termico(stato, T, t)
            # Temperatura vera della superficie interna nel punto osservato
            # dal sensore infrarosso.
            T_vera = sens.temperatura_osservata_C(interpolatore, T)
            T_letta = sens.aggiorna_inerzia_C(
                T_letta, T_vera, dt, sensore["costante_tempo_s"]
            )
            # I massimi delle temperature vere sono aggiornati a ogni passo
            # interno, non solo agli istanti di campionamento: servono come
            # ground truth e il picco è stretto rispetto al periodo di
            # campionamento CSV, che da solo lo taglierebbe di oltre un grado.
            T_max = float(T.max())
            T_vera_max = max(T_vera_max, T_vera)
            T_fascetta_max = max(T_fascetta_max, T_max)
            # Campionamento CSV.
            if t + 1e-12 >= prossimo_campione_t:
                misurata = T_letta + rng.gauss(0.0, sensore["rumore_std_C"])
                misurata = sens.quantizza(misurata, sensore["quantizzazione_C"])
                T_misurata_max = max(T_misurata_max, misurata)
                writer.writerow([
                    cfg_run["id_run"],
                    f"{t:.6f}",
                    f"{info['x_riferimento_m']:.9f}",
                    f"{sorgente['offset_y_percorso_m']:.9f}",
                    f"{info['flusso_al_sensore_W_m2']:.6f}",
                    f"{skin_depth:.9e}",
                    f"{T_vera:.6f}",
                    f"{misurata:.6f}",
                    f"{T_max:.6f}",
                    f"{float(T.mean()):.6f}",
                    f"{aria['temperatura_ambiente_C']:.6f}",
                    f"{sorgente['velocita_m_s']:.9f}",
                    f"{sorgente['sigma_punto_m']:.9f}",
                    f"{sorgente['flusso_termico_picco_W_m2']:.6f}",
                    nome_materiale,
                ])
                prossimo_campione_t += periodo_campionamento
            t += dt
    return {
        "id_run": cfg_run["id_run"],
        "file_csv": str(output_csv.name),
        "materiale": nome_materiale,
        "diametro_medio_m": fascetta["diametro_medio_mm"] / 1000.0,
        "lunghezza_m": fascetta["lunghezza_mm"] / 1000.0,
        "spessore_m": fascetta["spessore_mm"] / 1000.0,
        "n_elementi_x": mesh["n_elementi_x"],
        "n_elementi_theta": mesh["n_elementi_theta"],
        "n_nodi": mesh["n_nodi"],
        "lato_elemento_x_m": mesh["lato_x_m"],
        "lato_elemento_arco_m": mesh["lato_arco_m"],
        "durata_s": durata,
        "frequenza_campionamento_hz": SIMULAZIONE["frequenza_campionamento_hz"],
        "dt_interno_s": dt,
        "temperatura_ambiente_C": aria["temperatura_ambiente_C"],
        "h_esterno_W_m2K": aria["h_esterno_W_m2K"],
        "h_interno_W_m2K": aria["h_interno_W_m2K"],
        "h_bordi_W_m2K": aria["h_bordi_W_m2K"],
        "x_inizio_m": sorgente["x_inizio_m"],
        "x_fine_m": sorgente["x_fine_m"],
        "x_sensore_m": stato["x_sensore_m"],
        "distanza_sensore_parete_m": sensore["distanza_parete_mm"] / 1000.0,
        "offset_y_percorso_m": sorgente["offset_y_percorso_m"],
        "theta_percorso_rad": src.theta_percorso_rad(
            sorgente, stato["raggio_medio_m"]
        ),
        "velocita_m_s": sorgente["velocita_m_s"],
        "numero_sorgenti": sorgente.get("numero_sorgenti", 1),
        "distanza_sorgenti_m": sorgente.get("distanza_sorgenti_m", 0.0),
        "sigma_punto_m": sorgente["sigma_punto_m"],
        "flusso_termico_picco_W_m2": sorgente["flusso_termico_picco_W_m2"],
        "efficienza_riscaldamento": sorgente["efficienza_riscaldamento"],
        "frequenza_hz": sorgente["frequenza_hz"],
        "skin_depth_m": skin_depth,
        "numero_fourier_spessore": stato["numero_fourier_spessore"],
        "costante_tempo_sensore_s": sensore["costante_tempo_s"],
        "rumore_std_sensore_C": sensore["rumore_std_C"],
        "quantizzazione_sensore_C": sensore["quantizzazione_C"],
        "T_vera_max_lato_sensore_C": T_vera_max,
        "T_misurata_max_sensore_C": T_misurata_max,
        "T_max_fascetta_C": T_fascetta_max,
    }
 def _esegui_run(indice_e_seme: tuple[int, int]) -> dict:
    # Ogni run riceve un seme indipendente derivato dal seed globale, così
    # l'esecuzione in parallelo resta riproducibile indipendentemente
    # dall'ordine in cui i processi la completano.
    indice, seme = indice_e_seme
    rng = random.Random(seme)
    cfg_run = configurazione_randomizzata(indice, rng)
    cartella_output = Path(SIMULAZIONE["cartella_output"])
    percorso_csv = cartella_output / f"{cfg_run['id_run']}.csv"
    return simula_singolo(cfg_run, percorso_csv, rng)
 def main() -> None:
    cartella_output = Path(SIMULAZIONE["cartella_output"])
    if cartella_output.exists():
        shutil.rmtree(cartella_output)
    cartella_output.mkdir(parents=True, exist_ok=True)
    rng_semi = random.Random(SIMULAZIONE["seed"])
    num_run = SIMULAZIONE["num_run"]
    semi = [rng_semi.randrange(2**63) for _ in range(num_run)]
    num_processi = SIMULAZIONE["num_processi"] or os.cpu_count() or 1
    with ProcessPoolExecutor(max_workers=num_processi) as executor:
        righe_metadata = list(
            executor.map(_esegui_run, enumerate(semi, start=1))
        )
    percorso_metadata = cartella_output / "metadata.csv"
    with percorso_metadata.open("w", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=list(righe_metadata[0].keys()))
        writer.writeheader()
        writer.writerows(righe_metadata)
    print(f"Generati {len(righe_metadata)} run in: {cartella_output.resolve()}")
    print(f"Metadata: {percorso_metadata.resolve()}")
 if __name__ == "__main__":
    main()
@@ -0,0 +1,39 @@
 # Integrazione temporale implicita del sistema semidiscreto della shell.
 #
 # Dopo l'assemblaggio il problema è
 #
 #     C * dT/dt + K * T = f(t)        con K = K_cond + K_conv
 #
 # Con Eulero implicito:
 #
 #     (C/dt + K) * T^(n+1) = (C/dt) * T^n + f^(n+1)
 #
 # La matrice A = C/dt + K dipende solo da mesh, materiale, passo temporale e
 # condizioni convettive: nessuno di questi cambia durante un run, quindi si
 # fattorizza LU una volta sola e a ogni passo si risolvono solo i sistemi
 # triangolari. Se in futuro si introducessero proprietà dipendenti dalla
 # temperatura, A andrebbe ricostruita e rifattorizzata a ogni passo.
 import scipy.sparse as sp
 from scipy.sparse.linalg import splu
 def costruisci_solutore(C, K, dt_s: float) -> dict:
    # A è simmetrica e definita positiva: capacità, conduzione e convezione
    # lo sono tutte. L'ordinamento MMD sul grafo di A + A^T sfrutta la
    # simmetria e riduce il riempimento della fattorizzazione rispetto al
    # COLAMD di default, pensato per matrici non simmetriche.
    A = (C / dt_s + K).tocsc()
    return {
        "lu": splu(A, permc_spec="MMD_AT_PLUS_A"),
        "C_su_dt": (C / dt_s).tocsr(),
        "dt_s": dt_s,
    }
 def passo_implicito(solutore: dict, T, f_costante, f_sorgente):
    # Avanza il campo di un passo dt. f_costante raccoglie i termini noti
    # indipendenti dal tempo (convezione verso l'ambiente su facce e bordi),
    # f_sorgente il carico della sorgente all'istante di arrivo.
    rhs = solutore["C_su_dt"] @ T + f_costante + f_sorgente
    return solutore["lu"].solve(rhs)
@@ -0,0 +1,125 @@
 # Sorgenti a induzione in movimento sulla superficie esterna della fascetta.
 #
 # Il gruppo è formato da più sorgenti equidistanti che traslano insieme lungo
 # x con la stessa velocità. Tutte percorrono la stessa generatrice, cioè la
 # stessa posizione circonferenziale, definita dall'offset rispetto al punto
 # osservato dal sensore.
 #
 # Il flusso termico incidente su un punto della superficie è
 #
 #     q''(x, arco, t) = q_picco * somma_i exp(-((x - x_i(t))^2) / (2 sigma^2))
 #                                * exp(-(delta_arco^2) / (2 sigma^2))
 #
 # dove delta_arco è la distanza circonferenziale minima tra il punto e il
 # percorso delle sorgenti, misurata lungo la circonferenza e riavvolta su
 # 2*pi*R: passare da un lato o dall'altro della fascetta è equivalente.
 #
 # Il flusso è quindi il prodotto di un fattore che dipende solo da x e varia
 # nel tempo e di un fattore che dipende solo dall'arco ed è costante: questa
 # separazione è sfruttata dall'assemblaggio del vettore sorgente.
 import math
 import numpy as np
 def _spread_sorgenti_m(sorgente: dict) -> float:
    # Distanza lungo x tra la prima e l'ultima sorgente del gruppo.
    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
    return (numero_sorgenti - 1) * distanza
 def x_riferimento_iniziale_m(sorgente: dict, x_sensore_m: float) -> float:
    # Posizione a t=0 della sorgente di indice 0. Le sorgenti i sono a
    # x_i = riferimento + i * distanza, quindi per v >= 0 l'indice 0 è la
    # più arretrata nel verso di marcia, per v < 0 è la più avanzata.
    # x_inizio_m è la distanza dal sensore della sorgente più avanzata
    # (quella che lo raggiunge per prima).
    spread = _spread_sorgenti_m(sorgente)
    x_inizio = sorgente["x_inizio_m"]
    if sorgente["velocita_m_s"] >= 0:
        return (x_sensore_m - x_inizio) - spread
    return x_sensore_m + x_inizio
 def x_riferimento_finale_m(sorgente: dict, x_sensore_m: float) -> float:
    # Posizione di fine corsa della sorgente di indice 0. x_fine_m è la
    # distanza dal sensore della sorgente più arretrata nel verso di marcia
    # (quella che lo supera per ultima).
    spread = _spread_sorgenti_m(sorgente)
    x_fine = sorgente["x_fine_m"]
    if sorgente["velocita_m_s"] >= 0:
        return x_sensore_m + x_fine
    return (x_sensore_m - x_fine) - spread
 def _intervallo_attivo(inizio: float, fine: float, v: float, x_m: float) -> bool:
    if v >= 0:
        return inizio <= x_m <= fine
    return fine <= x_m <= inizio
 def theta_percorso_rad(sorgente: dict, raggio_medio_m: float) -> float:
    # Posizione circonferenziale del percorso delle sorgenti, ricavata
    # dall'offset in arco rispetto al sensore, che sta per convenzione a
    # theta = 0.
    return sorgente["offset_y_percorso_m"] / raggio_medio_m
 def fattore_circonferenziale(
    sorgente: dict, raggio_medio_m: float, arco_m: np.ndarray
 ) -> np.ndarray:
    # Attenuazione gaussiana del flusso in funzione della distanza
    # circonferenziale dal percorso delle sorgenti. La distanza è riavvolta
    # sulla circonferenza, così la sorgente scalda anche i punti raggiunti
    # "dall'altra parte" della fascetta.
    circonferenza = 2.0 * math.pi * raggio_medio_m
    delta = arco_m - sorgente["offset_y_percorso_m"]
    delta = (delta + circonferenza / 2.0) % circonferenza - circonferenza / 2.0
    sigma = sorgente["sigma_punto_m"]
    return np.exp(-0.5 * (delta * delta) / (sigma * sigma))
 def profilo_flusso_x_W_m2(
    sorgente: dict, x_sensore_m: float, t_s: float, x_m: np.ndarray
 ) -> tuple[float, np.ndarray]:
    # Restituisce la posizione della sorgente di riferimento e il fattore
    # assiale del flusso efficace [W/m²], somma dei contributi di tutte le
    # sorgenti attive valutata nei punti x richiesti.
    x_rif_iniziale = x_riferimento_iniziale_m(sorgente, x_sensore_m)
    x_rif_finale = x_riferimento_finale_m(sorgente, x_sensore_m)
    x_riferimento = x_rif_iniziale + sorgente["velocita_m_s"] * t_s
    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
    v = sorgente["velocita_m_s"]
    zero_dopo_fine = sorgente.get("zero_dopo_fine", True)
    sigma = sorgente["sigma_punto_m"]
    q_picco = (
        sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"]
    )
    q_x = np.zeros_like(x_m)
    for i in range(numero_sorgenti):
        x_i = x_riferimento + i * distanza
        if zero_dopo_fine:
            inizio_i = x_rif_iniziale + i * distanza
            fine_i = x_rif_finale + i * distanza
            if not _intervallo_attivo(inizio_i, fine_i, v, x_i):
                continue
        dx = x_m - x_i
        q_x += q_picco * np.exp(-0.5 * (dx * dx) / (sigma * sigma))
    return x_riferimento, q_x
 def posizioni_sorgenti_m(sorgente: dict, x_riferimento_m: float) -> np.ndarray:
    # Posizioni lungo x di tutte le sorgenti del gruppo a partire da quella
    # di riferimento.
    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza = sorgente.get("distanza_sorgenti_m", 0.0)
    return x_riferimento_m + np.arange(numero_sorgenti) * distanza
@@ -0,0 +1,252 @@
 # Validazione del solutore shell su casi con soluzione nota.
 #
 # Non è una suite di test automatica: stampa i risultati e i confronti perché
 # molti sono verifiche di convergenza, dove conta l'andamento dell'errore più
 # di una soglia di pass/fail. Va eseguita a mano dopo ogni modifica alla
 # fisica o all'assemblaggio:
 #
 #     python valida_shell.py
 import math
 import numpy as np
 import assemblaggio as asm
 import solutore as sol
 import sorgente as src
 from materials import MATERIALI
 from mesh import campo_su_griglia, costruisci_mesh_cilindrica
 MATERIALE = MATERIALI["banda_stagnata"]
 K_TERM = MATERIALE["conducibilita_termica_W_mK"]
 RHO = MATERIALE["densita_kg_m3"]
 CP = MATERIALE["calore_specifico_J_kgK"]
 ALPHA = K_TERM / (RHO * CP)
 LUNGHEZZA_M = 0.100
 RAGGIO_M = 0.035
 SPESSORE_M = 0.18e-3
 DT_S = 0.01
 def _esito(condizione: bool) -> str:
    return "OK" if condizione else "FALLITO"
 def _campo_nullo(mesh: dict) -> np.ndarray:
    return np.zeros(mesh["n_nodi"])
 def verifica_matrici(mesh: dict) -> None:
    print("=== 1. proprietà delle matrici assemblate ===")
    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
    attesa = RHO * CP * SPESSORE_M * mesh["area_totale_m2"]
    errore = abs(C.sum() / attesa - 1.0)
    print(f"  capacità totale {C.sum():.6f} vs rho*cp*s*A {attesa:.6f}"
          f"  err={errore:.1e} {_esito(errore < 1e-12)}")
    # A temperatura costante il flusso deve essere nullo, quindi ogni riga
    # della matrice di conduzione somma a zero.
    somma_righe = np.abs(np.asarray(K.sum(axis=1))).max()
    print(f"  max |somma righe K_cond| = {somma_righe:.1e}"
          f" {_esito(somma_righe < 1e-9 * K_TERM)}")
    asimmetria = abs(K - K.T).max()
    print(f"  simmetria di K_cond: scarto {asimmetria:.1e}"
          f" {_esito(asimmetria < 1e-12 * abs(K).max())}")
    autovalori = np.linalg.eigvalsh(K.toarray()[:400, :400])
    print(f"  blocco campione semidefinito positivo:"
          f" autovalore minimo {autovalori.min():.1e}"
          f" {_esito(autovalori.min() > -1e-9)}")
 def verifica_raffreddamento(mesh: dict) -> None:
    # Campo uniforme senza sorgente: deve restare uniforme e decadere verso la
    # temperatura ambiente con la costante di tempo rho*cp*s/(h_est+h_int).
    # I bordi assiali sono esclusi, altrimenti il decadimento non è uniforme.
    print("\n=== 2. raffreddamento convettivo uniforme ===")
    aria = {
        "temperatura_ambiente_C": 20.0,
        "h_esterno_W_m2K": 12.0,
        "h_interno_W_m2K": 8.0,
        "h_bordi_W_m2K": 0.0,
    }
    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
    K_conv, f_ambiente = asm.assembla_convezione(mesh, aria, SPESSORE_M)
    solutore = sol.costruisci_solutore(C, K + K_conv, DT_S)
    sovratemperatura = 200.0
    T = np.full(mesh["n_nodi"], aria["temperatura_ambiente_C"] + sovratemperatura)
    n_passi = 3000
    for _ in range(n_passi):
        T = sol.passo_implicito(solutore, T, f_ambiente, _campo_nullo(mesh))
    t = n_passi * DT_S
    tau = RHO * CP * SPESSORE_M / (aria["h_esterno_W_m2K"] + aria["h_interno_W_m2K"])
    attesa = aria["temperatura_ambiente_C"] + sovratemperatura * math.exp(-t / tau)
    errore = abs(T.mean() - attesa)
    # Eulero implicito è del primo ordine: l'errore atteso è dell'ordine di
    # dt/(2*tau) sulla sovratemperatura residua.
    tolleranza = 3.0 * (attesa - aria["temperatura_ambiente_C"]) * DT_S / (2 * tau)
    print(f"  tau = {tau:.3f} s, t = {t:.1f} s")
    print(f"  numerica {T.mean():.4f} °C vs analitica {attesa:.4f} °C"
          f"  err={errore:.4f} °C {_esito(errore < tolleranza)}")
    dispersione = T.max() - T.min()
    print(f"  campo ancora uniforme: max-min = {dispersione:.1e} °C"
          f" {_esito(dispersione < 1e-9)}")
 def _sorgente_prova() -> dict:
    return {
        "x_inizio_m": 0.5,
        "x_fine_m": 0.5,
        "offset_y_percorso_m": 0.0,
        "velocita_m_s": -1.0,
        "numero_sorgenti": 1,
        "distanza_sorgenti_m": 0.0,
        "sigma_punto_m": 0.012,
        "flusso_termico_picco_W_m2": 1e6,
        "efficienza_riscaldamento": 1.0,
        "zero_dopo_fine": True,
    }
 def verifica_sorgente(mesh: dict) -> None:
    print("\n=== 3. vettore della sorgente e conservazione dell'energia ===")
    sorgente = _sorgente_prova()
    assemblatore = asm.prepara_assemblatore_sorgente(mesh)
    fattore_arco = src.fattore_circonferenziale(
        sorgente, RAGGIO_M, assemblatore["arco_gauss_m"]
    )
    # Con la gaussiana interamente dentro la fascetta la potenza assemblata
    # deve valere l'integrale analitico 2*pi*sigma^2*q_max; con la sorgente
    # centrata su un bordo assiale deve valerne esattamente la metà.
    potenza_analitica = (
        sorgente["flusso_termico_picco_W_m2"]
        * 2.0 * math.pi * sorgente["sigma_punto_m"] ** 2
    )
    # x_riferimento(t) = x_sensore + x_inizio + v*t, con x_sensore = 0.05.
    for t, descrizione, atteso in [
        (0.50, "centrata a metà fascetta", 1.0),
        (0.55, "centrata sul bordo x = 0", 0.5),
    ]:
        _, flusso_x = src.profilo_flusso_x_W_m2(
            sorgente, 0.05, t, assemblatore["x_gauss_m"]
        )
        f = asm.assembla_sorgente(assemblatore, flusso_x, fattore_arco)
        rapporto = f.sum() / potenza_analitica
        print(f"  {descrizione:26s}: {f.sum():9.4f} W / {potenza_analitica:.4f} W"
              f" = {rapporto:.4f} (atteso {atteso:.1f})"
              f" {_esito(abs(rapporto - atteso) < 1e-3)}")
    # Senza convezione tutta l'energia iniettata deve finire nel campo.
    aria = {
        "temperatura_ambiente_C": 20.0,
        "h_esterno_W_m2K": 0.0,
        "h_interno_W_m2K": 0.0,
        "h_bordi_W_m2K": 0.0,
    }
    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
    _, f_ambiente = asm.assembla_convezione(mesh, aria, SPESSORE_M)
    solutore = sol.costruisci_solutore(C, K, DT_S)
    T = np.full(mesh["n_nodi"], 20.0)
    energia_iniziale = float((C @ T).sum())
    energia_iniettata = 0.0
    t = 0.0
    for _ in range(150):
        _, flusso_x = src.profilo_flusso_x_W_m2(
            sorgente, 0.05, t, assemblatore["x_gauss_m"]
        )
        f = asm.assembla_sorgente(assemblatore, flusso_x, fattore_arco)
        energia_iniettata += f.sum() * DT_S
        T = sol.passo_implicito(solutore, T, f_ambiente, f)
        t += DT_S
    incremento = float((C @ T).sum()) - energia_iniziale
    errore = abs(incremento / energia_iniettata - 1.0)
    print(f"  iniettata {energia_iniettata:.4f} J, accumulata {incremento:.4f} J"
          f"  err={errore:.1e} {_esito(errore < 1e-10)}")
 def verifica_modi_circonferenziali(mesh: dict) -> None:
    # Un modo sin(n*theta) su un cilindro adiabatico decade come
    # exp(-alpha * n^2 * t / R^2): confronto diretto con l'analitico.
    print("\n=== 4. modi sinusoidali circonferenziali ===")
    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
    solutore = sol.costruisci_solutore(C, K, DT_S)
    theta = np.tile(mesh["theta_nodi_rad"], mesh["n_nodi_x"])
    for n in (1, 2, 4):
        T = 100.0 * np.sin(n * theta)
        for _ in range(200):
            T = sol.passo_implicito(
                solutore, T, _campo_nullo(mesh), _campo_nullo(mesh)
            )
        t = 200 * DT_S
        attesa = 100.0 * math.exp(-ALPHA * n**2 / RAGGIO_M**2 * t)
        print(f"  n={n}: numerica {np.abs(T).max():8.4f} vs analitica {attesa:8.4f}"
              f"  err={abs(np.abs(T).max() / attesa - 1):.3%}")
 def verifica_convergenza() -> None:
    print("\n=== 5. convergenza al raffinamento della mesh (modo n=4) ===")
    for n_theta in (30, 60, 110, 220):
        mesh = costruisci_mesh_cilindrica(LUNGHEZZA_M, RAGGIO_M, 20, n_theta)
        C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
        K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
        solutore = sol.costruisci_solutore(C, K, DT_S)
        theta = np.tile(mesh["theta_nodi_rad"], mesh["n_nodi_x"])
        T = 100.0 * np.sin(4 * theta)
        for _ in range(200):
            T = sol.passo_implicito(
                solutore, T, _campo_nullo(mesh), _campo_nullo(mesh)
            )
        attesa = 100.0 * math.exp(-ALPHA * 16 / RAGGIO_M**2 * 200 * DT_S)
        print(f"  n_elementi_theta={n_theta:4d} (arco {mesh['lato_arco_m']*1e3:5.2f} mm):"
              f" err={abs(np.abs(T).max() / attesa - 1):.4%}")
 def verifica_periodicita(mesh: dict) -> None:
    # Un impulso su theta = 0 deve diffondere in modo identico nei due versi:
    # è la prova che la connettività riavvolge davvero la circonferenza.
    print("\n=== 6. periodicità circonferenziale ===")
    C = asm.assembla_capacita(mesh, MATERIALE, SPESSORE_M)
    K = asm.assembla_conduzione(mesh, MATERIALE, SPESSORE_M)
    solutore = sol.costruisci_solutore(C, K, DT_S)
    T = np.full(mesh["n_nodi"], 20.0)
    T[50 * mesh["n_elementi_theta"]] = 500.0
    for _ in range(50):
        T = sol.passo_implicito(solutore, T, _campo_nullo(mesh), _campo_nullo(mesh))
    griglia = campo_su_griglia(mesh, T)
    verso_positivo = griglia[50, 1:6]
    verso_negativo = griglia[50, -1:-6:-1]
    scarto = np.abs(verso_positivo - verso_negativo).max()
    print(f"  verso theta+ : {np.array2string(verso_positivo, precision=4)}")
    print(f"  verso theta- : {np.array2string(verso_negativo, precision=4)}")
    print(f"  scarto massimo = {scarto:.1e} {_esito(scarto < 1e-9)}")
 def main() -> None:
    mesh = costruisci_mesh_cilindrica(LUNGHEZZA_M, RAGGIO_M, 100, 110)
    verifica_matrici(mesh)
    verifica_raffreddamento(mesh)
    verifica_sorgente(mesh)
    verifica_modi_circonferenziali(mesh)
    verifica_convergenza()
    verifica_periodicita(mesh)
 if __name__ == "__main__":
    main()