def showSolution(self,dim=2):
     from mpl_toolkits.mplot3d import Axes3D
     from matplotlib import cm, pyplot
     x,y,tri,solution,grad = [],[],[],[],[]
     for n in self.node:
         x.append(n.x)
         y.append(n.y)
         solution.append(n.value)
     for e in self.element:
         tri.append([n.id for n in e.node])
         grad.append(e.grad())
     if dim==2:
         pyplot.figure(figsize=(17,6))
         pyplot.subplot(1,2,1)
         pyplot.title("Solution")
         pyplot.tripcolor(x, y, tri, solution, cmap=cm.jet,  edgecolors='black')
         pyplot.colorbar()
         pyplot.subplot(1,2,2)
         pyplot.title("Gradient")
         pyplot.tripcolor(x, y, tri, grad, cmap=cm.jet,  edgecolors='black')
         pyplot.colorbar()
     elif dim==3:
         fig = pyplot.figure()
         ax = fig.gca(projection='3d')
         ax.plot_trisurf(x, y, tri, z, cmap=cm.jet, linewidth=0.2)
     pyplot.show()

def el_plot(data, Map=False, show=True):
    """
    Plot the elevation for the region from the last time series
    
    :Parameters:
        **data** -- the standard python data dictionary
        
        **Map** -- {True, False} (optional): Optional argument.  If True,
            the elevation will be plotted on a map.  
    """
    trigrid = data['trigrid']
    plt.gca().set_aspect('equal')
    plt.tripcolor(trigrid, data['zeta'][-1,:])
    plt.colorbar()
    plt.title("Elevation")
    if Map:
        #we set the corners of where the map should show up
        llcrnrlon, urcrnrlon = plt.xlim()
        llcrnrlat, urcrnrlat = plt.ylim()
        #we construct the map.  Note that resolution serves to increase
        #or decrease the detail in the coastline.  Currently set to 
        #'i' for 'intermediate'
        m = Basemap(llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat, \
            resolution='i', suppress_ticks=False)
        #set color for continents.  Default is grey.
        m.fillcontinents(color='ForestGreen')
        m.drawmapboundary()
        m.drawcoastlines()
    if show:
        plt.show()

def field(domain, z, title, clim = None,  path=None, save=True, show =
          False, ics=1, ext='.png', cmap=plt.cm.jet, fmt=None):
    """
    Given a domain, plot the nodal value z
   
    :param domain: :class:`~polyadcirc.run_framework.domain`
    :param z: :class:`numpy.ndarray`
    :param string title: plot title
    :param clim: :class:`numpy.clim`
    :type path: string or None
    :param path: directory to store plots
    :type save: bool
    :param save: flag for whether or not to save plots
    :type show: bool
    :param show: flag for whether or not to show plots
    :param int ics:  polar coordinate option (1 = cart coords, 2 = polar
        coords)
    :param string ext: file extension
    :param callable fmt: formatter for color bar, takes ``(x, pos)`` returns
        ``string`` 

    """
    if path is None:
        path = os.getcwd()
    plt.figure()
    plt.tripcolor(domain.triangulation, z, shading='gouraud',
                  cmap=cmap)
    plt.gca().set_aspect('equal')
    plt.autoscale(tight=True)
    if clim:
        plt.clim(clim[0], clim[1])
    add_2d_axes_labels(ics=ics)    
    colorbar(fmt=fmt)
    #plt.title(title)
    save_show(os.path.join(path, 'figs', title), save, show, ext)

 def draw(self, x = (), y=()):
     #px,py = self.points.to(self.env_model.ROW_COL)
     plt.tripcolor(self.x, self.y, self.mesh.faces, facecolors=self.zfaces, edgecolors='k')
     if len(x) != 0:
         plt.plot(x, y)
     plt.axis('equal')
     plt.show()

def plotstate(Mesh, U, field, fname):
	V = Mesh['V']
	E = Mesh['E']
	BE = Mesh['BE']

	f = plt.figure(figsize=(12,6))

	F = pu.getField(U, field)
	plt.tripcolor(V[:,0], V[:,1], triangles=E, facecolors=F, shading='flat', vmin=1.55, vmax=1.9)

	for i in range(len(BE)):
		x = [V[BE[i,0],0], V[BE[i,1],0]]
		y = [V[BE[i,0],1], V[BE[i,1],1]]
		plt.plot(x, y, '-', linewidth=2, color='black')
	
	#dosave = (len(fname) != 0)

	plt.axis('equal')

	#plt.axis([-100, 100,-100, 100])
	plt.axis([-2, 10,-4, 4])
	plt.colorbar()
	#plt.clim(0, 0.8)
	plt.title(field, fontsize=16)

	f.tight_layout()
	plt.show()#block=(not dosave))
	#if (dosave):
	plt.savefig(fname)
	
	plt.close(f)

