# Cantilever example with beam element for "CG" order 2 space fails

Hi
I am solving simple linear elasticity problem like here . I am using 1D beam element in 3D space.
One question i have is the solution gives nan values for anything other than “CG” order 1 space.
I tried “CG” order 2 , “DG” order 1.
My MWE

from dolfin import *
from ufl import Jacobian, diag
#from plotting import plot

length = 10.0
N = 80 # number of elements

points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)
cells = np.array([[i, i+1] for i in range(N)])
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:

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

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))

# We now define the material and geometrical constants which will be used in the constitutive relation. We consider the case of a rectangular cross-section of width $b$ and height $h$ in directions $a_1$ and $a_2$. The bending inertia will therefore be $I_1 = bh^3/12$ and $I_2=hb^3/12$. The torsional inertia is $J=\beta hb^3$ with $\beta\approx 0.26$ for $h=3b$. Finally, the shear areas are approximated by $S_1=S_2=\kappa S$ with $\kappa=5/6$.

# In[3]:

thick = Constant(1)
width = thick#/3
E = Constant(1e5)
nu = Constant(0.3)
G = E/2/(1+nu)
rho = Constant(1e-3)
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

# We now consider a mixed $\mathbb{P}_1/\mathbb{P}_1$-Lagrange interpolation for the displacement and rotation fields. The variational form is built using a function generalized_strains giving the vector of six generalized strains as well as a function generalized_stresses which computes the dot product of the strains with the above-mentioned constitutive matrix (diagonal here). Note that since the 1D beams are embedded in an ambient 3D space, the gradient operator has shape (3,), we therefore define a tangential gradient operator tgrad by taking the dot product with the local tangent vector $t$.
#
# Finally, similarly to Reissner-Mindlin plates, shear-locking issues might arise in the thin beam limit. To avoid this, reduced integration is performed on the shear part $Q_1\widehat{\gamma}_1+Q_2\widehat{\gamma}_2$ of the variational form using a one-point rule.

# In[6]:

Ue = VectorElement("CG", mesh.ufl_cell(), 2, dim=3)
Te = VectorElement("CG", mesh.ufl_cell(), 1, dim=3)
W = FunctionSpace(mesh, MixedElement([Ue, Te]))

u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

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

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

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 = Constant(-rho*S*g)*w_[1]*dx

def bottom(x, on_boundary):
return near(x[0], 0.)
bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), bottom)

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

print("disp",max(abs(u.split()[0].vector()[:])))


Hi "DG",1 surely can’t work (without a proper DG formulation) as you need continuous interpolation of the displacement and rotation field. If you want higher degree of interpolation, you also need to increase the degree of the integration measure in dx_shear otherwise, reduced integration is “too reduced”. See for instance Reissner-Mindlin plate with Quadrilaterals — Numerical tours of continuum mechanics using FEniCS master documentation

1 Like

Thanks for your reply.  dx_shear, deg 4 works

Good evening @mohaiminul_islam_Mis,

I am just behind you on the cantilever 1D beam element in 3D space case.

Have you been able to work it through?

If so, could you share your final piece of code?

Regards,

The following is working

from dolfin import *
from ufl import Jacobian, diag
#from plotting import plot

length = 10.0
N = 80 # number of elements

points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)
cells = np.array([[i, i+1] for i in range(N)])
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:

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

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))

# We now define the material and geometrical constants which will be used in the constitutive relation. We consider the case of a rectangular cross-section of width $b$ and height $h$ in directions $a_1$ and $a_2$. The bending inertia will therefore be $I_1 = bh^3/12$ and $I_2=hb^3/12$. The torsional inertia is $J=\beta hb^3$ with $\beta\approx 0.26$ for $h=3b$. Finally, the shear areas are approximated by $S_1=S_2=\kappa S$ with $\kappa=5/6$.

# In[3]:

thick = Constant(1)
width = thick#/3
E = Constant(1e5)
nu = Constant(0.3)
G = E/2/(1+nu)
rho = Constant(1e-3)
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

# We now consider a mixed $\mathbb{P}_1/\mathbb{P}_1$-Lagrange interpolation for the displacement and rotation fields. The variational form is built using a function generalized_strains giving the vector of six generalized strains as well as a function generalized_stresses which computes the dot product of the strains with the above-mentioned constitutive matrix (diagonal here). Note that since the 1D beams are embedded in an ambient 3D space, the gradient operator has shape (3,), we therefore define a tangential gradient operator tgrad by taking the dot product with the local tangent vector $t$.
#
# Finally, similarly to Reissner-Mindlin plates, shear-locking issues might arise in the thin beam limit. To avoid this, reduced integration is performed on the shear part $Q_1\widehat{\gamma}_1+Q_2\widehat{\gamma}_2$ of the variational form using a one-point rule.

# In[6]:

Ue = VectorElement("CG", mesh.ufl_cell(), 2, dim=3)
Te = VectorElement("CG", mesh.ufl_cell(), 1, dim=3)
W = FunctionSpace(mesh, MixedElement([Ue, Te]))

u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

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

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

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 = Constant(-rho*S*g)*w_[1]*dx

def bottom(x, on_boundary):
return near(x[0], 0.)
bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), bottom)

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

print("disp",max(abs(u.split()[0].vector()[:])))

1 Like

Nice !

By any chance were you able to do concentrated loads too ?

I am trying to do something from this : How can point loads be applied to the end of the cantilever beam? - #7 by conpierce8

but for now it is just a mess …