hi, everyone.
I want to get the displacement solution and plotting it like Abaqus show. but i got arrow graph when plotting.
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
set_log_active(False)
L = 1.; H = 1.
Nx = 10; Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)
nodal_coordinates = mesh.coordinates()
Vm = FunctionSpace(mesh,'DG',1)
W = VectorFunctionSpace(mesh,'P',1)
F = TensorFunctionSpace(mesh, "Lagrange", 2)
u , v = TrialFunction(W), TestFunction(W)
E = 10e3
nu = 0.3
# plane problem type
model = "plane_stress"
# lame's parameters for plane strain
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
# lame's parameters for plane stress
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)
# displacement boundary condition
bottom = CompiledSubDomain("near(x[1], 0) && on_boundary")
left = CompiledSubDomain("near(x[0], 0) && on_boundary")
bc_u = [DirichletBC(W.sub(1), Constant(0.0), bottom),
DirichletBC(W.sub(0), Constant(0.0), left),]
# traction boundary condition
top = CompiledSubDomain("near(x[1], 1) && on_boundary")
load = Constant((0, 1e-3))
# index boundary no.
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
def epsilon(u):
return 0.5*(nabla_grad(u)+nabla_grad(u).T)
def sigma(u):
return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))
LHS = inner(sigma(u), epsilon(v))*dx
RHS = dot(load, v)*ds(1)
u = Function(W, name="Displacement")
K = assemble(LHS)
f = assemble(RHS)
[bc.apply(K, f) for bc in bc_u]
solve(K,u.vector(),f)
plot(u)
plt.show()