Index: /issm/trunk-jpl/src/m/plot/plot_unit.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 21445)
+++ /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 21446)
@@ -1,4 +1,5 @@
 from cmaptools import truncate_colormap
 from plot_quiver import plot_quiver
+from scipy.interpolate import griddata
 import numpy as  np
 try:
@@ -22,5 +23,7 @@
 	#if we are plotting 3d replace the current axis
 	if not is2d:
+		axgrid[gridindex].axis('off')
 		ax=inset_locator.inset_axes(axgrid[gridindex],width='100%',height='100%',loc=3,borderpad=0,axes_class=Axes3D)
+		ax.set_axis_bgcolor((0.7,0.7,0.7))
 	else:
 		ax=axgrid[gridindex]
@@ -91,4 +94,6 @@
 			tri=ax.tripcolor(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,edgecolors=edgecolor)
 		else:
+
+
 			#first deal with colormap
 			loccmap = plt.cm.ScalarMappable(cmap=cmap)
@@ -103,8 +108,8 @@
 			recface=np.vstack((recface,np.vstack((elements[:,2],elements[:,0],elements[:,3],elements[:,5])).T))
 			eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0])))
-			tmp = np.ascontiguousarray(recface).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1])))
-			_, idx = np.unique(tmp, return_index=True)
-			recel= recface[idx]
-			recindex=eltind[idx]
+			tmp = np.ascontiguousarray(np.sort(recface)).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1])))
+			_, idx, recur = np.unique(tmp, return_index=True, return_counts=True)
+			recel= recface[idx[np.where(recur==1)]]
+			recindex=eltind[idx[np.where(recur==1)]]
 			for i,rectangle in enumerate(recel):
 				rec=zip(x[rectangle],y[rectangle],z[rectangle])
@@ -120,7 +125,8 @@
 			eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0])))
 			tmp = np.ascontiguousarray(triface).view(np.dtype((np.void, triface.dtype.itemsize * triface.shape[1])))
-			_, idx = np.unique(tmp, return_index=True)
-			triel= triface[idx]
-			triindex=eltind[idx]
+			_, idx,recur = np.unique(tmp, return_index=True,return_counts=True)
+			#we keep only top and bottom elements
+			triel= triface[idx[np.where(recur==1)]]
+			triindex=eltind[idx[np.where(recur==1)]]
 			for i,triangle in enumerate(triel):
 				tri=zip(x[triangle],y[triangle],z[triangle])
@@ -154,12 +160,12 @@
 			loccmap.set_array([min(data),max(data)])
 			loccmap.set_clim(vmin=min(data),vmax=max(data))
-
+			
 			#deal with prism sides
 			recface=np.vstack((elements[:,0],elements[:,1],elements[:,4],elements[:,3])).T
 			recface=np.vstack((recface,np.vstack((elements[:,1],elements[:,2],elements[:,5],elements[:,4])).T))
 			recface=np.vstack((recface,np.vstack((elements[:,2],elements[:,0],elements[:,3],elements[:,5])).T))
-			tmp = np.ascontiguousarray(recface).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1])))
-			_, idx = np.unique(tmp, return_index=True)
-			recel= recface[idx]					
+			tmp = np.ascontiguousarray(np.sort(recface)).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1])))
+			_, idx, recur = np.unique(tmp, return_index=True, return_counts=True)
+			recel= recface[idx[np.where(recur==1)]]
 			for rectangle in recel:
 				rec=zip(x[rectangle],y[rectangle],z[rectangle])
@@ -169,10 +175,11 @@
 				pl3.set_color(color)
 				ax.add_collection3d(pl3)
-
+				
 			#deal with prism faces
 			triface=np.vstack((elements[:,0:3],elements[:,3:6]))
 			tmp = np.ascontiguousarray(triface).view(np.dtype((np.void, triface.dtype.itemsize * triface.shape[1])))
-			_, idx = np.unique(tmp, return_index=True)
-			triel= triface[idx]
+			_, idx,recur = np.unique(tmp, return_index=True,return_counts=True)
+			#we keep only top and bottom elements
+			triel= triface[idx[np.where(recur==1)]]
 			for triangle in triel:
 				tri=zip(x[triangle],y[triangle],z[triangle])
@@ -182,5 +189,5 @@
 				pl3.set_color(color)
 				ax.add_collection3d(pl3)
-	
+				
 			ax.set_xlim([min(x),max(x)])
 			ax.set_ylim([min(y),max(y)])