Plotting error: calculating principal stresses

1

Hey all - I am modeling simple 2D Stokes flow around a circular barrier:

stress_intensity
stress_intensity1280×960 33.3 KB

I am trying to calculate the maximum tensile strength (first principal stress) across the mesh. I solve for u, p using the residual method as found here and I calculate the maximum tensile stress like so:

def sigma(u,p): # Cauchy stress tensor
  return 2*eta*epsilon(u) - p*I

TS = TensorFunctionSpace(mesh, "CG", 2) 
cauchy = project(sigma(u,p), TS)
sigma_x = cauchy[0,0] 
sigma_y = cauchy[1, 1]
tau_xy = cauchy[0, 1]

def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
    return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy))

ps1 = sigma1(sigma_x, sigma_y, tau_xy)
plot(ps1, title="Maximum tensile stress")

I get the following error:

File "/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/common/plotting.py", line 438, in plot
    return _plot_matplotlib(object, mesh, kwargs)
  File "/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/common/plotting.py", line 304, in _plot_matplotlib
    return mplot_function(ax, obj, **kwargs)
  File "/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/common/plotting.py", line 152, in mplot_function
    return ax.tricontourf(mesh2triang(mesh), C, levels, **kwargs)
  File "/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/tri/tricontour.py", line 275, in tricontourf
    return TriContourSet(ax, *args, **kwargs)
  File "/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/tri/tricontour.py", line 35, in __init__
    ContourSet.__init__(self, ax, *args, **kwargs)
  File "/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/contour.py", line 822, in __init__
    kwargs = self._process_args(*args, **kwargs)
  File "/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/tri/tricontour.py", line 47, in _process_args
    tri, z = self._contour_args(args, kwargs)
  File "/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/tri/tricontour.py", line 94, in _contour_args
    raise ValueError('z array must not contain non-finite values '
ValueError: z array must not contain non-finite values within the triangulation

I believe this is because of the irrational values generated by the sqrt function required to calculate sigma1. I know that UFL does not provide any floor, ceil, or round operators. Is there any way to round irrational values to a finite values so that they can be plotted?

Thanks!

2

Shouldn’t tau_xy be squared? (With the current formulation, you could be getting negative values within the sqrt.)

3

LOL. I feel kind of dumb now. That was it!! Thank you so much.