From 624ea068b5d8183edcc015582f4d27a9444c73f9 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 9 May 2026 18:25:56 +0300 Subject: [PATCH] chore(test): remove THMEC wedge DOF regression Delete unmaintained penta-element coupling tests. - Drop `test/dofs/test_thmec_penta_physics.jl`. --- test/dofs/test_thmec_penta_physics.jl | 779 -------------------------- 1 file changed, 779 deletions(-) delete mode 100644 test/dofs/test_thmec_penta_physics.jl diff --git a/test/dofs/test_thmec_penta_physics.jl b/test/dofs/test_thmec_penta_physics.jl deleted file mode 100644 index 00c31e0..0000000 --- a/test/dofs/test_thmec_penta_physics.jl +++ /dev/null @@ -1,779 +0,0 @@ -""" -πŸš€ THE ULTIMATE: Thermo-Hydro-Mechanical-Electric-Chemical (THM-EC) - -This implements **FIVE-FIELD** coupled physics - the most complex multi-physics -system we've attempted! - -Field Variables: -- T: Temperature (Float64) at VERTICES - continuous HΒΉ field -- u: Displacement (Vec{3}) at VERTICES - continuous HΒΉ vector field -- p: Pore pressure (Float64) at CELLS - discontinuous LΒ² field -- Ο†: Electric potential (Float64) at EDGES - H(curl) field -- c: Chemical concentration (Float64) at VERTICES - continuous HΒΉ field - -═══════════════════════════════════════════════════════════════════════ -COMPLETE PHYSICS FORMULATION - FULLY COUPLED THM-EC SYSTEM -═══════════════════════════════════════════════════════════════════════ - -1️⃣ THERMAL (Heat Equation with ALL Couplings): - ρcβ‚š βˆ‚T/βˆ‚t - βˆ‡Β·(ΞΊβˆ‡T) = Q + Ξ±_TΒ·Tβ‚€Β·E/(1-2Ξ½) βˆ‡Β·βˆ‚u/βˆ‚t + Ξ²_TΒ·βˆ‚p/βˆ‚t + SΒ·βˆ‡Β·J + Q_chem(c,T) - - NEW: Q_chem = H_rxn Β· R(c,T) - Heat source from chemical reactions - - Coupling parameters: - - Ξ±_T: thermal expansion coefficient [1/K] - - Ξ²_T: thermal pressurization coefficient [K/Pa] - - S: Seebeck coefficient [V/K] - - H_rxn: heat of reaction [J/mol] - -2️⃣ MECHANICAL (Linear Elasticity with Multi-Physics): - ρ βˆ‚Β²u/βˆ‚tΒ² - βˆ‡Β·Οƒ = f - - Οƒ = C : Ξ΅(u) - Ξ±_TΒ·(T-Tβ‚€)Β·I - Ξ±_pΒ·pΒ·I - e^TΒ·E - Ξ±_cΒ·cΒ·I - - NEW: Ξ±_cΒ·cΒ·I - Chemomechanical coupling (swelling/shrinkage from concentration) - - Examples: - - Corrosion-induced expansion - - Polymer swelling in solvents - - Concrete alkali-silica reaction - -3️⃣ HYDRAULIC (Darcy Flow with Multi-Physics): - S_s βˆ‚p/βˆ‚t + Ξ±_p βˆ‚(βˆ‡Β·u)/βˆ‚t + Ξ²_T βˆ‚T/βˆ‚t - βˆ‡Β·(k/ΞΌ_f βˆ‡p) = q - ΞΆΒ·βˆ‡Β·J + q_chem(c) - - NEW: q_chem = Ξ½_f Β· R(c,T) - Fluid source from chemical reactions - - Examples: - - Dissolution creating pore space - - Precipitation clogging pores - - Gas generation from reactions - -4️⃣ ELECTRIC (Charge Conservation): - βˆ‡Β·D = ρ_e - D = Ρ·E + e:Ξ΅(u) - pΒ·βˆ‡ΞΆ - zΒ·FΒ·cΒ·βˆ‡ΞΌ_m - J = Οƒ_eΒ·E + SΒ·(-ΞΊβˆ‡T) + zΒ·FΒ·D_mΒ·βˆ‡c - - NEW: Migration current: zΒ·FΒ·D_mΒ·βˆ‡c - Charged species move in electric field - NEW: Electro-diffusion potential: zΒ·FΒ·cΒ·βˆ‡ΞΌ_m - - Examples: - - Ion transport in batteries - - Corrosion currents - - Electrochemical sensors - -5️⃣ CHEMICAL (Species Transport - NEW FIELD!): - βˆ‚c/βˆ‚t + βˆ‡Β·J_c = R(c,T) - - J_c = -D_eff(p,T)Β·βˆ‡c + uΜ‡Β·c + ΞΌ_mΒ·cΒ·E - D_TΒ·cΒ·βˆ‡T - - Transport mechanisms: - - Diffusion: D_effΒ·βˆ‡c (Fick's law) - - Advection: uΜ‡Β·c (carried by fluid) - - Migration: ΞΌ_mΒ·cΒ·E (in electric field) - - Thermal diffusion: D_TΒ·cΒ·βˆ‡T (Soret effect) - - Reaction: R(c,T) = kβ‚€Β·exp(-E_a/RT)Β·c^n - - Coupling dependencies: - - D_eff(p,T) = Dβ‚€Β·exp(Ξ±_DΒ·p + Ξ²_DΒ·T) - Pressure/temperature-dependent diffusivity - - uΜ‡ from mechanical deformation - - E from electric field - - T affects reaction rate exponentially - - Applications: - - CO2 sequestration (dissolution in brine) - - Nuclear waste (radionuclide transport) - - Geothermal (mineral dissolution/precipitation) - - Concrete (chloride ingress, ASR) - - Batteries (Li-ion transport) - - Corrosion (electrochemical reactions) - -═══════════════════════════════════════════════════════════════════════ -COUPLING MATRIX (20 OFF-DIAGONAL BLOCKS!): -═══════════════════════════════════════════════════════════════════════ - - β”‚ T u p Ο† c - ─────┼────────────────────────────────────────────────────── - T β”‚ K_TT K_Tu K_Tp K_TΟ† K_Tc - β”‚ (Ξ±_T) (Ξ²_T) (S) (H_rxn) - ─────┼────────────────────────────────────────────────────── - u β”‚ K_uT K_uu K_up K_uΟ† K_uc - β”‚ (Ξ±_T) (Ξ±_p) (e_kij) (Ξ±_c) - ─────┼────────────────────────────────────────────────────── - p β”‚ K_pT K_pu K_pp K_pΟ† K_pc - β”‚ (Ξ²_T) (Ξ±_p) (ΞΆ) (Ξ½_f) - ─────┼────────────────────────────────────────────────────── - Ο† β”‚ K_Ο†T K_Ο†u K_Ο†p K_φφ K_Ο†c - β”‚ (S) (e_kij) (ΞΆ) (zΒ·FΒ·D_m) - ─────┼────────────────────────────────────────────────────── - c β”‚ K_cT K_cu K_cp K_cΟ† K_cc - β”‚ (D_T) (adv) (D_eff) (ΞΌ_m) - - Legend: - - Ξ±_T: thermal expansion - - Ξ±_p: Biot coefficient - - Ξ±_c: chemomechanical expansion - - Ξ²_T: thermal pressurization - - ΞΆ: electro-osmotic coefficient - - S: Seebeck coefficient - - e_kij: piezoelectric tensor - - H_rxn: heat of reaction - - Ξ½_f: stoichiometric fluid coefficient - - zΒ·FΒ·D_m: ionic migration - - D_T: thermal diffusion (Soret) - - ΞΌ_m: electrophoretic mobility - -═══════════════════════════════════════════════════════════════════════ -PHYSICAL INTERPRETATION: -═══════════════════════════════════════════════════════════════════════ - -This is the MOTHER OF ALL COUPLING systems for porous media! - -Real-world scenarios: -1. **Geothermal reservoirs**: Fluid flow (p), heat (T), rock deformation (u), - mineral dissolution (c), electrokinetic effects (Ο†) - -2. **Nuclear waste disposal**: Radionuclide transport (c) in heated (T), - saturated (p), deforming (u) clay with electrochemical (Ο†) effects - -3. **CO2 sequestration**: Gas injection (p) causes cooling (T), formation - swelling (u), dissolution (c), pH changes affecting ΞΆ-potential (Ο†) - -4. **Concrete durability**: Chloride ingress (c) in heated (T), saturated (p), - cracking (u) concrete with corrosion currents (Ο†) - -5. **Battery electrodes**: Li-ion diffusion (c) with heat generation (T), - volume expansion (u), pore pressure (p), electric field (Ο†) - -6. **Corrosion**: Oxygen diffusion (c), galvanic currents (Ο†), crevice pressure (p), - stress corrosion (u), local heating (T) - -═══════════════════════════════════════════════════════════════════════ -""" - -using Test -using JuliaFEM -using LinearAlgebra -using SparseArrays -using StaticArrays -using Tensors - -@testset "πŸš€ THM-EC: PENTA-PHYSICS (5 Fields!) on All Entity Types" begin - - println("\n" * "="^70) - println("πŸš€ THM-EC: FIVE-FIELD COMPLETE PHYSICS ON ALL ENTITY TYPES") - println("="^70) - - # Create simple 3D mesh: 2 tetrahedra sharing a face - nodes = [ - Vec{3,Float64}((0.0, 0.0, 0.0)), # Node 1 - Vec{3,Float64}((1.0, 0.0, 0.0)), # Node 2 - Vec{3,Float64}((0.5, 1.0, 0.0)), # Node 3 - Vec{3,Float64}((0.5, 0.5, 1.0)), # Node 4 - Vec{3,Float64}((1.5, 0.5, 0.5)), # Node 5 - ] - connectivity = [ - (UInt32(1), UInt32(2), UInt32(3), UInt32(4)), # Tet 1 - (UInt32(2), UInt32(3), UInt32(4), UInt32(5)), # Tet 2 - ] - mesh = Mesh{Tetrahedron{4}}(nodes, connectivity) - - println("\n3D Mesh: 2 tetrahedra, 5 nodes") - - # Create ONE element type with ALL FIVE physics fields! - println("\nCreating multi-field elements with ALL FIVE physics fields...") - - # Define field spec with FIVE fields! - S = @DOFSet{T::DOF{Temperature, Vertex}, - u::DOF{Displacement{3}, Vertex}, - p::DOF{Pressure, Cell}, - Ο†::DOF{ElectricPotential, Edge}, - c::DOF{ChemicalConcentration, Vertex}} # NEW: Chemical concentration! - - # Step 1: Initialize DOF manager - dof_mgr = DOFManager(mesh) - - # Step 2: Register fields and create elements - register_fields!(dof_mgr, S) - elements = create_elements!(dof_mgr, Element{Tetrahedron{4}, Lagrange{1}, S}) - - n_total = dof_mgr.total_dofs - - # Count DOFs by field - elem1 = first(elements) - n_T = length(elem1.dof_indices.T) - n_u = length(elem1.dof_indices.u) - n_p = length(elem1.dof_indices.p) - n_Ο† = length(elem1.dof_indices.Ο†) - n_c = length(elem1.dof_indices.c) # NEW! - - # Total system DOFs - n_T_total = count_field_dofs(dof_mgr, :T) - n_u_total = count_field_dofs(dof_mgr, :u) - n_p_total = count_field_dofs(dof_mgr, :p) - n_Ο†_total = count_field_dofs(dof_mgr, :Ο†) - n_c_total = count_field_dofs(dof_mgr, :c) # NEW! - - # Actual field offsets in system (some fields share DOFs!) - offset_T = 0 - offset_u = n_T_total - offset_p = offset_u + n_u_total - offset_Ο† = offset_p + n_p_total - offset_c = n_total - n_c_total # c is at the end! - - println(" Temperature: $n_T DOFs per element (total: $n_T_total in system)") - println(" Displacement: $n_u DOFs per element (total: $n_u_total in system)") - println(" Pressure: $n_p DOFs per element (total: $n_p_total in system)") - println(" Electric: $n_Ο† DOFs per element (total: $n_Ο†_total in system)") - println(" Chemical: $n_c DOFs per element (total: $n_c_total in system)") # NEW! - println(" TOTAL SYSTEM DOFs: $n_total") - - @test n_T == 4 - @test n_u == 12 - @test n_p == 1 - @test n_Ο† == 6 - @test n_c == 4 # NEW: Same as temperature (both at vertices) - - # ======================================================================== - # ASSEMBLING REAL PHYSICS FROM MULTI-FIELD ELEMENTS - # ======================================================================== - - println("\n" * "="^70) - println("ASSEMBLING REAL PHYSICS FROM MULTI-FIELD ELEMENTS (NO MOCKS!)") - println("="^70) - - println("\nπŸ”₯ ONE ELEMENT LOOP - ALL FIVE PHYSICS FIELDS!") - println("="^70) - - # Allocate global system (5 fields now!) - K = spzeros(n_total, n_total) - F = zeros(n_total) - - # Material properties - ΞΊ = 50.0 # Thermal conductivity [W/(mΒ·K)] - E_young = 1e9 # Young's modulus [Pa] - Ξ½ = 0.3 # Poisson's ratio - k_perm = 1e-15 # Permeability [mΒ²] - ΞΌ_f = 1e-3 # Fluid viscosity [PaΒ·s] - Ξ΅_0 = 8.854e-12 # Vacuum permittivity [F/m] - Ξ΅_r = 80.0 # Relative permittivity (water) - Οƒ_e = 1e-2 # Electrical conductivity [S/m] - - # NEW: Chemical properties - D_0 = 1e-9 # Base diffusivity [mΒ²/s] - Ξ±_D = 1e-10 # Pressure dependence [1/Pa] - Ξ²_D = 0.01 # Temperature dependence [1/K] - k_rxn = 1e-6 # Reaction rate [1/s] - H_rxn = 5e4 # Heat of reaction [J/mol] - - # Coupling coefficients - Ξ±_T = 1e-5 # Thermal expansion [1/K] - Ξ±_p = 0.7 # Biot coefficient [-] - Ξ²_T = 1e-6 # Thermal pressurization [K/Pa] - ΞΆ = 1e-10 # Electro-osmotic [mΒ²/(VΒ·s)] - S_seebeck = 1e-6 # Seebeck coefficient [V/K] - Ξ±_c = 2e-4 # Chemomechanical expansion [1/(mol/mΒ³)] # NEW! - Ξ½_f = 1e-6 # Stoichiometric fluid coefficient [mΒ³/mol] # NEW! - z_F_Dm = 1e-11 # Ionic migration [mΒ²/(VΒ·s)] # NEW! - D_T = 1e-12 # Thermal diffusion (Soret) [mΒ²/(sΒ·K)] # NEW! - ΞΌ_m = 1e-10 # Electrophoretic mobility [mΒ²/(VΒ·s)] # NEW! - - # Piezoelectric tensor (3rd order) - e_piezo = Tensor{3,3}((k,i,j) -> k==i==j ? 1e-8 : 0.0) - - # Lame parameters - Ξ» = E_young * Ξ½ / ((1 + Ξ½) * (1 - 2*Ξ½)) - ΞΌ = E_young / (2 * (1 + Ξ½)) - - # Unit vectors for volumetric coupling - e_1 = Vec{3}((1.0, 0.0, 0.0)) - e_2 = Vec{3}((0.0, 1.0, 0.0)) - e_3 = Vec{3}((0.0, 0.0, 1.0)) - - # ======================================================================== - # MODULAR COUPLING FUNCTIONS - FULL PHYSICS (NO SIMPLIFICATIONS!) - # ======================================================================== - - # All coupling functions use proper tensor operations! - - @inline function thermal_expansion_coupling(Ξ±_T, E, Ξ½, βˆ‡N_T, βˆ‡N_u, e_Ξ±, vol) - # Full: Οƒ = C:Ξ΅ - Ξ±_TΒ·(T-Tβ‚€)Β·I - coupling_strength = Ξ±_T * E / (1 - 2*Ξ½) - return coupling_strength * (βˆ‡N_T β‹… e_Ξ±) * (e_Ξ± β‹… βˆ‡N_u) * vol - end - - @inline function biot_coupling(Ξ±_p, βˆ‡N_u, e_Ξ±, vol) - # Full: Οƒ_eff = Οƒ_total + Ξ±_pΒ·pΒ·I - return Ξ±_p * (e_Ξ± β‹… βˆ‡N_u) * vol - end - - @inline function seebeck_peltier_coupling(S, βˆ‡N_T, βˆ‡N_Ο†, vol) - # Full: J = Οƒ_eΒ·E + SΒ·(-ΞΊβˆ‡T) (Seebeck/Peltier thermoelectric) - return S * (βˆ‡N_T β‹… βˆ‡N_Ο†) * vol - end - - @inline function electroosmotic_coupling(ΞΆ, βˆ‡N_p, βˆ‡N_Ο†, vol) - # Full: v_f = -k/ΞΌ_fΒ·βˆ‡p + ΞΆΒ·E (fluid flow driven by electric field) - return ΞΆ * (βˆ‡N_p β‹… βˆ‡N_Ο†) * vol - end - - @inline function compute_strain_gradient_product(e::Tensor{3,3}, - βˆ‡N_u::Vec{3}, - βˆ‡N_Ο†::Vec{3}, - i_comp::Int, - vol::Float64) - # Contract: e_kij Β· (βˆ‚N_u^i/βˆ‚x_j) Β· (βˆ‚N_Ο†/βˆ‚x_k) - # This is the FULL piezoelectric coupling integral! - result = 0.0 - for k in 1:3, j in 1:3 - result += e[k,i_comp,j] * βˆ‡N_u[j] * βˆ‡N_Ο†[k] - end - return result * vol - end - - # ======================================================================== - # NEW CHEMICAL COUPLING FUNCTIONS! - # ======================================================================== - - @inline function chemomechanical_coupling(Ξ±_c, βˆ‡N_u, e_Ξ±, vol) - # Volumetric strain from concentration change: Ξ΅_vol = Ξ±_cΒ·c - # Couples to stress: Οƒ = C:Ξ΅ - Ξ±_cΒ·cΒ·I - return Ξ±_c * (e_Ξ± β‹… βˆ‡N_u) * vol - end - - @inline function chemical_reaction_heat(H_rxn, N_T, N_c, vol) - # Heat source from chemical reaction: Q = H_rxn Β· R(c,T) - # Simplified: R(c) β‰ˆ k_rxn Β· c - return H_rxn * k_rxn * N_T * N_c * vol - end - - @inline function chemical_fluid_source(Ξ½_f, N_p, N_c, vol) - # Fluid mass source from reaction: q = Ξ½_f Β· R(c,T) - return Ξ½_f * k_rxn * N_p * N_c * vol - end - - @inline function ionic_migration_coupling(z_F_Dm, βˆ‡N_c, βˆ‡N_Ο†, vol) - # Migration current: J = zΒ·FΒ·D_mΒ·βˆ‡c (charged species in electric field) - return z_F_Dm * (βˆ‡N_c β‹… βˆ‡N_Ο†) * vol - end - - @inline function thermal_diffusion_coupling(D_T, βˆ‡N_c, βˆ‡N_T, vol) - # Soret effect: J_c = -D_TΒ·cΒ·βˆ‡T (species move toward cold/hot) - return D_T * (βˆ‡N_c β‹… βˆ‡N_T) * vol - end - - @inline function advective_coupling(N_c, βˆ‡N_u, vol) - # Advection: J_c = uΜ‡Β·c (species carried by deformation) - # Simplified: ∫ N_c Β· (βˆ‡N_u) dV - return N_c * (βˆ‡N_u[1] + βˆ‡N_u[2] + βˆ‡N_u[3]) * vol - end - - @inline function pressure_dependent_diffusion_coupling(Ξ±_D, βˆ‡N_c, βˆ‡N_p, vol) - # D_eff(p) = Dβ‚€Β·exp(Ξ±_DΒ·p) β†’ linearized contribution - return Ξ±_D * D_0 * (βˆ‡N_c β‹… βˆ‡N_p) * vol - end - - # ======================================================================== - # ELEMENT ASSEMBLY - ONE LOOP FOR ALL PHYSICS! - # ======================================================================== - - for (elem_idx, elem) in enumerate(elements) - println("\nπŸ“¦ Element $elem_idx:") - - # Get LOCAL-GLOBAL mapping for coupled assembly - n_local = local_dof_count(elem) - dof_map = local_to_global_map(elem) - - # Get LOCAL DOF ranges for each field - T_local = field_dof_range(elem, :T) - u_local = field_dof_range(elem, :u) - p_local = field_dof_range(elem, :p) - Ο†_local = field_dof_range(elem, :Ο†) - c_local = field_dof_range(elem, :c) # NEW! - - println(" Total local DOFs: $n_local") - println(" T local range: $T_local ($(length(T_local)) DOFs)") - println(" u local range: $u_local ($(length(u_local)) DOFs)") - println(" p local range: $p_local ($(length(p_local)) DOFs)") - println(" Ο† local range: $Ο†_local ($(length(Ο†_local)) DOFs)") - println(" c local range: $c_local ($(length(c_local)) DOFs)") # NEW! - - # Allocate local matrices for FIVE FIELDS - K_local = zeros(n_local, n_local) - F_local = zeros(n_local) - - # Get element connectivity - conn = mesh.connectivity[elem_idx] - - # Integration over element - quad = Gauss{4}() # Order 4 for Tet4 - ips = integration_points(quad, Tetrahedron{4}()) - - for ip in ips - ΞΎ = Vec{3}(ip.ΞΎ) - w = ip.weight - - # Basis function derivatives - dN_dΞΎ = get_basis_derivatives(Tetrahedron{4}(), Lagrange{1}(), ΞΎ) - - # Compute Jacobian - X_nodes = [nodes[i] for i in conn] - J = X_nodes[1] βŠ— dN_dΞΎ[1] - @inbounds for i in 2:4 - J += X_nodes[i] βŠ— dN_dΞΎ[i] - end - - # Physical gradients - J_inv_T = transpose(inv(J)) - dN = ntuple(i -> J_inv_T β‹… dN_dΞΎ[i], 4) - - # Basis functions (for reaction terms) - N = get_basis_functions(Tetrahedron{4}(), Lagrange{1}(), ΞΎ) - - # Jacobian determinant and volume - J_det = det(J) - vol = w * J_det - - # ---------------------------------------------------------------------- - # DIAGONAL BLOCKS (Field self-interactions) - # ---------------------------------------------------------------------- - - # K_TT: Thermal diffusion βˆ«ΞΊβˆ‡TΒ·βˆ‡T' dV - for (i, T_i) in enumerate(T_local) - for (j, T_j) in enumerate(T_local) - K_local[T_i, T_j] += ΞΊ * (dN[i] β‹… dN[j]) * vol - end - end - - # K_uu: Elasticity ∫C:Ξ΅:Ξ΅ dV (simplified) - for i in 1:4, comp_i in 1:3 - u_i = u_local[(i-1)*3 + comp_i] - for j in 1:4, comp_j in 1:3 - u_j = u_local[(j-1)*3 + comp_j] - if comp_i == comp_j - K_local[u_i, u_j] += (Ξ» + 2*ΞΌ) * (dN[i] β‹… dN[j]) * vol - end - end - end - - # K_pp: Hydraulic diffusion ∫(k/ΞΌ_f)βˆ‡pΒ·βˆ‡p' dV (cell-local, 1Γ—1) - K_local[p_local[1], p_local[1]] += (k_perm / ΞΌ_f) * vol - - # K_φφ: Electric (edge-based, simplified 6Γ—6) - n_Ο† = length(Ο†_local) - for i in 1:n_Ο†, j in 1:n_Ο† - K_local[Ο†_local[i], Ο†_local[j]] += Οƒ_e * vol / (n_Ο† * n_Ο†) - end - - # K_cc: Chemical diffusion ∫D_effΒ·βˆ‡cΒ·βˆ‡c' dV (NEW!) - for (i, c_i) in enumerate(c_local) - for (j, c_j) in enumerate(c_local) - K_local[c_i, c_j] += D_0 * (dN[i] β‹… dN[j]) * vol - end - end - - # ---------------------------------------------------------------------- - # OFF-DIAGONAL BLOCKS (Coupling terms) - 20 BLOCKS NOW! - # ---------------------------------------------------------------------- - - # --- EXISTING 12 THM-E COUPLINGS --- - - # K_Tu & K_uT: Thermal expansion - for (T_i, i) in enumerate(T_local), comp_j in 1:3 - u_j_idx = (comp_j-1)*4 + 1 : (comp_j-1)*4 + 4 - e_Ξ± = [e_1, e_2, e_3][comp_j] - for (local_u_j, u_j) in enumerate(u_local[u_j_idx]) - val = thermal_expansion_coupling(Ξ±_T, E_young, Ξ½, dN[T_i], dN[local_u_j], e_Ξ±, vol) - K_local[T_i, u_j] += val - K_local[u_j, T_i] += val # Onsager reciprocity - end - end - - # K_up & K_pu: Biot poroelasticity - for comp_i in 1:3 - e_Ξ± = [e_1, e_2, e_3][comp_i] - u_i_idx = (comp_i-1)*4 + 1 : (comp_i-1)*4 + 4 - for (local_u_i, u_i) in enumerate(u_local[u_i_idx]) - val = biot_coupling(Ξ±_p, dN[local_u_i], e_Ξ±, vol) - K_local[u_i, p_local[1]] += val - K_local[p_local[1], u_i] += val # Reciprocity - end - end - - # K_TΟ† & K_Ο†T: Seebeck-Peltier thermoelectric (simplified - node-edge coupling) - for i_T in 1:4, i_Ο† in 1:length(Ο†_local) - T_idx = T_local[i_T] - Ο†_idx = Ο†_local[i_Ο†] - val = S_seebeck * sum(dN[i_T]) * vol / length(Ο†_local) # Simplified - K_local[T_idx, Ο†_idx] += val - K_local[Ο†_idx, T_idx] += val - end - - # K_pΟ† & K_Ο†p: Electro-osmotic (simplified - cell-edge coupling) - for i_Ο† in 1:length(Ο†_local) - Ο†_idx = Ο†_local[i_Ο†] - val = ΞΆ * vol / length(Ο†_local) # Simplified - K_local[p_local[1], Ο†_idx] += val - K_local[Ο†_idx, p_local[1]] += val - end - - # K_uΟ† & K_Ο†u: Piezoelectric (Tensor{3,3}! - node-edge coupling) - for comp_i in 1:3, node_i in 1:4, i_Ο† in 1:length(Ο†_local) - u_idx = u_local[(node_i-1)*3 + comp_i] - Ο†_idx = Ο†_local[i_Ο†] - # Simplified: use diagonal of 3rd-order tensor - val = e_piezo[comp_i,comp_i,comp_i] * dN[node_i][comp_i] * vol / length(Ο†_local) - K_local[u_idx, Ο†_idx] += val - K_local[Ο†_idx, u_idx] += val # Reciprocity - end - - # --- NEW 8 CHEMICAL COUPLINGS! --- - - # K_uc & K_cu: Chemomechanical expansion - for comp_i in 1:3 - e_Ξ± = [e_1, e_2, e_3][comp_i] - for node_i in 1:4, node_c in 1:4 - u_idx = u_local[(node_i-1)*3 + comp_i] - c_idx = c_local[node_c] - val = chemomechanical_coupling(Ξ±_c, dN[node_i], e_Ξ±, vol) - K_local[u_idx, c_idx] += val - K_local[c_idx, u_idx] += val # Reciprocity - end - end - - # K_Tc & K_cT: Chemical reaction heat + Thermal diffusion (Soret) - for i_T in 1:4, i_c in 1:4 - T_idx = T_local[i_T] - c_idx = c_local[i_c] - # Reaction heat - val_rxn = chemical_reaction_heat(H_rxn, N[i_T], N[i_c], vol) - K_local[T_idx, c_idx] += val_rxn - - # Soret effect (NOT symmetric!) - val_soret = thermal_diffusion_coupling(D_T, dN[i_c], dN[i_T], vol) - K_local[c_idx, T_idx] += val_soret - end - - # K_pc & K_cp: Chemical fluid source + Pressure-dependent diffusion - for i_c in 1:4 - c_idx = c_local[i_c] - # Fluid source from reaction (simplified) - val_src = Ξ½_f * k_rxn * N[i_c] * vol - K_local[p_local[1], c_idx] += val_src - - # Pressure-dependent diffusion (simplified) - val_diff = Ξ±_D * D_0 * sum(dN[i_c]) * vol - K_local[c_idx, p_local[1]] += val_diff - end - - # K_Ο†c & K_cΟ†: Ionic migration (node-edge coupling) - for i_Ο† in 1:length(Ο†_local), i_c in 1:4 - Ο†_idx = Ο†_local[i_Ο†] - c_idx = c_local[i_c] - val = z_F_Dm * sum(dN[i_c]) * vol / length(Ο†_local) # Simplified - K_local[Ο†_idx, c_idx] += val - K_local[c_idx, Ο†_idx] += val # Reciprocity - end - - end # Integration points - - # Apply body forces (small for demo) - F_local[T_local] .+= 0.01 # Heat source - - println(" πŸ“€ Scattering coupled local matrix ($(n_local)Γ—$(n_local)) to global") - - # Scatter to global (ONE operation for ALL physics!) - for i_local in 1:n_local, j_local in 1:n_local - i_global = dof_map[i_local] - j_global = dof_map[j_local] - K[i_global, j_global] += K_local[i_local, j_local] - end - for i_local in 1:n_local - i_global = dof_map[i_local] - F[i_global] += F_local[i_local] - end - end # Element loop - - println("\nβœ“ Assembly complete!") - println(" ONE coupled system matrix: $(size(K))") - println(" Total non-zeros: $(nnz(K))") - - # ======================================================================== - # BOUNDARY CONDITIONS AND SOLVE - # ======================================================================== - - println("\n" * "="^70) - println("APPLYING BOUNDARY CONDITIONS AND SOLVING") - println("="^70) - - println("\nBoundary conditions (FULL FIVE-FIELD MULTI-PHYSICS):") - println(" Thermal: Node 1 fixed at T=0 K (thermal ground)") - println(" Mechanical: Node 1 fully fixed u=(0,0,0) (mechanical ground)") - println(" Mechanical: Node 2 ux=0 (prevent rigid rotation)") - println(" Electric: Edge 1 fixed at Ο†=0 V (electric ground)") - println(" Chemical: Node 1 fixed at c=0 mol/mΒ³ (chemical ground)") # NEW! - println(" Hydraulic: Natural BCs (traction-free, no flow prescribed)") - - # Get global DOF indices for BCs - # Debug: Print field starting indices - println("\nDEBUG DOF layout:") - println(" T: $(offset_T+1):$(offset_T+n_T_total)") - println(" u: $(offset_u+1):$(offset_u+n_u_total)") - println(" p: $(offset_p+1):$(offset_p+n_p_total)") - println(" Ο†: $(offset_Ο†+1):$(offset_Ο†+(n_total-offset_c-n_c_total))") - println(" c: $(offset_c+1):$n_total") - println(" Total: $n_total DOFs") - - # Node 1: T, u, c all fixed - # Edge 1: Ο† fixed - bc_dofs = [ - offset_T+1, # T at node 1 - offset_u+1, offset_u+2, offset_u+3, # u at node 1 - offset_u+4, # ux at node 2 (prevent rotation) - offset_Ο†+1, # Ο† at edge 1 - offset_c+1 # c at node 1 (NEW!) - ] - bc_vals = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # 7 DOFs fixed - - for dof in bc_dofs - K[dof, :] .= 0.0 - K[:, dof] .= 0.0 - K[dof, dof] = 1.0 - F[dof] = 0.0 - end - - # Solve - println("\n🎯 Solving coupled system...") - println(" Matrix size: $(size(K))") - println(" Non-zeros: $(nnz(K))") - println(" Condition number estimate: checking...") - - # Add small regularization - Ξ΅_reg = 1e-12 - for i in 1:n_total - K[i,i] += Ξ΅_reg - end - - println(" Added regularization (Ξ΅=$Ξ΅_reg) for numerical stability") - - # Solve - sol = try - result = K \ F - println(" βœ“ Solution converged!") - result - catch e - println(" ERROR: System still singular!") - println(" This indicates physical model needs more constraints") - rethrow(e) - end - - # Extract fields (handle edge sharing) - T_sol = sol[offset_T+1:offset_T+n_T_total] - u_sol = sol[offset_u+1:offset_u+n_u_total] - p_sol = sol[offset_p+1:offset_p+n_p_total] - # Ο† and c may overlap in DOF numbering - just get last part - if offset_c > offset_Ο† - Ο†_sol = sol[offset_Ο†+1:offset_c] - c_sol = sol[offset_c+1:end] - else - # They overlap - extract what we can - Ο†_sol = Float64[] - c_sol = sol[offset_c+1:end] - end - - # ======================================================================== - # RESULTS - # ======================================================================== - - println("\n" * "="^70) - println("✨ SOLUTION (REAL FIVE-FIELD PHYSICS!)") - println("="^70) - - println("\n🌑️ Temperature field:") - for i in 1:n_T_total - println(" Node $i: T = $(T_sol[i]) K") - end - - println("\nπŸ—οΈ Displacement field:") - for i in 1:5 - u_i = u_sol[(i-1)*3+1:i*3] - println(" Node $i: u = ($(u_i[1]), $(u_i[2]), $(u_i[3])) m") - end - - println("\nπŸ’§ Pore pressure field:") - for i in 1:n_p_total - println(" Cell $i: p = $(p_sol[i]) Pa") - end - - println("\n⚑ Electric potential (edges):") - if !isempty(Ο†_sol) - for i in 1:length(Ο†_sol) - println(" Edge $i: Ο† = $(Ο†_sol[i]) V") - end - else - println(" (Edge DOFs overlap with other fields)") - end - - println("\nπŸ§ͺ Chemical concentration field (NEW!):") - for i in 1:n_c_total - println(" Node $i: c = $(c_sol[i]) mol/mΒ³") - end - - # ======================================================================== - # ACHIEVEMENTS - # ======================================================================== - - println("\n" * "="^70) - println("πŸŽ‰ ACHIEVEMENTS UNLOCKED:") - println("="^70) - println(" βœ… ONE element type with FIVE physics fields!") - println(" βœ… ONE local coupled matrix per element (27Γ—27)") - println(" βœ… ALL physics assembled together (true coupling!)") - println(" βœ… Thermo-mechanical coupling: K_Tu, K_uT (thermal expansion)") - println(" βœ… Hydro-mechanical coupling: K_up, K_pu (Biot poroelasticity)") - println(" βœ… Thermal-hydraulic coupling: K_Tp, K_pT (thermal pressurization)") - println(" βœ… Electro-osmotic coupling: K_Ο†p, K_pΟ† (electrokinetic flow)") - println(" βœ… Thermo-electric coupling: K_TΟ†, K_Ο†T (Seebeck/Peltier)") - println(" βœ… Piezoelectric coupling: K_uΟ†, K_Ο†u (Tensor{3,3} elegance!)") - println(" βœ… Chemomechanical coupling: K_uc, K_cu (swelling/shrinkage) πŸ†•") - println(" βœ… Chemical reaction heat: K_Tc (exothermic/endothermic) πŸ†•") - println(" βœ… Thermal diffusion: K_cT (Soret effect) πŸ†•") - println(" βœ… Chemical fluid source: K_pc (dissolution/precipitation) πŸ†•") - println(" βœ… Pressure-dependent diffusion: K_cp πŸ†•") - println(" βœ… Ionic migration: K_Ο†c, K_cΟ† (electrophoresis) πŸ†•") - println(" βœ… Total: 20 off-diagonal coupling blocks! (PENTA-PHYSICS!)") - println(" βœ… Modular coupling functions (inlined for zero overhead)") - println(" βœ… 3rd-order tensor formulation (e_kij via Tensor{3,3})") - println(" βœ… Onsager reciprocity respected (symmetric couplings)") - println(" βœ… Local-to-global mapping via type system") - println(" βœ… REAL thermal diffusion (βˆ«ΞΊβˆ‡TΒ·βˆ‡T' dV)") - println(" βœ… REAL 3D elasticity (∫C:Ξ΅:Ξ΅ dV)") - println(" βœ… REAL chemical diffusion (∫Dβˆ‡cΒ·βˆ‡c' dV) πŸ†•") - println(" βœ… Zero-allocation Tensors.jl operations") - println(" βœ… Cell-local pressure DOFs (discontinuous)") - println(" βœ… Edge-based electric DOFs") - println(" βœ… Vertex-based chemical DOFs (continuous) πŸ†•") - println(" βœ… Full $(n_total) Γ— $(n_total) coupled system solved") - println(" βœ… Type-safe field access: .T, .u, .p, .Ο†, .c") - println("="^70) - - println("\nπŸ’‘ THIS IS THE MOTHER OF ALL MULTI-PHYSICS SYSTEMS!") - println(" ONE element β†’ ONE local matrix β†’ ALL FIVE FIELDS COUPLED!") - println(" T ↔ u (thermal expansion), T ↔ p (thermal pressurization)") - println(" T ↔ Ο† (Seebeck/Peltier), T ↔ c (reaction heat + Soret)") - println(" u ↔ p (Biot poroelasticity), u ↔ Ο† (piezoelectric)") - println(" u ↔ c (chemomechanical swelling)") - println(" p ↔ Ο† (electro-osmotic), p ↔ c (fluid source + diff.)") - println(" Ο† ↔ c (ionic migration)") - println(" β†’ Complete system: 5 fields Γ— 10 pairs = 20 coupling blocks!") - println(" β†’ Applications: Geothermal, nuclear waste, CO2, batteries,") - println(" β†’ concrete durability, corrosion, electrochemistry!") - println(" β†’ Natural coupling, type-safe, composable, ULTIMATE! πŸš€") - - @test length(T_sol) == n_T_total - @test length(u_sol) == n_u_total - @test length(p_sol) == n_p_total - # Ο† and c DOFs may overlap in global numbering - skip test - @test length(c_sol) == n_c_total # NEW! - # Note: Total may be less than sum due to shared DOFs between fields - -end # testset