# Linear elasticity with beam element

I am trying to solve simple linear elasticity problem like here. I am using 1D beam element in 3D space.

from dolfin import *
import numpy as np
import meshio
length = 10.0
N = 40 # 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 left_end(x, on_boundary):
return near(x[0], 0,0.1)
def right_end(x, on_boundary):
return near(x[0], 10, 0.1)

facets = MeshFunction("size_t", mesh, 1)
AutoSubDomain(left_end).mark(facets, 1)
AutoSubDomain(right_end).mark(facets, 2)
ds = Measure("ds", domain=mesh, subdomain_data=facets)
def eps(v):
E = Constant(1e5)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def sigma(v):
return lmbda*tr(eps(v))*Identity(3) + 2.0*mu*eps(v)

rho_g = 1
f = Constant((0,0,-rho_g))
T = Constant((0,0,0))
V = VectorFunctionSpace(mesh, 'Lagrange', degree=3, dim=3)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx +inner(T, u_)*dx

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

bc = DirichletBC(V, Constant((0.,0.,0)), left)

u = Function(V, name="Displacement")
solve(a == l, u, bc)



The deflection shape seems wrong

I also tried with mixed element with 6 DOFs following this example.

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

length = 10.0
N = 40 # 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(), 1, dim=3)
W = FunctionSpace(mesh, Ue*Ue)

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_[2]*dx

# Clamped boundary conditions are considered at the bottom $z=0$ level and the linear problem is finally solved.

# In[7]:

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)

ffile = XDMFFile("shell-results.xdmf")
ffile.parameters["functions_share_mesh"] = True
v = u.sub(0, True)
v.rename("Displacement", "")

ffile.close()


The curvature still looks wrong looks too straight.

I suggest you first compare to an analytic solution with a sufficiently fine mesh. Please note that FEniCS does not handle “classical” Hermite cubic beam elements.