API Reference

load_mesh(filename; surface_dim=2, reverse_normals=false, verbose=true) -> MeshData

Load a mesh and extract all supported 1st- or 2nd-order elements belonging to named physical groups.

Any file format that the Gmsh API can open is supported — the format is inferred from the extension. This includes Gmsh .msh (all versions), and via Gmsh's importers .stl, .step/.stp, .bdf/.nas (Nastran), .med, .vtk (legacy ASCII unstructured grids), and others.

Supported element families (any mix in one mesh):

• Surfaces (surface_dim=2): Tri3, Quad4 (1st order); Tri6, Quad8/Quad9 (2nd order).
• Curves (surface_dim=1): Line2 (1st order); Line3 (2nd order).

surface_dim=2 (default) — surface mesh, 3D view factors. surface_dim=1 — planar curve mesh, 2D view factors.

reverse_normals=true flips every element's normal at load time (for a surface mesh wound back-to-front). For curve meshes it is applied after the automatic orientation described below. To flip only some groups, use reverse_group_normals on the loaded mesh instead.

For curve meshes read through Gmsh, each element's normal is oriented toward the interior of the adjacent surface, found from mesh connectivity (see the Mesh Requirements page). Surface-mesh normals follow node winding as written.

Named physical groups are used to partition the radiating surfaces. Formats that cannot carry physical groups (e.g. STL) have no named groups; in that case a single synthetic group named "default" covering every entity of surface_dim is created so the rest of the pipeline behaves uniformly.

XML VTK files (.vtu and XML-form .vtk) are detected automatically and read through ReadVTK.jl (a weak dependency — using ReadVTK must be in scope). Legacy ASCII/binary .vtk files are not handled by ReadVTK and fall through to the Gmsh importer. See load_vtu for VTK-specific options.

Nek5000/NekRS .re2 binary meshes are detected by extension and read by a dedicated in-tree parser (Gmsh cannot open them). The 3D hex volume mesh's boundary faces become radiating surfaces (Quad8 when the file carries curved-side records, else Quad4) grouped by Nek boundary-condition label. See load_re2.

load_vtu(filename; surface_dim=2, reverse_normals=false, verbose=true,
         group_field=nothing) -> MeshData

Load an XML VTK unstructured grid (.vtu) through ReadVTK.jl. Requires using ReadVTK to be in scope (ReadVTK is a weak dependency).

VTK has no native concept of physical groups. If group_field names a per-cell integer data array it is used to partition elements into groups; otherwise (or if the named array is absent) a single "default" group is created. Common region arrays (CellEntityIds, gmsh:physical, RegionId, MaterialIds, region, group) are tried automatically when group_field === nothing. Groups found this way are named "group_<value>".

Unlike the Gmsh path, no automatic normal orientation is applied to curve meshes (a warning is issued): normals follow node winding, so use reverse_normals if they come out backwards. Cells of a type not listed below are skipped.

Cell types are mapped to element families as: Line(3)→line2, QuadraticEdge(21)→line3, Triangle(5)→tri3, QuadraticTriangle(22)→tri6, Quad(9)→quad4, QuadraticQuad(23)→quad8.

load_re2(filename; surface_dim=2, reverse_normals=false, curved=true,
         verbose=true) -> MeshData

Load a Nek5000/NekRS .re2 binary mesh. The 3D hex volume mesh's boundary faces become radiating Quad8 surfaces when the file carries mid-side-node ('m') curved-side records and curved=true (the default), and Quad4 surfaces otherwise, grouped by their Nek boundary-condition label. Only genuinely internal faces (cbc = 'E' or blank) are skipped; periodic ('P') faces are kept as their own group, matching Nek5000/NekRS's own view-factor convention of including every non-internal boundary face in the radiation enclosure so that view factors close (row sums to 1) — solid walls vs. other boundaries is a runtime (cbc.eq.'W ') distinction made by the case's .usr file, not something this loader can determine from .re2 labels alone. Each element's Nek global element number and local face index (1-6) are recorded in SurfaceElement.eg/.iface for round-tripping view factors back into Nek5000's (element, face) bookkeeping — see write_nekrs_view_factors. Word size (4- or 8-byte reals) and byte order are auto-detected. Only surface_dim=2 (3D → surfaces) is supported.

