Add near-to-far-field transform and half-wave dipole validation

docs/src/maxwell.md

@@ -83,6 +83,44 @@ Z = FiniteElementMethod.z_parameters(model, S)
 The port couples all its edge dofs to each other, which widens the sparsity pattern, so keep the sheet a
 few elements across.
 
+## Far field, radiated power and directivity
+
+[`far_field`](@ref) transforms the solution on a closed facet set inside the mesh, the Huygens surface, into
+the far field: Love's equivalent currents `J = n × H` and `M = -n × E` on its outward normal radiate into
+vacuum, and with `exp(+iωt)` the result is the pair `E_θ, E_φ` of `r e^{ikr} E` as `r → ∞`. The surface must
+enclose every source, port and conductor, lie in vacuum, touch no PEC facet and be star-shaped about its
+`center` (default: its centroid), which decides the outward side; a sphere or box between the antenna and an
+absorbing boundary is the usual choice, and Gmsh physical groups name it. `H = (i/ωμ0) ∇ × E` is averaged
+over the cells on both sides of each facet. [`radiated_power`](@ref) integrates the Poynting flux through the
+surface (`method = :flux`) or the radiation intensity over the sphere (`method = :far_field`), and
+[`directivity`](@ref) is `4πU/P` with the far-field power, so it does not depend on the overall amplitude.
+
+The far field contains only the outgoing wave, while the flux and the port impedance also feel the wave the
+first-order absorbing sphere reflects back, about `1/(kR)²` in amplitude for a dipole field at radius `R`
+(5 % at `R = 0.75λ`): directivity and pattern are accurate on coarse meshes, the impedance needs the boundary
+a wavelength or more away.
+
+```julia
+using FiniteElementMethod, Ferrite
+# `grid` from Gmsh (FerriteGmsh.togrid) with facet sets "outer", "huygens", "strip" and "gap" (in metres).
+model = MaxwellModel(grid; pec = ["strip"], absorbing = ["outer"], ports = [LumpedPort("gap"; direction = (0, 0, 1))])
+sweep = solve_sweep(model, [2π * 3.0e8]; fields = true)
+Zin = only(FiniteElementMethod.z_parameters(model, sweep.S[:, :, 1]))
+sol = MaxwellSolution(model, 2π * 3.0e8, sweep.fields[1][:, 1])
+θs = range(0, π; length = 91)
+ff = far_field(sol, "huygens", [Vec(sin(θ), 0.0, cos(θ)) for θ in θs])  # E-plane cut
+D = directivity(sol, "huygens", [Vec(1.0, 0.0, 0.0)])
+```
+
+`examples/dipole` builds this model for a strip dipole (Gmsh OCC: an absorbing sphere, a Huygens sphere and a
+zero-thickness strip with a centre gap, graded towards the strip) and finds its resonance, impedance,
+directivity and beamwidth; `test/Maxwell/dipole.jl` runs a coarse version, and the full-size one with
+`FEM_MAXWELL_DIPOLE=1`. On the full mesh (1.4e5 edge dofs, a strip of equivalent radius 0.005λ, the absorbing
+sphere at 1λ) the dipole resonates at a length of 0.456λ with `R = 69.6 Ω`; at resonance `D = 1.636`
+(2.14 dBi, against 1.641 for the thin half-wave dipole), the half-power beamwidth is 78.5° (78.1°) and the
+E-plane pattern is within 0.06 dB of `cos(π/2 cos θ)/sin θ` over 30°–150°. The resonant length is still
+rising with refinement near the strip (0.437λ, 0.451λ, 0.456λ as the element size there halves).
+
 ```@docs
 MaxwellModel
 MaxwellMaterial
@@ -107,6 +145,9 @@ MaxwellSweep
 FiniteElementMethod.umfpack_solver
 FiniteElementMethod.tangential_mass
 FiniteElementMethod.port_excitations
+far_field
+radiated_power
+directivity
 FiniteElementMethod.s_parameters
 FiniteElementMethod.propagation_constant
 FiniteElementMethod.cutoff_frequency

examples/dipole/dipole.jl

@@ -0,0 +1,122 @@
+# A strip dipole fed by a lumped port, radiating inside a first-order absorbing sphere: input impedance against
+# frequency, the resonance, and the far field (pattern, directivity, half-power beamwidth) at resonance.
+#
+#     julia --project=<test environment> examples/dipole/dipole.jl
+#
+# Mesh units are metres and the reference frequency is f0 = c / (1 m), so lengths are in wavelengths at f0. The
+# dipole is a zero-thickness PEC strip along z in the plane y = 0, of width w = 4a (the equivalent radius of a
+# flat strip is a quarter of its width), with a gap of length `gap` at its centre carrying a 50 Ω LumpedPort.
+using FiniteElementMethod, Ferrite, LinearAlgebra
+include(joinpath(@__DIR__, "mesh.jl"))
+
+const DIPOLE_F0 = FiniteElementMethod.C0  # Hz: λ = 1 m
+
+"""
+    dipole_model(; length, radius, gap, outer_radius, huygens_radius, h, h_wire = 2radius, growth = 0.3,
+                 max_dofs = 3.0e5)
+
+The [`MaxwellModel`](@ref) of a strip dipole of `length` and equivalent `radius` (width `4radius`) with a
+centre `gap`, in a ball of radius `outer_radius` with the absorbing condition on `"outer"`, a Huygens sphere
+`"huygens"` of radius `huygens_radius`, and element size `h_wire` on the strip growing to `h`. A mesh above
+`max_dofs` edge dofs is refused before anything is assembled on it.
+"""
+function dipole_model(;
+        length, radius, gap, outer_radius, huygens_radius, h, h_wire = 2radius, growth = 0.3, max_dofs = 3.0e5
+    )
+    grid = antenna_grid(;
+        outer_radius, huygens_radius, h, h_fine = h_wire, growth, strip = (; length, width = 4radius, gap)
+    )
+    dh = DofHandler(grid)
+    add!(dh, :E, Nedelec{RefTetrahedron, 1}())
+    close!(dh)
+    ndofs(dh) <= max_dofs || error("the dipole mesh has $(ndofs(dh)) edge dofs, above max_dofs = $max_dofs")
+    return MaxwellModel(
+        grid; pec = ["strip"], absorbing = ["outer"], ports = [LumpedPort("gap"; direction = (0, 0, 1), Z0 = 50.0)]
+    )
+end
+
+# The sinusoidal-current pattern of a thin dipole of electrical half-length kl/2, |E_θ| up to a constant.
+thin_dipole_pattern(θ, half_kl) = abs((cos(half_kl * cos(θ)) - cos(half_kl)) / sin(θ))
+
+# The half-power beamwidth (degrees) of a power pattern `D` sampled at increasing `θs` (rad) on (0, π/2]
+# ending at broadside, by linear interpolation of the half-power crossing.
+function half_power_beamwidth(θs, D)
+    i = findfirst(>=(D[end] / 2), D)
+    (i === nothing || i == 1) && error("the pattern does not cross half power on the sampled range")
+    t = (D[end] / 2 - D[i - 1]) / (D[i] - D[i - 1])
+    return rad2deg(π - 2 * (θs[i - 1] + t * (θs[i] - θs[i - 1])))
+end
+
+"""
+    dipole_analysis(model; frequencies, length)
+
+Input impedance `Z(f)` of the dipole `model` of `length` (mesh units) at `frequencies` (Hz), the resonance (where `Im Z` changes sign,
+linearly interpolated and then solved again there), and at resonance: `Z`, the input and radiated powers, the
+directivity at broadside, the half-power beamwidth and the E-plane pattern (dB, normalised at broadside) on
+`θ = 30°:5°:150°` against `cos(π/2 cos θ)/sin θ` (the half-wave target) and against the sinusoidal-current
+pattern of the dipole's actual electrical length.
+"""
+function dipole_analysis(model; frequencies, length)
+    sweep = solve_sweep(model, 2π .* frequencies)
+    Z = [only(FiniteElementMethod.z_parameters(model, sweep.S[:, :, k])) for k in eachindex(frequencies)]
+    k = findfirst(i -> sign(imag(Z[i])) != sign(imag(Z[i + 1])), 1:(Base.length(Z) - 1))
+    k === nothing && error("no resonance in the sweep: Im Z = $(imag.(Z))")
+    t = imag(Z[k]) / (imag(Z[k]) - imag(Z[k + 1]))
+    f_res = frequencies[k] + t * (frequencies[k + 1] - frequencies[k])
+
+    at_res = solve_sweep(model, [2π * f_res]; fields = true)
+    Z_res = only(FiniteElementMethod.z_parameters(model, at_res.S[:, :, 1]))
+    sol = MaxwellSolution(model, 2π * f_res, at_res.fields[1][:, 1])
+    V = FiniteElementMethod.port_voltage(model, 1, sol.u)
+    I_in = FiniteElementMethod.port_current(model, 1, sol.u; source = 1)
+    P_in = real(V * conj(I_in)) / 2
+    P_flux = radiated_power(sol, "huygens")
+    P_far = radiated_power(sol, "huygens"; method = :far_field)
+
+    E_plane(θ) = Vec(sin(θ), 0.0, cos(θ))
+    D_broadside = directivity(sol, "huygens", [Vec(1.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0)])
+    θs = deg2rad.(30:5:150)
+    pattern = directivity(sol, "huygens", E_plane.(θs))
+    pattern_dB = 10 .* log10.(pattern ./ directivity(sol, "huygens", [E_plane(π / 2)]))
+    half_wave_dB = 20 .* log10.(thin_dipole_pattern.(θs, π / 2))
+    half_kl = π * length * f_res / DIPOLE_F0
+    sinusoidal_dB = 20 .* log10.(thin_dipole_pattern.(θs, half_kl) ./ thin_dipole_pattern(π / 2, half_kl))
+    θ_fine = deg2rad.(1:0.25:90)
+    hpbw = half_power_beamwidth(θ_fine, directivity(sol, "huygens", E_plane.(θ_fine)))
+    return (;
+        frequencies, Z, f_res, resonant_length = length * f_res / DIPOLE_F0, Z_res, P_in, P_flux, P_far,
+        D_broadside, hpbw, θ_deg = rad2deg.(θs), pattern_dB,
+        pattern_err_half_wave = maximum(abs, pattern_dB - half_wave_dB),
+        pattern_err_sinusoidal = maximum(abs, pattern_dB - sinusoidal_dB), ndofs = ndofs(model.dh),
+    )
+end
+
+"Strip dipole meshes (lengths in λ at f0): `coarse` (~2e4 dofs, for the Maxwell test group) and `full` (~1.4e5, 8 GB peak)."
+const DIPOLE_CASES = (
+    coarse = (;
+        length = 0.48, radius = 5.0e-3, gap = 0.02, outer_radius = 0.6, huygens_radius = 0.4, h = 1 / 7,
+        h_wire = 0.0075, frequencies = DIPOLE_F0 .* (0.85:0.05:1.0),
+    ),
+    full = (;
+        length = 0.48, radius = 5.0e-3, gap = 0.01, outer_radius = 1.0, huygens_radius = 0.5, h = 1 / 12,
+        h_wire = 0.00125, growth = 0.4, frequencies = DIPOLE_F0 .* (0.91:0.015:0.97),
+    ),
+)
+
+"""
+    run_dipole(case; max_dofs = 3.0e5)
+
+Build and analyse one of [`DIPOLE_CASES`](@ref) (`:coarse` or `:full`), with the build and analysis times.
+"""
+function run_dipole(case::Symbol; max_dofs = 3.0e5)
+    frequencies = DIPOLE_CASES[case].frequencies
+    mesh = Base.structdiff(DIPOLE_CASES[case], (; frequencies))
+    build = @elapsed model = dipole_model(; mesh..., max_dofs)
+    analysis = @elapsed result = dipole_analysis(model; frequencies, mesh.length)
+    return (; result..., build, analysis)
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    result = run_dipole(Symbol(get(ARGS, 1, "full")))
+    @info "strip dipole" result...
+end

examples/dipole/mesh.jl

@@ -0,0 +1,142 @@
+using Ferrite: Ferrite
+using FerriteGmsh: FerriteGmsh
+using Gmsh: Gmsh, gmsh
+
+"""
+    antenna_grid(; outer_radius, huygens_radius, h, h_fine = h, strip = nothing, source_radius = nothing,
+                 growth = 0.3)
+
+Graded tetrahedral mesh of the ball of radius `outer_radius` about the origin, built with the Gmsh OCC API.
+Physical groups, which become the facet and cell sets of the grid:
+
+  - `"outer"`: the outer sphere (for the absorbing boundary);
+  - `"huygens"`: the closed interior sphere of radius `huygens_radius`;
+  - with `strip = (; length, width, gap)`: a zero-thickness strip dipole in the plane `y = 0`, along `z`,
+    `|x| ≤ width/2`, `|z| ≤ length/2`, as the interior facet sets `"strip"` (the two arms, for `pec`) and
+    `"gap"` (the rectangle `|z| ≤ gap/2` between them, for a `LumpedPort`);
+  - with `source_radius = a`: the cell set `"source"`, the ball of radius `a`.
+
+The element size is `h_fine` on the strip and in the source ball, grows by `growth` per unit distance from
+them (a size `h_fine + growth d` at distance `d`) and is `h` elsewhere.
+"""
+function antenna_grid(;
+        outer_radius, huygens_radius, h, h_fine = h, strip = nothing, source_radius = nothing, growth = 0.3
+    )
+    Gmsh.initialize()
+    try
+        gmsh.option.setNumber("General.Terminal", 0)
+        gmsh.option.setNumber("General.NumThreads", 1)
+        gmsh.model.add("antenna")
+        occ = gmsh.model.occ
+        tools = [(3, occ.addSphere(0, 0, 0, huygens_radius))]
+        source_radius === nothing || push!(tools, (3, occ.addSphere(0, 0, 0, source_radius)))
+        if strip !== nothing
+            (; length, width, gap) = strip
+            arm = (length - gap) / 2
+            # Rectangles in the plane z = 0, turned into y = 0 by a quarter turn about x: (x, y, 0) -> (x, 0, y).
+            rects = [
+                occ.addRectangle(-width / 2, gap / 2, 0, width, arm),
+                occ.addRectangle(-width / 2, -length / 2, 0, width, arm),
+                occ.addRectangle(-width / 2, -gap / 2, 0, width, gap),
+            ]
+            occ.rotate([(2, r) for r in rects], 0, 0, 0, 1, 0, 0, π / 2)
+            append!(tools, [(2, r) for r in rects])
+        end
+        occ.fragment([(3, occ.addSphere(0, 0, 0, outer_radius))], tools)
+        occ.synchronize()
+
+        # Classify the surfaces by their bounding boxes rather than by OCC's tag numbering.
+        surfaces = [tag for (_, tag) in gmsh.model.getEntities(2)]
+        box(tag) = gmsh.model.getBoundingBox(2, tag)
+        extent(tag) = maximum(abs, box(tag))
+        tol = 1.0e-4 * outer_radius
+        sphere(r) = filter(t -> abs(extent(t) - r) < 1.0e-2 * r && abs(box(t)[2]) > tol, surfaces)
+        gmsh.model.addPhysicalGroup(2, sphere(outer_radius), -1, "outer")
+        gmsh.model.addPhysicalGroup(2, sphere(huygens_radius), -1, "huygens")
+        fine_curves = Int[]
+        if strip !== nothing
+            flat = filter(t -> abs(box(t)[2]) < tol && abs(box(t)[5]) < tol, surfaces)
+            gap_surface = filter(t -> box(t)[6] < strip.gap / 2 + tol && box(t)[3] > -strip.gap / 2 - tol, flat)
+            gmsh.model.addPhysicalGroup(2, setdiff(flat, gap_surface), -1, "strip")
+            gmsh.model.addPhysicalGroup(2, gap_surface, -1, "gap")
+            for t in flat
+                append!(fine_curves, abs.(last.(gmsh.model.getBoundary([(2, t)], false, false))))
+            end
+        end
+        if source_radius !== nothing
+            volumes = [tag for (_, tag) in gmsh.model.getEntities(3)]
+            ball = filter(t -> maximum(abs, gmsh.model.getBoundingBox(3, t)) < 1.1 * source_radius, volumes)
+            gmsh.model.addPhysicalGroup(3, ball, -1, "source")
+        end
+        gmsh.model.addPhysicalGroup(3, [tag for (_, tag) in gmsh.model.getEntities(3)], -1, "domain")
+
+        fields = Int[]
+        # Size h_fine within `near` of the field `distance`, growing linearly to h.
+        function refine(distance, near)
+            threshold = gmsh.model.mesh.field.add("Threshold")
+            gmsh.model.mesh.field.setNumber(threshold, "InField", distance)
+            gmsh.model.mesh.field.setNumber(threshold, "SizeMin", h_fine)
+            gmsh.model.mesh.field.setNumber(threshold, "SizeMax", h)
+            gmsh.model.mesh.field.setNumber(threshold, "DistMin", near)
+            gmsh.model.mesh.field.setNumber(threshold, "DistMax", near + (h - h_fine) / growth)
+            return push!(fields, threshold)
+        end
+        if !isempty(fine_curves)
+            distance = gmsh.model.mesh.field.add("Distance")
+            gmsh.model.mesh.field.setNumbers(distance, "CurvesList", unique(fine_curves))
+            gmsh.model.mesh.field.setNumber(distance, "Sampling", 400)
+            refine(distance, 0.0)
+        end
+        if source_radius !== nothing
+            radius = gmsh.model.mesh.field.add("MathEval")
+            gmsh.model.mesh.field.setString(radius, "F", "Sqrt(x^2 + y^2 + z^2)")
+            refine(radius, source_radius)
+        end
+        if isempty(fields)
+            constant = gmsh.model.mesh.field.add("MathEval")
+            gmsh.model.mesh.field.setString(constant, "F", string(h))
+            push!(fields, constant)
+        end
+        background = gmsh.model.mesh.field.add("Min")
+        gmsh.model.mesh.field.setNumbers(background, "FieldsList", fields)
+        gmsh.model.mesh.field.setAsBackgroundMesh(background)
+        gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
+        gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
+        gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
+        gmsh.model.mesh.generate(3)
+        return gmsh_grid()
+    finally
+        Gmsh.finalize()
+    end
+end
+
+"""
+    gmsh_grid()
+
+The `Ferrite.Grid` of the current Gmsh model, as `FerriteGmsh.togrid()` makes it (facet sets from the physical
+groups of dimension 2, listing every cell next to each facet; cell sets from those of dimension 3), but finding
+the facets through a hash of their nodes, in time linear in the mesh size: `togrid` compares every facet with
+every cell (about 17 s at 1.8 × 10⁵ cells, against 0.2 s here).
+"""
+function gmsh_grid()
+    gmsh.model.mesh.renumberNodes()
+    gmsh.model.mesh.renumberElements()
+    # getNodes lists the nodes by entity, not by tag.
+    tags, coordinates = gmsh.model.mesh.getNodes()
+    nodes = Vector{Ferrite.Node{3, Float64}}(undef, length(tags))
+    for (i, tag) in pairs(tags)
+        nodes[tag] = Ferrite.Node(Ferrite.Vec{3}(Tuple(coordinates[(3i - 2):(3i)])))
+    end
+    cells, element_tags = FerriteGmsh.toelements(3)
+    FerriteGmsh.reorient!(cells, nodes)
+    cellsets = FerriteGmsh.tocellsets(3, element_tags)
+    facets_of = Dict{Vector{Int}, Vector{Ferrite.FacetIndex}}()
+    for (i, cell) in pairs(cells), (k, facet) in pairs(Ferrite.facets(cell))
+        push!(get!(facets_of, sort!(collect(facet)), Ferrite.FacetIndex[]), Ferrite.FacetIndex(i, k))
+    end
+    facetsets = Dict{String, Set{Ferrite.FacetIndex}}()
+    for (name, triangles) in FerriteGmsh.toboundary(2)
+        facetsets[name] = Set(f for t in triangles for f in facets_of[sort!(collect(t))])
+    end
+    return Ferrite.Grid(cells, nodes; facetsets, cellsets)
+end

src/FiniteElementMethod.jl

@@ -58,9 +58,11 @@ include("maxwell/ports.jl")
 include("maxwell/model.jl")
 include("maxwell/driven.jl")
 include("maxwell/sweep.jl")
+include("maxwell/farfield.jl")
 export MaxwellMaterial, MaxwellModel, MaxwellTerm, maxwell_operator, maxwell_operator!
 export MaxwellSolution, driven_problem, solve_driven, field_values
 export WavePort, LumpedPort, MaxwellSweep, solve_sweep
+export far_field, radiated_power, directivity
 
 include("precompile.jl")