diff --git a/cxxfem/tests/billyjoe.py b/cxxfem/tests/billyjoe.py
new file mode 100644
index 0000000000000000000000000000000000000000..d4b2043e266f0ae5d0fad9f54809b45227e4213d
--- /dev/null
+++ b/cxxfem/tests/billyjoe.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Billy-Joe's test named "Fluid layer on a solid layer under an external pressure" (variant 2)
+# we check that the horizontal stress is constant and equal to the applied pressure.
+
+import cxxfem as fem
+import gmsh, math
+
+
+def build_trapezoid():
+
+    x1 = 2.5    
+    x2 = 3.5
+    w0 = 6.
+    h0 = 4.
+
+    gmsh.model.add("trapezoid")
+    lc = h0 / 10
+    gmsh.model.geo.addPoint(x1, 0, 0, lc, 1)
+    gmsh.model.geo.addPoint(w0, 0, 0, lc, 2)
+    gmsh.model.geo.addPoint(w0, h0, 0, lc, 3)
+    gmsh.model.geo.addPoint(x2, h0, 0, lc, 4)
+
+    gmsh.model.geo.addLine(1, 2, 1)
+    gmsh.model.geo.addLine(2, 3, 2)
+    gmsh.model.geo.addLine(4, 3, 3)
+    gmsh.model.geo.addLine(1, 4, 4)
+
+    gmsh.model.geo.addCurveLoop([1, 2, -3, -4], 1)
+
+    gmsh.model.geo.addPlaneSurface([1], 1)
+
+    gmsh.model.geo.synchronize() # sinon warning "unknown entity..." lors du addPhysicalGroup
+
+    # physical groups - lines
+    bottom = gmsh.model.addPhysicalGroup(1, [1], 1000)
+    right = gmsh.model.addPhysicalGroup(1, [2])
+    top = gmsh.model.addPhysicalGroup(1, [3])
+    left = gmsh.model.addPhysicalGroup(1, [4])
+    gmsh.model.setPhysicalName(1, bottom, "bottom")
+    gmsh.model.setPhysicalName(1, right, "right")
+    gmsh.model.setPhysicalName(1, top, "top")
+    gmsh.model.setPhysicalName(1, left, "left")
+
+    # physical groups - surfaces
+    interior = gmsh.model.addPhysicalGroup(2, [1])
+    gmsh.model.setPhysicalName(2, interior, "interior")
+
+    # finalize and mesh the geometry
+    gmsh.model.mesh.generate(2)
+
+    # calculate the angle of the interface
+    angle = math.atan2(h0, x2-x1)
+    return angle
+
+
+if __name__ == "__main__":
+
+    pbl = fem.Problem()
+    angle = build_trapezoid()
+
+    p = 100.0
+
+    m = fem.Medium(pbl, "interior", 60e3, 0.0)
+
+    d1 = fem.Dirichlet(pbl, "right", "X", 0.0)
+    d2 = fem.Dirichlet(pbl, "bottom", "Y", 0.0)
+
+    n1 = fem.Neumann(pbl, "top", "Y", -p)
+    n2 = fem.Neumann(pbl, "left", "Y", -p*math.cos(angle))
+    n3 = fem.Neumann(pbl, "left", "X", p*math.sin(angle))
+
+    d4 = fem.Dirichlet(pbl, "interior", "Z", 0.0)
+
+    solver = fem.Solver(pbl)
+    solver.solve()
+
+    post = fem.Post(solver)
+    post.build_views()
+    post.show_views( [ "stress_xx", "force_vector" ] )
+    post.write()
+    
+    args = fem.parseargs()
+    if not args.nogui:    
+        post.view()
+
+    
\ No newline at end of file