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:
infile.read(mesh)
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 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))
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 = 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()[:])))
