Cover the JuMP and Gmsh extensions; fix two JuMP-ext bugs the tests found

Extend test_jump.jl (87 tests) and test_gmsh.jl (60 tests) to cover the
previously untested surface: the variable/expression interface (identity,
zero/one/convert, broadcasting, show methods, MA promotion shims), the JuMP
model printer, objective plumbing and its error fallbacks, Uniform()
variables, region-restricted cones (a :power piece inside convex_piecewise),
MOI.SecondOrderCone entry, every affine-sugar operator (solved bit-identically
against the classical API), the convergence-failure status surface, lowering
rejections and the Dirichlet-never-differentiated warning; on the Gmsh side,
gmsh_import from .geo and .msh files (both init-lifecycle branches), subdomain
physical groups on the tensor path, the negative-orientation flip (clockwise
surface), and prism/line-only/empty/missing-file rejections. Measured
coverage: Gmsh ext 249/249 lines, JuMP ext 492/494 (the two rest are internal
invariant guards).

Fixes found by the new tests: vector comparisons with On crashed (the
GreaterThanZero-style marker methods assumed scalar functions; they now branch
vector/scalar like the Zeros/Nonnegatives/Nonpositives methods), and interval
syntax died with an internal JuMP error because variable_ref_type was not
defined for MGBExpr/Coef. Also delete the dead _model helpers.

Author: Sébastien Loisel

ext/MultiGridBarrierJuMPExt/MultiGridBarrierJuMPExt.jl

@@ -197,10 +197,6 @@ end
 
 # --- conversions and arithmetic ---------------------------------------------
 
-_model(x::MGBVarRef) = x.model
-_model(x::Coef) = x.model
-_model(x::MGBExpr) = x.model
-
 _to_expr(x::MGBExpr) = x
 function _to_expr(v::MGBVarRef{T}) where {T}
     e = MGBExpr{T}()
@@ -313,6 +309,12 @@ Base.convert(::Type{MGBExpr{T}}, x::MGBVarRef{T}) where {T} = _to_expr(x)
 Base.convert(::Type{MGBExpr{T}}, x::Coef{T}) where {T} = _to_expr(x)
 Base.convert(::Type{MGBExpr{T}}, r::Real) where {T} = _to_expr(T, r)
 
+# JuMP's macro plumbing asks expression types for their variable type (e.g. when
+# parsing interval constraints lb <= expr <= ub); without these the user gets an
+# internal JuMP error instead of our "unsupported set" rejection.
+JuMP.variable_ref_type(::Type{MGBExpr{T}}) where {T} = MGBVarRef{T}
+JuMP.variable_ref_type(::Type{Coef{T}}) where {T} = MGBVarRef{T}
+
 # ---------------------------------------------------------------------------
 # JuMP model interface: variables
 # ---------------------------------------------------------------------------
@@ -423,15 +425,21 @@ function JuMP.build_constraint(err::Function, f, set::JuMP.Nonpositives, on::On)
         JuMP.build_constraint(err, f, MOI.LessThan(0.0))
     RegionConstraint(inner, on.pairs)
 end
-# Newer JuMP versions route scalar comparisons with non-Number sides through
-# internal *Zero marker sets, again variadically; cover them when they exist.
-for (marker, moiset) in ((:GreaterThanZero, :(MOI.GreaterThan(0.0))),
-                         (:LessThanZero,    :(MOI.LessThan(0.0))),
-                         (:EqualToZero,     :(MOI.EqualTo(0.0))))
+# Newer JuMP versions route comparisons with non-Number sides through internal
+# *Zero marker sets, again variadically; cover them when they exist. Vector
+# comparisons arrive here too (f is then a Vector), so branch like the
+# Zeros/Nonnegatives/Nonpositives methods above.
+for (marker, moiset, vecset) in
+        ((:GreaterThanZero, :(MOI.GreaterThan(0.0)), :(MOI.Nonnegatives)),
+         (:LessThanZero,    :(MOI.LessThan(0.0)),    :(MOI.Nonpositives)),
+         (:EqualToZero,     :(MOI.EqualTo(0.0)),     :(MOI.Zeros)))
     if isdefined(JuMP, marker)
         @eval function JuMP.build_constraint(err::Function, f,
                                              set::JuMP.$marker, on::On)
