Refactor/moller trumbore ray triangle

[sbryngelson]

Summary

Refactors the STL ray-triangle intersection and inside/outside classification in m_model.fpp:

• Moller-Trumbore algorithm replaces the old cross-product sign test in f_intersects_triangle. The old approach was winding-order dependent. Moller-Trumbore computes barycentric coordinates directly from vertices, making it vertex-order agnostic.

• Generalized winding number replaces the 26-ray casting approach in f_model_is_inside_flat. Ray casting fails on small triangles because rays can miss or edge-graze tiny faces, corrupting the parity count. The winding number sums the solid angle (Van Oosterom-Strackee formula) subtended by each triangle. This is the standard approach in geometry processing (Jacobson et al., SIGGRAPH 2013).

• 2D branch uses the 2D winding number (signed angle sum over boundary edges) for p==0 simulations, since the 3D solid angle degenerates when all triangles are coplanar.

Known limitations

• Non-watertight or inconsistently-oriented meshes will produce fractional winding numbers. The threshold classification still works for mostly-consistent meshes.

[github-actions]

Claude Code Review

Head SHA: a1d3d9e
Files changed: 1 — src/common/m_model.fpp
+67 / -74

---

Summary

• Replaces multi-ray-casting point-in-surface test (f_model_is_inside_flat) with the generalized winding number method (Van Oosterom–Strackee, 1983), citing Jacobson et al. 2013.
• Replaces the old geometric (plane/sign) triangle–intersection test with the standard Möller–Trumbore barycentric-coordinate algorithm in f_intersects_triangle.
• Both algorithmic changes are well-known, numerically more robust than the previous approaches, and the math checks out.
• No GPU macro violations, no forbidden Fortran patterns, no precision-type issues observed.

---

Findings

1. 2D simulation case — winding number may degenerate (medium concern)
File: src/common/m_model.fpp, removed block (~line 613–634 in original)

The old code explicitly skipped z-direction rays when p == 0 (2D simulation) and used different normalization denominators (18 vs. 26). The new winding-number loop has no special p == 0 branch. For 2D simulations whose STL triangles are all coplanar (z-extent ≈ 0), the scalar triple product numerator will be near-zero for every triangle, so fraction will be ≈ 0 everywhere regardless of whether the query point is inside or outside. Please verify that at least one 2D STL test case still produces correct results, or document that winding-number mode requires a closed 3D surface.

2. pi computed via acos(-1.0_wp) — minor style note
File: src/common/m_model.fpp, new line (bottom of f_model_is_inside_flat):

fraction = fraction/(2.0_wp*acos(-1.0_wp))

If the codebase already defines a named pi or 4.0_wp*atan(1.0_wp) constant, prefer it for consistency. Otherwise this is fine; acos(-1.0_wp) is correct and portable.

3. Möller–Trumbore epsilon threshold — unchanged behaviour, noted for clarity
File: src/common/m_model.fpp, f_intersects_triangle (~line 671):

if (abs(a) < 1e-7_wp) return

The old implementation used 0.0000001_wp (same value). No regression here, just confirming the literal is written correctly with the _wp kind suffix.

4. triangle%n (precomputed normal) is now unused in f_intersects_triangle
File: src/common/m_model.fpp, updated f_intersects_triangle

Möller–Trumbore works from vertices alone and correctly drops the normal. However, check whether the t_triangle%n field is still populated/used elsewhere in the file (e.g. in f_model_is_inside or related routines), and whether gpu_trs_n still needs to be stored and uploaded to the GPU if no code path reads it anymore.

---

No other issues found

The winding-number normalization is mathematically correct: each triangle contributes 2*atan2(num,den) to the solid angle, so summing atan2 values and dividing by 2π yields the winding number. Degenerate triangles (numerator = denominator = 0) correctly contribute zero via IEEE 754 atan2(0,0) = 0.

[github-actions]

Claude Code Review

Head SHA: 28d8d500ffa8a0249db6ba61bb3023c6a3d55356

Files changed: 1

• src/common/m_model.fpp

Summary:

• Replaces 26-ray voting (ray casting) in f_model_is_inside_flat with the generalized winding number method (Jacobson et al. 2013), independently for 3D (Van Oosterom-Strackee solid angle formula) and 2D (signed angle sum over boundary edges).
• Replaces the geometric normal-based triangle intersection test with the Möller-Trumbore algorithm (barycentric coordinates) in f_intersects_triangle.
• Both algorithms are winding-order independent, which is the stated motivation.

---

Findings

1. [Critical] gpu_boundary_v and gpu_boundary_edge_count are not defined in this diff