Curvature (curved)

A Nek element face is curved, and .re2 records that curvature as a mid-side point per curved edge. Representing such a face by the flat quad through its 4 corners understates its area — by 1.5 % for a sphere meshed 6x6 per cubed-sphere block — and every view factor from it with it. Worse for a solver, the resulting F is reciprocal (AᵢFᵢⱼ = AⱼFⱼᵢ) with respect to those flat areas while Nek5000/NekRS integrates fluxes over the curved faces, so the enclosure energy balance Σᵢ Aᵢ qᵢ = 0 no longer holds: a 1.5 % area error became a 10 % error in the net radiative exchange in the single-pebble case that prompted this, because the net flux is a small difference of large radiosities.

With curved=true every boundary face therefore becomes a Quad8 whose mid-side nodes come from the file's 'm' records (edges without a record get the straight mid-point, so a mesh is all-Quad8 or all-Quad4 and the Duffy near-pair path stays available for every pair). curved=false restores the previous corners-only behaviour; files with no 'm' records load as Quad4 either way, unchanged.

Only the 'm' form is read. The analytic curved-side forms ('C' circle, 's' sphere) describe an edge by radius/centre rather than by a point and are ignored — gmsh2nek writes 'm' for every curved edge it exports.

MeshData

Fields

• coords : (3 × N_nodes) coordinate matrix
• surface_elems : vector of SurfaceElement (all radiating curves/surfaces)
• group_tags : Dict tag → name
• group_elems : Dict tag → element indices (into surface_elems)
• group_tri_soup : Dict tag → obstruction geometry for that group. For surface meshes (meshdim=2): (3, 3, Ntris) triangle soup. For curve meshes (meshdim=1): (3, 2, Nsegs) segment soup, axis 1 = xyz, axis 2 = endpoint (1 or 2), axis 3 = segment.
• mesh_dim : 1 for curve meshes, 2 for surface meshes

SurfaceElement(nodes, group, family, eg=0, iface=0, phys_tag=0)

eg and iface are the Nek5000/NekRS global element number and local face index (1-6) this element corresponds to, populated by load_re2 so that per-element view factors can be written back out in Nek's own (element, face) bookkeeping (see write_nekrs_view_factors). They default to 0 for elements loaded from any other mesh format, where the concept doesn't apply.

phys_tag is the original numeric Gmsh Physical Surface tag for this face (Nek's bc(5,ifc,iel,1), the same value usrdat2 reads), also populated by load_re2 (0 otherwise). It exists because group/ group_tags are keyed by the Nek boundary-condition code ('W', 'P', or a generic placeholder like "MSH" that gmsh2nek writes for every ordinary, non-periodic surface) — meshes with several differently-named Physical Surfaces that all end up with the same generic code (e.g. two concentric spheres, or a pebble bed's pebble/duct-wall surfaces) are otherwise indistinguishable by group alone, even though usrdat2 can tell them apart. Use split_groups_by_tag to regroup a mesh at (code, phys_tag) granularity when you need to select one such surface as an obstruction_groups entry.

reverse_group_normals(mesh::MeshData, groups) -> MeshData

Return a new MeshData with normals reversed for elements belonging to groups (a single tag or a collection of tags), leaving every other element's node order — and therefore normal — unchanged. Uses the same per-family swap as load_mesh's/load_re2's mesh-wide reverse_normals flag (see its docstring above for the corner/mid-side permutation table), just applied selectively.

This is the tool for meshes whose normal orientation needs differ between physical groups rather than being uniform across the whole mesh — most commonly, a CSG-cut body containing another (a smaller sphere or cylinder sitting inside a larger one, cut apart with a boolean cut/subtract): Gmsh/OCC's default orientation is typically already correct (outward from each body) for a body defined as its own solid, independent of any later cut, but a single global reverse_normals=true flips both surfaces — correcting the outer one (which does need flipping, from outward-from-the- whole-assembly to pointing into the gap) while wrongly over-correcting the inner one (which was already right), leaving it pointing into its own volume instead of away from it. Two symptoms distinguish this from a genuine modelling error: (1) the pair-area kernels (quadrature, monte_carlo) still recover roughly the correct aggregated view factor despite it — their point-pair cosine test doesn't care which of a pair's two normals is inward as long as the sign works out over the full double integral — while check_reciprocity/row-sum closure quietly does not (Σⱼ Fᵢⱼ stays far from 1, since the inner body's genuine zero self-view is being computed as if it weren't zero); and (2) raytrace=true fails outright rather than subtly, because a ray sampled from a cosine-weighted hemisphere around an inward-pointing normal is aimed into the body's own volume, and — starting exactly on that body's own surface — is geometrically guaranteed to exit back through the same body (a chord) rather than ever reaching anything else.