-            RegionConstraint(JuMP.build_constraint(err, f, $moiset), on.pairs)
+            inner = f isa AbstractVector ?
+                JuMP.build_constraint(err, f, $vecset(length(f))) :
+                JuMP.build_constraint(err, f, $moiset)
+            RegionConstraint(inner, on.pairs)
         end
     end
 end

test/test_gmsh.jl

@@ -158,6 +158,90 @@ end
             @test all(gm.geometry.x[v, e, 1] < 1e-9 for (v, e) in gm.regions["left_half"])
         end
 
+        @testset "subdomain physical group on quads (tensor path)" begin
+            # same split domain, all-quad: subdomain membership on the tensor
+            # families goes through the element-tag branch of _regions_tensor
+            # (every DOF of every element of the group), not node tags.
+            gmsh.clear()
+            r1 = gmsh.model.occ.addRectangle(-1.0, -1.0, 0.0, 1.0, 2.0)
+            r2 = gmsh.model.occ.addRectangle(0.0, -1.0, 0.0, 1.0, 2.0)
+            gmsh.model.occ.fragment([(2, r1)], [(2, r2)])
+            gmsh.model.occ.synchronize()
+            leftsurf = Int[]
+            for (d, t) in gmsh.model.getEntities(2)
+                xmin, _, _, xmax, _, _ = gmsh.model.getBoundingBox(2, t)
+                xmax < 1e-6 && push!(leftsurf, t)
+            end
+            gmsh.model.addPhysicalGroup(2, leftsurf, -1, "left_half")
+            gmsh.option.setNumber("Mesh.MeshSizeMax", 0.5)
+            gmsh.option.setNumber("Mesh.RecombineAll", 1)
+            gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)  # all-quad
+            gmsh.model.mesh.generate(2)
+            gm = gmsh_import()
+            gmsh.option.setNumber("Mesh.RecombineAll", 0)
+            gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 0)
+            @test typeof(gm.geometry.discretization) <: TensorFEM
+            left = gm.regions["left_half"]
+            @test !isempty(left)
+            @test all(gm.geometry.x[v, e, 1] < 1e-9 for (v, e) in left)
+            # whole-element semantics: every DOF of every member element
+            es = unique(e for (v, e) in left)
+            @test length(left) == size(gm.geometry.x, 1) * length(es)
+        end
+
+        @testset "import from a .msh file (gmsh already initialized)" begin
+            mktempdir() do dir
+                gmsh.clear()
+                gmsh.model.occ.addRectangle(-1.0, -1.0, 0.0, 2.0, 2.0)
+                gmsh.model.occ.synchronize()
+                gmsh.model.addPhysicalGroup(1,
+                    [t for (d, t) in gmsh.model.getEntities(1)], -1, "boundary")
+                # gmsh.write only saves elements in physical groups, so the
+                # surface needs one too (or the file would hold bare lines)
+                gmsh.model.addPhysicalGroup(2,
+                    [t for (d, t) in gmsh.model.getEntities(2)], -1, "domain")
+                gmsh.option.setNumber("Mesh.MeshSizeMax", 0.5)
+                gmsh.model.mesh.generate(2)
+                mshpath = joinpath(dir, "square.msh")
+                gmsh.write(mshpath)
+                gmsh.clear()
+                gm = gmsh_import(mshpath)   # meshed file: opened, not re-generated
+                @test typeof(gm.geometry.discretization) <: FEM2D_P1
+                @test length(gm.regions["boundary"]) > 8
+                V, N = size(gm.geometry.x, 1), size(gm.geometry.x, 2)
+                @test length(gm.regions["domain"]) == V * N
+            end
+        end
+
+        @testset "clockwise surface: negative orientation is flipped" begin
+            # geo-kernel curve loop wound clockwise -> surface normal -z ->
+            # negatively-oriented quads; the import must flip them and still
+            # reproduce linears exactly.
+            gmsh.clear()
+            p1 = gmsh.model.geo.addPoint(-1.0, -1.0, 0.0)
+            p2 = gmsh.model.geo.addPoint(1.0, -1.0, 0.0)
+            p3 = gmsh.model.geo.addPoint(1.0, 1.0, 0.0)
+            p4 = gmsh.model.geo.addPoint(-1.0, 1.0, 0.0)
+            l1 = gmsh.model.geo.addLine(p1, p4)
+            l2 = gmsh.model.geo.addLine(p4, p3)
+            l3 = gmsh.model.geo.addLine(p3, p2)
+            l4 = gmsh.model.geo.addLine(p2, p1)
+            loop = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
+            surf = gmsh.model.geo.addPlaneSurface([loop])
+            gmsh.model.geo.synchronize()
+            for l in (l1, l2, l3, l4)
+                gmsh.model.mesh.setTransfiniteCurve(l, 5)
+            end
+            gmsh.model.mesh.setTransfiniteSurface(surf)
+            gmsh.model.mesh.setRecombine(2, surf)
+            gmsh.model.addPhysicalGroup(1, [l1, l2, l3, l4], -1, "boundary")
+            gmsh.model.mesh.generate(2)
+            gm = gmsh_import()
+            @test typeof(gm.geometry.discretization) <: TensorFEM
+            @test abs(sum(gm.geometry.w) - 4.0) < 1e-9   # positive measure
+            @test _lin_gap(gm.geometry, gm.regions["boundary"]) < 1e-6
+        end
+
         @testset "unsupported elements are rejected with hints" begin
             gmsh.clear()
             gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0)
