# Wrong plot with plot_trisurf

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 =  * (3*nel[k])
YY =  * (3*nel[k])
ZZ =  * (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]
x1 = dof2coord[dof,0]
x2 = dof2coord[dof,0]
y0 = dof2coord[dof,1]
y1 = dof2coord[dof,1]
y2 = dof2coord[dof,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))