src/common/m_model.fpp (new 2D branch, ~line 620–636):

do q = 1, gpu_boundary_edge_count(pid)
    d1(1) = gpu_boundary_v(q, 1, 1, pid) - point(1)
    ...

These two arrays are referenced but do not appear anywhere in the diff. For this PR to be complete, they must be:

• Declared in m_global_parameters.fpp (all three targets that include m_model)
• Allocated/populated in the STL loading path (presumably somewhere in m_model.fpp)
• GPU-resident via @:ALLOCATE (or explicit GPU_ENTER_DATA) before the GPU_ROUTINE function executes on device
• Freed via @:DEALLOCATE in finalization

If these arrays are defined elsewhere in the same file outside this diff, please point to the line numbers. If they are not yet defined, the PR is incomplete and will fail to compile.

2. [Minor] atan2(0, 0) behavior is compiler-dependent

The comment at ~line 655:

! atan2(0,0) = 0 per IEEE 754, so degenerate
! triangles contribute nothing.

IEEE 754 does not define atan2(0, 0); it is implementation-defined. Fortran's intrinsic atan2(0.0, 0.0) returns 0 on most compilers (gfortran, nvfortran, ifx, CCE) in practice, but this is not guaranteed. Adding a guard (if (abs(r1_mag*r2_mag*r3_mag) < tiny(1.0_wp)) cycle) would make the degenerate-triangle handling explicit and portable.

3. [Minor] Hardcoded 1e-7_wp epsilon in Möller-Trumbore

src/common/m_model.fpp (~line 697):

if (abs(a) < 1e-7_wp) return

The old code used 0.0000001_wp (the same value). Neither is documented or tied to mesh scale. This is a pre-existing issue but worth a comment explaining the choice.

4. [Nit] PR description

The body is just "better stls". A one-paragraph description of the algorithmic change, why winding number is more robust than ray casting, and any known limitations (e.g., non-watertight meshes) would aid reviewers and future readers.

---

Correctness Assessment

The mathematics in both the 2D and 3D winding number implementations appears correct:

• 2D: atan2(d1×d2, d1·d2) gives the signed angle between edge endpoint vectors without requiring normalization (magnitudes cancel). Summing over boundary edges and dividing by 2π is standard.
• 3D (Van Oosterom-Strackee): Using non-unit r1, r2, r3 is valid — the magnitudes factor out of both numerator and denominator and cancel in the ratio passed to atan2. Division by 2π is correct given that each triangle contributes 2·atan2 to the total solid angle (Ω/2 = atan2(…)).
• Möller-Trumbore: Standard implementation, correctly winding-order independent.

The primary blocker is item 1 above: confirming gpu_boundary_v and gpu_boundary_edge_count are properly defined, populated, and GPU-resident.

[github-actions]

Claude Code Review

Head SHA: a819dd9
Files changed: 1 — src/common/m_model.fpp
+90 / -72

---

Summary

• Replaces multi-ray casting inside f_model_is_inside_flat with the generalized winding number method (Jacobson et al., 2013), which is robust to degenerate geometry and winding-order ambiguity.
• 3D path uses the Van Oosterom–Strackee solid-angle formula; 2D path sums signed angles over boundary edges.
• Rewrites f_intersects_triangle with the Möller–Trumbore barycentric algorithm — cleaner, winding-order independent, and avoids precomputed triangle normals.
• The math is correct in both cases.

---

Findings

1. Unverified 2D arrays gpu_boundary_edge_count / gpu_boundary_v (potentially blocking)

src/common/m_model.fpp, new 2D path (around line 621–634):

do q = 1, gpu_boundary_edge_count(pid)
    d1(1) = gpu_boundary_v(q, 1, 1, pid) - point(1)
    ...
    d2(1) = gpu_boundary_v(q, 2, 1, pid) - point(1)

These two arrays are referenced but are not declared, initialized, GPU-mapped, or deallocated anywhere in this diff. If they do not already exist in the file (outside the diff context), this will not compile. Please confirm that:

• They are declared (with @:ALLOCATE / @:DEALLOCATE pairs) elsewhere in m_model.fpp.
• They are populated before f_model_is_inside_flat is called.
• A $:GPU_UPDATE(device=[...]) or $:GPU_DECLARE(create=[...]) is in place so the GPU kernel can access them.

2. Near-parallel threshold 1e-7_wp in f_intersects_triangle — borderline in single-precision builds

src/common/m_model.fpp, Möller–Trumbore threshold:

if (abs(a) < 1e-7_wp) return