Example

mesh = load_mesh("concentric.msh"; surface_dim=2)   # reverse_normals=false: leave the default
mesh = reverse_group_normals(mesh, 2)               # flip only the outer body's group

restrict_to_radiating(mesh::MeshData, radiating_groups::Vector{Int}) -> MeshData

Return a MeshData whose radiating surface is only the elements of radiating_groups, while every original group is retained as obstruction geometry.

This is the lever for "I want one pair of surfaces, not the whole enclosure". compute_view_factors assembles a dense matrix over every element it is given, so loading a pin plus the two dozen pins that might shadow it costs O(N_total²) — and most of that work computes view factors between shadowing bodies that were never asked for. Restricting the radiating set first makes the cost O(N_radiating²) while leaving occlusion exactly as accurate, because group_tri_soup — the only thing the obstruction BVH reads — is carried over whole.

For a reactor fuel-assembly slice, computing FA_wall → Pin_11 with 19 shadowing pins loaded drops from 18400 radiating elements to 3040, a factor of 37 in pair count.

Fields are rebuilt as follows: coords and group_tri_soup are shared with the original (all groups, so obstruction is unchanged); surface_elems keeps only the radiating elements; group_tags and group_elems are restricted to the radiating groups and renumbered to the new element indices.

mesh = load_mesh("assembly.msh")
r = compute_view_factors(mesh;
                         radiating_groups   = [tag_wall, tag_pin11],
                         obstruction_groups = collect(keys(mesh.group_tags)))

split_groups_by_tag(mesh::MeshData) -> MeshData

Return a new MeshData with group/group_tags/group_elems/ group_tri_soup refined to (code, phys_tag) granularity instead of just code. load_re2's default grouping is by Nek boundary-condition code ('W', 'P', or a generic placeholder like "MSH"), which gmsh2nek writes identically for every ordinary, non-periodic surface regardless of which named Physical Surface it came from — so two concentric spheres, or a pebble bed's pebble/duct-wall surfaces, land in the same group and can't be told apart via group alone even though usrdat2 distinguishes them via bc(5,ifc,iel,1) (exposed here as SurfaceElement's phys_tag).

Call this once after load_re2 when you need to pass one particular named surface to compute_view_factors's obstruction_groups — e.g. the inner sphere blocking the outer sphere's concave self-view, or a pebble blocking view between two duct-wall patches. Elements with phys_tag == 0 (any mesh not loaded from a .re2 file) keep their original group unchanged, so this is a no-op for non-Nek meshes.

mesh  = load_re2("case.re2")
mesh2 = split_groups_by_tag(mesh)
# mesh2.group_tags now has one entry per (code, phys_tag) pair, e.g.
# 1 => "MSH", 2 => "MSH#2", ... instead of a single merged "MSH" group.
inner = findfirst(t -> occursin("#1", t), collect(values(mesh2.group_tags)))
compute_view_factors(mesh2; obstruction_groups=[inner], ...)

compute_view_factors(mesh; nquad=4, obstruction_groups=Int[],
                     radiating_groups=Int[], backend=CPU(),
                     self_vf=false, monte_carlo=false, n_samples=5000,
                     rng=Random.default_rng(), use_duffy=false,
                     factor=3.0, facing_cull=true,
                     raytrace=false, n_rays=10000, verbose=true)
                     -> ViewFactorResult

Assemble the full view factor matrix at element and physical-group level.

Arguments