def quick_plot(bmi, name, **kwds):
    gid = bmi.var_grid(name)
    gtype = bmi.grid_type(gid)
    grid = bmi.grid[gid]

    x, y = grid.node_x.values, grid.node_y.values
    z = bmi.get_value(name)

    x_label = "{name} ({units})".format(
        name=grid.node_x.standard_name, units=grid.node_x.units
    )
    y_label = "{name} ({units})".format(
        name=grid.node_y.standard_name, units=grid.node_y.units
    )

    if gtype in ("unstructured_triangular",):
        tris = bmi.grid_face_node_connectivity(gid).reshape((-1, 3))
        plt.tripcolor(x, y, tris, z, **kwds)
    elif gtype in ("uniform_rectilinear", "structured_quad"):
        shape = bmi.grid_shape(gid)
        spacing = bmi.grid_spacing(gid)
        origin = bmi.grid_origin(gid)
        x = np.arange(shape[-1]) * spacing[-1] + origin[-1]
        y = np.arange(shape[-2]) * spacing[-2] + origin[-2]
        plt.pcolormesh(x, y, z.reshape(shape), **kwds)
    else:
        raise ValueError("no plotter for {gtype}".format(gtype=gtype))

    plt.axis("tight")
    plt.gca().set_aspect("equal")
    plt.xlabel(x_label)
    plt.ylabel(y_label)

    cbar = plt.colorbar()
    cbar.ax.set_ylabel("{name} ({units})".format(name=name, units=bmi.var_units(name)))

def field(domain, z, title, clim = None,  path=None, save=True, show =
          False, ics=1, ext='.png', cmap=plt.cm.jet):
    """
    Given a domain, plot the nodal value z
   
    :param domain: :class:`~polyadcirc.run_framework.domain`
    :param z: :class:`np.array`
    :param string title: plot title
    :param clim: :class:`np.clim`
    :type path: string or None
    :param path: directory to store plots
    :type save: boolean
    :param save: flag for whether or not to save plots
    :type show: boolean
    :param show: flag for whether or not to show plots
    :param int ics:  polar coordinate option (1 = cart coords, 2 = polar
        coords)
    :param string ext: file extension

    """
    if path == None:
        path = os.getcwd()
    plt.figure()
    plt.tripcolor(domain.triangulation, z, shading='gouraud',
                  cmap=cmap)
    plt.gca().set_aspect('equal')
    plt.autoscale(tight=True)
    if clim:
        plt.clim(clim[0], clim[1])
    add_2d_axes_labels(ics)    
    colorbar()
    #plt.title(title)
    save_show(path+'/figs/'+title, save, show, ext)

def mplot_function(f, vmin, vmax, logscale):
    mesh = f.function_space().mesh()
    if (mesh.geometry().dim() != 2):
        raise AttributeError('Mesh must be 2D')
    # DG0 cellwise function
    if f.vector().size() == mesh.num_cells():
        C = f.vector().array()
        if logscale:
            return plt.tripcolor(mesh2triang(mesh), C, vmin=vmin, vmax=vmax, norm=cls.LogNorm() )
        else:
            return plt.tripcolor(mesh2triang(mesh), C, vmin=vmin, vmax=vmax)
    # Scalar function, interpolated to vertices
    elif f.value_rank() == 0:
        C = f.compute_vertex_values(mesh)
        if logscale:
            return plt.tripcolor(mesh2triang(mesh), C, vmin=vmin, vmax=vmax, norm=cls.LogNorm() )
        else:
            return plt.tripcolor(mesh2triang(mesh), C, shading='gouraud', vmin=vmin, vmax=vmax)
    # Vector function, interpolated to vertices
    elif f.value_rank() == 1:
        w0 = f.compute_vertex_values(mesh)
        if (len(w0) != 2*mesh.num_vertices()):
            raise AttributeError('Vector field must be 2D')
        X = mesh.coordinates()[:, 0]
        Y = mesh.coordinates()[:, 1]
        U = w0[:mesh.num_vertices()]
        V = w0[mesh.num_vertices():]
        C = np.sqrt(U*U+V*V)
        return plt.quiver(X,Y,U,V, C, units='x', headaxislength=7, headwidth=7, headlength=7, scale=4, pivot='middle')

