Plane Strain Linear Elasticity Program giving Error with Kernel Dead

dolfinx
1

I have a 2D geometry created in gmsh api and I import to fenicsx and try to perform a plane strain linear elastic simulation. I get this error of getting my kernel dead. Where I am going wrong?

gmsh.initialize()

gmsh.model.add("Winding Pack-1")

gdim = 3
mesh_comm = MPI.COMM_WORLD
model_rank = 0
h = 1

lc = 1e-3

# Added Points in the geometry

gmsh.model.geo.addPoint(0,0,0,lc,1)
gmsh.model.geo.addPoint(1,0,0,lc,2)
gmsh.model.geo.addPoint(1.5,1,0,lc,3)
gmsh.model.geo.addPoint(-0.5,1,0,lc,4)
#gmsh.model.geo.addPoint(0.3,0.3,0,lc,5)
#gmsh.model.geo.addPoint(0.7,0.3,0,lc,6)
#gmsh.model.geo.addPoint(0.7,0.6,0,lc,7)
#gmsh.model.geo.addPoint(0.3,0.6,0,lc,8)

# Added Lines in the geometry

gmsh.model.geo.addLine(1,2,1)
gmsh.model.geo.addLine(2,3,2)
gmsh.model.geo.addLine(3,4,3)
gmsh.model.geo.addLine(4,1,4)
#gmsh.model.geo.addLine(5,6,5)
#gmsh.model.geo.addLine(6,7,6)
#gmsh.model.geo.addLine(7,8,7)
#gmsh.model.geo.addLine(8,5,8)

gmsh.model.geo.addCurveLoop([1,2,3,4],1)

gmsh.model.geo.addPlaneSurface([1],1)

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(1, [1], name="bottom")
gmsh.model.addPhysicalGroup(1, [2], name="right")
gmsh.model.addPhysicalGroup(1, [3], name="top")
gmsh.model.addPhysicalGroup(1, [4], name="left")

gmsh.model.addPhysicalGroup(2, [1], name="my_surface")

gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)

gmsh.model.mesh.generate(gdim)

domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
#facet_markers.name = "Facet markers"

gmsh.write("Winding Pack-1.msh")
gmsh.fltk.run()
gmsh.finalize()

import ufl

def strain(u, repr ="vectorial"):
    eps_t = ufl.sym(ufl.grad(u))
    if repr =="vectorial":
        return ufl.as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
    elif repr =="tensorial":
        return eps_t

E = fem.Constant(domain, 210e9)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = ufl.as_matrix([[lmbda+2*nu, lmbda, 0.0],[lmbda, lmbda+2*nu,  0.0],[0.0, 0.0, nu]])

def stress(u, repr ="vectorial"):
    sigv = ufl.dot(C, strain(u))
    if repr =="vectorial":
        return sigv
    elif repr =="tensorial":
        return ufl.as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))

#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

# Applying a uniform pressure
rho = 7850
g = 9.81
fp = fem.Constant(domain, 10e4)
#fg = fem.Constant(domain, np.array([0, -rho*g]))

#Self-weight on the surface
n = FacetNormal(domain)
#L_form = dot(-fp*n,u_) * ds(5) + dot(fg*n,u_) * ds(5)

L_form = dot(-fp*n,u_) * ds(5)

#Boundary Conditions (simply supported on top and bottom edges)

V_uz, _ = V.sub(2).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(2), V_uz), gdim-1, facet_markers.find(1))
top_dofs = fem.locate_dofs_topological((V.sub(2), V_uz), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(2), V_uz), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(2), V_uz), gdim-1, facet_markers.find(2))


u0z = fem.Function(V_uz)
bcs = [fem.dirichletbc(u0z, bottom_dofs, V.sub(2)),fem.dirichletbc(u0z, top_dofs, V.sub(2)),fem.dirichletbc(u0z, left_dofs, V.sub(2)),fem.dirichletbc(u0z, right_dofs, V.sub(2))]


WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
WPsolve.solve()

V0 = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(stress(u), V0.element.interpolation_points())
stress_field = fem.Function(V0, name="Stress")
stress_field.interpolate(stress_expr)

with io.XDMFFile(domain.comm, "WPsolve_elasticity.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(u_sol)
    file.write_function(stress_field)
2

The code is not reproducible, since there are several missing imports.

Consider copying your code into a standalone python file, and then running it in a terminal rather than a jupyter notebook. That may give additional errors with more details that may be more helpful than “The kernel is dead” to debug your issue.

3

Hallo Francesco,

I am now getting the error “Segmentation fault (core dumped)”. What is the next step to this?

4

That may mean a lot of things, but there is no way to help you any further without an example we can run. It would be helpful to simplify the example, for instance (for starters) to see if that same error is obtained with a builtin mesh, say the unit square.

5

Thank you Francesco, I could overcome this problem as my boundary conditions provided were incorrect.