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)