• mesh : MeshData returned by load_mesh
• nquad : Gauss points per direction; nquad² points on each quadrilateral element (nquad on each curve element), so the pair integral costs nquad⁴ point-pairs per quadrilateral pair (nquad² per curve pair) — 256 point-pairs at the default nquad=4. Triangles use a Dunavant rule instead, with 1, 3, 7 or 13 points for nquad = 1, 2, 3 or ≥ 4. Compare n_samples below, which counts point-pairs directly. When monte_carlo=true, only used for the adjacent-pair Duffy patch (see monte_carlo below), not the bulk sampling. With raytrace=true, only used to compute element areas.
• obstruction_groups : physical group tags that may occlude rays. Source and destination groups are excluded automatically per pair. Also applied to the Monte Carlo near-pair Duffy patch. Obstruction geometry is triangulated from element corner nodes only, so a curved 2nd-order blocker is treated as flat.
• radiating_groups : restrict the radiating surface to these physical group tags; every other group still obstructs (if listed in obstruction_groups) but is not itself assembled. Empty (default) means every group radiates, the previous behaviour. Because assembly is dense over whatever it is handed, this turns an O(N_total²) problem into O(N_radiating²): asking for one pair of surfaces in a mesh full of shadowing bodies no longer computes the view factors between those shadowing bodies. On a reactor fuel-assembly slice, FA_wall → Pin_11 with 19 shadowing pins loaded drops from 18400 radiating elements to 3040 — 37x fewer pairs. The returned F_group covers only the radiating groups. The enclosure is deliberately open, so row sums no longer approach 1 and check_closure is not meaningful; reciprocity is unaffected. See restrict_to_radiating.
• backend : CPU() (default), CUDABackend(), or MetalBackend()
• self_vf : include self view factors (concave elements). CPU only — silently ignored on GPU backends — and an error with raytrace=true.
• monte_carlo : use stratified Monte Carlo area-sampling instead of Gauss–Legendre quadrature for the O(N²) bulk of element pairs — the fastest option for large meshes, especially on GPU. The 1/r² kernel has unbounded variance for pairs sharing a vertex or edge, so no amount of sampling fixes those; this function therefore always patches those O(N) pairs afterward with the deterministic Duffy transform (using nquad above; on the device on a GPU backend), on top of whichever backend the bulk ran on. use_duffy has no separate effect here — the patch is unconditional whenever monte_carlo=true.
• n_samples : MC sample pairs per element pair (default 5000). Ignored when monte_carlo=false. Variance decreases as O(1/N) for the (non-adjacent) pairs Monte Carlo handles. This counts point-pairs directly, so it is nquad⁴ — not nquad² — that it should be compared against: n_samples=5000 does ~20x the kernel evaluations of the default nquad=4, and with obstruction_groups set, ~20x the BVH ray casts, since both paths ray-cast once per kernel-positive point-pair. On CPU that made a 4320-element obstructed case take 1368 s at n_samples=5000 versus 60 s at nquad=4, for answers agreeing to 5 decimals. The default of 5000 is the value used throughout benchmarks/howell/ (115 catalog points, 96% within 1% of the published value, median error 0.008%) — at that sample count the MC estimator's own noise is typically 10-1000x smaller than the mesh's discretization error against the analytic answer, so raising it further mostly buys accuracy the mesh cannot use. Lower it for a quick look at a large mesh; raise it only if n_samples sweeps on your own geometry show the noise still dominates.
• rng : RNG for the CPU Monte Carlo paths (monte_carlo and raytrace). Pass a seeded RNG (e.g. MersenneTwister(42)) for reproducible results. Ignored on GPU, where a fresh random seed is drawn on each call.
• use_duffy : apply the Duffy singularity transformation (see DuffyKernel.jl) for same-order quad pairs (Quad4 or Quad8) sharing a vertex or edge, in the plain quadrature path, on the CPU or (for 3-D meshes) on a GPU backend, where the touching pairs are re-evaluated by a device kernel. Has no effect when monte_carlo=true (that path always Duffy-patches adjacent pairs regardless of this flag) or raytrace=true. Not applicable to curve meshes.
• factor : near-pair patch radius, in element diameters (see near_pairs in DuffyKernel.jl) — only used when monte_carlo=true, to decide which O(N) pairs get the deterministic Duffy patch instead of the raw MC estimate. Raising it does not, by itself, make those pairs more accurate — nquad (above) is what governs the patch's own resolution, and is the more effective lever for closure error on meshes with large or elongated elements.
• raytrace : use ray-shooting Monte Carlo (RayTraceKernel.jl) instead of either quadrature or the pair-area monte_carlo above. For each element, shoots n_rays cosine-weighted rays and tallies which element each first hits, via one BVH over the whole radiating mesh — O(N·nrays·log N) instead of the O(N²) pair loop the other two paths always pay. Obstruction is then a side effect of the same query, not a separate check: radiating elements obstruct each other automatically (a real behavioural difference from the other two paths, which only apply obstruction for groups explicitly listed in `obstructiongroups— see that argument above);obstructiongroupshere only adds *extra*, non-radiating blocker geometry to the scene. No adjacent-pair singularity exists for this estimator (it never evaluates 1/r²), souseduffy/the Duffy patch/factordo not apply.facingculldoes not apply either (there is no O(N²) pair loop to cull). Runs on CPU or GPU (backend=CUDABackend()/MetalBackend(), same asmontecarlo) with its own GPU kernel (GPURayTraceKernels.jl); one thread per *element* on GPU, not per pair. 3-D meshes only (mesh.meshdim == 2); incompatible withselfvf(a ray is never tested against its own origin element) and withmontecarlo=true. Reciprocity is enforced by construction (each unordered pair's two independent ray-based estimates — from each side — are averaged), so results are exactly symmetric despite being stochastic, unlike a raw hit-count would be. This same averaging is why row sums on a closed enclosure are only approximately 1 (ordinary MC noise, shrinking withnrays) rather than exact: each element's own rays alone would close exactly, but every off-diagonal entry blends in the *other* element's independent estimate too. See the module docstring inRayTraceKernel.jlfor the derivation andchangelog.md` for validation against analytic cases and against quadrature (including a case where quadrature itself was the less accurate one, under-resolving a sharp obstruction shadow boundary).
• n_rays : rays per element (default 10000) when raytrace=true; ignored otherwise. Unlike n_samples (which counts point-pairs per element pair), one ray-shooting pass from an element estimates its view factor to every other element at once, so this is not directly comparable to n_samples — tune per-case rather than assuming the same order of magnitude applies.
• facing_cull : reject element pairs that provably cannot see each other before integrating them (default true). Each element gets a conservative bounding box for its points and another for its normals (O(N)); a pair whose kernel is then provably zero at every point pair is skipped in O(1) instead of being discovered zero one quadrature point — or one MC sample, or one BVH ray cast — at a time. Results are unchanged; on closed convex bodies (tube bundles, pebble beds) most pairs face away and this is a large saving. Set false to disable, e.g. for a mesh whose 2nd-order elements are curved sharply enough that the sampled bounds might not contain them (see ElementBounds).
• verbose : print progress and row-sum diagnostics

Returns

A ViewFactorResult. Its F_group is ordered by the sorted group tags in result.group_tags; with radiating_groups set it covers only those groups. Curve meshes (mesh.mesh_dim == 1) give 2D view factors; they run on the CPU only and cannot use raytrace.

Examples

result = compute_view_factors(mesh; nquad=6)
result = compute_view_factors(mesh; nquad=6, use_duffy=true)
result = compute_view_factors(mesh; monte_carlo=true, n_samples=50000)               # CPU
result = compute_view_factors(mesh; monte_carlo=true, n_samples=50000,
                               backend=CUDABackend())                                 # GPU — fastest for large meshes
result = compute_view_factors(mesh; raytrace=true, n_rays=10000)                     # CPU, fastest for large/obstructed 3D meshes

# Just one pair of surfaces, with everything else still shadowing them:
all_tags = collect(keys(mesh.group_tags))
result = compute_view_factors(mesh; nquad=4,
                               radiating_groups   = [wall_tag, pin_tag],
                               obstruction_groups = all_tags)

ViewFactorResult

Stores element-level and group-level view factor matrices produced by compute_view_factors.

Fields

• F_elem : (Nelem × Nelem) element-level view factor matrix
• A_elem : N_elem-vector of element areas (or lengths for curve meshes)
• F_group : (Ngroup × Ngroup) group-level view factor matrix
• A_group : N_group-vector of group areas/lengths
• group_tags : sorted physical group tags (row/column labels for F_group)
• group_names : physical group names in the same order

aggregate_by_group(result, mesh) -> (F_group, A_group, tags, names)

Re-aggregate element-level results to group level (useful after modifying result.F_elem).

check_reciprocity(result; tol=1e-4) -> Bool

Verify Aᵢ Fᵢⱼ ≈ Aⱼ Fⱼᵢ for all element pairs. Prints the maximum relative error and returns true if it is below tol.

check_closure(result; tol=1e-3) -> Bool

Verify that no row of F_elem sums to more than 1 + tol. Prints the row-sum range and returns true if the check passes.

Only meaningful for a closed enclosure. A result computed with compute_view_factors(...; radiating_groups=...) is deliberately open — the non-radiating bodies absorb their share but are never assembled — so its rows sum to less than 1 by construction and a low row sum is not a defect. Use check_reciprocity, which is unaffected, to validate such a result.

enforce_closure(result, mesh; iters=200, tol=1e-12, verbose=true)
    -> ViewFactorResult

Return a copy of result whose F_elem satisfies both

Σⱼ Fᵢⱼ = 1            (closure, every row)
Aᵢ Fᵢⱼ = Aⱼ Fⱼᵢ        (reciprocity, with `result.A_elem`)

to machine precision, by alternating a symmetrisation of Aᵢ Fᵢⱼ with a rescaling of each row — each step perturbs the other, so they are iterated (at most iters times, stopping once both residuals are below tol).

Why a converged calculation still needs this

Net radiative flux is a difference of radiosities, qᵢ = Jᵢ - Σⱼ FᵢⱼJⱼ, and in a nearly-isothermal enclosure that difference is orders of magnitude smaller than J itself: in the single-pebble NekRS case this function was written for, J ≈ 5.9e4 W/m² and q ≈ 2.5e2 W/m². A row-sum error of ε therefore leaks ε·J straight into q — 0.3 % of closure becomes 70 % of the answer. Exact row sums remove that channel entirely: with Σⱼ Fᵢⱼ = 1 an isothermal enclosure gives q = 0 identically, whatever the remaining per-pair errors, so only genuine temperature differences drive a flux.

Reaching that by sampling alone is hopeless — Monte Carlo closure error falls as 1/√n_rays, so ε = 1e-5 would need about 1e10 rays per element — and quadrature on curved elements plateaus around 1e-3. Enforcing the two identities the view factors must satisfy anyway is the cheap route.

When not to use it

Only for a closed enclosure, where every row genuinely sums to 1. A result from compute_view_factors(...; radiating_groups=...), or any geometry open to the surroundings, has rows that legitimately sum to less than 1 (see check_closure), and forcing them to 1 would invent radiation that is not there.

write_nekrs_view_factors(path, result, mesh;
                          tol=0.0, max_visible_walls=10_000,
                          max_walls=4_000_000, verbose=true) -> Int

Write result.F_elem to path in the list-directed ASCII format read by Nek5000/NekRS's vf_read_view_factors (in the Nek-VF view_factors.f module):

nwalls
iw₁ ieg₁ ifc₁ nvwalls₁
  jw ieg ifc Fij        (nvwalls₁ lines)
iw₂ ieg₂ ifc₂ nvwalls₂
  ⋮

iw runs 1:N in mesh.surface_elems order (Nek5000 does not require any particular enclosure-face numbering — it rebuilds its own (eg,ifc) → iw map from this file's ieg/ifc columns at load time). ieg/ifc are read from SurfaceElement.eg/.iface, so mesh must come from load_re2 — an error is raised otherwise.

Only pairs with F_elem[i,j] > tol are written as visible faces per wall, matching how Nek5000's hemi-cube view-factor calculator only records genuinely visible (rasterized) pairs; tol=0.0 (default) keeps every strictly-positive entry, which for the quadrature/Monte Carlo/Duffy kernels here already excludes back-facing and obstructed pairs (vf_kernel returns exactly 0.0 for those). Raise tol to shrink the file for dense meshes at the cost of a small closure error.

The same file is used for both open_face_option=1 and =2 in vf_calculate_radiation_heat_flux — the option only changes which faces emit at runtime (cbc.eq.'W ', decided by the case's .usr file), not which faces appear in this file. This function therefore writes every element in mesh.surface_elems, walls and non-walls (e.g. periodic) alike — see the load_re2 docstring for why non-wall boundary faces must stay in the enclosure.

Prints (does not throw) a warning if any row's visible-face count would exceed max_visible_walls, or if N exceeds max_walls — these mirror the parameters of the same name in view_factors/VIEW_FACTORS and only warn, not abort, on the Fortran side too; increase those parameters if you see this warning and re-run the Nek5000/NekRS case.

Returns N, the number of walls written.

plot_mesh_normals(mesh; normal_scale=nothing, group_colors=nothing,
                  show_nodes=false, show_indices=false) -> Plots.Plot

Visualise mesh elements with normal arrows coloured by physical group.

Requires Plots.jl to be loaded first:

using Plots

Curve meshes (mesh.mesh_dim == 1) are drawn in the xy-plane. Surface meshes are drawn as an xy-projection of the element edges and normals — a quick sanity check, not a 3-D view.

Arguments

• mesh : MeshData from load_mesh
• normal_scale : arrow length in mesh units. Auto-estimated from the bounding box diagonal if omitted.
• group_colors : Dict{Int,Any} mapping physical group tag → any colour accepted by Plots.jl (e.g. :red, "#FF0000"), overriding the automatic palette for specified groups.
• show_nodes : scatter-plot all element nodes.
• show_indices : annotate each element with its index number.

Returns

The Plots.Plot object; save it with savefig(fig, "file.png").

Example

using Plots, RadiativeViewFactor
mesh = load_mesh("geometry.msh"; surface_dim=1)
fig  = plot_mesh_normals(mesh; normal_scale=0.05, show_nodes=true)

quad8_physical_point(coords, nodes, ξ, η) -> SVector{3,Float64}

Map (ξ,η) to physical space. nodes may be a Vector{Int} or SVector{8,Int}.

quad8_normal_and_area_element(coords, nodes, ξ, η) -> (n̂, dA)

Compute the unit normal and area element at (ξ,η). nodes may be a Vector{Int} or SVector{8,Int}.

quad4_shape(ξ, η) -> (N, dNdξ, dNdη)

Bilinear shape functions of the Quad4 element at (ξ,η) ∈ [-1,1]², and their derivatives, as SVector{4,Float64}s in Gmsh corner order (-1,-1), (1,-1), (1,1), (-1,1).

quad4_physical_point(coords, nodes, ξ, η) -> SVector{3,Float64}

Map (ξ,η) ∈ [-1,1]² to physical space for a Quad4 element. nodes is a length-4 index vector in Gmsh corner order.

quad4_normal_and_area_element(coords, nodes, ξ, η) -> (n̂, dA)

Compute the unit normal and area element at (ξ,η) for a Quad4 element. The normal follows the node winding (right-hand rule).

line3_physical_point(coords, nodes, ξ) -> SVector{3,Float64}

Map reference coordinate ξ ∈ [-1,1] to physical space for a Line3 element. nodes is a length-3 index vector [left corner, right corner, midpoint].

line3_normal_and_length_element(coords, nodes, ξ) -> (n̂, dL)

Compute the in-plane unit normal and arc-length element at ξ for a Line3 element lying in the xy-plane (z = const).

Returns

• n̂ : SVector{3,Float64} — unit normal, z-component is zero
• dL : Float64 — arc-length element |dx/dξ|

line2_shape(ξ) -> (N, dNdξ)

Linear shape functions of the Line2 element at ξ ∈ [-1,1], and their derivatives, as SVector{2,Float64}s ordered [first endpoint, second endpoint].

line2_physical_point(coords, nodes, ξ) -> SVector{3,Float64}

Map reference coordinate ξ ∈ [-1,1] to physical space for a Line2 element. nodes is a length-2 index vector [first endpoint, second endpoint].

line2_normal_and_length_element(coords, nodes, ξ) -> (n̂, dL)

Compute the in-plane unit normal and arc-length element for a Line2 element lying in the xy-plane. Both are constant along the element; the normal is the tangent rotated 90° counter-clockwise, and a zero-length element returns a zero normal and dL = 0.