No solution for simple cantilever beam with beam element

mohaiminul_islam_Mis · April 1, 2024, 4:51am

Hi,
I am trying to model a complex cantilever beam structure with beam elements (FEniCS 2019.1.0). I start with simple cases. I have 2 meshes (line1 and line2). the Line1 mesh works fine but Line2 doesn’t give any solution.

I use the same code . Here are my mesh files. I checked the marked domains it looks correct in both cases.

For line1.msh

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
11 10 0 0
1 0 0 0 0 
2 1 0 0 0 
3 2 0 0 0 
4 3 0 0 0 
5 4 0 0 0 
6 5 0 0 0 
7 6 0 0 0 
8 7 0 0 0 
9 8 0 0 0 
10 9 0 0 0 
11 10 0 0 0 
1 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 1 -2 
2 0.9999999000000001 -1e-07 -1e-07 2.0000001 1e-07 1e-07 0 2 2 -3 
3 1.9999999 -1e-07 -1e-07 3.0000001 1e-07 1e-07 0 2 3 -4 
4 2.9999999 -1e-07 -1e-07 4.0000001 1e-07 1e-07 0 2 4 -5 
5 3.9999999 -1e-07 -1e-07 5.0000001 1e-07 1e-07 0 2 5 -6 
6 4.9999999 -1e-07 -1e-07 6.0000001 1e-07 1e-07 0 2 6 -7 
7 5.9999999 -1e-07 -1e-07 7.0000001 1e-07 1e-07 0 2 7 -8 
8 6.9999999 -1e-07 -1e-07 8.000000099999999 1e-07 1e-07 0 2 8 -9 
9 7.9999999 -1e-07 -1e-07 9.000000099999999 1e-07 1e-07 0 2 9 -10 
10 8.999999900000001 -1e-07 -1e-07 10.0000001 1e-07 1e-07 0 2 10 -11 
$EndEntities
$Nodes
21 11 1 11
0 1 0 1
1
0 0 0
0 2 0 1
2
1 0 0
0 3 0 1
3
2 0 0
0 4 0 1
4
3 0 0
0 5 0 1
5
4 0 0
0 6 0 1
6
5 0 0
0 7 0 1
7
6 0 0
0 8 0 1
8
7 0 0
0 9 0 1
9
8 0 0
0 10 0 1
10
9 0 0
0 11 0 1
11
10 0 0
1 1 0 0
1 2 0 0
1 3 0 0
1 4 0 0
1 5 0 0
1 6 0 0
1 7 0 0
1 8 0 0
1 9 0 0
1 10 0 0
$EndNodes
$Elements
21 21 1 21
0 1 15 1
1 1 
0 2 15 1
2 2 
0 3 15 1
3 3 
0 4 15 1
4 4 
0 5 15 1
5 5 
0 6 15 1
6 6 
0 7 15 1
7 7 
0 8 15 1
8 8 
0 9 15 1
9 9 
0 10 15 1
10 10 
0 11 15 1
11 11 
1 1 1 1
12 1 2 
1 2 1 1
13 2 3 
1 3 1 1
14 3 4 
1 4 1 1
15 4 5 
1 5 1 1
16 5 6 
1 6 1 1
17 6 7 
1 7 1 1
18 7 8 
1 8 1 1
19 8 9 
1 9 1 1
20 9 10 
1 10 1 1
21 10 11 
$EndElements

and line2.msh

