V0.6.0

.github/workflows/CI.yml

@@ -1,18 +1,41 @@
 name: CI
-
 on:
   push:
-    branches: [main]
+    branches:
+      - main
+    tags: ['*']
   pull_request:
-    branches: [main]
-
+  workflow_dispatch:
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
 jobs:
   test:
-    runs-on: ubuntu-latest
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 60
+    permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created
+      actions: write
+      contents: read
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.10'
+          - '1.12'
+          - 'pre'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
     steps:
-      - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v2
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v3
         with:
-          version: '1'
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v2
+      - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1

Project.toml

@@ -1,6 +1,6 @@
 name = "RadiativeViewFactor"
 uuid = "27b7fd98-66cc-11f1-8ce2-791900f7bf38"
-version = "0.5.0"
+version = "0.6.0"
 authors = ["Alex Coxe <rot4te@gmail.com>"]
 
 [deps]
@@ -33,5 +33,5 @@ Plots              = "1"
 Metal              = "1.9.3"
 ReadVTK            = "0.2"
 StaticArrays       = "1"
-Test               = "1.11.0"
+Test               = "1"
 julia              = "1.9"

README.md

@@ -63,6 +63,11 @@ RadiativeViewFactor.jl/
 │   ├── RadiativeViewFactorMetalExt.jl    # Registers MetalBackend → MtlArray, Float32
 │   ├── RadiativeViewFactorPlotsExt.jl    # plot_mesh_normals (Plots.jl)
 │   └── RadiativeViewFactorReadVTKExt.jl  # XML VTK (.vtu) loading via ReadVTK.jl
+├── benchmarks/
+│   ├── common.jl                # Shared mesh generators and timing helpers
+│   ├── quadrature_bench.jl      # Deterministic assembly benchmark (sweeps N)
+│   ├── montecarlo_bench.jl      # Monte Carlo assembly benchmark (sweeps n_samples)
+│   └── RESULTS.md               # Before/after numbers for the pre-evaluation optimization
 ├── test/
 │   └── runtests.jl
 └── Project.toml
@@ -269,6 +274,19 @@ Physical Curve("obstruction") = {3};
 
 ## Performance Notes
 
+### Assembly cost and quadrature reuse
+
+Assembly is O(N²) in the element count. Each element's quadrature points and
+geometric quantities (deterministic path) or Monte Carlo samples (MC path) are
+**pre-evaluated once per element** and reused across every pair, rather than
+re-derived inside the pair loop. This keeps per-pair work to the kernel
+evaluation itself and avoids O(N²) shape-function and quadrature-rule
+reconstruction. On an 8-core CPU this is ~2.6–3.1× faster for the deterministic
+path and ~12–15× faster for Monte Carlo, with 60–450× fewer allocations,
+versus re-evaluating per pair — see [`benchmarks/RESULTS.md`](benchmarks/RESULTS.md)
+and the reproducible scripts in [`benchmarks/`](benchmarks/). Results are
+numerically identical (the change is evaluation order only).
+
 ### Integration method selection
 
 | Method | Best for | Avoid when |

benchmarks/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+RadiativeViewFactor = "27b7fd98-66cc-11f1-8ce2-791900f7bf38"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

benchmarks/RESULTS.md

