diff --git a/src/GenOpt.jl b/src/GenOpt.jl
index 52f76fb..bb8e003 100644
--- a/src/GenOpt.jl
+++ b/src/GenOpt.jl
@@ -6,6 +6,7 @@
 module GenOpt
 
 include("MOI_wrapper.jl")
+include("bridge.jl")
 include("JuMP_wrapper.jl")
 include("examodels.jl")
 
diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl
index c7081ba..a712cc7 100644
--- a/src/MOI_wrapper.jl
+++ b/src/MOI_wrapper.jl
@@ -13,6 +13,11 @@ end
 Base.copy(array::ContiguousArrayOfVariables) = array
 Base.size(array::ContiguousArrayOfVariables) = array.size
 
+function Base.getindex(A::ContiguousArrayOfVariables, I::Integer...)
+    index = A.offset + Base._to_linear_index(CartesianIndices(A.size), I...)
+    return MOI.VariableIndex(index)
+end
+
 """
     struct Iterator{T}
         values::Vector{T}
@@ -24,11 +29,16 @@ end
 
 Iterator(values::AbstractArray) = Iterator(vec(collect(values)))
 
+Base.length(it::Iterator) = length(it.values)
+
 struct IteratorIndex
     value::Int
 end
 
 Base.copy(i::IteratorIndex) = i
+function Base.isapprox(a::IteratorIndex, b::IteratorIndex; kwargs...)
+    return a.value == b.value
+end
 
 struct FunctionGenerator{F} <: MOI.AbstractVectorFunction
     func::MOI.ScalarNonlinearFunction
@@ -38,14 +48,30 @@ end
 function Base.copy(f::FunctionGenerator{F}) where {F}
     return FunctionGenerator{F}(copy(f.func), f.iterators)
 end
+
+function Base.isapprox(a::FunctionGenerator, b::FunctionGenerator; kwargs...)
+    return isapprox(a.func, b.func; kwargs...) &&
+           length(a.iterators) == length(b.iterators) &&
+           all(
+               isapprox(ai.values, bi.values; kwargs...) for
+               (ai, bi) in zip(a.iterators, b.iterators)
+           )
+end
 function MOI.Utilities.is_canonical(f::FunctionGenerator)
     return MOI.Utilities.is_canonical(f.func)
 end
+
+function MOI.output_dimension(f::FunctionGenerator)
+    return prod(length, f.iterators)
+end
+
 function MOI.Utilities.is_coefficient_type(
-    ::Type{FunctionGenerator{E}},
+    ::Type{<:FunctionGenerator},
     ::Type{T},
-) where {E,T}
-    return MOI.Utilities.is_coefficient_type(E, T)
+) where {T}
+    # Return false so standard MOI bridges (ScalarizeBridge, FlipSignBridge, etc.)
+    # don't try to handle FunctionGenerator. Only FunctionGeneratorBridge should.
+    return false
 end
 
 struct SumGenerator{F} <: MOI.AbstractScalarFunction
diff --git a/src/bridge.jl b/src/bridge.jl
new file mode 100644
index 0000000..c0d5611
--- /dev/null
+++ b/src/bridge.jl
@@ -0,0 +1,247 @@
+# Copyright (c) 2024: Benoît Legat and contributors
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+"""
+    FunctionGeneratorBridge{T,F,S}
+
+Bridge that expands a `FunctionGenerator{F}` constraint into individual
+scalar `F`-in-`S` constraints.
+
+The `FunctionGenerator` contains a template `MOI.ScalarNonlinearFunction`
+with `IteratorIndex` placeholders. This bridge substitutes concrete values
+from the iterators for each combination and converts each expanded expression
+to the type `F` declared by the `FunctionGenerator{F}` type parameter.
+If the expanded expression cannot be converted to `F`, an error is thrown.
+"""
+struct FunctionGeneratorBridge{T,F,S,VS<:MOI.Utilities.VectorLinearSet} <:
+       MOI.Bridges.Constraint.AbstractBridge
+    constraints::Vector{MOI.ConstraintIndex{F,S}}
+    func::FunctionGenerator{F}
+    set::VS
+end
+
+function MOI.Bridges.Constraint.bridge_constraint(
+    ::Type{FunctionGeneratorBridge{T,F,S,VS}},
+    model::MOI.ModelLike,
+    func::FunctionGenerator{F},
+    set::VS,
+) where {T,F,S,VS<:MOI.Utilities.VectorLinearSet}
+    scalar_set = S(zero(T))
+    constraints = MOI.ConstraintIndex{F,S}[]
+    sizes = Tuple(length.(func.iterators))
+    for idx in CartesianIndices(sizes)
+        values = [
+            func.iterators[k].values[idx[k]] for k in eachindex(func.iterators)
+        ]
+        expanded = _expand(func.func, values)
+        scalar_func = _convert(F, expanded)
+        ci = MOI.Utilities.normalize_and_add_constraint(
+            model,
+            scalar_func,
+            scalar_set,
+        )
+        push!(constraints, ci)
+    end
+    return FunctionGeneratorBridge{T,F,S,VS}(constraints, func, set)
+end
+
+function MOI.supports_constraint(
+    ::Type{FunctionGeneratorBridge{T}},
+    ::Type{<:FunctionGenerator},
+    ::Type{<:MOI.Utilities.VectorLinearSet},
+) where {T}
+    return true
+end
+
+function MOI.Bridges.Constraint.concrete_bridge_type(
+    ::Type{FunctionGeneratorBridge{T}},
+    ::Type{FunctionGenerator{F}},
+    VS::Type{<:MOI.Utilities.VectorLinearSet},
+) where {T,F}
+    S = MOI.Utilities.scalar_set_type(VS, T)
+    return FunctionGeneratorBridge{T,F,S,VS}
+end
+
+function MOI.Bridges.added_constraint_types(
+    ::Type{<:FunctionGeneratorBridge{T,F,S}},
+) where {T,F,S}
+    return Tuple{Type,Type}[(F, S)]
+end
+
+function MOI.Bridges.added_constrained_variable_types(
+    ::Type{<:FunctionGeneratorBridge},
+)
+    return Tuple{Type}[]
+end
+
+function MOI.get(
+    bridge::FunctionGeneratorBridge{T,F,S},
+    ::MOI.NumberOfConstraints{F,S},
+) where {T,F,S}
+    return length(bridge.constraints)
+end
+
+function MOI.get(
+    bridge::FunctionGeneratorBridge{T,F,S},
+    ::MOI.ListOfConstraintIndices{F,S},
+) where {T,F,S}
+    return copy(bridge.constraints)
+end
+
+function MOI.get(
+    ::MOI.ModelLike,
+    ::MOI.ConstraintFunction,
+    bridge::FunctionGeneratorBridge,
+)
+    return copy(bridge.func)
+end
+
+function MOI.get(
+    ::MOI.ModelLike,
+    ::MOI.ConstraintSet,
+    bridge::FunctionGeneratorBridge,
+)
+    return bridge.set
+end
+
+function MOI.delete(model::MOI.ModelLike, bridge::FunctionGeneratorBridge)
+    for ci in bridge.constraints
+        MOI.delete(model, ci)
+    end
+    return
+end
+
+# --- Expression expansion ---
+
+"""
+    _expand(func, values)
+
+Recursively expand an `MOI.ScalarNonlinearFunction` template by substituting
+`IteratorIndex` placeholders with concrete `values`.
+
+Evaluates `:getindex` nodes (for variable and data lookups) and simplifies
+pure-numeric subexpressions.
+"""
+function _expand(func::MOI.ScalarNonlinearFunction, values)
+    args = [_expand(arg, values) for arg in func.args]
+    if func.head == :getindex
+        collection = args[1]
+        indices = [_to_index(a) for a in args[2:end]]
+        return getindex(collection, indices...)
+    end
+    if all(_is_numeric, args)
+        return _eval_op(func.head, args)
+    end
+    return MOI.ScalarNonlinearFunction(func.head, args)
+end
+
+_expand(idx::IteratorIndex, values) = values[idx.value]
+_expand(x, _) = x  # constants, MOI.VariableIndex, etc.
+
+_is_numeric(::Number) = true
+_is_numeric(_) = false
+
+_to_index(x::Integer) = Int(x)
+_to_index(x::AbstractFloat) = Int(x)
+_to_index(x) = x
+
+function _eval_op(head::Symbol, args::Vector)
+    registry = MOI.Nonlinear.OperatorRegistry()
+    float_args = Float64.(args)
+    if length(float_args) == 1
+        return MOI.Nonlinear.eval_univariate_function(
+            registry,
+            head,
+            float_args[1],
+        )
+    else
+        return MOI.Nonlinear.eval_multivariate_function(
+            registry,
+            head,
+            float_args,
+        )
+    end
+end
+
+# --- Conversion from expanded ScalarNonlinearFunction to target type F ---
+
+function _convert(
+    ::Type{MOI.ScalarNonlinearFunction},
+    expr::MOI.ScalarNonlinearFunction,
+)
+    return expr
+end
+_convert(::Type{F}, expr::F) where {F} = expr
+
+function _convert(
+    ::Type{MOI.ScalarAffineFunction{T}},
+    expr::MOI.ScalarNonlinearFunction,
+) where {T}
+    terms, constant = _collect_affine_terms(T, expr)
+    if terms === nothing
+        throw(InexactError(:convert, MOI.ScalarAffineFunction{T}, expr))
+    end
+    return MOI.ScalarAffineFunction(terms, T(constant))
+end
+
+function _collect_affine_terms(
+    ::Type{T},
+    expr::MOI.ScalarNonlinearFunction,
+) where {T}
+    if expr.head == :+ && length(expr.args) == 2
+        t1, c1 = _collect_affine_terms(T, expr.args[1])
+        t2, c2 = _collect_affine_terms(T, expr.args[2])
+        (t1 === nothing || t2 === nothing) && return (nothing, zero(T))
+        return (vcat(t1, t2), c1 + c2)
+    elseif expr.head == :- && length(expr.args) == 2
+        t1, c1 = _collect_affine_terms(T, expr.args[1])
+        t2, c2 = _collect_affine_terms(T, expr.args[2])
+        (t1 === nothing || t2 === nothing) && return (nothing, zero(T))
+        neg_t2 = [MOI.ScalarAffineTerm(-t.coefficient, t.variable) for t in t2]
+        return (vcat(t1, neg_t2), c1 - c2)
+    elseif expr.head == :- && length(expr.args) == 1
+        t1, c1 = _collect_affine_terms(T, expr.args[1])
+        t1 === nothing && return (nothing, zero(T))
+        neg_t1 = [MOI.ScalarAffineTerm(-t.coefficient, t.variable) for t in t1]
+        return (neg_t1, -c1)
+    elseif expr.head == :* && length(expr.args) == 2
+        a1, a2 = expr.args
+        if a1 isa Number && a2 isa MOI.VariableIndex
+            return ([MOI.ScalarAffineTerm(T(a1), a2)], zero(T))
+        elseif a2 isa Number && a1 isa MOI.VariableIndex
+            return ([MOI.ScalarAffineTerm(T(a2), a1)], zero(T))
+        elseif a1 isa Number && a2 isa MOI.ScalarNonlinearFunction
+            t2, c2 = _collect_affine_terms(T, a2)
+            t2 === nothing && return (nothing, zero(T))
+            scaled = [
+                MOI.ScalarAffineTerm(T(a1) * t.coefficient, t.variable) for
+                t in t2
+            ]
+            return (scaled, T(a1) * c2)
+        elseif a2 isa Number && a1 isa MOI.ScalarNonlinearFunction
+            t1, c1 = _collect_affine_terms(T, a1)
+            t1 === nothing && return (nothing, zero(T))
+            scaled = [
+                MOI.ScalarAffineTerm(T(a2) * t.coefficient, t.variable) for
+                t in t1
+            ]
+            return (scaled, T(a2) * c1)
+        elseif a1 isa Number && a2 isa Number
+            return (MOI.ScalarAffineTerm{T}[], T(a1 * a2))
+        else
+            return (nothing, zero(T))
+        end
+    else
+        return (nothing, zero(T))
+    end
+end
+
+function _collect_affine_terms(::Type{T}, x::MOI.VariableIndex) where {T}
+    return ([MOI.ScalarAffineTerm(one(T), x)], zero(T))
+end
+
+function _collect_affine_terms(::Type{T}, x::Number) where {T}
+    return (MOI.ScalarAffineTerm{T}[], T(x))
+end
diff --git a/test/Project.toml b/test/Project.toml
index 3b980f2..8351306 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,6 +1,8 @@
 [deps]
 GenOpt = "f2c049d8-7489-4223-990c-4f1c121a4cde"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [sources]
diff --git a/test/bridge.jl b/test/bridge.jl
new file mode 100644
index 0000000..0edadb8
--- /dev/null
+++ b/test/bridge.jl
@@ -0,0 +1,429 @@
+# Copyright (c) 2024: Benoît Legat and contributors
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+module TestBridge
+
+using Test
+import MathOptInterface as MOI
+import GenOpt
+import HiGHS
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$(name)", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+    return
+end
+
+function _create_optimizer()
+    inner = HiGHS.Optimizer()
+    MOI.set(inner, MOI.RawOptimizerAttribute("output_flag"), false)
+    optimizer = MOI.Bridges.full_bridge_optimizer(inner, Float64)
+    MOI.Bridges.add_bridge(optimizer, GenOpt.FunctionGeneratorBridge{Float64})
+    return optimizer
+end
+
+function _affine_sum(vars)
+    return MOI.ScalarAffineFunction(
+        [MOI.ScalarAffineTerm(1.0, xi) for xi in vars],
+        0.0,
+    )
+end
+
+function _minimize!(optimizer, obj)
+    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    return MOI.set(optimizer, MOI.ObjectiveFunction{typeof(obj)}(), obj)
+end
+
+function test_runtests_simple()
+    # x[i] - 1 >= 0 for i in 1..3
+    # normalize_and_add moves constant to set: 1.0*x[i] + 0.0 >= 1.0
+    return MOI.Bridges.runtests(
+        GenOpt.FunctionGeneratorBridge,
+        model -> begin
+            x = MOI.add_variables(model, 3)
+            x_block = GenOpt.ContiguousArrayOfVariables(0, (3,))
+            template = MOI.ScalarNonlinearFunction(
+                :-,
+                Any[
+                    MOI.ScalarNonlinearFunction(
+                        :getindex,
+                        Any[x_block, GenOpt.IteratorIndex(1)],
+                    ),
+                    1.0,
+                ],
+            )
+            iterators = [GenOpt.Iterator([1, 2, 3])]
+            func_gen =
+                GenOpt.FunctionGenerator{MOI.ScalarAffineFunction{Float64}}(
+                    template,
+                    iterators,
+                )
+            MOI.add_constraint(model, func_gen, MOI.Nonnegatives(3))
+        end,
+        model -> begin
+            x = MOI.add_variables(model, 3)
+            for xi in x
+                MOI.add_constraint(
+                    model,
+                    MOI.ScalarAffineFunction(
+                        [MOI.ScalarAffineTerm(1.0, xi)],
+                        0.0,
+                    ),
+                    MOI.GreaterThan(1.0),
+                )
+            end
+        end,
+    )
+end
+
+function test_runtests_equality()
+    # x[i] - 5 == 0 for i in 1..2
+    # normalize_and_add moves constant to set: 1.0*x[i] + 0.0 == 5.0
+    return MOI.Bridges.runtests(
+        GenOpt.FunctionGeneratorBridge,
+        model -> begin
+            x = MOI.add_variables(model, 2)
+            x_block = GenOpt.ContiguousArrayOfVariables(0, (2,))
+            template = MOI.ScalarNonlinearFunction(
+                :-,
+                Any[
+                    MOI.ScalarNonlinearFunction(
+                        :getindex,
+                        Any[x_block, GenOpt.IteratorIndex(1)],
+                    ),
+                    5.0,
+                ],
+            )
+            iterators = [GenOpt.Iterator([1, 2])]
+            func_gen =
+                GenOpt.FunctionGenerator{MOI.ScalarAffineFunction{Float64}}(
+                    template,
+                    iterators,
+                )
+            MOI.add_constraint(model, func_gen, MOI.Zeros(2))
+        end,
+        model -> begin
+            x = MOI.add_variables(model, 2)
+            for xi in x
+                MOI.add_constraint(
+                    model,
+                    MOI.ScalarAffineFunction(
+                        [MOI.ScalarAffineTerm(1.0, xi)],
+                        0.0,
+                    ),
+                    MOI.EqualTo(5.0),
+                )
+            end
+        end,
+    )
+end
+
+function test_runtests_consecutive()
+    # x[i] + x[i+1] - 2 >= 0 for i in 1..2 (3 variables)
+    # normalize_and_add: 1.0*x[i] + 1.0*x[i+1] + 0.0 >= 2.0
+    return MOI.Bridges.runtests(
+        GenOpt.FunctionGeneratorBridge,
+        model -> begin
+            x = MOI.add_variables(model, 3)
+            x_block = GenOpt.ContiguousArrayOfVariables(0, (3,))
+            idx = GenOpt.IteratorIndex(1)
+            idx_plus_1 = MOI.ScalarNonlinearFunction(:+, Any[idx, 1])
+            template = MOI.ScalarNonlinearFunction(
+                :-,
+                Any[
+                    MOI.ScalarNonlinearFunction(
+                        :+,
+                        Any[
+                            MOI.ScalarNonlinearFunction(
+                                :getindex,
+                                Any[x_block, idx],
+                            ),
+                            MOI.ScalarNonlinearFunction(
+                                :getindex,
+                                Any[x_block, idx_plus_1],
+                            ),
+                        ],
+                    ),
+                    2.0,
+                ],
+            )
+            iterators = [GenOpt.Iterator([1, 2])]
+            func_gen =
+                GenOpt.FunctionGenerator{MOI.ScalarAffineFunction{Float64}}(
+                    template,
+                    iterators,
+                )
+            MOI.add_constraint(model, func_gen, MOI.Nonnegatives(2))
+        end,
+        model -> begin
+            x = MOI.add_variables(model, 3)
+            MOI.add_constraint(
+                model,
+                MOI.ScalarAffineFunction(
+                    [
+                        MOI.ScalarAffineTerm(1.0, x[1]),
+                        MOI.ScalarAffineTerm(1.0, x[2]),
+                    ],
+                    0.0,
+                ),
+                MOI.GreaterThan(2.0),
+            )
+            MOI.add_constraint(
+                model,
+                MOI.ScalarAffineFunction(
+                    [
+                        MOI.ScalarAffineTerm(1.0, x[2]),
+                        MOI.ScalarAffineTerm(1.0, x[3]),
+                    ],
+                    0.0,
+                ),
+                MOI.GreaterThan(2.0),
+            )
+        end,
+    )
+end
+
+function test_expand_constant()
+    func = MOI.ScalarNonlinearFunction(:+, Any[GenOpt.IteratorIndex(1), 1.0])
+    result = GenOpt._expand(func, [5.0])
+    @test result == 6.0
+end
+
+function test_expand_variable()
+    x = GenOpt.ContiguousArrayOfVariables(0, (3,))
+    index_expr =
+        MOI.ScalarNonlinearFunction(:+, Any[GenOpt.IteratorIndex(1), 1])
+    func = MOI.ScalarNonlinearFunction(:getindex, Any[x, index_expr])
+    result = GenOpt._expand(func, [1])
+    @test result == MOI.VariableIndex(2)
+end
+
+function test_expand_with_variable_in_expr()
+    func = MOI.ScalarNonlinearFunction(
+        :+,
+        Any[MOI.VariableIndex(1), GenOpt.IteratorIndex(1)],
+    )
+    result = GenOpt._expand(func, [3.0])
+    @test result isa MOI.ScalarNonlinearFunction
+    @test result.head == :+
+    @test result.args[1] == MOI.VariableIndex(1)
+    @test result.args[2] == 3.0
+end
+
+function test_convert_to_affine()
+    x1 = MOI.VariableIndex(1)
+    # x1 - 1.0 should convert to ScalarAffineFunction
+    func = MOI.ScalarNonlinearFunction(:-, Any[x1, 1.0])
+    result = GenOpt._convert(MOI.ScalarAffineFunction{Float64}, func)
+    @test result isa MOI.ScalarAffineFunction{Float64}
+    @test length(result.terms) == 1
+    @test result.terms[1].coefficient == 1.0
+    @test result.terms[1].variable == x1
+    @test result.constant == -1.0
+end
+
+function test_simple_constraint_group()
+    # min sum(x)  s.t.  x[i] >= 1 for i in 1..3
+    optimizer = _create_optimizer()
+    x = MOI.add_variables(optimizer, 3)
+    for xi in x
+        MOI.add_constraint(optimizer, xi, MOI.GreaterThan(0.0))
+    end
+
+    x_block = GenOpt.ContiguousArrayOfVariables(0, (3,))
+    template = MOI.ScalarNonlinearFunction(
+        :-,
+        Any[
+            MOI.ScalarNonlinearFunction(
+                :getindex,
+                Any[x_block, GenOpt.IteratorIndex(1)],
+            ),
+            1.0,
+        ],
+    )
+    iterators = [GenOpt.Iterator([1, 2, 3])]
+    func_gen = GenOpt.FunctionGenerator{MOI.ScalarAffineFunction{Float64}}(
+        template,
+        iterators,
+    )
+    MOI.add_constraint(optimizer, func_gen, MOI.Nonnegatives(3))
+
+    _minimize!(optimizer, _affine_sum(x))
+    MOI.optimize!(optimizer)
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+
+    for xi in x
+        @test MOI.get(optimizer, MOI.VariablePrimal(), xi) ≈ 1.0 atol = 1e-6
+    end
+end
+
+function test_consecutive_constraint_group()
+    # min sum(x)  s.t.  x[i] + x[i+1] >= 2 for i in 1..4, x[i] >= 0
+    optimizer = _create_optimizer()
+    x = MOI.add_variables(optimizer, 5)
+    for xi in x
+        MOI.add_constraint(optimizer, xi, MOI.GreaterThan(0.0))
+    end
+
+    x_block = GenOpt.ContiguousArrayOfVariables(0, (5,))
+    idx = GenOpt.IteratorIndex(1)
+    idx_plus_1 = MOI.ScalarNonlinearFunction(:+, Any[idx, 1])
+    template = MOI.ScalarNonlinearFunction(
+        :-,
+        Any[
+            MOI.ScalarNonlinearFunction(
+                :+,
+                Any[
+                    MOI.ScalarNonlinearFunction(:getindex, Any[x_block, idx]),
+                    MOI.ScalarNonlinearFunction(
+                        :getindex,
+                        Any[x_block, idx_plus_1],
+                    ),
+                ],
+            ),
+            2.0,
+        ],
+    )
+    iterators = [GenOpt.Iterator([1, 2, 3, 4])]
+    func_gen = GenOpt.FunctionGenerator{MOI.ScalarAffineFunction{Float64}}(
+        template,
+        iterators,
+    )
+    MOI.add_constraint(optimizer, func_gen, MOI.Nonnegatives(4))
+
+    _minimize!(optimizer, _affine_sum(x))
+    MOI.optimize!(optimizer)
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+
+    for i in 1:4
+        xi = MOI.get(optimizer, MOI.VariablePrimal(), x[i])
+        xi1 = MOI.get(optimizer, MOI.VariablePrimal(), x[i+1])
+        @test xi + xi1 >= 2.0 - 1e-6
+    end
+end
+
+function test_parameter_data_lookup()
+    # min sum(x)  s.t.  x[i] >= demand[i] for i in 1..3
+    # demand = [1.0, 2.0, 3.0]
+    optimizer = _create_optimizer()
+    x = MOI.add_variables(optimizer, 3)
+    for xi in x
+        MOI.add_constraint(optimizer, xi, MOI.GreaterThan(0.0))
+    end
+
+    demand = [1.0, 2.0, 3.0]
+    x_block = GenOpt.ContiguousArrayOfVariables(0, (3,))
+    idx = GenOpt.IteratorIndex(1)
+    template = MOI.ScalarNonlinearFunction(
+        :-,
+        Any[
+            MOI.ScalarNonlinearFunction(:getindex, Any[x_block, idx]),
+            MOI.ScalarNonlinearFunction(:getindex, Any[demand, idx]),
+        ],
+    )
+    iterators = [GenOpt.Iterator([1, 2, 3])]
+    func_gen = GenOpt.FunctionGenerator{MOI.ScalarAffineFunction{Float64}}(
+        template,
+        iterators,
+    )
+    MOI.add_constraint(optimizer, func_gen, MOI.Nonnegatives(3))
+
+    _minimize!(optimizer, _affine_sum(x))
+    MOI.optimize!(optimizer)
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+
+    for (i, d) in enumerate(demand)
+        @test MOI.get(optimizer, MOI.VariablePrimal(), x[i]) ≈ d atol = 1e-6
+    end
+end
+
+function test_multidim_constraint_group()
+    # x[2*(i-1) + j] >= 1 for i in 1..2, j in 1..2  (4 constraints, 4 variables)
+    optimizer = _create_optimizer()
+    x = MOI.add_variables(optimizer, 4)
+    for xi in x
+        MOI.add_constraint(optimizer, xi, MOI.GreaterThan(0.0))
+    end
+
+    x_block = GenOpt.ContiguousArrayOfVariables(0, (4,))
+    idx_expr = MOI.ScalarNonlinearFunction(
+        :+,
+        Any[
+            MOI.ScalarNonlinearFunction(
+                :*,
+                Any[
+                    2,
+                    MOI.ScalarNonlinearFunction(
+                        :-,
+                        Any[GenOpt.IteratorIndex(1), 1],
+                    ),
+                ],
+            ),
+            GenOpt.IteratorIndex(2),
+        ],
+    )
+    template = MOI.ScalarNonlinearFunction(
+        :-,
+        Any[
+            MOI.ScalarNonlinearFunction(:getindex, Any[x_block, idx_expr]),
+            1.0,
+        ],
+    )
+    iterators = [GenOpt.Iterator([1, 2]), GenOpt.Iterator([1, 2])]
+    func_gen = GenOpt.FunctionGenerator{MOI.ScalarAffineFunction{Float64}}(
+        template,
+        iterators,
+    )
+    MOI.add_constraint(optimizer, func_gen, MOI.Nonnegatives(4))
+
+    _minimize!(optimizer, _affine_sum(x))
+    MOI.optimize!(optimizer)
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+
+    for xi in x
+        @test MOI.get(optimizer, MOI.VariablePrimal(), xi) ≈ 1.0 atol = 1e-6
+    end
+end
+
+function test_equality_constraint_group()
+    # x[i] == 5 for i in 1..3
+    optimizer = _create_optimizer()
+    x = MOI.add_variables(optimizer, 3)
+
+    x_block = GenOpt.ContiguousArrayOfVariables(0, (3,))
+    template = MOI.ScalarNonlinearFunction(
+        :-,
+        Any[
+            MOI.ScalarNonlinearFunction(
+                :getindex,
+                Any[x_block, GenOpt.IteratorIndex(1)],
+            ),
+            5.0,
+        ],
+    )
+    iterators = [GenOpt.Iterator([1, 2, 3])]
+    func_gen = GenOpt.FunctionGenerator{MOI.ScalarAffineFunction{Float64}}(
+        template,
+        iterators,
+    )
+    MOI.add_constraint(optimizer, func_gen, MOI.Zeros(3))
+
+    _minimize!(optimizer, _affine_sum(x))
+    MOI.optimize!(optimizer)
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+
+    for xi in x
+        @test MOI.get(optimizer, MOI.VariablePrimal(), xi) ≈ 5.0 atol = 1e-6
+    end
+end
+
+end  # module
+
+TestBridge.runtests()