diff --git a/CLAUDE.md b/CLAUDE.md index fb590d1..ed2a9e0 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -24,6 +24,7 @@ python main.py simula # genera dataset/run_XXXX.csv + dataset/metadata.csv python main.py grafico # grafici del primo run python main.py anima # animazione 2D della sezione python main.py anima3d # animazione 3D isometrica del barattolo +python main.py fem # animazione 3D del campo FEM sulla mesh shell ``` Ogni modulo resta eseguibile anche direttamente (`python simulate.py`, `python plot_csv.py`, ...). @@ -45,10 +46,26 @@ Generatore di dataset per misurazioni termiche pseudo-realistiche di una fascett 2. `materials.py` — dizionario MATERIALI con proprietà termofisiche ed elettriche per materiale 3. `simulate.py` — motore principale: genera N run randomizzati, scrive i CSV, scrive `metadata.csv` 4. `plot_csv.py` — visualizzazione autonoma per un singolo run -5. `mesh.py` — mesh a elementi shell quadrangolari della superficie media del cilindro (base per la futura analisi FEM), visualizzata da `plot_mesh.py` +5. `mesh.py` — mesh a elementi shell quadrangolari della superficie media del cilindro, visualizzata da `plot_mesh.py` +6. `fem.py` — solutore termico transitorio a elementi finiti sulla mesh shell, animato da `plot_animazione_fem.py` **Mesh shell (`mesh.py`):** griglia strutturata `n_elementi_x × n_elementi_circonferenza` di quadrilateri a 4 nodi sulla superficie media (raggio = (diametro − spessore)/2); lo spessore non è discretizzato, è un attributo della shell. La mesh è chiusa lungo la circonferenza (nessun nodo duplicato sulla cucitura) e usa lo stesso sistema di coordinate globale di `plot_animazione_3d.py`: x = asse, y = R·sin(theta), z = R·cos(theta), con theta = 0 sul piano del sensore. +## Due modelli termici distinti + +Il progetto contiene due discretizzazioni indipendenti della stessa fisica; scambiano solo la cinematica delle sorgenti (`_x_riferimento_iniziale_m`, `_x_riferimento_finale_m`, `_intervallo_attivo` in [simulate.py](simulate.py)), che deve restare unica. + +| | volumi finiti ([simulate.py](simulate.py)) | elementi finiti ([fem.py](fem.py)) | +|---|---|---| +| dominio | sezione x-z (spessore risolto) | superficie media x-s (circonferenza risolta) | +| circonferenza y | attenuazione gaussiana + sink di aletta | direzione risolta, conduzione reale | +| spessore | `n_nodi_z` celle, deposizione esponenziale da skin depth | collassato: T uniforme nella parete | +| griglia | `FASCETTA["n_nodi_x"]`, `n_nodi_z` | `MESH["n_elementi_x"]`, `n_elementi_circonferenza` | +| passo temporale | `SIMULAZIONE["dt_interno_s"]` | `FEM["dt_s"]` | +| output | dataset CSV + animazioni 2D/3D | animazione 3D sulla mesh | + +**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 skin depth non entra nel modello shell, conta solo il flusso assorbito. 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 run (`prepara_stato_fem`), poi `passo_implicito_fem` assembla solo il termine noto. + **Pipeline fisica dentro `simula_singolo()` in [simulate.py](simulate.py):** - La skin depth è calcolata dalla resistività elettrica del materiale e dalla frequenza di induzione (`calcola_skin_depth_m`) @@ -71,7 +88,7 @@ Ogni parametro in [config.py](config.py) ha un commento che spiega solo cos'è ( ## Vincoli progettuali chiave -- Il modello è 2D nella sezione (x = lunghezza, z = spessore). La coordinata circonferenziale y non è risolta spazialmente — l'offset y del percorso delle sorgenti è collassato in un'attenuazione gaussiana del flusso, e la conduzione lungo y è un termine di scambio lineare verso la temperatura ambiente; il diametro è registrato solo come geometria del setup. +- Il modello ai volumi finiti di `simulate.py` è 2D nella sezione (x = lunghezza, z = spessore). La coordinata circonferenziale y non è risolta spazialmente — l'offset y del percorso delle sorgenti è collassato in un'attenuazione gaussiana del flusso, e la conduzione lungo y è un termine di scambio lineare verso la temperatura ambiente; il diametro è registrato solo come geometria del setup. Il modello FEM di `fem.py` fa il contrario (circonferenza risolta, spessore collassato): non modificare uno assumendo che valgano le ipotesi dell'altro. - Le posizioni di inizio/fine corsa delle sorgenti (`x_inizio_m`, `x_fine_m`) sono distanze dal punto x del sensore lungo il verso di marcia; il segno di `velocita_m_s` determina il verso (negativo = -x). - La matrice implicita è costruita e fattorizzata una volta per run (proprietà del materiale costanti, nessun coefficiente dipendente dalla temperatura). Se si aggiungono proprietà dipendenti dalla temperatura, la matrice deve essere ricostruita e rifattorizzata ad ogni passo temporale. - `simulate.py` cancella e ricrea l'intera cartella di output ad ogni esecuzione (`shutil.rmtree`). diff --git a/config.py b/config.py index 9d72336..f3b6b55 100644 --- a/config.py +++ b/config.py @@ -89,6 +89,11 @@ MESH = { "n_elementi_circonferenza": 48, } +FEM = { + # Passo di integrazione temporale del solutore a elementi finiti. + "dt_s": 0.001, +} + ARIA = { # Temperatura dell'aria ambiente. "temperatura_ambiente_C": 25.0, diff --git a/fem.py b/fem.py new file mode 100644 index 0000000..f1301bf --- /dev/null +++ b/fem.py @@ -0,0 +1,313 @@ +# Solutore termico transitorio a elementi finiti sulla mesh shell della +# fascetta (vedi mesh.py per la geometria). +# +# Differenza rispetto al modello ai volumi finiti di simulate.py: +# - lì il dominio è la sezione x-z (spessore risolto, circonferenza y +# collassata in un'attenuazione gaussiana del flusso e in un sink lineare +# di tipo aletta); +# - qui 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 skin depth non entra nel modello shell: conta solo il flusso +# totale assorbito per unità di superficie. +# +# 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 numpy as np +import scipy.sparse as sp +from scipy.sparse.linalg import splu + +from config import FEM +from materials import MATERIALI +from mesh import genera_mesh +# Helper della cinematica delle sorgenti condivisi con il modello ai volumi +# finiti: la corsa del gruppo di induttori deve essere identica nei due modelli. +from simulate import ( + _intervallo_attivo, + _x_riferimento_finale_m, + _x_riferimento_iniziale_m, +) + +# 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 _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, + 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 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"]] + 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 + + # Coordinate nodali nel piano parametrico (x, s) della superficie media. + x_nodi_m = np.repeat(mesh_dati["x_nodi_m"], n_theta) + s_nodi_m = np.tile(raggio_m * mesh_dati["theta_nodi_rad"], n_x + 1) + + # 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": x_nodi_m, + "s_nodi_m": s_nodi_m, + "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, + } + + +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); la distanza circonferenziale + è valutata sull'immagine più vicina, perché la superficie è chiusa. + """ + 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"] + + circonferenza = stato["circonferenza_m"] + ds = stato["s_nodi_m"] - sorgente["offset_y_percorso_m"] + ds -= circonferenza * np.round(ds / circonferenza) + peso_circonferenziale = np.exp(-0.5 * (ds * ds) / (sigma * sigma)) + + q_nodale = np.zeros(stato["n_nodi"]) + 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 + q_nodale += ( + q_picco + * peso_circonferenziale + * np.exp(-0.5 * (dx * dx) / (sigma * sigma)) + ) + + return x_riferimento, q_nodale + + +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( + cfg_run: dict, + durata_s: float, + dt_frame_s: float, + mesh_dati: dict | None = None, +) -> dict: + """Integra il campo FEM fino a durata_s, salvando un campo ogni dt_frame_s.""" + fascetta = cfg_run["fascetta"] + aria = cfg_run["aria"] + sorgente = cfg_run["sorgente"] + sensore = cfg_run["sensore"] + + stato = prepara_stato_fem(fascetta, aria, 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 = [], [] + + 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: + tempi.append(t) + campi.append(T.copy()) + x_riferimenti.append(x_rif) + T_vere.append(float(T[indice_sensore])) + T_lette.append(T_sensore) + 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), + "sorgente": sorgente, + "T_ambiente_C": stato["T_ambiente_C"], + } diff --git a/main.py b/main.py index 2ca39c6..e8b9f2d 100644 --- a/main.py +++ b/main.py @@ -6,6 +6,7 @@ # python main.py grafico # grafici del primo run generato # python main.py anima # animazione 2D della sezione # python main.py anima3d # animazione 3D isometrica del barattolo +# python main.py fem # animazione 3D del campo FEM sulla mesh shell # # Senza argomenti stampa l'elenco delle azioni disponibili. @@ -43,6 +44,12 @@ def _azione_anima3d() -> None: plot_animazione_3d.main() +def _azione_fem() -> None: + import plot_animazione_fem + + plot_animazione_fem.main() + + # Chiave da riga di comando -> (funzione, descrizione mostrata nell'help). AZIONI = { "mesh": (_azione_mesh, "Disegna la mesh a elementi shell della fascetta"), @@ -50,6 +57,7 @@ AZIONI = { "grafico": (_azione_grafico, "Grafici temperatura e flusso del primo run"), "anima": (_azione_anima, "Animazione 2D del campo di temperatura nella sezione"), "anima3d": (_azione_anima3d, "Animazione 3D isometrica del barattolo"), + "fem": (_azione_fem, "Animazione 3D del campo FEM sulla mesh a elementi shell"), } diff --git a/plot_animazione_fem.py b/plot_animazione_fem.py new file mode 100644 index 0000000..678e021 --- /dev/null +++ b/plot_animazione_fem.py @@ -0,0 +1,170 @@ +# Animazione 3D isometrica del barattolo con il campo di temperatura calcolato +# a elementi finiti sulla mesh shell (fem.py). +# +# È l'equivalente di plot_animazione_3d.py ma senza ricostruzioni: lì il campo +# circonferenziale era estruso riusando la gaussiana della sorgente, qui la +# coordinata circonferenziale è una direzione risolta del modello FEM, quindi +# la temperatura dipinta su ogni elemento è quella effettivamente calcolata. +# +# 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. + +import random +from pathlib import Path + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from matplotlib import cm +from matplotlib.animation import FuncAnimation, PillowWriter +from mpl_toolkits.mplot3d.art3d import Poly3DCollection + +from config import SIMULAZIONE +from fem import simula_campo_fem +from mesh import genera_mesh, riepilogo_mesh +from plot_animazione import T_FINE_ANIMAZIONE_S, T_INIZIO_ANIMAZIONE_S +from simulate import configurazione_randomizzata + +# 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)) + + rng = random.Random(SIMULAZIONE["seed"]) + cfg_run = configurazione_randomizzata(1, rng) + + print("Integrazione FEM in corso...") + dati = simula_campo_fem( + cfg_run, + durata_s=T_FINE_ANIMAZIONE_S, + 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" + ) + + 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, 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() + 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([]) + barra = fig.colorbar(mappabile, ax=ax, shrink=0.6, pad=0.05) + barra.set_label("T [°C]") + + titolo = ax.set_title("") + + 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), + ) + + titolo.set_text( + f"FEM shell — t = {tempi[k]:.2f} s — " + f"T sensore = {dati['T_lette'][k]:.1f} °C" + ) + return collezione, marker_sorgenti, titolo + + animazione = FuncAnimation( + fig, disegna_frame, frames=indici_frame, + interval=INTERVALLO_RIPRODUZIONE_MS, blit=False, + ) + + if matplotlib.get_backend().lower() == "agg": + cartella = Path("dataset") + 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()