feat: implement Robin boundary conditions (convection) for thermal optimization

Project.toml

@@ -1,49 +1,11 @@
-authors = ["Mohamed Tarek <mohdibntarek@gmail.com>", "Yijiang Huang <yijiang94817@gmail.com>"]
 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 <mohdibntarek@gmail.com>", "Yijiang Huang <yijiang94817@gmail.com>"]
 
 [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"

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