Hi, all. As I’ve mentioned elsewhere, I have been (slowly) working on implementing this linear shell model from @bleyerj in Dolfinx. I’m working via docker with Dolfinx version 0.3.1.0, git commit fdb1d1871945e308db07f94175fe3dd08033be5f. I have a code that runs, but the results are NaN at present (I will share the mesh geometry if there is a way to upload these non-image files that that I am unaware of…):
from dolfinx.fem import Constant, Function, FunctionSpace, locate_dofs_geometrical, dirichletbc
from dolfinx.fem.petsc import LinearProblem
from dolfinx import io
import ufl
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
# read mesh from file
filename = "I_beam.xdmf"
infile = io.XDMFFile(MPI.COMM_WORLD, filename, "r")
mesh = infile.read_mesh(name="Grid")
# define material parameters
thick = Constant(mesh, 1e-3)
E = Constant(mesh, 210e9)
nu = Constant(mesh, 0.3)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)
lmbda_ps = 2 * lmbda * mu / (lmbda + 2 * mu)
# loading (self-weight)
f = Constant(mesh, (0, 0.2, -1))
def local_frame(mesh):
t = ufl.Jacobian(mesh)
if mesh.topology.dim == 2:
t1 = ufl.as_vector([t[0, 0], t[1, 0], 0])
t2 = ufl.as_vector([t[0, 1], t[1, 1], 0])
else:
t1 = ufl.as_vector([t[0, 0], t[1, 0], t[2, 0]])
t2 = ufl.as_vector([t[0, 1], t[1, 1], t[2, 1]])
e3 = ufl.cross(t1, t2)
e3 /= ufl.sqrt(ufl.dot(e3, e3))
ey = ufl.as_vector([0, 1, 0])
ez = ufl.as_vector([0, 0, 1])
e1 = ufl.cross(ey, e3)
norm_e1 = ufl.sqrt(ufl.dot(e1, e1))
e1 = ufl.conditional(ufl.lt(norm_e1, 0.5), ez, e1 / norm_e1)
e2 = ufl.cross(e3, e1)
e2 /= ufl.sqrt(ufl.dot(e2, e2))
return e1, e2, e3
frame = local_frame(mesh)
e1, e2, e3 = frame
# create test and trial function spaces
Ue = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2, dim=3)
Te = ufl.VectorElement("CR", mesh.ufl_cell(), 1, dim=3)
V = FunctionSpace(mesh, ufl.MixedElement([Ue, Te]))
v = Function(V)
u, theta = ufl.split(v)
v_ = ufl.TestFunction(V)
u_, theta_ = ufl.split(v_)
dv = ufl.TrialFunction(V)
# boundary conditions
def clamped_boundary(x):
xMin1 = -0.1250
xMin2 = 0.1250
return np.logical_or(np.isclose(x[0], xMin1), np.isclose(x[0], xMin2))
u_D = Function(V)
with u_D.vector.localForm() as loc:
loc.set(0.0)
# locate dofs
Vs = [V.sub(i).collapse()[0] for i in range(V.num_sub_spaces)]
dofs = [locate_dofs_geometrical((V.sub(i), Vs[i]), clamped_boundary) for i in range(V.num_sub_spaces)]
bc = [dirichletbc(u_D, dofs[i]) for i in range(V.num_sub_spaces)]
# define projectors
def vstack(vectors):
"""Stack a list of vectors vertically."""
return ufl.as_matrix([[v[i] for i in range(len(v))] for v in vectors])
def hstack(vectors):
"""Stack a list of vectors horizontally."""
return vstack(vectors).T
# In-plane projection
P_plane = hstack([e1, e2])
# define tangential gradient operator
def t_grad(u):
"""Tangential gradient operator"""
g = ufl.grad(u)
return ufl.dot(g, P_plane)
# extract gradient of in-plane displacement, define membrane strain, bending curvature, shear strain
t_gu = ufl.dot(P_plane.T, t_grad(u))
eps = ufl.sym(t_gu)
kappa = ufl.sym(ufl.dot(P_plane.T, t_grad(ufl.cross(e3, theta))))
gamma = t_grad(ufl.dot(u, e3)) - ufl.dot(P_plane.T, ufl.cross(e3, theta))
eps_ = ufl.replace(eps, {v: v_})
kappa_ = ufl.replace(kappa, {v: v_})
gamma_ = ufl.replace(gamma, {v: v_})
# stress measures
def plane_stress_elasticity(e):
return lmbda_ps * ufl.tr(e) * ufl.Identity(2) + 2 * mu * e
N = thick * plane_stress_elasticity(eps)
M = thick ** 3 / 12 * plane_stress_elasticity(kappa)
Q = mu * thick * gamma
# drilling rotation stabilization
drilling_strain = (t_gu[0, 1] - t_gu[1, 0]) / 2 - ufl.dot(theta, e3)
drilling_strain_ = ufl.replace(drilling_strain, {v: v_})
drilling_stress = E * thick ** 3 * drilling_strain
# build equations
Wdef = (
ufl.inner(N, eps_)
+ ufl.inner(M, kappa_)
+ ufl.dot(Q, gamma_)
+ drilling_stress * drilling_strain_
) * ufl.dx
a = ufl.derivative(Wdef, v, dv)
Wext = ufl.dot(f, u_) * ufl.dx
# solve
problem = LinearProblem(a, Wext, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
u = problem.solve()
If, for example, I do
print(u.sub(0).x.array)
the output is [nan nan nan ... nan nan nan].
So there are problems. The are numerous changes from the original implementation, some of which may well be ill-advised(?). I will of course accept any help whatsoever, but what would be particularly useful to assist in my debugging would be an explanation of how I can map the DOF result vector the nodes of the input mesh? Specifically, the sizes of the mesh nodal geometry and one of the DOF vectors are:
print(mesh.geometry.x.shape)
(306, 3)
print(u.sub(0).x.array.shape)
(5748,)
How can I pull out, for example, only the X-direction displacements (in which case I would expect a size of (306,))? I see various instances of dofmap properties but it isn’t clear to me out these work.
Any help is appreciated, thanks in advance!
, but this approach will definitely be helpful as I continue to debug. Thanks!