Why i plot displacement arrow graph rather than cloud graph?

1

hi, everyone.
I want to get the displacement solution and plotting it like Abaqus show. but i got arrow graph when plotting.

image
image832×817 107 KB

image
image1320×737 39.1 KB

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()
2

Hi,
You need to plot the norm since u is a vector field. \| u\| = \sqrt{\sum_{i}u_i^2}

Use

norm_u = sqrt(inner(u,u))
plot(norm_u)

That should give you the result to compare with Abaqus.

3

thanks, it works but there is no scale bar (showing which colour represent which range data). how ti fix it?

4 
surf = plot(norm_u)
plt.colorbar(surf)

should add the scalar bar.

5

In general, you can save your solution as a pvd or xdmf-file for post processing with for instance Paraview (this is especially recommended in cases of 3D geometries).

6

yes, paraview is referable for 3D cases, here, I just want to look around what happen when running code.

7

hi, I want to change the colour of map like Abaqus. for doing this, i write:

norm_u = sqrt(inner(u,u))
plot(norm_u)
surf = plot(norm_u)
plt.colorbar(surf,cmap='jet')
plt.show()

but i can not get the correct colour.

image
image966×815 62.6 KB

8

You need to pass the cmap into the plotting instruction, as in this MWE:

from dolfin import *

mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, "CG", 1)
v = Function(V)
v.vector()[0] = 1

import matplotlib.pyplot as plt

surf = plot(v, cmap=plt.cm.jet)
plt.colorbar(surf)
plt.savefig("u.png")

This is also easy to achieve in Paraview, which offers a large variety of colormaps.

9

thanks dokken. it works