When --single is active, wp becomes single precision and epsilon(1.0_wp) ≈ 1.19e-7. The threshold is then nearly indistinguishable from machine epsilon, so near-parallel rays may either pass through or be spuriously rejected depending on rounding. The old code used the same literal (0.0000001_wp), so this is not a regression, but worth noting. Consider max(1e-7_wp, 10*epsilon(1.0_wp)*norm2(ray%d)*norm2(edge1)) for a scale-relative guard, or at minimum document that the value is intentionally kept for compatibility.

3. Degeneracy guard is effectively an exact-zero check

src/common/m_model.fpp, 3D winding number loop:

if (r1_mag*r2_mag*r3_mag < tiny(1.0_wp)) cycle

tiny(1.0_wp) is the smallest normal floating-point number (~2.2e-308 for double). The product of three magnitudes that are each, say, 1e-100_wp would be 1e-300_wp, which is greater than tiny and would not be skipped — yet atan2(num, denom) would still be numerically valid because atan2 handles small arguments gracefully. The guard only fires on subnormal/zero inputs. This is not a bug (the atan2 is robust), but the comment "Skip degenerate triangles" is slightly over-broad. The actual guard only prevents a divide-by-zero or NaN in the magnitude computation itself, which only arises when the query point is exactly coincident with a triangle vertex. Low risk, but the comment could be more precise (e.g., "Query point coincident with a vertex; atan2 undefined").

4. Minor: normalization comment in 3D path is imprecise

! Winding number = total solid angle / (4*pi)
! Each triangle contributes 2*atan2, so sum / (2*pi)

The Van Oosterom–Strackee formula gives tan(Ω/2) = num/denom, so each atan2 call returns Ω/2 (half the solid angle). Dividing by 2π therefore gives sum(Ω/2) / (2π) = sum(Ω) / (4π) — numerically correct, but the comment "each triangle contributes 2*atan2" reads as though the factor of 2 is already applied. Suggest: "Each atan2 returns Ω/2 per triangle; divide by 2π to get winding number = sum(Ω)/(4π)."

---

Verdict

The algorithmic replacement is a clear improvement. The Möller–Trumbore refactor is straightforward and correct. The winding-number math checks out. The main concern before merging is item 1: confirm gpu_boundary_edge_count and gpu_boundary_v exist, are GPU-resident, and are populated for the 2D path. Items 2–4 are low-severity.

[github-actions]

Claude Code Review

Head SHA: 58273b1
Files changed: 1 — `src/common/m_model.fpp`
+90 / -72

---

Summary

• Replaces 26-ray casting in `f_model_is_inside_flat` with the generalized winding number (Van Oosterom-Strackee formula) for 3D, and a signed-angle sum for 2D.
• Replaces the cross-product sign test in `f_intersects_triangle` with the Möller-Trumbore barycentric algorithm.
• The Van Oosterom-Strackee un-normalized form (raw vectors, with per-magnitude denominators) is mathematically correct and verified against the formula.
• Möller-Trumbore implementation is correct and winding-order independent.
• Since `src/common/` is shared by all three executables, the full test suite (`./mfc.sh test -j 8`) is required — not just a subset.

---

Findings

1. [Critical] Undeclared arrays gpu_boundary_edge_count and gpu_boundary_v in 2D branch

! m_model.fpp, 2D branch (lines ~630–643)
do q = 1, gpu_boundary_edge_count(pid)
    d1(1) = gpu_boundary_v(q, 1, 1, pid) - point(1)
    d1(2) = gpu_boundary_v(q, 1, 2, pid) - point(2)
    d2(1) = gpu_boundary_v(q, 2, 1, pid) - point(1)
    d2(2) = gpu_boundary_v(q, 2, 2, pid) - point(2)

Neither `gpu_boundary_edge_count` nor `gpu_boundary_v` appear anywhere in this diff. If they are not already declared as module-level GPU-resident arrays with correct dimensions, this will produce a compile error. Please confirm:

• Where they are declared (module-level in `m_model.fpp` or an imported module)?
• Where they are allocated, populated, and GPU-updated?
• Whether GPU data entry/exit (`@:ALLOCATE` / `@:DEALLOCATE`) is already handled?

This is in a `GPU_ROUTINE(parallelism='[seq]')`, so the data must be device-resident when called from a parallel kernel.

---

2. [Minor] Product underflow in coincidence guard

! m_model.fpp, 3D winding number loop
if (r1_mag*r2_mag*r3_mag < tiny(1.0_wp)) cycle

