I am trying to make a surface plot of a solution on posed on a domain which is a nonconvex polygon, specifically, the square [-1,1]x[-1,1] minus the square [0,1]x[0,1]. I am using this code (which seems to be too long!)
#####
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#####
#####
# graphics
coordinates = mesh2d.coordinates()
XX = [0] * (3*nel[k])
YY = [0] * (3*nel[k])
ZZ = [0] * (3*nel[k])
dofmap = V.dofmap()
nodal_values = u.vector()
dof2coord= coordinates[dof_to_vertex_map(V)]
for cell in cells(mesh2d):
dof = dofmap.cell_dofs(cell.index())
x0 = dof2coord[dof[0],0]
x1 = dof2coord[dof[1],0]
x2 = dof2coord[dof[2],0]
y0 = dof2coord[dof[0],1]
y1 = dof2coord[dof[1],1]
y2 = dof2coord[dof[2],1]
#check orientation
if (x1-x0)*(y2-y0)-(x2-x0)*(y1-y0)>0:
i0 = 0
i1 = 1
i2 = 2
else:
i0 = 0
i1 = 2
i2 = 1
XX[3*cell.index()] = dof2coord[dof[i0],0]
XX[3*cell.index()+1] = dof2coord[dof[i1],0]
XX[3*cell.index()+2] = dof2coord[dof[i2],0]
YY[3*cell.index()] = dof2coord[dof[i0],1]
YY[3*cell.index()+1] = dof2coord[dof[i1],1]
YY[3*cell.index()+2] = dof2coord[dof[i2],1]
ZZ[3*cell.index()] = nodal_values[dof[i0]]
ZZ[3*cell.index()+1] = nodal_values[dof[i1]]
ZZ[3*cell.index()+2] = nodal_values[dof[i2]]
fig = plt.figure(figsize =(16, 9))
ax = fig.add_subplot(projection='3d')
# Creating color map
my_cmap = plt.get_cmap('summer')
trisurf = ax.plot_trisurf(XX, YY, ZZ, linewidth=0.2, antialiased=False, cmap=my_cmap)
fig.colorbar(trisurf, ax = ax, shrink = 0.5, aspect = 5)
plt.show()
It results that it draws some triangles which do not correspond to the plot, I mean, close to the nonconvex vertex there appear wrong triangles, like the ones showed in this figure:
I wanted to know if I am using in a wrong way the command plot_trisurf. Thanks!