@@ -166,8 +250,56 @@ end
             @test_throws ArgumentError gmsh_import()
             err = try gmsh_import(); "" catch e; sprint(showerror, e) end
             @test occursin("SubdivisionAlgorithm", err)
+
+            # prisms (recombined extrusion of a triangle mesh): generic rejection
+            gmsh.clear()
+            r = gmsh.model.occ.addRectangle(0.0, 0.0, 0.0, 1.0, 1.0)
+            gmsh.model.occ.extrude([(2, r)], 0.0, 0.0, 1.0, [1], [], true)
+            gmsh.model.occ.synchronize()
+            gmsh.option.setNumber("Mesh.MeshSizeMax", 0.6)
+            gmsh.model.mesh.generate(3)
+            @test_throws ArgumentError gmsh_import()
+
+            # a mesh with only line elements, and no mesh at all
+            gmsh.clear()
+            gmsh.model.occ.addRectangle(-1.0, -1.0, 0.0, 2.0, 2.0)
+            gmsh.model.occ.synchronize()
+            gmsh.model.mesh.generate(1)
+            @test_throws ArgumentError gmsh_import()
+            gmsh.clear()
+            @test_throws ArgumentError gmsh_import()
         end
     finally
         gmsh.finalize()
     end
 end
+
+@testset "Gmsh import from a .geo file (gmsh not initialized)" begin
+    # gmsh was finalized above, so gmsh_import(path) owns the
+    # initialize/finalize lifecycle here, and a .geo file forces the
+    # mesh-generation branch.
+    @test_throws ArgumentError gmsh_import("no_such_file_xyz.msh")
+    mktempdir() do dir
+        geopath = joinpath(dir, "square.geo")
+        write(geopath, """
+            Point(1) = {-1, -1, 0, 0.6};
+            Point(2) = {1, -1, 0, 0.6};
+            Point(3) = {1, 1, 0, 0.6};
+            Point(4) = {-1, 1, 0, 0.6};
+            Line(1) = {1, 2};
+            Line(2) = {2, 3};
+            Line(3) = {3, 4};
+            Line(4) = {4, 1};
+            Curve Loop(5) = {1, 2, 3, 4};
+            Plane Surface(6) = {5};
+            Physical Curve("boundary") = {1, 2, 3, 4};
+            Physical Surface("domain") = {6};
+            """)
+        gm = gmsh_import(geopath)
+        @test typeof(gm.geometry.discretization) <: FEM2D_P1
+        V, N = size(gm.geometry.x, 1), size(gm.geometry.x, 2)
+        @test length(gm.regions["boundary"]) > 4
+        # the whole-domain group contains every (vertex, element) pair
+        @test sort(gm.regions["domain"]) == sort([(v, e) for e in 1:N for v in 1:V])
+    end
+end

test/test_jump.jl

