Hi, I have been trying to bend a simple cantilever beam under self load and an applied distributed load on the top face with a very simple gmsh and fenics script. For some reason it just doesn't work and the output I get is displayed below. I would truly appreciate if someone could help me out here.
The complete .geo file is noted below:
Point(1) = {0, 0, 0, 0.125};
Point(2) = {8, 0, 0, 0.125};
Point(3) = {8, 0, 1, 0.125};
Point(4) = {0, 0, 1, 0.125};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(6) = {5};
Extrude {0, 1, 0} {
Surface{6};
}
Physical Surface(1) = {28}; // top face (force to bend in -ve 'y' direction)
Physical Surface(2) = {15};
Physical Surface(3) = {23};
Physical Surface(4) = {6}; // bottom face
Physical Surface(5) = {27}; // fixed left side
Physical Surface(6) = {19}; // right end
Physical Volume(1) = {1};
The complete fenics script along with .msh, mesh.xml, facet.xml, physical.xml are posted at THIS LINK HERE
This is how I am implying the elasticity equation:
mesh = Mesh("beam.xml")
boundaries = MeshFunction("size_t", mesh, "beam_facet_region.xml")
f= Constant((0.0, 10.0, 0.0))
ds = Measure("ds")[boundaries]
E = 100
nu = 0.3
mu, lmbda = E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 - 2.0*nu))
# Stress
sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(v.cell().d)
a = inner(sigma, grad(v))*dx
L = dot(f, v)*dx + dot(f, v)*ds(1) # self load + load on top face
clamp = Constant((0.0, 0.0, 0.0))
bc = [DirichletBC(V, clamp, boundaries, 5)]
This is what I get on running the script and I would highly appreciate some help/ guidance.
EDIT: I just ran a similar python script based on a beam mesh created in Fenics. It gives the same result when plotted with "mode=displacement" or even in paraview warp by vector. However, when I use "glyph" filter for fenics or gmsh mesh, it shows the foll. result:
So, I am not sure but I suppose there may be a problem with the code rather than the mesh.
