add billy-joe's test in 3d

cxxfem/tests/billyjoe.py

@@ -3,6 +3,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.
+# This is done here in 2d - plane stress.
+
 
 import cxxfem as fem
 import gmsh, math

cxxfem/tests/billyjoe3d.py

+#!/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)
+
+    outDimTags = gmsh.model.geo.extrude([(2, 1)], 0, 0, 1)
+    t_base2, t_volume, t_bottom, t_right, t_top, t_left = outDimTags
+
+    print(outDimTags)
+
+    gmsh.model.geo.synchronize() # sinon warning "unknown entity..." lors du addPhysicalGroup
+
+    #gmsh.fltk.run()
+
+    # physical groups - lines
+    bottom = gmsh.model.addPhysicalGroup(2, [ t_bottom[1] ], 1000)
+    right = gmsh.model.addPhysicalGroup(2, [ t_right[1] ])
+    top = gmsh.model.addPhysicalGroup(2, [ t_top[1] ])
+    left = gmsh.model.addPhysicalGroup(2, [ t_left[1] ])
+    gmsh.model.setPhysicalName(2, bottom, "bottom")
+    gmsh.model.setPhysicalName(2, right, "right")
+    gmsh.model.setPhysicalName(2, top, "top")
+    gmsh.model.setPhysicalName(2, left, "left")
+
+
+    # physical groups - surfaces
+    base1 = gmsh.model.addPhysicalGroup(2, [1])
+    gmsh.model.setPhysicalName(2, base1, "base1")
+    base2 = gmsh.model.addPhysicalGroup(2, [ t_base2[1] ])
+    gmsh.model.setPhysicalName(2, base2, "base2")
+
+    # physical groups - volumes
+    interior = gmsh.model.addPhysicalGroup(3, [1])
+    gmsh.model.setPhysicalName(3, interior, "interior")
+
+    # finalize and mesh the geometry
+    gmsh.model.mesh.generate(3)
+
+
+    # gmsh.fltk.run()
+
+
+    # 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.3)
+
+    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))
+
+    d3 = fem.Dirichlet(pbl, "base1", "Z", 0.0)
+    d4 = fem.Dirichlet(pbl, "base2", "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