The triple product `r1_magr2_magr3_mag` can silently underflow to `0.0` when each magnitude is individually valid but small (e.g., each ~1e-103 in double, ~1e-13 in single). This would cause valid triangles to be skipped, producing incorrect winding numbers for very fine meshes near the origin. A more robust guard compares magnitudes independently:

if (min(r1_mag, r2_mag, r3_mag) < epsilon(1.0_wp)) cycle

Or equivalently check before computing magnitudes' product:

if (r1_mag < epsilon(1.0_wp) .or. r2_mag < epsilon(1.0_wp) .or. r3_mag < epsilon(1.0_wp)) cycle

---

3. [Minor] Hardcoded 1e-7_wp floor in f_intersects_triangle effectively disables double-precision tolerance

if (abs(a) < max(1e-7_wp, 10.0_wp*epsilon(1.0_wp))) return

In double precision: `10epsilon(1.0_wp) ≈ 2.2e-15`, so `max(1e-7, 2.2e-15) = 1e-7`. The guard is effectively hardcoded to `1e-7` in double-precision builds — the same threshold as single precision — so very nearly parallel rays are treated as parallel even when double precision could resolve the intersection. The comment says "use 10epsilon as a floor" but the `max` inverts that intent for double. Consider:

if (abs(a) < 10.0_wp*epsilon(1.0_wp)) return  ! precision-adaptive threshold

---

4. [Info] src/common/ blast radius — all three targets require testing

This PR modifies `src/common/m_model.fpp`, which is compiled into `pre_process`, `simulation`, and `post_process`. Per project convention, the full test suite must pass:

./mfc.sh test -j 8

Not just STL-specific or 2D/3D subsets.

---

Overall

The mathematical approach is sound and well-motivated. The Van Oosterom-Strackee and Möller-Trumbore implementations are both correct. Items 1 (undeclared arrays) and 2 (underflow) should be resolved before merging; item 3 is a robustness improvement worth considering.

[github-actions]

Claude Code Review

Head SHA: 2078c83a404995ef03fd3efeca904eea45e6431f
Files changed: 1 — src/common/m_model.fpp
+90 / -72

---

Summary

• Replaces 26-ray casting in f_model_is_inside_flat with the generalized winding number (Van Oosterom-Strackee / Jacobson 2013). This is a well-known, theoretically sounder approach for STL inside/outside queries.
• Replaces the winding-order-dependent cross-product sign test in f_intersects_triangle with the Möller-Trumbore barycentric algorithm. Correct and standard.
• 2D branch now sums signed angles over boundary edges instead of a ray-cast parity count. Numerically cleaner.
• Code is significantly shorter, better documented, and more robust to degenerate geometry.
• Precision discipline is good: uses wp literals throughout; no forbidden intrinsics or raw pragmas.

---

Findings

1. Product underflow in coincidence guard — medium

! src/common/m_model.fpp (new code, ~line 658)
if (r1_mag*r2_mag*r3_mag < tiny(1.0_wp)) cycle

Multiplying three magnitudes can silently underflow to 0.0 even when each individual magnitude is much larger than zero. In single-precision (--single build), tiny(1.0_wp) ≈ 1.18e-38; three magnitudes of ~2.3e-13 each yield a product ≈ 1e-38 that may underflow before the comparison. A safer guard checks each magnitude individually:

if (r1_mag < tiny(1.0_wp)**0.34_wp .or. &
    r2_mag < tiny(1.0_wp)**0.34_wp .or. &
    r3_mag < tiny(1.0_wp)**0.34_wp) cycle

or simply:

if (min(r1_mag, r2_mag, r3_mag) < sqrt(tiny(1.0_wp))) cycle

2. Misleading inline comment — minor

! Each atan2 returns Omega/2 per triangle; divide
! by 2*pi to get winding number = sum(Omega)/(4*pi).
fraction = fraction/(2.0_wp*acos(-1.0_wp))

The comment is partially correct: each atan2 call does return Ω/2. But the expression before division is Σ(Ω/2), not Σ(Ω). The division by 2π therefore gives Σ(Ω/2)/(2π) = Σ(Ω)/(4π) — the math is right, but the comment implies the accumulator holds Σ(Ω) which it doesn't. A cleaner note:

! accumulator = sum(Omega/2); winding number = sum(Omega)/(4*pi) = accumulator/(2*pi)

3. ntrs parameter unused in 2D branch — minor

In the 2D path, ntrs (the triangle count) is never referenced; gpu_boundary_edge_count(pid) is used instead. The parameter still has a clear purpose for the 3D path so this isn't a bug, but a brief comment like ! ntrs unused in 2D; gpu_boundary_edge_count drives the loop near the 2D branch would prevent confusion.