@@ -0,0 +1,89 @@
+# Benchmark Results — Quadrature-Point / Sample Pre-evaluation
+
+This directory benchmarks the CPU assembly path before and after the
+**pre-evaluation** optimization, in which each element's quadrature points
+(deterministic path) or Monte Carlo samples (MC path) are computed **once**
+per element — O(N) — and reused across every element pair, instead of being
+re-derived inside the O(N²) pair loop.
+
+## What changed
+
+The assembly matrix is dense and O(N²) in the element count `N`. Previously,
+`element_pair_view_factor(coords, elem_i, elem_j, …)` re-evaluated the shape
+functions, the quadrature rule, and the physical points/normals for **both**
+elements on **every** call. Element *i*'s data depends only on *i*, yet it was
+rebuilt for all `N − i` partners; for `nquad > 5` this even re-ran a
+Golub–Welsch eigensolve per pair.
+
+The optimization:
+
+- **Deterministic path** — `precompute_quad` builds one `ElementQuad` per
+  element; the pair integrator consumes cached points
+  (`src/ViewFactorKernel.jl`, `src/Assembly.jl`).
+- **Monte Carlo path** — `sample_element_mc` draws one `ElementSamples` set per
+  element, reused across the row/column; the diagonal self-pair draws a fresh
+  independent set so `self_vf` stays correct (`src/MCKernel.jl`,
+  `src/Assembly.jl`). Each per-entry estimate stays unbiased.
+- **Rule memoization** — `gauss_legendre_1d` caches the Golub–Welsch rule for
+  `nquad > 5` behind a lock (`src/Quadrature.jl`).
+
+Because the math is unchanged (only the evaluation order), results are
+numerically identical; reciprocity holds to machine precision and the
+parallel-plate view factors match the analytic value.
+
+## How to reproduce
+
+```bash
+julia --project=benchmarks -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'
+julia --project=benchmarks --threads=auto benchmarks/quadrature_bench.jl
+julia --project=benchmarks --threads=auto benchmarks/montecarlo_bench.jl
+```
+
+The "before" columns below were produced by running the same scripts against
+the pre-optimization sources (`git stash` of the `src/` changes).
+
+## Environment
+
+| | |
+|---|---|
+| CPU | Apple M1 (8 logical cores) |
+| Threads | 8 (`--threads=auto`) |
+| Julia | 1.12.6 |
+| Timing | minimum of 3 runs (warm) |
+
+## Deterministic quadrature — two facing unit plates, `nquad=6`
+
+Analytic F(bottom→top) ≈ 0.19982; every configuration reproduced it to a
+relative error of 2.4×10⁻⁵.
+
+| N (elements) | before (s) | after (s) | speedup | before alloc | after alloc | alloc ↓ |
+|---:|---:|---:|---:|---:|---:|---:|
+| 240  | 0.0100 | 0.0032 | 3.1× | 75.9 MiB   | 1.2 MiB  | 63× |
+| 396  | 0.0243 | 0.0089 | 2.7× | 206.3 MiB  | 2.9 MiB  | 71× |
+| 692  | 0.0708 | 0.0265 | 2.7× | 629.3 MiB  | 8.2 MiB  | 77× |
+| 1088 | 0.1716 | 0.0658 | 2.6× | 1554.8 MiB | 19.5 MiB | 80× |
+
+After the change, allocations are dominated by the two `N×N` result matrices
+(constant per `N`) rather than per-pair temporaries.
+
+## Monte Carlo — two facing unit plates, N = 484
+
+| n_samples | before (s) | after (s) | speedup | before alloc | after alloc | alloc ↓ |
+|---:|---:|---:|---:|---:|---:|---:|
+| 1000 | 1.189 | 0.096 | 12.4× | 16.2 GiB | 47.4 MiB  | 350× |
+| 2000 | 2.475 | 0.205 | 12.1× | 32.3 GiB | 81.4 MiB  | 406× |
+| 5000 | 6.865 | 0.464 | 14.8× | 86.0 GiB | 194.8 MiB | 452× |
+
+The MC path re-sampled both elements on every pair, so the redundancy — and
+therefore the speedup and allocation reduction — is even larger than for the
+deterministic path. Estimated view factors matched the analytic value to a
+relative error ≈ 1.5×10⁻⁵ across sample counts.
+
+## Takeaways
+
+- The optimization is a pure implementation change: identical results, no new
+  approximations.
+- Speedups are **~2.6–3.1× (deterministic)** and **~12–15× (Monte Carlo)** with
+  **60–450× fewer allocations**, growing with `N` and `n_samples`.
+- The remaining cost is the genuine O(N²) kernel work plus the dense result
+  matrices; both are inherent to full-matrix assembly.

benchmarks/common.jl

