Wrong result when multiplying integrand with scalar when using `InterfaceTriangulation`

[simonsticko]

When integrating over an InterfaceTriangulation (see the run()-function below) the value of an integral changes from 1 to 0 if I multiply the integrand with a scalar 1. However, this only happens on one side of the interface. Can someone confirm if this is a bug?

(Sorry about the somewhat long example. I was not able to set up an InterfaceTriangulation without gmsh)

Example

using Gridap

import Gmsh: gmsh

using GridapGmsh

"""
Set up a mesh, as illustrated below, and write it to file with the incoming
filename.
            * ------ * -------- *
            | "Left"  | "Right" |
            * ------- * ------- *
"""
function create_mesh(filename)

    gmsh.initialize()

    dim2 = 2

    # Define the dimensions of the rectangle
    x_min = -1
    x_interface = 0
    x_max = 1.0
    y_min = 0.0
    y_max = 1.0
    z = 0.0  # 2D geometry, so z-coordinate is constant

    # Define the mesh size
    mesh_size = 1

    # Numbering of mesh entities:
    #
    #  p4       p3       p6
    #  * --l3-- * --l7-- *
    #  |        |        |
    #  l4       l2       l6
    #  |        |        |
    #  * --l1-- * --l5-- *
    #  p1       p2       p5

    # Add points for the rectangle corners
    p1 = gmsh.model.geo.addPoint(x_min, y_min, z, mesh_size)
    p2 = gmsh.model.geo.addPoint(x_interface, y_min, z, mesh_size)
    p3 = gmsh.model.geo.addPoint(x_interface, y_max, z, mesh_size)
    p4 = gmsh.model.geo.addPoint(x_min, y_max, z, mesh_size)
    p5 = gmsh.model.geo.addPoint(x_max, y_min, z, mesh_size)
    p6 = gmsh.model.geo.addPoint(x_max, y_max, z, mesh_size)

    # Add lines connecting the points to form the rectangle
    l1 = gmsh.model.geo.addLine(p1, p2)
    l2 = gmsh.model.geo.addLine(p2, p3)
    l3 = gmsh.model.geo.addLine(p3, p4)
    l4 = gmsh.model.geo.addLine(p4, p1)
    l5 = gmsh.model.geo.addLine(p5, p2)
    l6 = gmsh.model.geo.addLine(p6, p5)
    l7 = gmsh.model.geo.addLine(p3, p6)

    # Add curve loop to define the boundary of the rectangle
    left_loop = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
    right_loop = gmsh.model.geo.addCurveLoop([-l5, -l6, -l7, -l2])

    # Add plane surface using the curve loop
    left_surface = gmsh.model.geo.addPlaneSurface([left_loop])
    right_surface = gmsh.model.geo.addPlaneSurface([right_loop])

    gmsh.model.geo.synchronize()

    no_tag = -1

    gmsh.model.addPhysicalGroup(dim2, [left_surface], no_tag, "Left")
    gmsh.model.addPhysicalGroup(dim2, [right_surface], no_tag, "Right")

    gmsh.model.mesh.generate(dim2)
    gmsh.write(filename)
    gmsh.finalize()
end

"""
Create the mesh above, write it to file and read it back in.
"""
function create_model()
    filename = "rectangle_mesh.msh"
    create_mesh(filename)
    model = GmshDiscreteModel(filename)
    return model
end

"""
Create the following model:
            * ------ * -------- *
            | "Left"  | "Right" |
            * ------- * ------- *
Set up a space over the left and right side, and interpolate
constant one into each space: u_L, u_R. Integrate u_L or u_R over
the interface between "Left" and "Right". This should give 1. Multiplying the
the integrand u_R with a scalar makes the integral 0. This looks like a bug.
"""
function run()
    model = create_model()
    Ω_L = Triangulation(model, tags="Left")
    Ω_R = Triangulation(model, tags="Right")

    element_order = 1
    element = ReferenceFE(lagrangian, Float64, element_order)
    V_L = TestFESpace(Ω_L,
        element,
        conformity=:H1)
    V_R = TestFESpace(Ω_R,
        element,
        conformity=:H1)

    U_L = TrialFESpace(V_L)
    U_R = TrialFESpace(V_R)

    one(x) = 1
    u_L = interpolate_everywhere(one, U_L)
    u_R = interpolate_everywhere(one, U_R)

    Γ = InterfaceTriangulation(Ω_L, Ω_R)

    degree = 2 * element_order
    dΓ = Measure(Γ, degree)

    println()
    println("I expect all these integrals to be equal to one. Why is the last one zero?")
    println("left side =", sum(∫(u_L.⁺)dΓ))
    println("right side =", sum(∫(u_R.⁻)dΓ))
    println("1*(left side) =", sum(∫(1.0 * u_L.⁺)dΓ))
    println("1*(right side) =", sum(∫(1.0 * u_R.⁻)dΓ))
end

run()

Output

I expect all these integrals to be equal to one. Why is the last one zero?
left side =0.9999999999999998
right side =0.9999999999999998
1*(left side) =0.9999999999999998
1*(right side) =0.0