From 678162084b61fd1a9f44c9eeeff00dbebe70d521 Mon Sep 17 00:00:00 2001
From: Romain Boman <romain.boman@gmail.com>
Date: Wed, 22 Nov 2023 16:54:21 +0100
Subject: [PATCH] add billy-joe's test in 3d
---
cxxfem/tests/billyjoe.py | 2 +
cxxfem/tests/billyjoe3d.py | 107 +++++++++++++++++++++++++++++++++++++
2 files changed, 109 insertions(+)
create mode 100644 cxxfem/tests/billyjoe3d.py
diff --git a/cxxfem/tests/billyjoe.py b/cxxfem/tests/billyjoe.py
index d4b2043..5f3f2c6 100644
--- a/cxxfem/tests/billyjoe.py
+++ b/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
diff --git a/cxxfem/tests/billyjoe3d.py b/cxxfem/tests/billyjoe3d.py
new file mode 100644
index 0000000..b8e8439
--- /dev/null
+++ b/cxxfem/tests/billyjoe3d.py
@@ -0,0 +1,107 @@
+#!/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
--
GitLab