def test_tripcolor():
    x = np.asarray([0, 0.5, 1, 0,   0.5, 1,   0, 0.5, 1, 0.75])
    y = np.asarray([0, 0,   0, 0.5, 0.5, 0.5, 1, 1,   1, 0.75])
    triangles = np.asarray([
        [0, 1, 3], [1, 4, 3],
        [1, 2, 4], [2, 5, 4],
        [3, 4, 6], [4, 7, 6],
        [4, 5, 9], [7, 4, 9], [8, 7, 9], [5, 8, 9]])

    # Triangulation with same number of points and triangles.
    triang = mtri.Triangulation(x, y, triangles)

    Cpoints = x + 0.5*y

    xmid = x[triang.triangles].mean(axis=1)
    ymid = y[triang.triangles].mean(axis=1)
    Cfaces = 0.5*xmid + ymid

    plt.subplot(121)
    plt.tripcolor(triang, Cpoints, edgecolors='k')
    plt.title('point colors')

    plt.subplot(122)
    plt.tripcolor(triang, facecolors=Cfaces, edgecolors='k')
    plt.title('facecolors')

def triangle_vert_plot(filename, args):
    dset = gusto_dataset.GustoDataset(filename)
    x = dset.get_vert_variable('x1')
    z = dset.get_vert_variable('x3')
    f = dset.get_vert_variable(args.data)
    plt.tripcolor(x, z, f)
    plt.axis('equal')

def draw_map(triangulation, options):
    '''
    get the triangle tuple (concentration, triangle] prepared before
    and draw the map of triangles
    options :
    "map_format": "svg",
    "map_file": "../../mapa"
    '''
    lab_x = options['xlabel'] if value_set(options, 'xlabel') else 'mesh X coord'
    lab_y = options['ylabel'] if value_set(options, 'ylabel') else 'mesh Y coord'
    lab_tit = options['title'] if value_set(options, 'title') else 'Map of concentrations'
    
    
    plt.figure()
    plt.gca().set_aspect('equal')
    plt.tripcolor(triangulation['x_np'],
                  triangulation['y_np'],
                  triangulation['triangles'],
                  facecolors=triangulation['zfaces'],
                  edgecolors='k')
    plt.colorbar()
    plt.title(lab_tit)
    plt.xlabel(lab_x)
    plt.ylabel(lab_y)
    
    plt.savefig(options["map_file"], format=options["map_format"])               

def plot_latlon_tri(lon=None, lat=None, data=None, title='Title', 
	vmin_in=CBAR_MINT, vmax_in=CBAR_MAXT):
	triang = tri.Triangulation(lon, lat)
	fig, ax = plt.subplots()
	plt.gca().set_aspect('equal')
	plt.tripcolor(triang, data, cmap=cm.jet, vmin=vmin_in, vmax=vmax_in)
	label_plot(fig, ax, title)
	return fig

def triangle_variable_plot(filename, args):
    dset = gusto_dataset.GustoDataset(filename)
    x = dset.get_cell_variable('x1')
    z = dset.get_cell_variable('x3')
    f = dset.get_cell_variable('dg')
    plt.tripcolor(x, z, f)
    plt.axis('equal')
    plt.colorbar()

def plot_faces(nodes,faces,fn,cmap=None):
    fn = np.ma.array(fn,mask=~np.isfinite(fn))
    assert(fn.size == faces.shape[0])
    if cmap is None:
        cmap = plt.get_cmap('jet')
    cmap.set_bad('w',1.)
    plt.gca()
    plt.tripcolor(nodes[:,0],nodes[:,1],faces,facecolors=fn,
                  edgecolor='k',cmap=cmap)
    plt.colorbar()

def postprocessor(nodes, val):
    x = nodes[:, 0]
    y = nodes[:, 1]

    fig = plt.figure()
    plt.gca().set_aspect('equal')
    plt.tripcolor(x,y,elements,facecolors=val,edgecolors='k')
    plt.colorbar()
    plt.title("Poisson's equation")
    plt.xlabel('x')
    plt.ylabel('y')
    #plt.savefig("FEM.png")
    plt.show()

def highlightTriangle(index): #give this function the index of element and it will highlight
    clickedTriangleX = [x[elements[index][0]],x[elements[index][1]],x[elements[index][2]]]
    clickedTriangleY = [y[elements[index][0]],y[elements[index][1]],y[elements[index][2]]]
    plt.tripcolor(clickedTriangleX,clickedTriangleY,[0,1,2],cmap=plt.cm.Accent)
    a = elements[index]
    x1 = x[a[0]]
    x2 = x[a[1]]
    x3 = x[a[2]]
    y1 = y[a[0]]
    y2 = y[a[1]]
    y3 = y[a[2]]
    centroid = tuple(((x1+x2+x3)/3,(y1+y2+y3)/3))
    if centroid not in visible:
        visible[centroid] = plt.text(centroid[0], centroid[1], '.' + str(index))