@@ -0,0 +1,89 @@
+# benchmarks/common.jl
+# ---------------------------------------------------------------------------
+# Shared mesh generators and helpers for the benchmark scripts.
+# Requires Gmsh (see benchmarks/Project.toml).
+# ---------------------------------------------------------------------------
+
+import Gmsh: gmsh
+
+"""
+    make_box_msh(path; lc=0.09, order=1)
+
+Write a closed unit-cube surface mesh (6 physical surfaces, one group) to
+`path`. Smaller `lc` → more elements. `order=1` gives Quad4/Tri3, `order=2`
+gives Quad8/Tri6.
+"""
+function make_box_msh(path; lc=0.09, order=1)
+    gmsh.initialize(); gmsh.option.setNumber("General.Verbosity", 0)
+    gmsh.model.add("box")
+    gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)
+    gmsh.model.occ.synchronize()
+    ptag = gmsh.model.addPhysicalGroup(2, collect(1:6))
+    gmsh.model.setPhysicalName(2, ptag, "box")
+    gmsh.option.setNumber("Mesh.MeshSizeMax", lc)
+    gmsh.option.setNumber("Mesh.ElementOrder", order)
+    gmsh.model.mesh.generate(2)
+    gmsh.write(path)
+    gmsh.finalize()
+    return path
+end
+
+"""
+    make_two_plates_msh(path; lc=0.1, order=1)
+
+Two coaxial unit-square plates a distance 1 apart, wound to face each other so
+the analytic bottom→top view factor is ≈ 0.19982. Physical groups "bottom" and
+"top". Uses Gmsh's default unstructured triangulation.
+"""
+function make_two_plates_msh(path; lc=0.1, order=1)
+    gmsh.initialize(); gmsh.option.setNumber("General.Verbosity", 0)
+    gmsh.model.add("plates")
+    # Bottom plate (z=0), outward normal +z (toward the top plate)
+    gmsh.model.geo.addPoint(0,0,0, lc, 1); gmsh.model.geo.addPoint(1,0,0, lc, 2)
+    gmsh.model.geo.addPoint(1,1,0, lc, 3); gmsh.model.geo.addPoint(0,1,0, lc, 4)
+    for (i,(a,b)) in enumerate([(1,2),(2,3),(3,4),(4,1)]); gmsh.model.geo.addLine(a,b,i); end
+    gmsh.model.geo.addCurveLoop([1,2,3,4], 1); gmsh.model.geo.addPlaneSurface([1], 1)
+    # Top plate (z=1), wound clockwise so its outward normal points -z (downward)
+    gmsh.model.geo.addPoint(0,0,1, lc, 5); gmsh.model.geo.addPoint(1,0,1, lc, 6)
+    gmsh.model.geo.addPoint(1,1,1, lc, 7); gmsh.model.geo.addPoint(0,1,1, lc, 8)
+    for (i,(a,b)) in enumerate([(5,6),(6,7),(7,8),(8,5)]); gmsh.model.geo.addLine(a,b,i+4); end
+    gmsh.model.geo.addCurveLoop([-8,-7,-6,-5], 2); gmsh.model.geo.addPlaneSurface([2], 2)
+    gmsh.model.geo.synchronize()
+    gmsh.model.setPhysicalName(2, gmsh.model.addPhysicalGroup(2, [1]), "bottom")
+    gmsh.model.setPhysicalName(2, gmsh.model.addPhysicalGroup(2, [2]), "top")
+    gmsh.option.setNumber("Mesh.ElementOrder", order)
+    gmsh.model.mesh.generate(2)
+    gmsh.write(path)
+    gmsh.finalize()
+    return path
+end
+
+"""
+    best_time(f; samples=3)
+
+Run `f()` once to warm up, then return the minimum wall-clock time (seconds)
+and the allocated bytes over `samples` runs. Minimum time is the standard
+choice for microbenchmarks — it is the run least perturbed by the OS / GC.
+"""
+function best_time(f; samples=3)
+    f()                                   # warm up / compile
+    GC.gc()
+    t = Inf; a = 0
+    for _ in 1:samples
+        GC.gc()
+        b = @allocated f()
+        s = @elapsed f()
+        t = min(t, s)
+        a = b
+    end
+    return t, a
+end
+
+"Print the machine / thread configuration used for a benchmark run."
+function print_env()
+    println("─"^64)
+    println("CPU     : ", Sys.CPU_NAME, "  (", Sys.CPU_THREADS, " logical cores)")
+    println("Machine : ", Sys.MACHINE)
+    println("Julia   : ", VERSION, "   threads=", Threads.nthreads())
+    println("─"^64)
+end

benchmarks/montecarlo_bench.jl

@@ -0,0 +1,45 @@
+# benchmarks/montecarlo_bench.jl
+# ---------------------------------------------------------------------------
+# Benchmark the stratified Monte Carlo assembly path.
+#
+# Measures full `compute_view_factors(monte_carlo=true)` wall-clock time and
+# allocations on two facing unit plates, sweeping `n_samples`, and reports the
+# estimated view factor against the analytic value (≈ 0.19982) so the accuracy
+# of the sample-reuse optimization stays visible alongside the timing.
+#
+# Run with threads:
+#     julia --project=benchmarks --threads=auto benchmarks/montecarlo_bench.jl
+# ---------------------------------------------------------------------------
+
+using RadiativeViewFactor
+using Random
+using Printf
+include(joinpath(@__DIR__, "common.jl"))
+
+print_env()
+
+const ANALYTIC = 0.19982   # two directly-opposed unit squares, separation 1
+
+f = tempname() * ".msh"
+make_two_plates_msh(f; lc=0.1, order=1)
+mesh = load_mesh(f; verbose=false)
+N = length(mesh.surface_elems)
+println("Monte Carlo assembly — two facing unit plates, N = $N elements")
+println("Analytic F(bottom→top) ≈ $ANALYTIC\n")
+
+bi(r) = findfirst(==("bottom"), r.group_names)
+ti(r) = findfirst(==("top"),    r.group_names)
+
+@printf("%-11s %12s %14s %12s %12s\n",
+        "n_samples", "time (s)", "alloc (MiB)", "F_bot→top", "rel.err")
+for ns in (1_000, 2_000, 5_000, 20_000)
+    run() = compute_view_factors(mesh; monte_carlo=true, n_samples=ns,
+                                 rng=MersenneTwister(1), verbose=false)
+    t, a = best_time(run)
+    r    = run()
+    F    = r.F_group[bi(r), ti(r)]
+    @printf("%-11d %12.4f %14.1f %12.6f %12.2e\n",
+            ns, t, a/2^20, F, abs(F - ANALYTIC)/ANALYTIC)
+end
+
+rm(f)

