Finding stress values at each vertex

mesh
1

I’m trying to find the stress values at each vertex so that I can plot a graph having a stress distribution over the entire mesh. But it seems there are 4 stress values at each vertex and I’m not sure if I need to do further calculations or not.

I’m trying to plot the graph with plotly. Here is the code:

# --------------------
# Functions and classes
# --------------------
def bottom(x, on_boundary):
    return (on_boundary and near(x[1], 0.0))

# Strain function
def epsilon(u):
    return sym(grad(u))

# Stress function
def sigma(u):
    return lambda_*div(u)*Identity(2) + 2*mu*epsilon(u)

# --------------------
# Parameters
# --------------------

# Density (kg/m3)
rho = Constant(1240.0)

# Young's modulus (Pa) and Poisson's ratio
E = 33e6 
nu = 0.39

# Lame's constants
lambda_ = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

l_x, l_y = 5.0, 5.0 # Domain dimensions
n_x, n_y = 20, 20 # Number of elements

# Load
g_z = -10000
b_z = -10
g = Constant((0.0, g_z))
b = Constant((0.0, b_z))

# Model type
model = "plane_strain"

if model == "plane_stress":
  lambda_ = 2*mu*lambda_/(lambda_+2*mu)

# --------------------
# Geometry
# --------------------
mesh = RectangleMesh(Point(0.0, 0.0), Point(l_x, l_y), n_x, n_y)

# --------------------
# Function spaces
# --------------------

V = VectorFunctionSpace(mesh, "CG", 1) # displacement function space
u_tr = TrialFunction(V)
u_test = TestFunction(V)

# --------------------
# Boundary conditions
# --------------------

bc = DirichletBC(V.sub(1), 0.0, bottom)

# Neumann boundary conditions (top)
top = AutoSubDomain(lambda x: near(x[1], l_y))

# Definition of Neumann boundary condition domain
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries, 1)

ds = ds(subdomain_data=boundaries)

# --------------------
# Weak form
# --------------------
a = inner(sigma(u_tr), epsilon(u_test))*dx
l = rho*dot(b, u_test)*dx + inner(g, u_test)*ds(1)

# --------------------
# Solver
# --------------------

u = Function(V)
A, L = assemble_system(a, l, bc)

# --------------------
# Problem Definition & Solution
# --------------------

problem = LinearVariationalProblem(a,l,u,bc)
solver = LinearVariationalSolver(problem)
solve(A, u.vector(), L)

# --------------------
# Calculate stress
# --------------------

stress = sigma(u)

# --------------------
# Project Stress
# --------------------
def local_project(fce, space):
    lp_trial, lp_test = TrialFunction(space), TestFunction(space)
    lp_a = inner(lp_trial, lp_test)*dx
    lp_L = inner(fce, lp_test)*dx
    local_solver = LocalSolver(lp_a, lp_L)
    local_solver.factorize() # caches the lhs
    lp_f = Function(space)
    local_solver.solve_local_rhs(lp_f)
    return lp_f

V0 = TensorFunctionSpace(mesh, "DG", 0)
stress_2 = local_project(stress, V0)

#Plot deformed body
import plotly.graph_objects as go
def mesh2triang(mesh):
  import matplotlib.tri as tri
  xy = mesh.coordinates() #mesh.points
  return tri.Triangulation(xy[:, 0], xy[:, 1], mesh.cells())

gdim = mesh.topology().dim() 
nv = mesh.num_vertices() 
w0 = u.compute_vertex_values(mesh) 
U = [w0[i*nv: (i+1) *nv] for i in range(gdim)]
X = mesh.coordinates()
X = [X[:, i] for i in range(gdim)]

xdef = [X[i] + U[i] for i in range(gdim)]

stress_values = stress_2.compute_vertex_values(mesh)

xy = mesh.coordinates()
x, y, z = xdef[0], xdef[1], np.zeros(xy.shape[0])
ijk = mesh2triang(mesh).triangles
i, j, k = ijk[:, 0], ijk[:, 1], ijk[:, 2] 

fig = go.Figure(data=[
    go.Mesh3d(
        x=x, y=y, z=z,
        # Intensity of each vertex, which will be interpolated and color-coded
        # intensity=disp,
        intensity = stress_values,
        flatshading=True,
        # i, j and k give the vertices of triangles
        # here we represent the 4 triangles of the tetrahedron surface
        i=i, j=j, k=k,
        name='y', showscale=True,
        opacity=0.5  # make the mesh translucent
    )
])

fig.show()
2

This is a piecewise discontinuous function. I am not familiar with plotly.graph_objects but if you want to plot a discontinuous function at nodes, you are likely to encounter multiple values from the neighboring elements. The best way to visualize it would be using write_checkpoint and visualizing the result in Paraview

1 Like