Plot my Stress Tensor Components

Question

Hi,
I am working on learning how to use FEniCS and I produce the code bellow. However I don't understand how I am suppose to plot or write into a file my stress tensor (here sigma) what format it is? I tried plotting for instance sigma[0,0] for sigma_xx but it doesn't work. Any help would be welcome.
Regards
David

    from dolfin import *
# Create mesh
x0, y0, x1, y1, nx, ny = 0, 0, 10, 20, 20, 40
mesh = Rectangle(x0, y0, x1, y1, nx, ny, 'crossed')
# Create function space
V = VectorFunctionSpace(mesh, "Lagrange", 2)
# Create test and trial functions, and source term
u, w = TrialFunction(V), TestFunction(V)
b = Constant((0.0, 0.0))

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 - 2.0*nu))

# Stress
sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(w.cell().d)

# Governing balance equation
F = inner(sigma, grad(w))*dx - dot(b, w)*dx

# Extract bilinear and linear forms from F
a, L = lhs(F), rhs(F)

# Mark boundary subdomians
left, right = compile_subdomains(["(std::abs(x[0])       < DOLFIN_EPS) && on_boundary",
                                  "(std::abs(x[0] - 10) < DOLFIN_EPS) && on_boundary"])

bottom, top = compile_subdomains(["(std::abs(x[1])       < DOLFIN_EPS) && on_boundary",
                                  "(std::abs(x[1] - 20) < DOLFIN_EPS) && on_boundary"])

# Define boundary condiition                                  
zero = Constant(0.0)
zero2 = Constant((0.0, 0.0))
bc1 = DirichletBC(V.sub(0), zero, left)     # u.n = 0
bc2 = DirichletBC(V.sub(1), zero, bottom)   # u.n = 0
bc3 = DirichletBC(V.sub(0), zero, right)    # u.n = 0
c = Constant((0.0, -2.0))
bc4 = DirichletBC(V, c, top)
bc = [bc1, bc2, bc3, bc4]

# Set up PDE and solve
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bcs=bc)
solver = LinearVariationalSolver(problem)
solver.parameters["symmetric"] = True
solver.solve()


# Plot Results
plot(mesh, axes=True, title = "Meshed Geometry" )
plot(u, mode="displacement", axes=True, interactive = True,title = "Displacement Vectors" )

Comments

Please improve formatting of code. (Mark code and use click the curly braces).

---

Is this better?

Answer 1

Hi, adding these 3 lines to your code should make plotting components of sigma (you access them
as sigma[i, j]) possible :

W = TensorFunctionSpace(mesh, "Lagrange", 2)
sigma_w = project(2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(w.cell().d), W)
plot(sigma_w[0, 0], interactive=True)