$MeshFormat
4.1 0 8
$EndMeshFormat
$Entities
21 20 0 0
1 0 0 0 0 
2 1 0 0 0 
3 2 0 0 0 
4 3 0 0 0 
5 4 0 0 0 
6 5 0 0 0 
7 6 0 0 0 
8 7 0 0 0 
9 8 0 0 0 
10 9 0 0 0 
13 9 1 0 0 
14 8 2 0 0 
15 7 3 0 0 
16 6 4 0 0 
17 5 5 0 0 
18 4 6 0 0 
19 3 7 0 0 
20 2 8 0 0 
21 1 9 0 0 
22 0 10 0 0 
23 10 0 0 0 
1 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 1 -2 
2 0.9999999000000001 -1e-07 -1e-07 2.0000001 1e-07 1e-07 0 2 2 -3 
3 1.9999999 -1e-07 -1e-07 3.0000001 1e-07 1e-07 0 2 3 -4 
4 2.9999999 -1e-07 -1e-07 4.0000001 1e-07 1e-07 0 2 4 -5 
5 3.9999999 -1e-07 -1e-07 5.0000001 1e-07 1e-07 0 2 5 -6 
6 4.9999999 -1e-07 -1e-07 6.0000001 1e-07 1e-07 0 2 6 -7 
7 5.9999999 -1e-07 -1e-07 7.0000001 1e-07 1e-07 0 2 7 -8 
8 6.9999999 -1e-07 -1e-07 8.000000099999999 1e-07 1e-07 0 2 8 -9 
9 7.9999999 -1e-07 -1e-07 9.000000099999999 1e-07 1e-07 0 2 9 -10 
10 8.999999900000001 -1e-07 -1e-07 10.0000001 1e-07 1e-07 0 2 10 -23 
11 8.999999900000001 -9.999999994736442e-08 -1e-07 10.0000001 1.0000001 1e-07 0 2 23 -13 
12 7.9999999 0.9999999000000001 -1e-07 9.000000099999999 2.0000001 1e-07 0 2 13 -14 
13 6.9999999 1.9999999 -1e-07 8.000000099999999 3.0000001 1e-07 0 2 14 -15 
14 5.9999999 2.9999999 -1e-07 7.0000001 4.0000001 1e-07 0 2 15 -16 
15 4.9999999 3.9999999 -1e-07 6.0000001 5.0000001 1e-07 0 2 16 -17 
16 3.9999999 4.9999999 -1e-07 5.0000001 6.0000001 1e-07 0 2 17 -18 
17 2.9999999 5.9999999 -1e-07 4.0000001 7.0000001 1e-07 0 2 18 -19 
18 1.9999999 6.9999999 -1e-07 3.0000001 8.000000099999999 1e-07 0 2 19 -20 
19 0.9999999000000001 7.9999999 -1e-07 2.0000001 9.000000099999999 1e-07 0 2 20 -21 
20 -9.999999994736442e-08 8.999999900000001 -1e-07 1.0000001 10.0000001 1e-07 0 2 21 -22 
$EndEntities
$Nodes
41 31 1 31
0 1 0 1
1
0 0 0
0 2 0 1
2
1 0 0
0 3 0 1
3
2 0 0
0 4 0 1
4
3 0 0
0 5 0 1
5
4 0 0
0 6 0 1
6
5 0 0
0 7 0 1
7
6 0 0
0 8 0 1
8
7 0 0
0 9 0 1
9
8 0 0
0 10 0 1
10
9 0 0
0 13 0 1
11
9 1 0
0 14 0 1
12
8 2 0
0 15 0 1
13
7 3 0
0 16 0 1
14
6 4 0
0 17 0 1
15
5 5 0
0 18 0 1
16
4 6 0
0 19 0 1
17
3 7 0
0 20 0 1
18
2 8 0
0 21 0 1
19
1 9 0
0 22 0 1
20
0 10 0
0 23 0 1
21
10 0 0
1 1 0 0
1 2 0 0
1 3 0 0
1 4 0 0
1 5 0 0
1 6 0 0
1 7 0 0
1 8 0 0
1 9 0 0
1 10 0 0
1 11 0 1
22
9.5 0.5 0
1 12 0 1
23
8.5 1.5 0
1 13 0 1
24
7.5 2.5 0
1 14 0 1
25
6.5 3.5 0
1 15 0 1
26
5.5 4.5 0
1 16 0 1
27
4.5 5.5 0
1 17 0 1
28
3.5 6.5 0
1 18 0 1
29
2.5 7.5 0
1 19 0 1
30
1.5 8.5 0
1 20 0 1
31
0.5 9.5 0
$EndNodes
$Elements
41 51 1 51
0 1 15 1
1 1 
0 2 15 1
2 2 
0 3 15 1
3 3 
0 4 15 1
4 4 
0 5 15 1
5 5 
0 6 15 1
6 6 
0 7 15 1
7 7 
0 8 15 1
8 8 
0 9 15 1
9 9 
0 10 15 1
10 10 
0 13 15 1
11 11 
0 14 15 1
12 12 
0 15 15 1
13 13 
0 16 15 1
14 14 
0 17 15 1
15 15 
0 18 15 1
16 16 
0 19 15 1
17 17 
0 20 15 1
18 18 
0 21 15 1
19 19 
0 22 15 1
20 20 
0 23 15 1
21 21 
1 1 1 1
22 1 2 
1 2 1 1
23 2 3 
1 3 1 1
24 3 4 
1 4 1 1
25 4 5 
1 5 1 1
26 5 6 
1 6 1 1
27 6 7 
1 7 1 1
28 7 8 
1 8 1 1
29 8 9 
1 9 1 1
30 9 10 
1 10 1 1
31 10 21 
1 11 1 2
32 21 22 
33 22 11 
1 12 1 2
34 11 23 
35 23 12 
1 13 1 2
36 12 24 
37 24 13 
1 14 1 2
38 13 25 
39 25 14 
1 15 1 2
40 14 26 
41 26 15 
1 16 1 2
42 15 27 
43 27 16 
1 17 1 2
44 16 28 
45 28 17 
1 18 1 2
46 17 29 
47 29 18 
1 19 1 2
48 18 30 
49 30 19 
1 20 1 2
50 19 31 
51 31 20 
$EndElements

