Example 6 — 3D Surface Enclosure

Last updated

30 September 2026

This example solves radiative equilibrium in a unit cube with transparent (non-participating) media. Two opposing faces have prescribed temperatures (1000 K and 0 K); the four side walls are in radiative equilibrium (unknown temperature, zero net heat flux). All surfaces are black (ε = 1). View factors are computed semi-analytically using the formulation of Narayanaswamy (2015), which means no ray tracing is needed.

Step 1: Define the cube geometry

using RayTraceHeatTransfer
using GLMakie

# Cube vertices
points = [
    0.0 0.0 0.0;  # 1
    0.0 0.0 1.0;  # 2
    0.0 1.0 0.0;  # 3
    0.0 1.0 1.0;  # 4
    1.0 0.0 0.0;  # 5
    1.0 0.0 1.0;  # 6
    1.0 1.0 0.0;  # 7
    1.0 1.0 1.0   # 8
]

# Six faces. Winding does not matter: the constructor orients all faces
# consistently and fixes the global sign from the enclosed volume, so the
# same face list works for convex and non-convex enclosures alike.
faces = [
    1 2 4 3;  # Face 1 (x = 0) — hot
    5 6 8 7;  # Face 2 (x = 1) — cold
    1 5 7 3;  # Face 3 (z = 0)
    2 6 8 4;  # Face 4 (z = 1)
    3 4 8 7;  # Face 5 (y = 1)
    1 2 6 5   # Face 6 (y = 0)
]

6×4 Matrix{Int64}:
 1  2  4  3
 5  6  8  7
 1  5  7  3
 2  6  8  4
 3  4  8  7
 1  2  6  5

true

Step 2: Set boundary conditions and mesh the domain

Ndim = 11  # 11 × 11 subdivisions per face

epsilon = ones(size(faces, 1))                     # black surfaces
q_in_w  = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]           # zero net heat flux on sides
T_in_w  = [1000.0, 0.0, -1.0, -1.0, -1.0, -1.0]    # hot, cold, four unknown

domain3D = ViewFactorDomain3D(points, faces, Ndim, q_in_w, T_in_w, epsilon); # mesh domain

Faces 1 and 2 are the hot and cold walls at opposing ends of the cube. The four side faces have T_in_w = -1.0 (unknown) and q_in_w = 0.0 (radiative equilibrium), so their temperature distributions emerge from the solution.

Step 3: Visualise the domain and identify elements

Each subface is one element: a row and column of the exchange factor matrix and one entry of the solution vectors. The mesh is drawn as a single surface built in exactly that order, so elements can be identified by hovering over them. Pass inspect = true and add a DataInspector; the hovered element is outlined and named in a tooltip, and if a Label is supplied its full property list is written there. Hovering works in an interactive Julia session; on this page the figure is a static image.

fig  = Figure(size = (1200, 700))
ax   = LScene(fig[1, 1], scenekw = (camera = cam3d!, show_axis = true))
info = Label(fig[1, 2], ""; tellheight = false, tellwidth = true,
             halign = :left, justification = :left,
             fontsize = 14, font = "DejaVu Sans Mono")
colsize!(fig.layout, 2, Relative(0.5))

plotMesh(ax, domain3D; inspect = true, label = info)
DataInspector(fig)
fig

The panel reports the element’s global index and its superface, its area and midpoint, its boundary conditions, and — once the domain has been solved — its temperature, net heat flux and radiosity. Before the solve those fields read —.

Elements are numbered superface by superface. With Ndim = 11 each quadrilateral face contributes 121 subfaces, so the hot face is elements 1–121 and the cold face 122–242. Hovering across the edge between two faces shows the index jumping between blocks. Triangular superfaces contribute blocks of a different length (Example 7), but the ordering rule is the same. These are the indices used in domain3D.F_smooth and in the system matrices.

The scene can be zoomed and entered, so elements can be inspected from inside the enclosure — which can be useful in a non-convex geometry (Example 9).

Step 4: Compute view factors

View factors are computed directly on the mesh object (requires a convex domain). Rows are normalised on assembly, so F_raw conserves energy exactly:

domain3D(; parallel = true, verbose = false)

Step 5: Smooth

smooth!(domain3D; verbose = false);  # enforce reciprocity and energy conservation

Energy conservation is exact by construction, so what smoothing recovers here is reciprocity. How much work that is depends on the geometry. The contour integral is ill-conditioned for polygons sharing an edge, and the implementation evaluates such pairs on a slightly perturbed geometry. On the cube most adjacent pairs are the coplanar ones within a single face, whose true view factor is zero and which the perturbation therefore cannot corrupt; only the subcells meeting along the twelve cube edges are affected. The raw reciprocity defect is correspondingly small, around 10⁻¹³. Example 7 shows the opposite case.

Step 6: Solve and visualise

Passing inspect = true to plotField gives the same hover on the solved field, now with every property populated:

solveEquilibrium!(domain3D, domain3D.F_smooth; verbose = false)
fig  = Figure(size = (1200, 700))
ax   = LScene(fig[1, 1], scenekw = (camera = cam3d!, show_axis = true))
info = Label(fig[1, 2], ""; tellheight = false, tellwidth = true,
             halign = :left, justification = :left,
             fontsize = 14, font = "DejaVu Sans Mono")
colsize!(fig.layout, 2, Relative(0.5))

plotField(ax, domain3D; field = :T, inspect = true, label = info)
DataInspector(fig)
fig

This is the quickest check that a solution is what it should be: on a prescribed-temperature face, T equals T_in and q is whatever flux maintains it; on an equilibrium face, T_in reads unknown and q sits at the noise floor.

The temperature field shows a smooth gradient from the hot face (1000 K) to the cold face (0 K), with the side walls at intermediate temperatures determined by radiative equilibrium. The analytical view factors ensure exact geometric accuracy without statistical noise.

Step 7: Cross-validation against ray tracing

The same enclosure can be solved by tracing rays instead of evaluating view factors analytically. Both produce an exchange factor matrix, so everything downstream is identical:

domainMC = RayTracingDomain3D_surfaces(points, faces, Ndim, q_in_w, T_in_w, epsilon)
domainMC(10^8; verbose = false)                  # total rays, split across emitters
smooth!(domainMC; verbose = false)
solveEquilibrium!(domainMC, domainMC.F_smooth; verbose = false)

T_MC = [sf.T_w for f in domainMC.facesMesh for sf in f.subFaces]
T_VF = [sf.T_w for f in domain3D.facesMesh  for sf in f.subFaces]
free = findall(sf -> sf.T_in_w < 0, [sf for f in domainMC.facesMesh for sf in f.subFaces])

println("rms |T_MC - T_VF| = ", sqrt(sum(abs2, T_MC[free] - T_VF[free]) / length(free)))

rms |T_MC - T_VF| = 0.08806408823560741

The two agree to an rms of 0.09 K on the side walls at this ray count, about 1 × 10⁻⁴ of their temperature, and the deviation keeps falling with the number of rays. That is the essential trade: view factors are machine precision and cost nothing to converge, while ray tracing pays for every digit.

Ray tracing earns its cost where view factors cannot go at all: enclosures that are not convex, where surfaces shadow one another. See Example 9.

Reference

Narayanaswamy, A. (2015). An analytic expression for radiation view factor between two arbitrarily oriented planar polygons. International Journal of Heat and Mass Transfer, 91, 841–847.