@@ -216,4 +216,194 @@ _maxdiff(ref, cols...) =
         @test termination_status(ma) == JuMP.MOI.LOCALLY_SOLVED
         @test maximum(abs.(value(ua) .- value(ub))) < 1e-6
     end
+
+    @testset "variable/expression interface + model printing" begin
+        # The small JuMP/Base interface surface: identity, hashing, zero/one,
+        # convert, broadcasting, string forms, objective plumbing, and the
+        # native-JuMP-style model printer. No solve — this is all modeling-time.
+        mp = MGBModel(geom); _quiet!(mp)
+        @variable(mp, u); @variable(mp, s, Broken())
+        E = typeof(zero(u))                      # MGBExpr{Float64}
+
+        # identity / accessors
+        @test JuMP.num_variables(mp) == 2
+        @test JuMP.owner_model(u) === mp
+        @test JuMP.name(u) == "u"
+        @test JuMP.name(deriv(u, :dx)) == "deriv(u, :dx)"
+        @test deriv(u, :id) == u && u != s
+        @test copy(u) == u
+        @test JuMP.isequal_canonical(u, u)
+        @test Dict(u => 1)[deriv(u, :id)] == 1   # hash + == in real use
+        @test sprint(show, deriv(u, :dy)) == "deriv(u, :dy)"
+        @test occursin("MGBModel over", sprint(show, mp))
+
+        # zero / one / convert / broadcastable
+        @test sprint(show, zero(u)) == "0.0"
+        @test sprint(show, zero(typeof(u))) == "0.0"
+        @test sprint(show, zero(E)) == "0.0" && sprint(show, zero(zero(u))) == "0.0"
+        @test sprint(show, one(typeof(u))) == "1.0"
+        @test sprint(show, convert(E, Coef(mp, 2.5))) == "2.5"
+        @test sprint(show, convert(E, 1.5)) == "1.5"
+        @test (2.0 .* u) isa E && ((u + s) .+ 1.0) isa E && (Coef(mp, 1.0) .* 2.0) isa E
+        es = sprint(show, 2.0 * u + Coef(mp, x -> x[1]) * s + Coef(mp, x -> x[1]))
+        @test occursin("2.0*", es) && occursin("⟨coef⟩", es)
+
+        # MutableArithmetics promotion (the Real variants; the macros need these
+        # on some JuMP versions, so exercise them through the ext's own alias)
+        MA = JuMP._MA
+        @test MA.promote_operation(+, typeof(u), Float64) == E
+        @test MA.promote_operation(*, Float64, typeof(u)) == E
+        @test MA.promote_operation(-, typeof(Coef(mp, 1.0)), Float64) == E
+        @test MA.promote_operation(+, E, Float64) == E
+
+        # objective plumbing: sense accessors, split sense/function route,
+        # and the not-an-integral error fallbacks
+        JuMP.set_objective_sense(mp, JuMP.MOI.MAX_SENSE)
+        @test JuMP.objective_sense(mp) == JuMP.MOI.MAX_SENSE
+        JuMP.set_objective_sense(mp, JuMP.MOI.MIN_SENSE)
+        JuMP.set_objective_function(mp, integral(1.0 * s + 1.0))
+        @test_throws ArgumentError JuMP.set_objective_function(mp, s)
+        @test_throws ArgumentError @objective(mp, Min, s)
+        @objective(mp, Min, integral(1.0 * s + 1.0))
+
+        # printing with no constraints yet, then with every constraint flavor
+        @test occursin("∫(", sprint(print, mp))
+        cr = @constraint(mp, u == Coef(mp, 0.0), On(bd))
+        @test sprint(show, cr) isa String
+        @test_throws ArgumentError JuMP.dual(cr)
+        @constraint(mp, u >= -10.0)
+        @constraint(mp, [deriv(u, :dx); deriv(u, :dy); s] in EpiPower(x -> 1.5 + 0.1 * x[1]^2))
+        @constraint(mp, [s; deriv(u, :dx); deriv(u, :dy)] in JuMP.MOI.SecondOrderCone(3))
+        str = sprint(print, mp)
+        @test occursin("== ⟨data⟩", str) && occursin("on ", str)   # Dirichlet + region
+        @test occursin("≥ 0", str)                                 # nonneg row
+        @test occursin("^⟨coef⟩", str) && occursin("^1.0", str)    # spatial + uniform exponent
+        @test occursin("∫(", str) && occursin("+ 1.0", str)        # objective with constant
+
+        # a fresh model prints "0" for the missing objective
+        m0 = MGBModel(geom)
+        @variable(m0, w0)
+        @test occursin("0", sprint(print, m0))
+
+        # RegionConstraint is a JuMP.AbstractConstraint: function/set accessors
+        rc = JuMP.build_constraint(error, 1.0 * u, JuMP.MOI.GreaterThan(0.0), On(bd))
+        @test JuMP.jump_function(rc) isa E
+        @test JuMP.moi_set(rc) isa JuMP.MOI.GreaterThan{Float64}
+
+        # vector comparisons (with and without On) must land in the same
+        # region/cone machinery as their scalar spellings
+        @constraint(mp, [u, u] >= [Coef(mp, -5.0), Coef(mp, -6.0)], On(bd))
+        @constraint(mp, [u, u] <= [Coef(mp, 5.0), Coef(mp, 6.0)], On(bd))
+        @constraint(mp, [u, u] >= [Coef(mp, -7.0), Coef(mp, -8.0)])
+        @test_throws ArgumentError @constraint(mp, [u, u] == [Coef(mp, 0.0), Coef(mp, 0.0)], On(bd))
+        @test JuMP.num_constraints(mp) == 7
+
+        # rejections with explanatory errors
+        @test_throws ArgumentError u * s                              # nonlinear
+        @test_throws ArgumentError @constraint(mp, -1.0 <= u <= 1.0)  # Interval set
+        @test_throws ArgumentError @constraint(mp, [u, s, u] in JuMP.MOI.ExponentialCone())
+    end
+
+    @testset "lowering guards" begin
+        # a Broken variable cannot carry a Dirichlet condition
+        mb = MGBModel(geom); _quiet!(mb)
+        @variable(mb, ub, Broken()); @variable(mb, sb, Broken())
+        @constraint(mb, ub == Coef(mb, 0.0), On(bd))
+        @constraint(mb, [ub; sb] in EpiPower(2.0))
+        @objective(mb, Min, integral(1.0 * sb))
+        @test_throws ArgumentError optimize!(mb)
+
+        # Dirichlet-but-never-differentiated is legal but warns (likely accident)
+        mw = MGBModel(geom); _quiet!(mw)
+        @variable(mw, uw); @variable(mw, sw, Broken())
+        set_start(sw, 10.0)
+        @constraint(mw, uw == Coef(mw, 0.0), On(bd))
+        @constraint(mw, [uw; sw] in EpiPower(2.0))
+        @objective(mw, Min, integral(1.0 * sw))
+        @test_logs (:warn, r"never differentiated") match_mode = :any optimize!(mw)
+        @test termination_status(mw) == JuMP.MOI.LOCALLY_SOLVED
+    end
+
+    @testset "scalar-arithmetic sugar lowers to the classical problem" begin
+        # Objective spelled with every affine-sugar form (unary +/-, Real on
+        # either side, /); algebraically it is 0.5u + s, i.e. the package's
+        # default p-Laplacian, so the solve must match the classical reference.
+        g1 = x -> x[1]^2 + x[2]^2
+        m = MGBModel(geom); _quiet!(m)
+        @variable(m, u); @variable(m, s, Broken())
+        set_start(u, g1)
+        JuMP.set_start_value(s, 100.0)           # JuMP-native start API delegates
+        @constraint(m, u == Coef(m, g1), On(bd))
+        @constraint(m, [deriv(u, :dx); deriv(u, :dy); s] in EpiPower(1.5))
+        @objective(m, Min,
+            integral(+(u / 2) - (-s) + (u + 1.0) + (1.0 + u) + (1.0 - u) - u - 3.0))
+        optimize!(m)
+        sol_ref = mgb_solve(assemble(amg(geom); p = 1.5); verbose = false)
+        @test maximum(abs.(value(u) .- sol_ref.z[:, 1])) < _JUMP_MATCH_TOL
+        # success-side status accessors
+        @test occursin("converged", JuMP.raw_status(m))
+        @test JuMP.solve_time(m) > 0
+        @test JuMP.primal_status(m) == JuMP.MOI.FEASIBLE_POINT
+        @test JuMP.dual_status(m) == JuMP.MOI.NO_SOLUTION
+    end
+
+    @testset "region-restricted cone + Uniform variable" begin
+        # ROF-style fidelity cone imposed only on the left half (a :power piece
+        # inside convex_piecewise), with a global r >= 0 keeping the slack
+        # bounded off the region.
+        fdata = x -> 0.5 * tanh(5 * x[1])
+        mask = reshape(geom.x, :, 2)[:, 1] .< 0
+        mr = MGBModel(geom); _quiet!(mr)
+        @variable(mr, w); @variable(mr, sw, Broken()); @variable(mr, r, Broken())
+        set_start(w, fdata); set_start(sw, 10.0); set_start(r, 10.0)
+        fd = Coef(mr, fdata)
+        @constraint(mr, w == fd, On(bd))
+        @constraint(mr, [deriv(w, :dx); deriv(w, :dy); sw] in EpiPower(1.0))
+        @constraint(mr, [w - fd; r] in EpiPower(2.0), On(geom, mask))
+        @constraint(mr, r >= 0.0)
+        @objective(mr, Min, integral(sw + Coef(mr, 0.5) * r))
+        optimize!(mr)
+        @test termination_status(mr) == JuMP.MOI.LOCALLY_SOLVED
+        rv = value(r); gap2 = (value(w) .- value(fd)) .^ 2
+        @test maximum(rv[.!mask]) < 1e-3               # cone absent off-region
+        @test all(rv[mask] .>= gap2[mask] .- 1e-8)     # epigraph holds on-region
+
+        # Uniform(): one global constant c with c >= u pointwise; minimizing a
+        # small weight on c drives it to max(u), and its column is constant.
+        g0 = x -> x[1]^2 + x[2]^2
+        mu = MGBModel(geom); _quiet!(mu)
+        @variable(mu, uu); @variable(mu, su, Broken()); @variable(mu, c, Uniform())
+        set_start(uu, g0); set_start(su, 100.0); set_start(c, 10.0)
+        @constraint(mu, uu == Coef(mu, g0), On(bd))
+        @constraint(mu, [deriv(uu, :dx); deriv(uu, :dy); su] in EpiPower(2.0))
+        @constraint(mu, c >= uu)
+        @objective(mu, Min, integral(Coef(mu, 0.5) * uu + su + Coef(mu, 0.1) * c))
+        optimize!(mu)
+        @test termination_status(mu) == JuMP.MOI.LOCALLY_SOLVED
+        cv = value(c)
+        @test maximum(abs.(cv .- cv[1])) < 1e-8        # genuinely one global dof
+        @test cv[1] > maximum(value(uu)) - 1e-6        # dominates u ...
+        @test cv[1] < maximum(value(uu)) + 1e-2        # ... and binds
+    end
+
+    @testset "convergence failure surfaces as OTHER_ERROR" begin
+        # maxit too small for the barrier to reach its target t: optimize! must
+        # catch MGBConvergenceFailure, report it, and refuse to hand out values.
+        g1 = x -> x[1]^2 + x[2]^2
+        mf = MGBModel(geom); _quiet!(mf)
+        set_attribute(mf, "maxit", 2)
+        @variable(mf, u); @variable(mf, s, Broken())
+        set_start(u, g1); set_start(s, 100.0)
+        @constraint(mf, u == Coef(mf, g1), On(bd))
+        @constraint(mf, [deriv(u, :dx); deriv(u, :dy); s] in EpiPower(1.5))
+        @objective(mf, Min, integral(Coef(mf, 0.5) * u + s))
+        optimize!(mf)                                   # must not throw
+        @test termination_status(mf) == JuMP.MOI.OTHER_ERROR
+        @test occursin("onvergence failure", JuMP.raw_status(mf))
+        @test JuMP.primal_status(mf) == JuMP.MOI.NO_SOLUTION
+        @test_throws ArgumentError value(u)
+        @test_throws ArgumentError objective_value(mf)
+        @test_throws ArgumentError mgb_solution(mf)
+        @test_throws ArgumentError solver_log(mf)
+    end
 end