and MWE

from dolfin import *
#from dolfin_adjoint import *
import meshio
from ufl import Jacobian, diag
parameters['allow_extrapolation'] = True
msh = meshio.read("gmsh/t1/line1.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points = msh.points, cells={"line": msh.cells_dict['line']}))
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
thick = Constant(1)
width = thick#/3
E = Constant(1)
nu = Constant(0.3)
G = E/2/(1+nu)
rho = Constant(1)
g = Constant(1)

S = thick*width
ES = E*S
EI1 = E*width*thick**3/12
EI2 = E*width**3*thick/12
GJ = G*0.26*thick*width**3
kappa = Constant(5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S
f = Constant(-rho*g*S)
def tangent(mesh):
    t = Jacobian(mesh)
    return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))
def tgrad(u):
    return dot(grad(u), t)
def generalized_strains(u):
    (w, theta) = split(u)
    return as_vector([dot(tgrad(w), t),
                      dot(tgrad(w), a1)-dot(theta, a2),
                      dot(tgrad(w), a2)+dot(theta, a1),
                      dot(tgrad(theta), t),
                      dot(tgrad(theta), a1),
                      dot(tgrad(theta), a2)])
def generalized_stresses(u):
    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])), generalized_strains(u))

t = tangent(mesh)



ez = as_vector([0, 0, 1])
a1 = cross(t, ez)
a1 /= sqrt(dot(a1, a1))
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))

def left_end(x, on_boundary):
    return near(x[0], 0)
def right_end(x, on_boundary):
    return near(x[0], 10,1e-5) and near(x[1],0,1e-5)


Ue = VectorElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, Ue*Ue)
M = VectorFunctionSpace(mesh, "CG", 1, dim=3)
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), left_end)
facets = MeshFunction("size_t", mesh, 0)
AutoSubDomain(left_end).mark(facets, 1)
AutoSubDomain(right_end).mark(facets, 2)
# Right = CompiledSubDomain("x[1]< 0.2-tol ", l=length,tol=1e-8) ## cantilever
# Right.mark(facets, 2)

ds = Measure("ds", domain=mesh, subdomain_data=facets, metadata = {"quadrature_scheme":"default","quadrature_degree": 1})


Sig = generalized_stresses(du)
Eps =  generalized_strains(u_)