4. Near-parallel threshold is independent of a's scale — low

if (abs(a) < max(1e-7_wp, 10.0_wp*epsilon(1.0_wp))) return

a = dot_product(edge1, h) is a dot product of edge and cross-product vectors, so its magnitude scales with the triangle's area times the sine of the ray-to-normal angle. A fixed threshold of 1e-7 can incorrectly reject large-area triangles with nearly-parallel rays, or accept tiny triangles with non-parallel rays. The old code used the same fixed threshold (0.0000001_wp), so this is a pre-existing limitation and acceptable to leave as-is, but noting it for the record.

---

Not an issue

The 3D normalization fraction / (2π) is mathematically correct: Van Oosterom-Strackee returns Ω/2 per triangle via atan2, so Σ atan2 / (2π) = Σ(Ω) / (4π) as required. The 2D normalization by 2π is also correct for a simple closed polygon. Both normalizations yield ≈1.0 for interior points on a watertight mesh.

[github-actions]

Claude Code Review

Head SHA: e8aff26

Files changed: 1

• src/common/m_model.fpp (+90 / -72)

---

Summary

• Replaces 26-ray casting in f_model_is_inside_flat with the generalized winding number (Jacobson et al., 2013), correctly using the Van Oosterom-Strackee formula.
• The 2D branch now uses the signed-angle winding number over boundary edges rather than over triangles; gpu_boundary_v and gpu_boundary_edge_count are confirmed to exist in master and are indexed consistently with their existing usage elsewhere in the file.
• Replaces the winding-order-dependent cross-product sign test in f_intersects_triangle with standard Möller–Trumbore; implementation is correct.
• Precision discipline is respected (_wp literals, generic intrinsics, no d-exponent literals).
• GPU annotations (GPU_ROUTINE) are preserved on both modified functions; no raw !/! pragmas introduced.

---

Findings

[Minor – numerical stability] src/common/m_model.fpp — f_intersects_triangle, new parallelism check:

if (abs(a) < max(1e-7_wp, 10.0_wp*epsilon(1.0_wp))) return

a is the dot product of an unnormalized ray direction and an unnormalized edge vector, so its magnitude depends on geometry scale. The fixed floor of 1e-7_wp is scale-dependent and may be too aggressive (skipping valid intersections) for unit-scale geometry, or too lenient for very small triangles. Consider normalizing ray%d before entering the function, or document the assumption that directions are unit-length. The old code had the same literal (0.0000001_wp), so this is not a regression, but the refactor is a good opportunity to address it.

[Minor – edge case, 3D winding number] src/common/m_model.fpp — new f_model_is_inside_flat 3D branch:

if (r1_mag*r2_mag*r3_mag < tiny(1.0_wp)) cycle

The guard catches exact coincidence with a vertex, but not near-coincidence (e.g. query point ~sqrt(tiny) from a vertex). In that regime atan2(numerator, denominator) can be ill-conditioned and return ±π/2, contributing a spurious solid-angle term. A slightly larger guard (e.g. sqrt(tiny(1.0_wp))) would improve robustness at negligible cost. This is a pre-existing limitation of the algorithm that is acknowledged in the PR description.

[Correctness – verified ✓] The 3D division:

fraction = fraction/(2.0_wp*acos(-1.0_wp))

Each atan2 call returns Ω/2 per triangle; dividing the accumulated sum by 2π correctly yields Σ(Ω)/(4π) → 1.0 inside a closed mesh. Math checks out.

[Correctness – verified ✓] 2D indexing matches master convention: gpu_boundary_v(edge, vertex, coord, pid) with vertex∈{1,2} and coord∈{1=x,2=y}. Consistent with lines ~1000–1003 and ~1022–1023 in the base file.

---

Opportunities (non-blocking)

• The acos(-1.0_wp) idiom for π appears twice. A named parameter (real(wp), parameter :: pi = acos(-1.0_wp)) at module scope would make intent clearer and is consistent with what other CFD codes in the vicinity do.
• The comment "divide by 2*pi to get winding number = sum(Omega)/(4*pi)" is mathematically correct but slightly cryptic. A one-line derivation (e.g. ! sum(Ω_i/2) / (2π) = sum(Ω_i) / (4π) → 1.0 inside) would help future readers.

---

Overall the algorithmic replacement is well-motivated, mathematically sound, and cleaner than the old approach. The two minor numerical points above are both pre-existing in spirit and do not block merging.