Stress contour in frequency domain

1

Hi all,
I am doing dynamic analysis in frequency domain. I have computed the displacement response in frequency domain. I have to extend it to calculate Von mises stress.
The MWE is given below-

from dolfin import *
import numpy
from ufl import indices
import matplotlib.pyplot as plt
frequency=50000


L = 0.2
H = 0.01
Nx = 200
Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)


Element1 = VectorElement("CG", mesh.ufl_cell(), 1)
Element2 = VectorElement("CG", mesh.ufl_cell(), 1)

V_elem = MixedElement([Element1, Element2])

Vcomplex = FunctionSpace(mesh, V_elem)


E, rho, nu = 300e8, 8000, 0.3

alpha = Constant(0)
beta = Constant(1e-6)

mu = E/2./(1+nu)
lmbda = E * nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))

def sigma(v):
    dim = v.geometric_dimension()
    return mu*eps(v)*2 + lmbda*tr(eps(v))*Identity(dim)

zero=Constant((0.0, 0.0))
one=Constant((0.0, -1.0))

def left(x, on_boundary):
    return near(x[0], 0.)

def rigth(x, on_boundary):
    return near(x[0], 0.2)

boundary = [DirichletBC(Vcomplex.sub(0), zero, left), DirichletBC(Vcomplex.sub(1), zero, left)]
source = [DirichletBC(Vcomplex.sub(0), one, rigth), DirichletBC(Vcomplex.sub(1), one, rigth)]
bc=boundary+source


uR, uI = TrialFunctions(Vcomplex)
wR, wI = TestFunctions(Vcomplex)
Nsteps  = 100
import numpy as np
frequency = np.linspace(2000, 4000, Nsteps+1)
u_tip = np.zeros((Nsteps+1,))
for (i, dw) in enumerate(np.diff(frequency)):
    w = frequency[i+1]
    M = (frequency)** 2 * rho * (dot(uR,wR) + dot(uI,wI)) *dx
    K = (inner(sigma(uR), grad(wR))+inner(sigma(uI), grad(wI)))*dx
    C = alpha * M + beta * K

    F = M - K + C
    a, L = lhs(F), rhs(F)

    solnU = Function(Vcomplex)
    problem = LinearVariationalProblem(a, L, solnU, bcs=bc)
    solver = LinearVariationalSolver(problem)
    solver.solve()
    solnUR, solnUI=solnU.split()

    file_results = XDMFFile("Amplitude_response.xdmf")
    file_results.parameters["flush_output"] = True
    file_results.parameters["functions_share_mesh"]= True
    file_results.write(solnUR, w)
    file_results.write(solnUI, w)

For static case, demo for von mises stress computation can be found out here.