Theory
View factor definition
The view factor $F_{ij}$ from surface $i$ to surface $j$ is the fraction of diffuse radiation leaving surface $i$ that arrives at surface $j$:
3D (surface meshes):
\[F_{ij} = \frac{1}{A_i} \iint_{A_i} \iint_{A_j} \frac{\cos\theta_i \cos\theta_j}{\pi r^2} \, H_{ij} \, dA_j \, dA_i\]
2D (curve meshes):
\[F_{ij} = \frac{1}{L_i} \int_{L_i} \int_{L_j} \frac{\cos\theta_i \cos\theta_j}{2r} \, H_{ij} \, dL_j \, dL_i\]
where:
- $\theta_i$ is the angle between the outward normal at $dA_i$ and the line of sight $r_{ij}$
- $\theta_j$ is the angle between the outward normal at $dA_j$ and the reverse line of sight
- $r$ is the distance between the two differential elements
- $H_{ij} \in \{0, 1\}$ is the visibility function (0 = obstructed)
The factor $2$ rather than $\pi$ in the 2D kernel follows from integrating the 2D radiation intensity over the hemisphere, which gives $\pi/2$ rather than $\pi$.
Reciprocity
The reciprocity relation:
\[A_i F_{ij} = A_j F_{ji}\]
is a consequence of the symmetry of the kernel and holds exactly in the continuous case. At the discrete (element) level, reciprocity holds to within quadrature error; verify with check_reciprocity.
Gauss–Legendre quadrature
The isoparametric map $\mathbf{x}(\xi, \eta)$ transforms the reference element $[-1,1]^2$ to physical space. The integral over one Quad8 element becomes:
\[\int_{A} f \, dA = \int_{-1}^{1} \int_{-1}^{1} f(\mathbf{x}(\xi,\eta)) \left|\frac{\partial\mathbf{x}}{\partial\xi} \times \frac{\partial\mathbf{x}}{\partial\eta}\right| d\xi \, d\eta \approx \sum_{p=1}^{n^2} w_p \, f(\mathbf{x}(\xi_p,\eta_p)) \, J_p\]
where $J_p = |\partial_\xi \mathbf{x} \times \partial_\eta \mathbf{x}|$ is the area Jacobian. A tensor-product $n \times n$ Gauss–Legendre rule is used on quadrilateral elements (Quad4, Quad8) and an $n$-point rule on curve elements (Line2, Line3). Triangles (Tri3, Tri6) use a Dunavant rule on the reference triangle instead, with 1, 3, 7, or 13 points for $n$ = 1, 2, 3, or $\geq 4$. Rules for $n > 5$ are generated by the Golub–Welsch algorithm and memoized.
The double integral over an element pair is then evaluated by applying the kernel at every pair of quadrature points, and the physical-group view factor is the area-weighted sum of element rows, $A_g F_{g \to h} = \sum_{i \in g} A_i \sum_{j \in h} F_{ij}$.
Duffy transformation
For element pairs sharing a vertex at $\mathbf{u}_0 = (u_0, v_0)$ in element $i$'s unit square and $\mathbf{s}_0 = (s_0, t_0)$ in element $j$'s, the kernel diverges as $r \to 0$. Near the singularity:
\[K \sim \frac{1}{r^2}, \quad r \sim \sqrt{(u-u_0)^2+(v-v_0)^2+(s-s_0)^2+(t-t_0)^2}\]
The implementation uses an elementary "biggest-coordinate" Duffy decomposition (a generalization of Duffy's original vertex-singularity transformation), not the specific region formulas of Sauter & Schwab's boundary element method — a related but different decomposition of the same singularity, and still the reference for the underlying idea.
Common vertex (one shared corner node): shift coordinates to $\tilde u = u-u_0$, etc., so the singularity sits at the origin of $[0,1]^4$. Split the hypercube into 4 regions by which of $(\tilde u, \tilde v, \tilde s, \tilde t)$ has the largest magnitude; in each region that coordinate is set to a radial variable $\rho \in [0,1]$ and the other three to $\rho \eta_k$ ($\eta_k \in [0,1]$), giving a lower-triangular Jacobian with determinant $\rho^3$:
\[K \cdot |\text{Jac}| \sim \frac{1}{\rho^2} \cdot \rho^3 = \rho \to 0 \quad \text{as } \rho \to 0\]
The 4 regions exactly tile $[0,1]^4$ (ties between coordinates have measure zero); each is evaluated with an $nquad^4$ tensor-product rule, so a common-vertex pair costs $4 \times nquad^4$ evaluation points in total.
Common edge (two shared corner nodes): in edge-local coordinates the singularity is the line $u=s,\, v=t=0$ rather than a single point. Split the $(u,s)$ square into the two triangles $u \geq s$ and $u < s$ (Jacobian $u$ or $s$ respectively, a standard ratio parametrization); within each triangle the remaining three-variable point singularity $(w, v, t) \to 0$ (where $w = (u-s)/u$ or $(s-u)/s$) is resolved the same "biggest-coordinate" way, now with 3 regions and Jacobian $\rho^2$. This gives 6 regions total (2 triangles × 3 sub-regions), with combined Jacobian $(u\text{ or }s) \cdot \rho^2$; each is evaluated with an $nquad^4$ rule, so a common-edge pair costs $6 \times nquad^4$ evaluation points in total.
The transformed integrand is smooth at $\rho = 0$ in both cases and is efficiently resolved by standard Gauss–Legendre quadrature.
Monte Carlo estimator
The unbiased area-sampling estimator for each element pair $(i,j)$ is:
\[\iint K \, dA_j \, dA_i \approx \frac{A_i \cdot A_j}{N} \sum_{k=1}^{N} K\!\left(x_i^{(k)}, n_i^{(k)}, x_j^{(k)}, n_j^{(k)}\right) \cdot H_{ij}^{(k)}\]
where $(x_i^{(k)}, x_j^{(k)})$ are i.i.d. uniform samples on $A_i \times A_j$.
Stratified sampling subdivides the reference square into $\lfloor\sqrt{N}\rfloor \times \lfloor\sqrt{N}\rfloor$ strata and draws one point per stratum. For smooth integrands, stratified sampling achieves $O(1/N)$ variance convergence rather than $O(1/\sqrt{N})$ for plain Monte Carlo, matching the rate of a 1D Gauss rule.
Variance near singularities: for the $1/r^2$ kernel the variance of the MC estimator is proportional to $\iint K^2 \, dA$, which diverges at shared edges. Infinite variance means no amount of increasing $N$ gives reliable convergence — use the Duffy transformation instead for such pairs. The pair-area Monte Carlo path does this automatically: every pair within factor element diameters of each other (a superset of the touching pairs) is replaced by the deterministic Duffy/quadrature value.
Ray-shooting Monte Carlo estimator
A second, unrelated Monte Carlo estimator (raytrace=true) rewrites the double-area integral as a single-area integral over the (occlusion-limited) solid angle $\Omega_j$ that surface $j$ subtends from each point $x \in A_i$:
\[F_{i \to j} = \frac{1}{A_i} \int_{A_i} \left[ \int_{\Omega_j \text{ visible from } x} \frac{\cos\theta_i}{\pi} \, d\Omega \right] dA_i\]
which follows from the solid-angle identity $dA_j \cos\theta_j / r^2 = d\Omega$. Drawing a ray direction $\omega$ from the cosine-weighted hemisphere pdf $\cos\theta_i / \pi$ gives
\[\mathbb{E}\big[\mathbb{1}(\text{ray's first hit is } j)\big] = \int_{\Omega_j \text{ visible}} \frac{\cos\theta_i}{\pi} \, d\Omega = F_{x \to j}\]
— the hit indicator's expectation is exactly the point-to-area form factor, occlusion included, with no $\cos\theta_j$, no $1/r^2$, and no separate visibility test: a ray hitting the wrong surface first is the visibility test. Averaging over points $x$ drawn uniformly by area on $A_i$ gives $F_{i \to j}$ directly.
Because this estimator never evaluates $1/r^2$, it has ordinary bounded (binomial) variance for every element pair, including adjacent ones — no Duffy patch is applied. Each unordered pair $\{i,j\}$ gets two independent estimates (one from each element's own rays); Assembly.jl averages them before dividing back out to $F$, which enforces reciprocity by construction rather than as an emergent property of enough samples. Element areas $A_i$ come from quadrature, and no self view factor is estimated (a ray is never tested against its own origin element). See src/RayTraceKernel.jl's module docstring for the full derivation and a front/back-face subtlety found while implementing it.
Closure enforcement
Net radiative flux is a difference of radiosities, $q_i = J_i - \sum_j F_{ij} J_j$, and in a nearly isothermal enclosure it is orders of magnitude smaller than $J$ itself. A row-sum error $\epsilon$ therefore leaks $\epsilon J$ straight into $q$. Sampling cannot remove it (Monte Carlo closure error falls only as $1/\sqrt{n_\text{rays}}$), and quadrature on curved elements plateaus at a small but nonzero level.
enforce_closure instead imposes the two identities the view factors must satisfy anyway, on a closed enclosure. It alternates
\[M \leftarrow \tfrac{1}{2}\,(M + M^{\mathsf T}),\quad M_{ij} = A_i F_{ij}; \qquad F_{ij} \leftarrow \frac{M_{ij}/A_i}{\sum_k M_{ik}/A_i}\]
— symmetrising $A_i F_{ij}$ (reciprocity) and rescaling each row to sum to 1 (closure). Each step slightly perturbs the other, so the pair is iterated until both residuals fall below tol (or iters is reached). With exact row sums an isothermal enclosure gives $q = 0$ identically, so only genuine temperature differences drive a flux. On an open enclosure, where rows legitimately sum to less than 1, it would invent radiation and must not be used.