Hai,

I am running simple elasticity problem. But facing the some solver issue leading to termination.

from dolfin import *

E=3300 #in mpa

nu=0.37

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

# Create mesh and define function space

mesh = Mesh(‘plate.xml’)

V = VectorFunctionSpace(mesh, “CG”, 1)

WV=TensorFunctionSpace(mesh, “CG”, 2)

# Define boundary condition

tol = 1E-14

bot = CompiledSubDomain(“near(x[1],-20) && on_boundary”) #bottom fixed

top=CompiledSubDomain(“near(x[1],20) && on_boundary”) #bottom fixed

bc1 = DirichletBC(V, Constant((0, 0, 0)), bot)

#bc2= DirichletBC(V.sub(0), Constant((0)), top)

bc=[bc1]

boundary = MeshFunction(“size_t”, mesh,mesh.topology().dim() - 1)

boundary.set_all(0)

top.mark(boundary, 1)

ds = Measure(“ds”)(subdomain_data=boundary)

Trc=Constant((0,10,0))

# Define strain and stress

def sigma(u):

return 2.0*mu*sym(grad(u)) + lmbda*tr(sym(grad(u)))*Identity(len(u))

# Define variational problem

u = TrialFunction(V)

v = TestFunction(V)

a = inner(grad(v),sigma(u))*dx

L = dot(v, Trc)*ds(1)

# Compute solution

u = Function(V)

solve(a==L,u,bc,solver_parameters={“linear_solver”: “mumps”})

m= File ("./ResultsDir-Von/fringeplot.pvd")

umap = File ("./ResultsDir-Von/u.pvd")

# Plot stress

sigmaT = project(sigma(u),WV)

sigma_max = 0.5*(sigmaT[0,0] + sigmaT[1,1]) + sqrt( (0.5*(sigmaT[0,0] - sigmaT[1,1]))**2 + (sigmaT[0,1])* 2 )*(sigmaT[0,0] - sigmaT[1,1]))**2 + (sigmaT[0,1])**2 )

sigma_min = 0.5(sigmaT[0,0] + sigmaT[1,1]) - sqrt( (0.5

fringes = project((sigma_max-sigma_min))

m<<fringes

umap<<u