H(curl)/H(div) fields cannot be plotted: transfer_solution uses PointValues, which rejects Piola mappings

[koehlerson]

FEData cannot resolve fields whose interpolation uses a non-identity mapping, so no H(curl) (Nedelec) or H(div) (RaviartThomas, BrezziDouglasMarini) field can be plotted by name:

ds = FEData(dh, u)          # dh has a Nedelec field :E
FerriteViz.point_data(ds, :E)
# ERROR: ArgumentError: The cell::AbstractCell input is required to reinit!
#        non-identity function mappings

which means solutionplot, arrowplot, surfaceplot and WarpByVector all fail on such a field.

Affected

Checked against Ferrite 1.5.0 (also reproduces on the 1.3 this package pins), 2D and 3D, all orders:

Interpolation	Mapping	point_data(ds, :u)
Lagrange{RefTriangle,1}()^2	Identity	OK
DiscontinuousLagrange{RefTriangle,1}()^2	Identity	OK
Nedelec{RefTriangle,1}()	Covariant Piola	FAIL
Nedelec{RefTriangle,2}()	Covariant Piola	FAIL
Nedelec{RefTetrahedron,1}()	Covariant Piola	FAIL
RaviartThomas{RefTriangle,1}()	Contravariant Piola	FAIL
RaviartThomas{RefTriangle,2}()	Contravariant Piola	FAIL
RaviartThomas{RefTetrahedron,1}()	Contravariant Piola	FAIL
BrezziDouglasMarini{RefTriangle,1}()	Contravariant Piola	FAIL

What still works

On a DofHandler that contains a Nédélec field:

FEData construction / geometry (meshplot)   OK
point_data of a Lagrange field :p           OK     <- other fields unaffected
point_data of the Nedelec field :E          FAIL
set_point_data! + filters (Magnitude, ...)  OK     <- a manual transfer works
set_cell_data! (cellplot)                   OK
WarpByVector(:E)                            FAIL   <- same root cause
Gradient(:E)                                FAIL   <- different cause, see below

Root cause

transfer_solution evaluates the field with Ferrite.PointValues:

https://github.com/Ferrite-FEM/FerriteViz.jl/blob/master/src/dataset.jl

pv = Ferrite.PointValues(ip_field, ip_geo; update_gradients=false)
...
Ferrite.reinit!(pv, local_coords, ξ)

The covariant/contravariant Piola pullback needs the cell to fix edge/facet orientations and signs, and PointValues has no reinit! overload taking one — Ferrite deliberately rejects it. CellValues does accept the cell, via reinit!(cv, getcells(grid, cellid), coords).

Suggested fix

Branch in transfer_solution on Ferrite.mapping_type(ip_field) !== Ferrite.IdentityMapping() and take a CellValues path for the mapped case: build a QuadratureRule whose points are the tessellation's reference coordinates (they are identical for every cell of a given reference shape) and evaluate cell by cell. One fix covers point_data, solutionplot, arrowplot, surfaceplot and WarpByVector.

A working prototype is transfer_hcurl in the reproducer below.

Separate gap: Gradient

Gradient(:E) fails with a different error — MethodError: no method matching get_gradient_interpolation(...) — because src/gradient.jl has no gradient interpolation for these spaces. Fixing transfer_solution will not make Gradient work on them. Arguably these spaces want curl/div filters rather than a full gradient, so this may be worth documenting as a limitation rather than implementing.

Reproducer

Adapted from the Ferrite tutorial Maxwell: the good, the bad and the ugly — the curl-curl problem on a rotated L-shape, solved with both Nédélec and nodal Lagrange, then visualised with FerriteViz + GLMakie. It reproduces the tutorial's result (Nédélec converges, Lagrange stalls at the corner singularity):

 n |    Nedelec |   Lagrange     L2 error of E
 4 |    0.06307 |      0.156
 8 |    0.04058 |    0.12462
16 |    0.02592 |    0.10444
32 |    0.01648 |    0.09274

The Nédélec field is rendered through the manual transfer_hcurl workaround, not through transfer_solution — it is registered under :E_hcurl because set_point_data! refuses names that shadow a dof field. That is the workaround this issue asks to make unnecessary.