--- title: "Example 6 — 3D Surface Enclosure" execute: warning: false --- This example solves radiative equilibrium in a unit cube with transparent (non-participating) media. Two opposing faces have prescribed temperatures (1000 K and 0 K); the four side walls are in radiative equilibrium (unknown temperature, zero net heat flux). All surfaces are black (ε = 1). View factors are computed semi-analytically using the formulation of Narayanaswamy (2015), which means no ray tracing is needed. ## Step 1: Define the cube geometry ```{julia} using RayTraceHeatTransfer using GLMakie # Cube vertices points = [ 0.0 0.0 0.0; # 1 0.0 0.0 1.0; # 2 0.0 1.0 0.0; # 3 0.0 1.0 1.0; # 4 1.0 0.0 0.0; # 5 1.0 0.0 1.0; # 6 1.0 1.0 0.0; # 7 1.0 1.0 1.0 # 8 ] # Six faces. Winding does not matter: the constructor orients all faces # consistently and fixes the global sign from the enclosed volume, so the # same face list works for convex and non-convex enclosures alike. faces = [ 1 2 4 3; # Face 1 (x = 0) — hot 5 6 8 7; # Face 2 (x = 1) — cold 1 5 7 3; # Face 3 (z = 0) 2 6 8 4; # Face 4 (z = 1) 3 4 8 7; # Face 5 (y = 1) 1 2 6 5 # Face 6 (y = 0) ] ``` ```{julia} #| echo: false Makie.inline!(true) # draw figures into the page rather than a window ``` ## Step 2: Set boundary conditions and mesh the domain ```{julia} Ndim = 11 # 11 × 11 subdivisions per face epsilon = ones(size(faces, 1)) # black surfaces q_in_w = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # zero net heat flux on sides T_in_w = [1000.0, 0.0, -1.0, -1.0, -1.0, -1.0] # hot, cold, four unknown domain3D = ViewFactorDomain3D(points, faces, Ndim, q_in_w, T_in_w, epsilon); # mesh domain ``` Faces 1 and 2 are the hot and cold walls at opposing ends of the cube. The four side faces have `T_in_w = -1.0` (unknown) and `q_in_w = 0.0` (radiative equilibrium), so their temperature distributions emerge from the solution. ## Step 3: Visualise the domain and identify elements Each subface is one element: a row and column of the exchange factor matrix and one entry of the solution vectors. The mesh is drawn as a single surface built in exactly that order, so elements can be identified by hovering over them. Pass `inspect = true` and add a `DataInspector`; the hovered element is outlined and named in a tooltip, and if a `Label` is supplied its full property list is written there. Hovering works in an interactive Julia session; on this page the figure is a static image. ```{julia} fig = Figure(size = (1200, 700)) ax = LScene(fig[1, 1], scenekw = (camera = cam3d!, show_axis = true)) info = Label(fig[1, 2], ""; tellheight = false, tellwidth = true, halign = :left, justification = :left, fontsize = 14, font = "DejaVu Sans Mono") colsize!(fig.layout, 2, Relative(0.5)) plotMesh(ax, domain3D; inspect = true, label = info) DataInspector(fig) fig ``` The panel reports the element's global index and its superface, its area and midpoint, its boundary conditions, and — once the domain has been solved — its temperature, net heat flux and radiosity. Before the solve those fields read `—`. Elements are numbered superface by superface. With `Ndim = 11` each quadrilateral face contributes 121 subfaces, so the hot face is elements 1–121 and the cold face 122–242. Hovering across the edge between two faces shows the index jumping between blocks. Triangular superfaces contribute blocks of a different length ([Example 7](example-07.qmd)), but the ordering rule is the same. These are the indices used in `domain3D.F_smooth` and in the system matrices. The scene can be zoomed and entered, so elements can be inspected from inside the enclosure — which can be useful in a non-convex geometry ([Example 9](example-09.qmd)). ## Step 4: Compute view factors View factors are computed directly on the mesh object (requires a convex domain). Rows are normalised on assembly, so `F_raw` conserves energy exactly: ```{julia} domain3D(; parallel = true, verbose = false) ``` ## Step 5: Smooth ```{julia} smooth!(domain3D; verbose = false); # enforce reciprocity and energy conservation ``` Energy conservation is exact by construction, so what smoothing recovers here is reciprocity. How much work that is depends on the geometry. The contour integral is ill-conditioned for polygons sharing an edge, and the implementation evaluates such pairs on a slightly perturbed geometry. On the cube most adjacent pairs are the coplanar ones within a single face, whose true view factor is zero and which the perturbation therefore cannot corrupt; only the subcells meeting along the twelve cube edges are affected. The raw reciprocity defect is correspondingly small, around 10⁻¹³. [Example 7](example-07.qmd) shows the opposite case. ## Step 6: Solve and visualise Passing `inspect = true` to `plotField` gives the same hover on the solved field, now with every property populated: ```{julia} solveEquilibrium!(domain3D, domain3D.F_smooth; verbose = false) fig = Figure(size = (1200, 700)) ax = LScene(fig[1, 1], scenekw = (camera = cam3d!, show_axis = true)) info = Label(fig[1, 2], ""; tellheight = false, tellwidth = true, halign = :left, justification = :left, fontsize = 14, font = "DejaVu Sans Mono") colsize!(fig.layout, 2, Relative(0.5)) plotField(ax, domain3D; field = :T, inspect = true, label = info) DataInspector(fig) fig ``` This is the quickest check that a solution is what it should be: on a prescribed-temperature face, `T` equals `T_in` and `q` is whatever flux maintains it; on an equilibrium face, `T_in` reads `unknown` and `q` sits at the noise floor. The temperature field shows a smooth gradient from the hot face (1000 K) to the cold face (0 K), with the side walls at intermediate temperatures determined by radiative equilibrium. The analytical view factors ensure exact geometric accuracy without statistical noise. ## Step 7: Cross-validation against ray tracing The same enclosure can be solved by tracing rays instead of evaluating view factors analytically. Both produce an exchange factor matrix, so everything downstream is identical: ```{julia} domainMC = RayTracingDomain3D_surfaces(points, faces, Ndim, q_in_w, T_in_w, epsilon) domainMC(10^8; verbose = false) # total rays, split across emitters smooth!(domainMC; verbose = false) solveEquilibrium!(domainMC, domainMC.F_smooth; verbose = false) T_MC = [sf.T_w for f in domainMC.facesMesh for sf in f.subFaces] T_VF = [sf.T_w for f in domain3D.facesMesh for sf in f.subFaces] free = findall(sf -> sf.T_in_w < 0, [sf for f in domainMC.facesMesh for sf in f.subFaces]) println("rms |T_MC - T_VF| = ", sqrt(sum(abs2, T_MC[free] - T_VF[free]) / length(free))) ``` The two agree to an rms of 0.09 K on the side walls at this ray count, about 1 × 10⁻⁴ of their temperature, and the deviation keeps falling with the number of rays. That is the essential trade: view factors are machine precision and cost nothing to converge, while ray tracing pays for every digit. Ray tracing earns its cost where view factors cannot go at all: enclosures that are not convex, where surfaces shadow one another. See [Example 9](example-09.qmd). ## Reference Narayanaswamy, A. (2015). An analytic expression for radiation view factor between two arbitrarily oriented planar polygons. *International Journal of Heat and Mass Transfer*, 91, 841–847.