From b6a5fbce37205c6d9849a2eeeda63d6a7b8dd564 Mon Sep 17 00:00:00 2001 From: Rainer Rodrigues Date: Sat, 30 May 2026 07:12:26 +0100 Subject: [PATCH 1/2] feat: implement Robin boundary conditions (convection) for thermal optimization --- Project.toml | 86 ++++++++-------- demo.jl | 79 +++++++++++++++ src/FEA/solvers_api.jl | 30 ++++++ src/TopOptProblems/TopOptProblems.jl | 1 + src/TopOptProblems/matrices_and_vectors.jl | 108 ++++++++++++++++++++- src/TopOptProblems/problem_types.jl | 14 ++- src/thermal_quad4_element.jl | 30 ++++++ test/examples/heat_sink.jl | 21 +++- 8 files changed, 320 insertions(+), 49 deletions(-) create mode 100644 demo.jl create mode 100644 src/thermal_quad4_element.jl diff --git a/Project.toml b/Project.toml index 44961617..8db61e75 100644 --- a/Project.toml +++ b/Project.toml @@ -1,49 +1,11 @@ -authors = ["Mohamed Tarek ", "Yijiang Huang "] name = "TopOpt" uuid = "53a1e1a5-51bb-58a9-8a02-02056cc81109" version = "0.11.0" - -[compat] -AbstractDifferentiation = "0.6" -ChainRulesCore = "1" -ColorSchemes = "3" -Distributions = "0.24, 0.25" -Einsum = "0.4" -Ferrite = "=0.3.0" -FileIO = "1" -Flux = "0.11, 0.12, 0.13, 0.14" -ForwardDiff = "<0.10.35, 0.10, 1" -GeometryBasics = "0.5" -GeometryTypes = "0.8" -IterativeSolvers = "0.9" -JSON = "0.21" -JuliaFormatter = "1" -MacroTools = "0.5" -Makie = "0.23" -MappedArrays = "0.4" -NearestNeighbors = "0.4" -Nonconvex = "2" -NonconvexMMA = "1" -NonconvexPercival = "0.1.5" -NonconvexSemidefinite = "0.1.10" -Parameters = "0.12, 0.13" -Preconditioners = "0.3, 0.4, 0.5, 0.6" -Preferences = "1" -Reexport = "1" -Revise = "3" -Setfield = "0.7, 0.8, 1" -StaticArrays = "1" -Statistics = "1" -StatsFuns = "0.9, 1" -Tensors = "1.4" -TimerOutputs = "0.5" -VTKDataTypes = "0.1" -WriteVTK = "1.9" -Zygote = "0.6" -julia = "1.9" +authors = ["Mohamed Tarek ", "Yijiang Huang "] [deps] AbstractDifferentiation = "c29ec348-61ec-40c8-8164-b8c60e9d9f3d" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" @@ -59,6 +21,7 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" MappedArrays = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" Nonconvex = "01bcebdf-4d21-426d-b5c4-6132c1619978" @@ -81,11 +44,50 @@ Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" VTKDataTypes = "10d27dd1-1d0f-5a4c-b178-bd2d0045a217" +WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [extensions] TopOptMakieExt = "Makie" -[weakdeps] -Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +[compat] +AbstractDifferentiation = "0.6" +CairoMakie = "0.14.0" +ChainRulesCore = "1" +ColorSchemes = "3" +Distributions = "0.24, 0.25" +Einsum = "0.4" +Ferrite = "=0.3.0" +FileIO = "1" +Flux = "0.11, 0.12, 0.13, 0.14" +ForwardDiff = "<0.10.35, 0.10, 1" +GeometryBasics = "0.5" +GeometryTypes = "0.8" +IterativeSolvers = "0.9" +JSON = "0.21" +JuliaFormatter = "1" +MacroTools = "0.5" +Makie = "0.23" +MappedArrays = "0.4" +NearestNeighbors = "0.4" +Nonconvex = "2" +NonconvexMMA = "1" +NonconvexPercival = "0.1.5" +NonconvexSemidefinite = "0.1.10" +Parameters = "0.12, 0.13" +Preconditioners = "0.3, 0.4, 0.5, 0.6" +Preferences = "1" +Reexport = "1" +Revise = "3" +Setfield = "0.7, 0.8, 1" +StaticArrays = "1" +Statistics = "1" +StatsFuns = "0.9, 1" +Tensors = "1.4" +TimerOutputs = "0.5" +VTKDataTypes = "0.1" +WGLMakie = "0.12.0" +WriteVTK = "1.9" +Zygote = "0.6" +julia = "1.9" diff --git a/demo.jl b/demo.jl new file mode 100644 index 00000000..6fc6739c --- /dev/null +++ b/demo.jl @@ -0,0 +1,79 @@ +using StaticArrays +using LinearAlgebra + +# 1. Define the Gauss Quadrature points and weights for a 2x2 rule +# The points are at ±1/√3 +const GAUSS_PTS = SVector(-1/sqrt(3), 1/sqrt(3)) +const GAUSS_WTS = SVector(1.0, 1.0) + +""" + quad4_thermal_matrix(nodes, k) + +Calculates the 4x4 local thermal conductivity matrix for a Quad4 element. +`nodes` is a 4x2 matrix (or vector of vectors) of the [x, y] coordinates of the 4 nodes. +`k` is the thermal conductivity of the element. +""" +function quad4_thermal_matrix(nodes::SMatrix{4, 2, Float64}, k::Float64) + # Initialize an empty 4x4 mutable static matrix for performance + Ke = @MMatrix zeros(4, 4) + + # Loop over the 2x2 Gauss integration points + for i in 1:2 + for j in 1:2 + xi = GAUSS_PTS[i] + eta = GAUSS_PTS[j] + weight = GAUSS_WTS[i] * GAUSS_WTS[j] + + # 2. Derivatives of shape functions with respect to local coords (ξ, η) + # dN/dξ + dN_dxi = SVector( + -0.25 * (1 - eta), + 0.25 * (1 - eta), + 0.25 * (1 + eta), + -0.25 * (1 + eta) + ) + # dN/dη + dN_deta = SVector( + -0.25 * (1 - xi), + -0.25 * (1 + xi), + 0.25 * (1 + xi), + 0.25 * (1 - xi) + ) + + # Combine into a 2x4 matrix: [dN/dξ ; dN/dη] + dN_dxi_eta = SMatrix{2, 4}( + dN_dxi[1], dN_deta[1], + dN_dxi[2], dN_deta[2], + dN_dxi[3], dN_deta[3], + dN_dxi[4], dN_deta[4] + ) + + # 3. Calculate the Jacobian matrix (J = dN_dxi_eta * coordinates) + J = dN_dxi_eta * nodes + detJ = det(J) + + # 4. Calculate the B matrix (derivatives w.r.t real coords x, y) + # B = J \ dN_dxi_eta (using left division for inverse multiplication) + B = J \ dN_dxi_eta + + # 5. Add the contribution of this Gauss point to the element matrix + # Integral evaluation: B^T * k * B * det(J) * weight + Ke += B' * k * B * detJ * weight + end + end + + return SMatrix(Ke) # Return as immutable for safety and speed +end + +""" + penalized_thermal_matrix(nodes, x_e, p, k0, kmin) + +Returns the local conductivity matrix for an element, scaled by its SIMP density. +""" +function penalized_thermal_matrix(nodes::SMatrix{4, 2, Float64}, x_e::Float64, p::Float64, k0::Float64, kmin::Float64) + # 1. Calculate the penalized thermal conductivity + k_penalized = kmin + (x_e^p) * (k0 - kmin) + + # 2. Generate the standard local matrix using that conductivity + return quad4_thermal_matrix(nodes, k_penalized) +end \ No newline at end of file diff --git a/src/FEA/solvers_api.jl b/src/FEA/solvers_api.jl index 69b36511..a59c8388 100644 --- a/src/FEA/solvers_api.jl +++ b/src/FEA/solvers_api.jl @@ -212,7 +212,22 @@ function (s::GenericFEASolver{T,Physics,Solver})( # Multiple RHS columns - solve each one # Assemble stiffness matrix (and force vector if assemble_f=true) assemble!(s.globalinfo, s.problem, s.elementinfo, s.vars, getpenalty(s), s.xmin; assemble_f=assemble_f) + + if Physics === HeatTransfer + K_conv = assemble_convection_matrix(s.problem) + # Mutate the pre-allocated matrix in-place so the Ferrite assembler doesn't break + K_data = s.globalinfo.K.data + for col in 1:size(K_conv, 2) + for ptr in K_conv.colptr[col]:(K_conv.colptr[col+1]-1) + row = K_conv.rowval[ptr] + val = K_conv.nzval[ptr] + # Add the convection values directly into the existing memory slots + K_data[row, col] += val + end + end + end + # Solve for each column of the matrix RHS for j in 1:size(rhs, 2) # Get the RHS for this column - use the provided matrix columns @@ -242,6 +257,21 @@ function (s::GenericFEASolver{T,Physics,Solver})( # Single RHS case (original behavior) assemble!(s.globalinfo, s.problem, s.elementinfo, s.vars, getpenalty(s), s.xmin; assemble_f=assemble_f) + if Physics === HeatTransfer + K_conv = assemble_convection_matrix(s.problem) + + # Mutate the pre-allocated matrix in-place so the Ferrite assembler doesn't break + K_data = s.globalinfo.K.data + for col in 1:size(K_conv, 2) + for ptr in K_conv.colptr[col]:(K_conv.colptr[col+1]-1) + row = K_conv.rowval[ptr] + val = K_conv.nzval[ptr] + # Add the convection values directly into the existing memory slots + K_data[row, col] += val + end + end + end + # Apply boundary conditions to rhs if needed (only for vectors) if !assemble_f && rhs !== s.globalinfo.f && ndims(rhs) == 1 rhs = copy(rhs) diff --git a/src/TopOptProblems/TopOptProblems.jl b/src/TopOptProblems/TopOptProblems.jl index bf3e897c..fb90dac3 100644 --- a/src/TopOptProblems/TopOptProblems.jl +++ b/src/TopOptProblems/TopOptProblems.jl @@ -46,6 +46,7 @@ export PointLoadCantilever, YoungsModulus, assemble, assemble_f!, + assemble_convection_matrix, buckling, get_Kσs, RaggedArray, diff --git a/src/TopOptProblems/matrices_and_vectors.jl b/src/TopOptProblems/matrices_and_vectors.jl index 3609effb..24419390 100644 --- a/src/TopOptProblems/matrices_and_vectors.jl +++ b/src/TopOptProblems/matrices_and_vectors.jl @@ -319,7 +319,38 @@ function _make_dloads(fes, problem::HeatTransferTopOptProblem, facevalues) dloads[cellid] = fe end end - + convectiondict = getconvectiondict(problem) + if !isempty(convectiondict) + for (faceset_name, (h, T_inf)) in convectiondict + faceset = getfacesets(problem)[faceset_name] + for (cellid, faceid) in faceset + boundary_matrix[faceid, cellid] || + throw("Face $((cellid, faceid)) not on boundary.") + + fe = dloads[cellid] + getcoordinates!(cell_coords, grid, cellid) + reinit!(facevalues, cell_coords, faceid) + + for q_point in 1:getnquadpoints(facevalues) + dΓ = getdetJdV(facevalues, q_point) + for i in 1:n_basefuncs + ϕ = shape_value(facevalues, q_point, i) + + # Convection load: fe[i] += ∫ ϕi * h * T_inf * dΓ + # (Note: h * T_inf has the same units as heat flux q) + convection_flux = h * T_inf + + if fe isa SArray + fe = @set fe[i] += ϕ * convection_flux * dΓ + else + fe[i] += ϕ * convection_flux * dΓ + end + end + end + dloads[cellid] = fe + end + end + end return dloads end @@ -459,6 +490,81 @@ function _make_Kes_and_weights_heat( return Kes, weights end +using SparseArrays + +""" + assemble_convection_matrix(problem::HeatConductionProblem) + +Builds the constant global convection matrix for Robin boundary conditions. +Evaluates the surface integral: ∫ h * N^T * N dΓ +""" +function assemble_convection_matrix(problem::HeatConductionProblem{dim, T}) where {dim, T} + convectiondict = getconvectiondict(problem) + dh = getdh(problem) + + # If no convection is defined, return an empty sparse matrix + if isempty(convectiondict) + return spzeros(T, ndofs(dh), ndofs(dh)) + end + + grid = dh.grid + boundary_matrix = grid.boundary_matrix + + # We need FaceScalarValues to evaluate the surface integrals + ip = dh.field_interpolations[1] + qr = Ferrite.QuadratureRule{dim - 1, Ferrite.RefCube}(2) + facevalues = Ferrite.FaceScalarValues(qr, ip) + + n_basefuncs = getnbasefunctions(facevalues) + cell_coords = zeros(Ferrite.Vec{dim,T}, nnodespercell(problem)) + + I = Int[] + J = Int[] + V = T[] + + global_dofs = zeros(Int, ndofs_per_cell(dh)) + + for (faceset_name, (h, T_inf)) in convectiondict + faceset = getfacesets(problem)[faceset_name] + + for (cellid, faceid) in faceset + boundary_matrix[faceid, cellid] || throw("Face not on boundary.") + + getcoordinates!(cell_coords, grid, cellid) + reinit!(facevalues, cell_coords, faceid) + celldofs!(global_dofs, dh, cellid) + + # Creating a local face matrix + K_face = zeros(T, n_basefuncs, n_basefuncs) + + # Integrate: ∫ h * N^T * N dΓ + for q_point in 1:getnquadpoints(facevalues) + dΓ = getdetJdV(facevalues, q_point) + + for i in 1:n_basefuncs + ϕ_i = shape_value(facevalues, q_point, i) + for j in 1:n_basefuncs + ϕ_j = shape_value(facevalues, q_point, j) + + K_face[i, j] += ϕ_i * ϕ_j * h * dΓ + end + end + end + + # Scatter local face matrix to global I, J, V vectors + for i in 1:n_basefuncs + for j in 1:n_basefuncs + push!(I, global_dofs[i]) + push!(J, global_dofs[j]) + push!(V, K_face[i, j]) + end + end + end + end + + return sparse(I, J, V, ndofs(dh), ndofs(dh)) +end + # Fallback for non-static arrays # function _make_Kes_and_weights_heat( # dh::DofHandler{dim, N, T}, diff --git a/src/TopOptProblems/problem_types.jl b/src/TopOptProblems/problem_types.jl index 2d6aa0c0..e271bea5 100644 --- a/src/TopOptProblems/problem_types.jl +++ b/src/TopOptProblems/problem_types.jl @@ -640,6 +640,7 @@ getmetadata(p::HeatTransferTopOptProblem) = p.metadata getdh(p::HeatTransferTopOptProblem) = p.ch.dh getpressuredict(p::HeatTransferTopOptProblem{dim,T}) where {dim,T} = Dict{String,T}() getheatfluxdict(p::HeatTransferTopOptProblem{dim,T}) where {dim,T} = Dict{String,T}() +getconvectiondict(p::HeatTransferTopOptProblem{dim,T}) where {dim,T} = Dict{String,Tuple{T,T}}() getfacesets(p::HeatTransferTopOptProblem) = getdh(p).grid.facesets Ferrite.getncells(problem::HeatTransferTopOptProblem) = Ferrite.getncells(getdh(problem).grid) getgeomorder(p::HeatTransferTopOptProblem) = nnodespercell(p) in (9, 27) ? 2 : 1 @@ -693,12 +694,14 @@ struct HeatConductionProblem{ Tr<:RectilinearGrid{dim, T, N, M}, Tc<:ConstraintHandler{<:DofHandler{dim, <:Cell{dim, N, M}, T}, T}, Th<:AbstractDict{String,T}, + Tconv<:AbstractDict{String,Tuple{T,T}}, Tm<:Metadata, } <: HeatTransferTopOptProblem{dim, T} rect_grid::Tr k::T ch::Tc heatfluxdict::Th + convectiondict::Tconv metadata::Tm end @@ -708,6 +711,8 @@ end getheatfluxdict(p::HeatConductionProblem) = p.heatfluxdict +getconvectiondict(p::HeatConductionProblem) = p.convectiondict + """ HeatConductionProblem(::Type{Val{CellType}}, nels, sizes, k=1.0; Tleft=0.0, Tright=0.0, heatflux=Dict{String,Float64}()) @@ -734,6 +739,7 @@ function HeatConductionProblem( Tleft=0.0, Tright=0.0, heatflux=Dict{String,Float64}(), + convection=Dict{String,Tuple{Float64,Float64}}() ) where {dim, CellType} T = float(promote_type(eltype(sizes), typeof(k), typeof(Tleft), typeof(Tright))) @@ -781,9 +787,15 @@ function HeatConductionProblem( for (key, val) in heatflux heatfluxdict[key] = T(val) end + + convectiondict = Dict{String, Tuple{T, T}}() + for (key, val) in convection + # val[1] is h, val[2] is T_ambient + convectiondict[key] = (T(val[1]), T(val[2])) + end return HeatConductionProblem( - rect_grid, T(k), ch, heatfluxdict, metadata + rect_grid, T(k), ch, heatfluxdict, convectiondict, metadata ) end diff --git a/src/thermal_quad4_element.jl b/src/thermal_quad4_element.jl new file mode 100644 index 00000000..97fc38c7 --- /dev/null +++ b/src/thermal_quad4_element.jl @@ -0,0 +1,30 @@ +using StaticArrays +using LinearAlgebra +using Ferrite + +function thermal_conductivity_matrix(cellvalues::CellScalarValues{dim, T}, k_prop::Float64) where {dim, T} + + n_basefuncs = getnbasefunctions(cellvalues) + + # Initialize an empty mutable static matrix for speed + Ke = @MMatrix zeros(T, n_basefuncs, n_basefuncs) + + for q_point in 1:getnquadpoints(cellvalues) + # dΩ is the determinant of the Jacobian multiplied by the quadrature weight + dΩ = getdetJdV(cellvalues, q_point) + + # Integrate B^T * k * B + for i in 1:n_basefuncs + ∇N_i = shape_gradient(cellvalues, q_point, i) + for j in 1:n_basefuncs + ∇N_j = shape_gradient(cellvalues, q_point, j) + + # The dot product of the shape gradients scales the conductivity + Ke[i, j] += (∇N_i ⋅ ∇N_j) * k_prop * dΩ + end + end + end + + # Return as an immutable SMatrix for thread safety + return SMatrix(Ke) +end diff --git a/test/examples/heat_sink.jl b/test/examples/heat_sink.jl index 362d0005..ac11f4f3 100644 --- a/test/examples/heat_sink.jl +++ b/test/examples/heat_sink.jl @@ -7,19 +7,30 @@ println("Heat Sink Topology Optimization Example") println("=" ^ 50) # Problem parameters -nels = (60, 30) # Mesh resolution -sizes = (1.0, 1.0) # Element sizes +nels = (60, 30, 10) +# nels = (60, 30) # Mesh resolution +sizes = (1.0, 1.0, 1.0) +# sizes = (1.0, 1.0) # Element sizes k = 1.0 # Thermal conductivity # Heat flux on top boundary (W/m²) - heat entering the domain -heatflux = Dict{String,Float64}("top" => 100.0) +heatflux = Dict{String,Float64}("bottom" => 100.0) V = 0.5 # Volume fraction constraint +convection = Dict{String,Tuple{Float64,Float64}}( + "top" => (10.0, 20.0), + "left" => (10.0, 20.0), + "right" => (10.0, 20.0) +) + # Create heat conduction problem # Temperature is fixed at left and right edges (heat sink fins) # Heat flux enters from top boundary problem = HeatConductionProblem( Val{:Linear}, nels, sizes, k; - Tleft=100.0, Tright=0.0, heatflux=heatflux + Tleft=100.0, + Tright=0.0, + heatflux=heatflux, + convection=convection ) println("Created heat conduction problem with $(Ferrite.getncells(problem)) elements") @@ -82,7 +93,7 @@ println(" Gradient norm: $(norm(grad))") println(" All gradients negative (expected): $(all(grad .< 0))") using Makie -using CairoMakie +using WGLMakie # alternatively, `using GLMakie` fig = visualize(problem; topology=result.minimizer) Makie.display(fig) From 2b169e8016cbaaf837e89fa27c1474d02f622f13 Mon Sep 17 00:00:00 2001 From: Rainer Rodrigues Date: Mon, 1 Jun 2026 10:18:06 +0100 Subject: [PATCH 2/2] feat: implement Robin boundary conditions (convection) for thermal optimization --- src/FEA/solvers_api.jl | 211 ++++++++++++++++----- src/TopOptProblems/TopOptProblems.jl | 4 +- src/TopOptProblems/matrices_and_vectors.jl | 69 +++---- src/TopOptProblems/problem_types.jl | 62 +++--- test/topopt_problems/test_convection.jl | 35 ++++ 5 files changed, 265 insertions(+), 116 deletions(-) create mode 100644 test/topopt_problems/test_convection.jl diff --git a/src/FEA/solvers_api.jl b/src/FEA/solvers_api.jl index a59c8388..6065bfbf 100644 --- a/src/FEA/solvers_api.jl +++ b/src/FEA/solvers_api.jl @@ -72,8 +72,15 @@ export GenericFEASolver # These functions dispatch on the linear solver type to solve the system # Direct solver (factorization-based) -function solve_system!(::Type{DirectSolver}, solver::GenericFEASolver{T,Physics,DirectSolver}, K, f, lhs; - reuse_fact=false, safe=false) where {T,Physics} +function solve_system!( + ::Type{DirectSolver}, + solver::GenericFEASolver{T,Physics,DirectSolver}, + K, + f, + lhs; + reuse_fact=false, + safe=false, +) where {T,Physics} if safe m = meandiag(K) for i in 1:size(K, 1) @@ -102,8 +109,16 @@ function solve_system!(::Type{DirectSolver}, solver::GenericFEASolver{T,Physics, end # CG with assembled matrix -function solve_system!(::Type{CGAssemblySolver}, solver::GenericFEASolver{T,Physics,CGAssemblySolver}, K, f, lhs; - safe=false, initially_zero=true, kwargs...) where {T,Physics} +function solve_system!( + ::Type{CGAssemblySolver}, + solver::GenericFEASolver{T,Physics,CGAssemblySolver}, + K, + f, + lhs; + safe=false, + initially_zero=true, + kwargs..., +) where {T,Physics} if safe m = meandiag(K) for i in 1:size(K, 1) @@ -151,8 +166,15 @@ function solve_system!(::Type{CGAssemblySolver}, solver::GenericFEASolver{T,Phys end # Matrix-free CG -function solve_system!(::Type{CGMatrixFreeSolver}, solver::GenericFEASolver{T,Physics,CGMatrixFreeSolver}, K, f, lhs; - initially_zero=true, kwargs...) where {T,Physics} +function solve_system!( + ::Type{CGMatrixFreeSolver}, + solver::GenericFEASolver{T,Physics,CGMatrixFreeSolver}, + K, + f, + lhs; + initially_zero=true, + kwargs..., +) where {T,Physics} @unpack cg_max_iter, abstol, cg_statevars = solver @unpack preconditioner, preconditioner_initialized = solver @unpack elementinfo, meandiag, vars, xmin, fixed_dofs, free_dofs, xes = solver @@ -162,8 +184,16 @@ function solve_system!(::Type{CGMatrixFreeSolver}, solver::GenericFEASolver{T,Ph # Build matrix-free operator penalty = getpenalty(solver) operator = MatrixFreeOperator( - f, elementinfo, meandiag, vars, xes, - fixed_dofs, free_dofs, xmin, penalty, solver.conv + f, + elementinfo, + meandiag, + vars, + xes, + fixed_dofs, + free_dofs, + xmin, + penalty, + solver.conv, ) if !(preconditioner === identity) @@ -205,21 +235,29 @@ function (s::GenericFEASolver{T,Physics,Solver})( assemble_f=true, rhs=assemble_f ? s.globalinfo.f : s.rhs, lhs=assemble_f ? s.u : s.lhs, - kwargs... + kwargs..., ) where {T,Physics,Solver,safe} # Handle matrix RHS by solving for each column if ndims(rhs) == 2 && size(rhs, 2) > 1 # Multiple RHS columns - solve each one # Assemble stiffness matrix (and force vector if assemble_f=true) - assemble!(s.globalinfo, s.problem, s.elementinfo, s.vars, getpenalty(s), s.xmin; assemble_f=assemble_f) + assemble!( + s.globalinfo, + s.problem, + s.elementinfo, + s.vars, + getpenalty(s), + s.xmin; + assemble_f=assemble_f, + ) if Physics === HeatTransfer K_conv = assemble_convection_matrix(s.problem) - + # Mutate the pre-allocated matrix in-place so the Ferrite assembler doesn't break K_data = s.globalinfo.K.data for col in 1:size(K_conv, 2) - for ptr in K_conv.colptr[col]:(K_conv.colptr[col+1]-1) + for ptr in K_conv.colptr[col]:(K_conv.colptr[col + 1] - 1) row = K_conv.rowval[ptr] val = K_conv.nzval[ptr] # Add the convection values directly into the existing memory slots @@ -234,43 +272,70 @@ function (s::GenericFEASolver{T,Physics,Solver})( # Apply boundary conditions to the column rhs_j = copy(rhs[:, j]) apply_zero!(rhs_j, s.problem.ch) - + # Get the view of lhs for this column using @view macro lhs_j = @view lhs[:, j] - + # Pass initially_zero only for CG solvers if Solver === DirectSolver # Filter out kwargs that DirectSolver doesn't accept - filtered_kwargs = filter(p -> p.first ∉ (:initially_zero, :solver), collect(kwargs)) - solve_system!(Solver, s, s.globalinfo.K, rhs_j, lhs_j; - reuse_fact=(j > 1 || reuse_fact), safe=safe, filtered_kwargs...) + filtered_kwargs = filter( + p -> p.first ∉ (:initially_zero, :solver), collect(kwargs) + ) + solve_system!( + Solver, + s, + s.globalinfo.K, + rhs_j, + lhs_j; + reuse_fact=(j > 1 || reuse_fact), + safe=safe, + filtered_kwargs..., + ) else # For CG solvers, start from zero initial guess for each column lhs_j .= zero(T) - solve_system!(Solver, s, s.globalinfo.K, rhs_j, lhs_j; - reuse_fact=(j > 1 || reuse_fact), safe=safe, initially_zero=true, kwargs...) + solve_system!( + Solver, + s, + s.globalinfo.K, + rhs_j, + lhs_j; + reuse_fact=(j > 1 || reuse_fact), + safe=safe, + initially_zero=true, + kwargs..., + ) end end return nothing end # Single RHS case (original behavior) - assemble!(s.globalinfo, s.problem, s.elementinfo, s.vars, getpenalty(s), s.xmin; assemble_f=assemble_f) + assemble!( + s.globalinfo, + s.problem, + s.elementinfo, + s.vars, + getpenalty(s), + s.xmin; + assemble_f=assemble_f, + ) if Physics === HeatTransfer - K_conv = assemble_convection_matrix(s.problem) - - # Mutate the pre-allocated matrix in-place so the Ferrite assembler doesn't break - K_data = s.globalinfo.K.data - for col in 1:size(K_conv, 2) - for ptr in K_conv.colptr[col]:(K_conv.colptr[col+1]-1) - row = K_conv.rowval[ptr] - val = K_conv.nzval[ptr] - # Add the convection values directly into the existing memory slots - K_data[row, col] += val - end + K_conv = assemble_convection_matrix(s.problem) + + # Mutate the pre-allocated matrix in-place so the Ferrite assembler doesn't break + K_data = s.globalinfo.K.data + for col in 1:size(K_conv, 2) + for ptr in K_conv.colptr[col]:(K_conv.colptr[col + 1] - 1) + row = K_conv.rowval[ptr] + val = K_conv.nzval[ptr] + # Add the convection values directly into the existing memory slots + K_data[row, col] += val end end + end # Apply boundary conditions to rhs if needed (only for vectors) if !assemble_f && rhs !== s.globalinfo.f && ndims(rhs) == 1 @@ -279,27 +344,51 @@ function (s::GenericFEASolver{T,Physics,Solver})( end # Solve system (physics-independent, solver-algorithm dependent) - solve_system!(Solver, s, s.globalinfo.K, rhs, lhs; reuse_fact=reuse_fact, safe=safe, kwargs...) + solve_system!( + Solver, s, s.globalinfo.K, rhs, lhs; reuse_fact=reuse_fact, safe=safe, kwargs... + ) return nothing end # Show methods -function Base.show(io::IO, ::MIME{Symbol("text/plain")}, ::GenericFEASolver{T,LinearElasticity,DirectSolver}) where {T} +function Base.show( + io::IO, + ::MIME{Symbol("text/plain")}, + ::GenericFEASolver{T,LinearElasticity,DirectSolver}, +) where {T} return println(io, "TopOpt direct structural solver (GenericFEASolver)") end -function Base.show(io::IO, ::MIME{Symbol("text/plain")}, ::GenericFEASolver{T,HeatTransfer,DirectSolver}) where {T} +function Base.show( + io::IO, ::MIME{Symbol("text/plain")}, ::GenericFEASolver{T,HeatTransfer,DirectSolver} +) where {T} return println(io, "TopOpt direct heat transfer solver (GenericFEASolver)") end -function Base.show(io::IO, ::MIME{Symbol("text/plain")}, ::GenericFEASolver{T,LinearElasticity,CGAssemblySolver}) where {T} +function Base.show( + io::IO, + ::MIME{Symbol("text/plain")}, + ::GenericFEASolver{T,LinearElasticity,CGAssemblySolver}, +) where {T} return println(io, "TopOpt CG with assembly structural solver (GenericFEASolver)") end -function Base.show(io::IO, ::MIME{Symbol("text/plain")}, ::GenericFEASolver{T,HeatTransfer,CGAssemblySolver}) where {T} +function Base.show( + io::IO, + ::MIME{Symbol("text/plain")}, + ::GenericFEASolver{T,HeatTransfer,CGAssemblySolver}, +) where {T} return println(io, "TopOpt CG with assembly heat transfer solver (GenericFEASolver)") end -function Base.show(io::IO, ::MIME{Symbol("text/plain")}, ::GenericFEASolver{T,LinearElasticity,CGMatrixFreeSolver}) where {T} +function Base.show( + io::IO, + ::MIME{Symbol("text/plain")}, + ::GenericFEASolver{T,LinearElasticity,CGMatrixFreeSolver}, +) where {T} return println(io, "TopOpt matrix-free CG structural solver (GenericFEASolver)") end -function Base.show(io::IO, ::MIME{Symbol("text/plain")}, ::GenericFEASolver{T,HeatTransfer,CGMatrixFreeSolver}) where {T} +function Base.show( + io::IO, + ::MIME{Symbol("text/plain")}, + ::GenericFEASolver{T,HeatTransfer,CGMatrixFreeSolver}, +) where {T} return println(io, "TopOpt matrix-free CG heat transfer solver (GenericFEASolver)") end @@ -359,10 +448,10 @@ function FEASolver( preconditioner=identity, # Matrix-free options conv=DefaultCriteria(), - kwargs... + kwargs..., ) where {Physics<:AbstractPhysics,Solver<:AbstractLinearSolver} T = TopOptProblems.floattype(problem) - _xmin = xmin === nothing ? T(1)/1000 : T(xmin) + _xmin = xmin === nothing ? T(1) / 1000 : T(xmin) _penalty = penalty === nothing ? PowerPenalty{T}(1) : penalty _prev_penalty = prev_penalty === nothing ? deepcopy(_penalty) : prev_penalty _abstol = abstol === nothing ? T(1e-7) : T(abstol) @@ -396,13 +485,41 @@ function FEASolver( free_dofs = Int[] end - return GenericFEASolver{T,Physics,Solver,typeof(_penalty),typeof(problem), - typeof(globalinfo),typeof(elementinfo),typeof(u), - typeof(cg_max_iter),typeof(cg_statevars),typeof(preconditioner),typeof(conv)}( - problem, globalinfo, elementinfo, u, lhs, rhs, vars, - _penalty, _prev_penalty, _xmin, qr, - cg_max_iter, _abstol, cg_statevars, preconditioner, Ref(false), conv, - meandiag, fixed_dofs, free_dofs, xes + return GenericFEASolver{ + T, + Physics, + Solver, + typeof(_penalty), + typeof(problem), + typeof(globalinfo), + typeof(elementinfo), + typeof(u), + typeof(cg_max_iter), + typeof(cg_statevars), + typeof(preconditioner), + typeof(conv), + }( + problem, + globalinfo, + elementinfo, + u, + lhs, + rhs, + vars, + _penalty, + _prev_penalty, + _xmin, + qr, + cg_max_iter, + _abstol, + cg_statevars, + preconditioner, + Ref(false), + conv, + meandiag, + fixed_dofs, + free_dofs, + xes, ) end diff --git a/src/TopOptProblems/TopOptProblems.jl b/src/TopOptProblems/TopOptProblems.jl index fb90dac3..a9948f6c 100644 --- a/src/TopOptProblems/TopOptProblems.jl +++ b/src/TopOptProblems/TopOptProblems.jl @@ -55,6 +55,8 @@ export PointLoadCantilever, bcmatrix, save_mesh, RandomMagnitude, - MultiLoad + MultiLoad, + getconvectiondict, # Add your getter + assemble_convection_matrix end # module diff --git a/src/TopOptProblems/matrices_and_vectors.jl b/src/TopOptProblems/matrices_and_vectors.jl index 24419390..a4bba851 100644 --- a/src/TopOptProblems/matrices_and_vectors.jl +++ b/src/TopOptProblems/matrices_and_vectors.jl @@ -222,7 +222,8 @@ function _make_dloads(fes, problem::StiffnessTopOptProblem, facevalues) N = nnodespercell(problem) T = floattype(problem) dloads = deepcopy(fes) - eltype(dloads) <: StaticArray || throw("Expected dloads to be StaticArrays for stiffness problems.") + eltype(dloads) <: StaticArray || + throw("Expected dloads to be StaticArrays for stiffness problems.") for i in 1:length(dloads) if eltype(dloads) <: SArray dloads[i] = zero(eltype(dloads)) @@ -326,20 +327,20 @@ function _make_dloads(fes, problem::HeatTransferTopOptProblem, facevalues) for (cellid, faceid) in faceset boundary_matrix[faceid, cellid] || throw("Face $((cellid, faceid)) not on boundary.") - + fe = dloads[cellid] getcoordinates!(cell_coords, grid, cellid) reinit!(facevalues, cell_coords, faceid) - + for q_point in 1:getnquadpoints(facevalues) dΓ = getdetJdV(facevalues, q_point) for i in 1:n_basefuncs ϕ = shape_value(facevalues, q_point, i) - + # Convection load: fe[i] += ∫ ϕi * h * T_inf * dΓ # (Note: h * T_inf has the same units as heat flux q) convection_flux = h * T_inf - + if fe isa SArray fe = @set fe[i] += ϕ * convection_flux * dΓ else @@ -411,19 +412,19 @@ For thermal compliance minimization: - dJ/dx_e = -T_e^T Ke T_e · dρ_e/dx_e """ function make_Kes_and_fes( - problem::HeatTransferTopOptProblem{dim, T}, quad_order, ::Type{Val{mat_type}} -) where {dim, T, mat_type} + problem::HeatTransferTopOptProblem{dim,T}, quad_order, ::Type{Val{mat_type}} +) where {dim,T,mat_type} dh = getdh(problem) k = getk(problem) refshape = Ferrite.getrefshape(dh.field_interpolations[1]) # Shape functions for scalar field (temperature) - interpolation_space = Lagrange{dim, refshape, 1}() - quadrature_rule = QuadratureRule{dim, refshape}(quad_order) + interpolation_space = Lagrange{dim,refshape,1}() + quadrature_rule = QuadratureRule{dim,refshape}(quad_order) cellvalues = CellScalarValues(quadrature_rule, interpolation_space) facevalues = FaceScalarValues( - QuadratureRule{dim - 1, refshape}(quad_order), interpolation_space + QuadratureRule{dim - 1,refshape}(quad_order), interpolation_space ) # Calculate element conductivity matrices @@ -433,7 +434,7 @@ function make_Kes_and_fes( MatrixType, VectorType = gettypes(T, Val{mat_type}, Val{Kesize}) Kes, weights = _make_Kes_and_weights_heat( dh, - Tuple{MatrixType, VectorType}, + Tuple{MatrixType,VectorType}, Val{n_basefuncs}, Val{Kesize}, k, @@ -451,17 +452,18 @@ end # No body forces in heat transfer - weights should be zeros # Surface heat flux is computed separately via _make_dloads function _make_Kes_and_weights_heat( - dh::DofHandler{dim, N, T}, - ::Type{Tuple{MatrixType, VectorType}}, + dh::DofHandler{dim,N,T}, + ::Type{Tuple{MatrixType,VectorType}}, ::Type{Val{n_basefuncs}}, ::Type{Val{Kesize}}, k::T, quadrature_rule, cellvalues, -) where {dim, N, T, MatrixType <: StaticArray, VectorType, n_basefuncs, Kesize} - MatrixType <: SizedMatrix && throw("SizedMatrix not supported for heat transfer problems with StaticArrays.") +) where {dim,N,T,MatrixType<:StaticArray,VectorType,n_basefuncs,Kesize} + MatrixType <: SizedMatrix && + throw("SizedMatrix not supported for heat transfer problems with StaticArrays.") nel = getncells(dh.grid) - Kes = Symmetric{T, MatrixType}[] + Kes = Symmetric{T,MatrixType}[] sizehint!(Kes, nel) # No body forces in heat transfer - weights should be zeros weights = [zeros(VectorType) for i in 1:nel] @@ -496,12 +498,13 @@ using SparseArrays assemble_convection_matrix(problem::HeatConductionProblem) Builds the constant global convection matrix for Robin boundary conditions. -Evaluates the surface integral: ∫ h * N^T * N dΓ +Evaluates the surface integral `∫ h * N^T * N dΓ` over all faces defined in `problem.convectiondict`. +Returns a sparse matrix of size `(ndofs, ndofs)`. """ -function assemble_convection_matrix(problem::HeatConductionProblem{dim, T}) where {dim, T} +function assemble_convection_matrix(problem::HeatConductionProblem{dim,T}) where {dim,T} convectiondict = getconvectiondict(problem) dh = getdh(problem) - + # If no convection is defined, return an empty sparse matrix if isempty(convectiondict) return spzeros(T, ndofs(dh), ndofs(dh)) @@ -509,48 +512,48 @@ function assemble_convection_matrix(problem::HeatConductionProblem{dim, T}) wher grid = dh.grid boundary_matrix = grid.boundary_matrix - + # We need FaceScalarValues to evaluate the surface integrals - ip = dh.field_interpolations[1] - qr = Ferrite.QuadratureRule{dim - 1, Ferrite.RefCube}(2) + ip = dh.field_interpolations[1] + qr = Ferrite.QuadratureRule{dim - 1,Ferrite.RefCube}(2) facevalues = Ferrite.FaceScalarValues(qr, ip) - + n_basefuncs = getnbasefunctions(facevalues) cell_coords = zeros(Ferrite.Vec{dim,T}, nnodespercell(problem)) - + I = Int[] J = Int[] V = T[] - + global_dofs = zeros(Int, ndofs_per_cell(dh)) for (faceset_name, (h, T_inf)) in convectiondict faceset = getfacesets(problem)[faceset_name] - + for (cellid, faceid) in faceset boundary_matrix[faceid, cellid] || throw("Face not on boundary.") - + getcoordinates!(cell_coords, grid, cellid) reinit!(facevalues, cell_coords, faceid) celldofs!(global_dofs, dh, cellid) - + # Creating a local face matrix K_face = zeros(T, n_basefuncs, n_basefuncs) - + # Integrate: ∫ h * N^T * N dΓ for q_point in 1:getnquadpoints(facevalues) dΓ = getdetJdV(facevalues, q_point) - + for i in 1:n_basefuncs ϕ_i = shape_value(facevalues, q_point, i) for j in 1:n_basefuncs ϕ_j = shape_value(facevalues, q_point, j) - + K_face[i, j] += ϕ_i * ϕ_j * h * dΓ end end end - + # Scatter local face matrix to global I, J, V vectors for i in 1:n_basefuncs for j in 1:n_basefuncs @@ -561,7 +564,7 @@ function assemble_convection_matrix(problem::HeatConductionProblem{dim, T}) wher end end end - + return sparse(I, J, V, ndofs(dh), ndofs(dh)) end diff --git a/src/TopOptProblems/problem_types.jl b/src/TopOptProblems/problem_types.jl index e271bea5..f2dbd60a 100644 --- a/src/TopOptProblems/problem_types.jl +++ b/src/TopOptProblems/problem_types.jl @@ -166,9 +166,7 @@ function PointLoadCantilever( N = nnodespercell(rect_grid) M = nfacespercell(rect_grid) - return PointLoadCantilever( - rect_grid, E, ν, ch, force, force_dof, metadata - ) + return PointLoadCantilever(rect_grid, E, ν, ch, force, force_dof, metadata) end """ @@ -366,12 +364,7 @@ end - `metadata`: Metadata having various cell-node-dof relationships """ struct LBeam{ - T, - N, - M, - Tc<:ConstraintHandler{<:DofHandler{2,<:Cell{2,N,M},T},T}, - Tf<:Integer, - Tm<:Metadata, + T,N,M,Tc<:ConstraintHandler{<:DofHandler{2,<:Cell{2,N,M},T},T},Tf<:Integer,Tm<:Metadata } <: StiffnessTopOptProblem{2,T} E::T ν::T @@ -380,7 +373,9 @@ struct LBeam{ force_dof::Tf metadata::Tm end -Base.show(io::IO, ::MIME{Symbol("text/plain")}, ::LBeam) = println(io, "TopOpt L-beam problem") +function Base.show(io::IO, ::MIME{Symbol("text/plain")}, ::LBeam) + return println(io, "TopOpt L-beam problem") +end """ LBeam(::Type{Val{CellType}}, ::Type{T}=Float64; length = 100, height = 100, upperslab = 50, lowerslab = 50, E = 1.0, ν = 0.3, force = 1.0) where {T, CellType} @@ -539,11 +534,7 @@ end - `metadata`: Metadata having various cell-node-dof relationships """ struct TieBeam{ - T, - N, - M, - Tc<:ConstraintHandler{<:DofHandler{2,<:Cell{2,N,M},T},T}, - Tm<:Metadata, + T,N,M,Tc<:ConstraintHandler{<:DofHandler{2,<:Cell{2,N,M},T},T},Tm<:Metadata } <: StiffnessTopOptProblem{2,T} E::T ν::T @@ -566,12 +557,7 @@ end - `CellType`: can be either `:Linear` or `:Quadratic` to determine the order of the geometric and field basis functions and element type. Only isoparametric elements are supported for now. """ function TieBeam( - ::Type{Val{CellType}}, - (::Type{T})=Float64; - refine=1, - force=T(1), - E=T(1), - ν=T(0.3), + ::Type{Val{CellType}}, (::Type{T})=Float64; refine=1, force=T(1), E=T(1), ν=T(0.3) ) where {T,CellType} grid = TieBeamGrid(Val{CellType}, T; refine=refine) @@ -640,9 +626,13 @@ getmetadata(p::HeatTransferTopOptProblem) = p.metadata getdh(p::HeatTransferTopOptProblem) = p.ch.dh getpressuredict(p::HeatTransferTopOptProblem{dim,T}) where {dim,T} = Dict{String,T}() getheatfluxdict(p::HeatTransferTopOptProblem{dim,T}) where {dim,T} = Dict{String,T}() -getconvectiondict(p::HeatTransferTopOptProblem{dim,T}) where {dim,T} = Dict{String,Tuple{T,T}}() +function getconvectiondict(p::HeatTransferTopOptProblem{dim,T}) where {dim,T} + return Dict{String,Tuple{T,T}}() +end getfacesets(p::HeatTransferTopOptProblem) = getdh(p).grid.facesets -Ferrite.getncells(problem::HeatTransferTopOptProblem) = Ferrite.getncells(getdh(problem).grid) +function Ferrite.getncells(problem::HeatTransferTopOptProblem) + return Ferrite.getncells(getdh(problem).grid) +end getgeomorder(p::HeatTransferTopOptProblem) = nnodespercell(p) in (9, 27) ? 2 : 1 getcloaddict(p::HeatTransferTopOptProblem{dim,T}) where {dim,T} = Dict{String,Vector{T}}() @@ -691,12 +681,12 @@ struct HeatConductionProblem{ T, N, M, - Tr<:RectilinearGrid{dim, T, N, M}, - Tc<:ConstraintHandler{<:DofHandler{dim, <:Cell{dim, N, M}, T}, T}, + Tr<:RectilinearGrid{dim,T,N,M}, + Tc<:ConstraintHandler{<:DofHandler{dim,<:Cell{dim,N,M},T},T}, Th<:AbstractDict{String,T}, Tconv<:AbstractDict{String,Tuple{T,T}}, Tm<:Metadata, -} <: HeatTransferTopOptProblem{dim, T} +} <: HeatTransferTopOptProblem{dim,T} rect_grid::Tr k::T ch::Tc @@ -733,14 +723,14 @@ problem = HeatConductionProblem(Val{:Linear}, nels, sizes, k; Tleft=0.0, Tright= """ function HeatConductionProblem( ::Type{Val{CellType}}, - nels::NTuple{dim, Int}, + nels::NTuple{dim,Int}, sizes::NTuple{dim}, k=1.0; Tleft=0.0, Tright=0.0, heatflux=Dict{String,Float64}(), - convection=Dict{String,Tuple{Float64,Float64}}() -) where {dim, CellType} + convection=Dict{String,Tuple{Float64,Float64}}(), +) where {dim,CellType} T = float(promote_type(eltype(sizes), typeof(k), typeof(Tleft), typeof(Tright))) if CellType === :Linear @@ -765,7 +755,7 @@ function HeatConductionProblem( if CellType === :Linear push!(dh, :T, 1) # Temperature is a scalar field else - ip = Lagrange{dim, RefCube, 2}() + ip = Lagrange{dim,RefCube,2}() push!(dh, :T, 1, ip) end close!(dh) @@ -773,7 +763,9 @@ function HeatConductionProblem( # Apply temperature boundary conditions ch = ConstraintHandler(dh) dbc_left = Dirichlet(:T, getnodeset(rect_grid.grid, "left_boundary"), (x, t) -> Tleft) - dbc_right = Dirichlet(:T, getnodeset(rect_grid.grid, "right_boundary"), (x, t) -> Tright) + dbc_right = Dirichlet( + :T, getnodeset(rect_grid.grid, "right_boundary"), (x, t) -> Tright + ) add!(ch, dbc_left) add!(ch, dbc_right) close!(ch) @@ -787,11 +779,11 @@ function HeatConductionProblem( for (key, val) in heatflux heatfluxdict[key] = T(val) end - - convectiondict = Dict{String, Tuple{T, T}}() + + convectiondict = Dict{String,Tuple{T,T}}() for (key, val) in convection # val[1] is h, val[2] is T_ambient - convectiondict[key] = (T(val[1]), T(val[2])) + convectiondict[key] = (T(val[1]), T(val[2])) end return HeatConductionProblem( @@ -799,4 +791,4 @@ function HeatConductionProblem( ) end -nnodespercell(p::HeatConductionProblem) = nnodespercell(p.rect_grid) \ No newline at end of file +nnodespercell(p::HeatConductionProblem) = nnodespercell(p.rect_grid) diff --git a/test/topopt_problems/test_convection.jl b/test/topopt_problems/test_convection.jl new file mode 100644 index 00000000..66ae4703 --- /dev/null +++ b/test/topopt_problems/test_convection.jl @@ -0,0 +1,35 @@ +using TopOpt +using TopOpt.TopOptProblems +using Test +using LinearAlgebra +using SparseArrays + +@testset "Thermal Convection Matrix Assembly" begin + # Creates a tiny 2x2 dummy grid + nels = (2, 2) + sizes = (1.0, 1.0) + k = 1.0 + + # Applying convection to the top edge + convection = Dict{String,Tuple{Float64,Float64}}("top" => (10.0, 20.0)) + problem = TopOpt.TopOptProblems.HeatConductionProblem( + Val{:Linear}, nels, sizes, k; convection=convection + ) + + # Building the matrix using explicit path to avoid export issues + K_conv = TopOpt.TopOptProblems.assemble_convection_matrix(problem) + + # --- THE TESTS --- + # A 2x2 linear grid has 3x3 = 9 nodes. 1 DOF per node = 9x9 matrix. + @test size(K_conv) == (9, 9) + # It must be sparse for performance + @test issparse(K_conv) + # The math dictates it must be perfectly symmetric + @test issymmetric(K_conv) + # Because we applied convection, there should be non-zero elements + @test nnz(K_conv) > 0 + # Check the fallback (Empty convection should return an empty sparse matrix) + prob_empty = TopOpt.TopOptProblems.HeatConductionProblem(Val{:Linear}, nels, sizes, k) + K_empty = TopOpt.TopOptProblems.assemble_convection_matrix(prob_empty) + @test nnz(K_empty) == 0 +end