benchmarks/quadrature_bench.jl

@@ -0,0 +1,43 @@
+# benchmarks/quadrature_bench.jl
+# ---------------------------------------------------------------------------
+# Benchmark the deterministic Gauss–Legendre / Dunavant assembly path.
+#
+# The assembly is O(N²) in the element count N, and the optimization being
+# measured pre-evaluates each element's quadrature once (O(N)) instead of
+# re-deriving it inside every pair. So we sweep the mesh refinement (N) at a
+# fixed `nquad`, reporting wall-clock time, allocations, and the estimated
+# view factor against the analytic value (≈ 0.19982) to keep accuracy visible.
+#
+# Run with threads for a representative number:
+#     julia --project=benchmarks --threads=auto benchmarks/quadrature_bench.jl
+# ---------------------------------------------------------------------------
+
+using RadiativeViewFactor
+using Printf
+include(joinpath(@__DIR__, "common.jl"))
+
+print_env()
+
+const ANALYTIC = 0.19982   # two directly-opposed unit squares, separation 1
+const NQUAD    = 6
+
+println("Quadrature assembly — two facing unit plates, nquad = $NQUAD")
+println("Analytic F(bottom→top) ≈ $ANALYTIC\n")
+
+bi(r) = findfirst(==("bottom"), r.group_names)
+ti(r) = findfirst(==("top"),    r.group_names)
+
+@printf("%-8s %10s %12s %14s %12s %12s\n",
+        "lc", "N", "time (s)", "alloc (MiB)", "F_bot→top", "rel.err")
+for lc in (0.16, 0.12, 0.09, 0.07)
+    f = tempname() * ".msh"
+    make_two_plates_msh(f; lc=lc, order=1)
+    mesh = load_mesh(f; verbose=false)
+    N = length(mesh.surface_elems)
+    t, a = best_time(() -> compute_view_factors(mesh; nquad=NQUAD, verbose=false))
+    r    = compute_view_factors(mesh; nquad=NQUAD, verbose=false)
+    F    = r.F_group[bi(r), ti(r)]
+    @printf("%-8.2f %10d %12.4f %14.1f %12.6f %12.2e\n",
+            lc, N, t, a/2^20, F, abs(F - ANALYTIC)/ANALYTIC)
+    rm(f)
+end

docs/src/manual/integration_methods.md

@@ -57,10 +57,18 @@ standard quadrature for inclined-plate geometries.
 result = compute_view_factors(mesh; monte_carlo=true, n_samples=50000)
 ```
 
-Each element pair draws `n_samples` stratified random sample pairs. Samples are
-placed on a ⌊√N⌋ × ⌊√N⌋ grid of strata within the reference element, giving
+Each element pair is estimated from `n_samples` stratified sample pairs. Samples
+are placed on a ⌊√N⌋ × ⌊√N⌋ grid of strata within the reference element, giving
 O(1/N) variance convergence rather than O(1/√N) for plain Monte Carlo.
 
+For efficiency, one independent stratified sample set is drawn **once per
+element** and reused across that element's pairings (the diagonal self-pair,
+with `self_vf=true`, draws a fresh second set). Within any pair the two
+elements' samples are still independent, so each entry's estimate remains
+unbiased with the stated variance; only estimates in the same row/column become
+correlated. This reuse is what makes the MC path ~12–15× faster than
+re-sampling both elements per pair — see the `benchmarks/` directory.
+
 **When to use:**
 - Many obstructions (MC pays the BVH cost only for kernel-positive pairs)
 - Near-singular pairs where MC variance is still finite (unlike the `1/r²`