maxwell-nedelec.jl (self-contained, no Gmsh / FerriteAssembly)

# Nédélec (H(curl)) elements — visual check with FerriteViz + GLMakie
#
# Adapted from the Ferrite tutorial "Maxwell: the good, the bad and the ugly"
# https://ferrite-fem.github.io/Ferrite.jl/previews/PR798/tutorials/maxwell_good_bad_ugly/
#
# The problem is the curl-curl system on a rotated L-shaped domain,
#
#     curl(curl(E)) = 0   in Ω,      div(E) = 0   in Ω,      E ⋅ t = g   on Γ,
#
# whose exact solution is the gradient of the classic reentrant-corner potential
# r^(2/3) sin(2θ/3). The point of the tutorial: a nodal Lagrange discretization
# satisfies the equations but converges to the *wrong* function near the corner,
# while the H(curl)-conforming Nédélec space converges properly.
#
# Deviations from the upstream tutorial, so this runs on the Ferrite version
# this package is tested against (1.3) with no extra dependencies:
#
#   * the L-shape is built from a structured triangle grid with one quadrant
#     removed, instead of Gmsh + FerriteGmsh;
#   * assembly is written out by hand instead of using FerriteAssembly;
#   * `WeakDirichlet` does not exist in Ferrite 1.3, and `Dirichlet` /
#     `apply_analytical!` do not support Nédélec (no `reference_coordinates`),
#     so the tangential boundary condition is imposed with a penalty term.
#
# GLMakie is not a dependency of this package or of its docs environment, so run
# this from an environment that has it, e.g. from the repository root:
#
#   julia -e 'using Pkg; Pkg.activate(; temp=true); Pkg.develop(path="."); \
#             Pkg.add(["GLMakie", "Ferrite", "Tensors"]); \
#             include("docs/src/ferrite-examples/maxwell-nedelec.jl")'

using Ferrite, Tensors, LinearAlgebra, SparseArrays
using FerriteViz
using FerriteViz: FEData, Component, Magnitude
import GLMakie

#########
# Grid  #
#########

# The rotated L: the square [-1,1]² minus the quadrant x>0, y<0, so the
# reentrant corner sits at the origin.
function lshape_grid(n::Int)
    base = generate_grid(Triangle, (2n, 2n), Vec(-1.0, -1.0), Vec(1.0, 1.0))
    keep = Int[]
    for (cid, cell) in enumerate(Ferrite.getcells(base))
        c = sum(Ferrite.get_node_coordinate(Ferrite.getnodes(base)[i]) for i in cell.nodes) / 3
        (c[1] > 0 && c[2] < 0) || push!(keep, cid)
    end
    used = falses(Ferrite.getnnodes(base))
    for cid in keep, i in Ferrite.getcells(base, cid).nodes
        used[i] = true
    end
    newid = zeros(Int, Ferrite.getnnodes(base))
    nodes = Ferrite.Node{2,Float64}[]
    for i in 1:Ferrite.getnnodes(base)
        used[i] || continue
        push!(nodes, Ferrite.getnodes(base)[i])
        newid[i] = length(nodes)
    end
    cells = [Triangle(ntuple(k -> newid[Ferrite.getcells(base, cid).nodes[k]], 3)) for cid in keep]
    grid = Grid(cells, nodes)
    top = ExclusiveTopology(grid)
    addboundaryfacetset!(grid, top, "vertical",   x -> abs((x[1] - 1) * x[1] * (x[1] + 1)) ≤ 1e-6)
    addboundaryfacetset!(grid, top, "horizontal", x -> abs((x[2] - 1) * x[2] * (x[2] + 1)) ≤ 1e-6)
    addfacetset!(grid, "boundary", union(getfacetset(grid, "vertical"), getfacetset(grid, "horizontal")))
    return grid
end

##################
# Exact solution #
##################