dx_shear = dx(scheme="default",metadata={"quadrature_scheme":"default", "quadrature_degree": 1})
k_form = sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]]) + (Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear
#
l_form = inner(f, u_[2])*ds(2)


u = Function(W)
solve(k_form == l_form, u, bc)

print(0.5*assemble(inner(generalized_strains(u), generalized_stresses(u))*dx))

mohaiminul_islam_Mis · April 1, 2024, 10:07pm

One additional detail. The code works only for line1.msh if I apply load at the tip. If I apply the load at the middle say at (5,0,0). It gives no solution.

Any suggestions is much appreciated.

bleyerj · April 2, 2024, 9:24am

First check that the points corresponding to your boundaries or load application are correctly captured by the facets and boundary conditions.

mohaiminul_islam_Mis · April 2, 2024, 3:25pm

Thanks for the reply.
I checked the markers in Paraview. they look correct for both cases

dokken · April 2, 2024, 6:28pm

The issue is that you are using the ds measure, which is used for exterior facets, those facets only connected to a single cell.
Here is a minimal working code

from dolfin import *
import meshio
from ufl import Jacobian, diag

msh = meshio.read("line2.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points,
             cells={"line": msh.cells_dict['line']}))
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
thick = Constant(1)
width = thick  # /3
E = Constant(1)
nu = Constant(0.3)
G = E/2/(1+nu)
rho = Constant(1)
g = Constant(1)

S = thick*width
ES = E*S
EI1 = E*width*thick**3/12
EI2 = E*width**3*thick/12
GJ = G*0.26*thick*width**3
kappa = Constant(5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S
f = Constant(-rho*g*S)


def tangent(mesh):
    t = Jacobian(mesh)
    return as_vector([t[0, 0], t[1, 0], t[2, 0]])/sqrt(inner(t, t))


def tgrad(u):
    return dot(grad(u), t)


def generalized_strains(u):
    (w, theta) = split(u)
    return as_vector([dot(tgrad(w), t),
                      dot(tgrad(w), a1)-dot(theta, a2),
                      dot(tgrad(w), a2)+dot(theta, a1),
                      dot(tgrad(theta), t),
                      dot(tgrad(theta), a1),
                      dot(tgrad(theta), a2)])


def generalized_stresses(u):
    return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])), generalized_strains(u))


t = tangent(mesh)


ez = as_vector([0, 0, 1])
a1 = cross(t, ez)
a1 /= sqrt(dot(a1, a1))
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))


def left_end(x, on_boundary):
    return near(x[0], 0)


def right_end(x, on_boundary):
    return near(x[0], 10, 1e-5) and near(x[1], 0, 1e-5)


Ue = VectorElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, Ue*Ue)
M = VectorFunctionSpace(mesh, "CG", 1, dim=3)
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), left_end)
facets = MeshFunction("size_t", mesh, 0, 0)
AutoSubDomain(left_end).mark(facets, 1)
AutoSubDomain(right_end).mark(facets, 2)

metadata = {"quadrature_scheme": "default", "quadrature_degree": 2}
dS = Measure("dS", domain=mesh, subdomain_data=facets, metadata=metadata)


Sig = generalized_stresses(du)
Eps = generalized_strains(u_)

dx_shear = dx(scheme="default", metadata=metadata)
k_form = sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]]) + \
    (Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear
#
l_form = inner(f, 2*avg(u_[2]))*dS(2)

u = Function(W)
solve(k_form == l_form, u, bc)

print(0.5*assemble(inner(generalized_strains(u), generalized_stresses(u))*dx))


XDMFFile("u.xdmf").write(u)

mohaiminul_islam_Mis · April 2, 2024, 6:53pm

Thanks for the reply.
However, I can’t run the code.
It gives KeyError: 'FE11_C0_D1_Q1'

dokken · April 2, 2024, 7:01pm

What version of dolfin are your running and how did you install it?

mohaiminul_islam_Mis · April 2, 2024, 7:04pm

I use it in Lab computer. Can’t tell the exact details.
my conda list shows fenics-dolfin 2019.1.0 py38h4d3837b_26 conda-forge

mohaiminul_islam_Mis · April 3, 2024, 4:35am

Thanks, the code is working now.
Don’t know what happened. I just restarted jupyter and it started working

Topic		Replies	Views
Cantilever exact solution General	0	414	May 23, 2023
Static analysis of cantilever beam using beam element General	7	534	May 8, 2024
3D Elastic Beam: Kinematic & Dynamics mathematics	0	513	January 2, 2021
1D Cantilever Beam problem	2	999	March 22, 2020
Converting legacy FEniCS code to FEniCSx causes issue General	19	377	May 29, 2024

No solution for simple cantilever beam with beam element

Related topics