def onclick(event): #print('button=%d, x=%d, y=%d, xdata=%f, ydata=%f'%(event.button, event.x, event.y, event.xdata, event.ydata))
    #print(event.button, event.x, event.y, event.xdata, event.ydata)
    if event.xdata < maxX and event.xdata > minX and event.ydata < maxY and event.ydata > minY:
        if event.key in ['alt+q','q'] or event.key in ['alt+w','w']:
            # TODO: put clicked triangles into an array
            # after a certain number of clicks, we want to remake the plot, to prevent recursion depth error
            if event.xdata < maxX or event.xdata > minX or event.ydata < maxY or event.ydata > minY:
                jump = index = len(orderedCentroids) // 2
                
                # Complexity of log(n) - ordered binary search to get to the nearest X
                while True:
                    jump = jump //2
                    if jump < 1 or event.xdata == orderedCentroids[index][0]:
                        break
                    index += jump * (int(event.xdata > orderedCentroids[index][0]) - int(event.xdata < orderedCentroids[index][0]))
                
                # Bounce outward checking if point is inside triangle who's centroid is centroid of list (to skip Y checking)
                jump = 1
                while True:
                    if inTriangle(event.xdata,event.ydata,allTriangles[orderedCentroids[index]]):
                        break
                    index += jump * (1-2*(jump % 2 == 0))
                    jump += 1

                # Find which 3 points make up the triangle
                clickedTriangle = allTriangles[orderedCentroids[index]]
                clickedTriangleX = [clickedTriangle[0][0],clickedTriangle[1][0],clickedTriangle[2][0]]
                clickedTriangleY = [clickedTriangle[0][1],clickedTriangle[1][1],clickedTriangle[2][1]]

            if event.key in ['alt+q','q']:
                if orderedCentroids[index] not in visible:
                    plt.tripcolor(clickedTriangleX,clickedTriangleY,[0,1,2],cmap=plt.cm.Accent)
                    visible[orderedCentroids[index]] = plt.text(orderedCentroids[index][0], orderedCentroids[index][1], '.' + str(clickedTriangle[3]))
            elif event.key in ['alt+w','w']:
                if orderedCentroids[index] in visible:
                    plt.tripcolor(clickedTriangleX,clickedTriangleY,[0,1,2],cmap=plt.cm.Greys)
                    visible[orderedCentroids[index]].remove()
                    del visible[orderedCentroids[index]]
            
        elif event.key in ['alt+e','e']:
            typedInput = input('Enter triangle index to highlight:  ')
            try:
                typedInput = int(typedInput)
                highlightTriangle(typedInput)
            except:
                print('Invalid Input')

        plt.draw()

def plot_f_m3d_hdf5_2D(filename,time_id,data_id,plane_id):
	import matplotlib.pyplot as plt
	import matplotlib.tri as triangle
	from matplotlib.pyplot import figure, axes, plot, xlabel, ylabel, title, \
		grid, savefig, show
	(timeframe,nnode_1_plane,R,Z,trianlist,data_name,data_2D)=readhdf5_f_time_data(filename,time_id,data_id,plane_id)
	#plt.gca().set_aspect('equal')
	tri=triangle.Triangulation(R,Z)
	fig = figure(figsize=(6,8))
	plt.tripcolor(tri,data_2D)
	plt.xlabel('R')
	plt.ylabel('Z');
	title_str=format('%s, t=%f'%(data_name,timeframe))
	plt.title(title_str)
	plt.colorbar()
	plt.show()

def plot_fill(value_array,title,sym=0):
    x = nodearray[:,0]
    y = nodearray[:,1]
    triangles = np.array(elements)
    
    if sym:
        x = np.append(x,-x)
        y = np.append(y,y)
        triangles = np.vstack((triangles, triangles+np.amax(triangles)+1))
        value_array = np.append(value_array,value_array) 
    
    plt.figure()
    plt.title(title)
    plt.axes().set_aspect('equal')        
    plt.tripcolor(x,y,triangles,value_array, edgecolors='k',cmap=plt.get_cmap('jet'))
    plt.colorbar()

  def plot(self, savename, astr='.*_fcomp(.*)_focean(.*)_CPL', xlabel='fcomp', ylabel='focean'):
    """ astr is used to build indices from files"""
    matcher = re.compile(astr)
    matches = [(matcher.findall(af)[0][0], matcher.findall(af)[0][1], val) for af, val in zip(self.files, self.vals) if matcher.findall(af) != []]
    points = np.squeeze(np.asarray(matches, dtype='f8'))

    plt.figure()
    tri = Triangulation(points[:,0], points[:,1])
    plt.tripcolor(tri, points[:,2], cmap='viridis')
    plt.colorbar()
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(self.units)
    plt.savefig(savename + '.png')

    return #}}}