function exact_potential(x::Vec{2})
    Δθ = -3π / 4
    xp = rotate(x, Δθ)
    r = sqrt(x ⋅ x + eps())
    θ = r ≤ 1e-6 ? zero(eltype(x)) : (atan(xp[2], xp[1]) - Δθ)
    return r^(2 // 3) * sin(2θ / 3)
end
exact_solution(x::Vec{2}) = Tensors.gradient(exact_potential, x)

##########
# Solves #
##########

# Mixed (E, ϕ): curl-curl for E with a Lagrange multiplier ϕ enforcing div(E) = 0.
function solve_nedelec(grid; γ = 1e3)
    ipE = Ferrite.Nedelec{RefTriangle,1}()
    ipϕ = Lagrange{RefTriangle,1}()
    dh = DofHandler(grid)
    add!(dh, :E, ipE); add!(dh, :ϕ, ipϕ); close!(dh)

    ipg = Ferrite.geometric_interpolation(Triangle)
    qr = QuadratureRule{RefTriangle}(2)
    cvE, cvϕ = CellValues(qr, ipE, ipg), CellValues(qr, ipϕ, ipg)
    fvE = FacetValues(FacetQuadratureRule{RefTriangle}(2), ipE, ipg)

    K = allocate_matrix(dh); f = zeros(ndofs(dh))
    asm = start_assemble(K, f)
    nd = ndofs_per_cell(dh)
    Ke, fe = zeros(nd, nd), zeros(nd)
    sdh = dh.subdofhandlers[1]
    rE, rϕ = dof_range(sdh, :E), dof_range(sdh, :ϕ)
    bset = getfacetset(grid, "boundary")

    for cell in CellIterator(dh)
        fill!(Ke, 0); fill!(fe, 0)
        reinit!(cvE, cell); reinit!(cvϕ, cell)
        for q in 1:getnquadpoints(cvE)
            dΩ = getdetJdV(cvE, q)
            for (i, I) in pairs(rE)
                cδE, δE = shape_curl(cvE, q, i), shape_value(cvE, q, i)
                for (j, J) in pairs(rE)
                    Ke[I, J] += (cδE ⋅ shape_curl(cvE, q, j)) * dΩ
                end
                for (j, J) in pairs(rϕ)
                    g = shape_gradient(cvϕ, q, j)
                    Ke[I, J] += (δE ⋅ g) * dΩ
                    Ke[J, I] += (δE ⋅ g) * dΩ
                end
            end
        end
        # Penalty enforcement of the tangential trace E ⋅ t = g ⋅ t. Ferrite 1.3
        # has no WeakDirichlet, and the Nédélec edge dofs are tangential moments
        # that `Dirichlet` cannot interpolate — but FacetValues applies the
        # covariant Piola mapping for us, so a penalty term needs no knowledge
        # of the dof convention.
        for fi in 1:nfacets(getcells(grid, cellid(cell)))
            (cellid(cell), fi) in bset || continue
            reinit!(fvE, cell, fi)
            for q in 1:getnquadpoints(fvE)
                dΓ = getdetJdV(fvE, q)
                n = getnormal(fvE, q)
                t = Vec(-n[2], n[1])
                h = sqrt(dΓ)
                gt = exact_solution(spatial_coordinate(fvE, q, getcoordinates(cell))) ⋅ t
                for (i, I) in pairs(rE)
                    δt = shape_value(fvE, q, i) ⋅ t
                    fe[I] += γ / h * gt * δt * dΓ
                    for (j, J) in pairs(rE)
                        Ke[I, J] += γ / h * (shape_value(fvE, q, j) ⋅ t) * δt * dΓ
                    end
                end
            end
        end
        assemble!(asm, celldofs(cell), Ke, fe)
    end

    ch = ConstraintHandler(dh)
    add!(ch, Dirichlet(:ϕ, getfacetset(grid, "boundary"), Returns(0.0)))
    close!(ch); apply!(K, f, ch)
    return dh, K \ f
end

# The "bad": nodal Lagrange, curl-curl + div-div, components fixed on the boundary.
function solve_lagrange(grid)
    ip = Lagrange{RefTriangle,1}()^2
    dh = close!(add!(DofHandler(grid), :E, ip))
    cv = CellValues(QuadratureRule{RefTriangle}(2), ip, Ferrite.geometric_interpolation(Triangle))
    K = allocate_matrix(dh); f = zeros(ndofs(dh))
    asm = start_assemble(K, f)
    nd = ndofs_per_cell(dh); Ke = zeros(nd, nd)
    for cell in CellIterator(dh)
        fill!(Ke, 0); reinit!(cv, cell)
        for q in 1:getnquadpoints(cv)
            dΩ = getdetJdV(cv, q)
            for i in 1:getnbasefunctions(cv)
                dδ, cδ = shape_divergence(cv, q, i), shape_curl(cv, q, i)
                for j in 1:getnbasefunctions(cv)
                    Ke[i, j] += (cδ ⋅ shape_curl(cv, q, j) + dδ * shape_divergence(cv, q, j)) * dΩ
                end
            end
        end
        assemble!(asm, celldofs(cell), Ke, zeros(nd))
    end
    ch = ConstraintHandler(dh)
    add!(ch, Dirichlet(:E, getfacetset(grid, "horizontal"), (x, _) -> exact_solution(x)[2], [2]))
    add!(ch, Dirichlet(:E, getfacetset(grid, "vertical"),   (x, _) -> exact_solution(x)[1], [1]))
    close!(ch); apply!(K, f, ch)
    return dh, K \ f
end

function l2error(dh, a, fieldname, ip)
    cv = CellValues(QuadratureRule{RefTriangle}(4), ip, Ferrite.geometric_interpolation(Triangle))
    r = dof_range(dh.subdofhandlers[1], fieldname)
    err = vol = 0.0
    for cell in CellIterator(dh)
        reinit!(cv, cell); ae = a[celldofs(cell)]
        for q in 1:getnquadpoints(cv)
            dΩ = getdetJdV(cv, q)
            x = spatial_coordinate(cv, q, getcoordinates(cell))
            err += norm(function_value(cv, q, ae, r) - exact_solution(x))^2 * dΩ
            vol += dΩ
        end
    end
    return sqrt(err) / vol
end

###########################################
# Transferring an H(curl) field to FEData #
###########################################

# FerriteViz cannot resolve a Nédélec field by name: `point_data(ds, :E)` goes
# through `transfer_solution`, which evaluates the field with
# `Ferrite.PointValues` and `reinit!(pv, coords, ξ)` — and Ferrite refuses that
# for any non-identity mapping:
#
#     ArgumentError: The cell::AbstractCell input is required to reinit!
#                    non-identity function mappings
#
# The covariant Piola mapping needs the cell itself to fix the edge
# orientations, and `PointValues` has no `reinit!` overload taking one.
# `CellValues` does, so we build a quadrature rule whose points *are* the
# tessellation's reference coordinates and evaluate cell by cell, then register
# the result as an ordinary named point-data array.
function transfer_hcurl(ds::FEData, dh, a, fieldname::Symbol, ip)
    grid = Ferrite.get_grid(dh)
    r = dof_range(dh.subdofhandlers[1], fieldname)
    # every cell of a single-refshape grid carries the same reference points
    ξs = [Vec{2}(ntuple(d -> ds.reference_coords[v, d], 2))
          for v in FerriteViz.vertices_on_cell(ds, 1)]
    qr = QuadratureRule{RefTriangle}(zeros(length(ξs)), ξs)
    cv = CellValues(qr, ip, Ferrite.geometric_interpolation(Triangle))
    out = fill(NaN, FerriteViz.num_vertices(ds), 2)
    dofs = zeros(Int, ndofs_per_cell(dh))
    for cid in 1:getncells(grid)
        reinit!(cv, getcells(grid, cid), getcoordinates(grid, cid))
        celldofs!(dofs, dh, cid)
        ae = a[dofs]
        for (q, v) in enumerate(FerriteViz.vertices_on_cell(ds, cid))
            E = function_value(cv, q, ae, r)
            out[v, 1], out[v, 2] = E[1], E[2]
        end
    end
    return out
end

##########
# Driver #
##########

grid = lshape_grid(16)
@info "L-shaped grid" cells = getncells(grid) nodes = getnnodes(grid)

dh_ned, a_ned = solve_nedelec(grid)
dh_lag, a_lag = solve_lagrange(grid)

# reference: the exact field sampled into a discontinuous Lagrange space
dh_exact = close!(add!(DofHandler(grid), :E, DiscontinuousLagrange{RefTriangle,1}()^2))
a_exact = zeros(ndofs(dh_exact))
apply_analytical!(a_exact, dh_exact, :E, exact_solution)

@info "L2 error of E" nedelec = l2error(dh_ned, a_ned, :E, Ferrite.Nedelec{RefTriangle,1}()) lagrange = l2error(dh_lag, a_lag, :E, Lagrange{RefTriangle,1}()^2)

# --- datasets -------------------------------------------------------------
ds_exact = FEData(dh_exact, a_exact)
ds_lag   = FEData(dh_lag, a_lag)
ds_ned   = FEData(dh_ned, a_ned)

# The Lagrange and exact fields resolve natively; the Nédélec one is registered
# under its own name because `:E` would shadow the (unusable) dof field.
FerriteViz.set_point_data!(ds_ned, :E_hcurl, transfer_hcurl(ds_ned, dh_ned, a_ned, :E, Ferrite.Nedelec{RefTriangle,1}()))

panels = [
    ("Exact",    ds_exact |> Component(1; input=:E,       output=:E1)),
    ("Lagrange", ds_lag   |> Component(1; input=:E,       output=:E1)),
    ("Nédélec",  ds_ned   |> Component(1; input=:E_hcurl, output=:E1)),
]

# Same colormap and colorrange as the upstream tutorial, so the top row can be
# compared against its figure directly.
cmap = GLMakie.Makie.wong_colors()
crange = (-2.0, 0.0)

# --- figure ---------------------------------------------------------------
fig = GLMakie.Figure(size = (1400, 780))

for (j, (name, d)) in enumerate(panels)
    ax = GLMakie.Axis(fig[1, j]; aspect = GLMakie.DataAspect(), title = "$name — E₁", xlabel = "x₁", ylabel = "x₂")
    FerriteViz.solutionplot!(ax, d; color = :E1, colormap = cmap, colorrange = crange)
    FerriteViz.meshplot!(ax, d; plotnodes = false, linewidth = 0.4, color = (:black, 0.25))
end
GLMakie.Colorbar(fig[1, 4]; colormap = cmap, colorrange = crange, label = "E₁")

# second row: magnitude of the Nédélec field, and the field itself as arrows
ax_mag = GLMakie.Axis(fig[2, 1]; aspect = GLMakie.DataAspect(), title = "Nédélec — |E|", xlabel = "x₁", ylabel = "x₂")
mag = ds_ned |> Magnitude(input = :E_hcurl)
FerriteViz.solutionplot!(ax_mag, mag; color = :magnitude, colormap = :inferno)

ax_arr = GLMakie.Axis(fig[2, 2]; aspect = GLMakie.DataAspect(), title = "Nédélec — E", xlabel = "x₁", ylabel = "x₂")
FerriteViz.meshplot!(ax_arr, ds_ned; plotnodes = false, linewidth = 0.4, color = (:black, 0.2))
FerriteViz.arrowplot!(ax_arr, ds_ned; field = :E_hcurl, lengthscale = 0.06, normalize = true)

ax_ex = GLMakie.Axis(fig[2, 3]; aspect = GLMakie.DataAspect(), title = "Exact — E", xlabel = "x₁", ylabel = "x₂")
FerriteViz.meshplot!(ax_ex, ds_exact; plotnodes = false, linewidth = 0.4, color = (:black, 0.2))
FerriteViz.arrowplot!(ax_ex, ds_exact; field = :E, lengthscale = 0.06, normalize = true)

GLMakie.Label(fig[0, 1:4],
    "curl-curl on the L-shape: nodal Lagrange misses the corner singularity, Nédélec captures it";
    fontsize = 18, font = :bold)

display(fig)
fig