Clean Plot grid code and add to documentation

This commit is contained in:
rowanc1
2014-03-03 15:50:06 -08:00
parent eeae3ec783
commit 48e506f4bb
6 changed files with 128 additions and 120 deletions
+37 -39
View File
@@ -334,7 +334,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
# Plotting Functions #
#############################################
def plotGrid(self, length=0.05, showIt=False):
def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False):
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
@@ -352,55 +352,54 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
from mpl_toolkits.mplot3d import Axes3D
mkvc = Utils.mkvc
axOpts = {'projection':'3d'} if self.dim == 3 else {}
if ax is None: ax = plt.subplot(111, **axOpts)
NN = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
fig = plt.figure(2)
fig.clf()
ax = plt.subplot(111)
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
if lines:
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
X = np.r_[X1, X2]
Y = np.r_[Y1, Y2]
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
plt.plot(X, Y)
X = np.r_[X1, X2]
Y = np.r_[Y1, Y2]
plt.hold(True)
Nx = self.r(self.normals, 'F', 'Fx', 'V')
Ny = self.r(self.normals, 'F', 'Fy', 'V')
Tx = self.r(self.tangents, 'E', 'Ex', 'V')
Ty = self.r(self.tangents, 'E', 'Ey', 'V')
ax.plot(X, Y, 'b-')
if centers:
ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
# Nx = self.r(self.normals, 'F', 'Fx', 'V')
# Ny = self.r(self.normals, 'F', 'Fy', 'V')
# Tx = self.r(self.tangents, 'E', 'Ex', 'V')
# Ty = self.r(self.tangents, 'E', 'Ey', 'V')
nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
plt.plot(nX, nY, 'r-')
# ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
#plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
plt.plot(nX, nY, 'g-')
# nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
# nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
# ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
# ax.plot(nX, nY, 'r-')
tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
plt.plot(tX, tY, 'r-')
# nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
# nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
# #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
# ax.plot(nX, nY, 'g-')
nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
#plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
plt.plot(nX, nY, 'g-')
plt.axis('equal')
# tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
# tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
# ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
# ax.plot(tX, tY, 'r-')
# nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
# nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
# #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
# ax.plot(nX, nY, 'g-')
elif self.dim == 3:
fig = plt.figure(3)
fig.clf()
ax = fig.add_subplot(111, projection='3d')
X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
@@ -417,11 +416,10 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
Y = np.r_[Y1, Y2, Y3]
Z = np.r_[Z1, Z2, Z3]
plt.plot(X, Y, 'b', zs=Z)
ax.plot(X, Y, 'b', zs=Z)
ax.set_zlabel('x3')
ax.grid(True)
ax.hold(False)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
+35 -59
View File
@@ -369,7 +369,7 @@ class TensorView(object):
if showIt: plt.show()
return out
def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False):
def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False):
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
:param bool nodes: plot nodes
@@ -399,35 +399,26 @@ class TensorView(object):
mesh.plotGrid(nodes=True, faces=True, centers=True, lines=True, showIt=True)
"""
if self.dim == 1:
fig = plt.figure(1)
fig.clf()
ax = plt.subplot(111)
xn = self.gridN
xc = self.gridCC
ax.hold(True)
ax.plot(xn, np.ones(np.shape(xn)), 'bs')
ax.plot(xc, np.ones(np.shape(xc)), 'ro')
ax.plot(xn, np.ones(np.shape(xn)), 'k--')
ax.grid(True)
ax.hold(False)
ax.set_xlabel('x1')
if showIt: plt.show()
elif self.dim == 2:
fig = plt.figure(2)
fig.clf()
ax = plt.subplot(111)
xn = self.gridN
xc = self.gridCC
xs1 = self.gridFx
xs2 = self.gridFy
ax.hold(True)
if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs')
if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro')
axOpts = {'projection':'3d'} if self.dim == 3 else {}
if ax is None: ax = plt.subplot(111, **axOpts)
if self.dim == 1:
if nodes:
ax.plot(xn, np.ones(self.nN), 'bs')
if centers:
ax.plot(xc, np.ones(self.nC), 'ro')
if lines:
ax.plot(xn, np.ones(self.nN), 'b-')
ax.set_xlabel('x1')
elif self.dim == 2:
if nodes:
ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bs')
if centers:
ax.plot(self.gridCC[:, 0], self.gridCC[:, 1], 'ro')
if faces:
ax.plot(xs1[:, 0], xs1[:, 1], 'g>')
ax.plot(xs2[:, 0], xs2[:, 1], 'g^')
ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'g>')
ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'g^')
if edges:
ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'c>')
ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'c^')
@@ -441,38 +432,23 @@ class TensorView(object):
Y2 = np.c_[mkvc(NN[1][:, 0]), mkvc(NN[1][:, self.nCy]), mkvc(NN[1][:, 0])*np.nan].flatten()
X = np.r_[X1, X2]
Y = np.r_[Y1, Y2]
plt.plot(X, Y)
ax.plot(X, Y, 'b-')
ax.grid(True)
ax.hold(False)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
if showIt: plt.show()
elif self.dim == 3:
fig = plt.figure(3)
fig.clf()
ax = fig.add_subplot(111, projection='3d')
xn = self.gridN
xc = self.gridCC
xfs1 = self.gridFx
xfs2 = self.gridFy
xfs3 = self.gridFz
xes1 = self.gridEx
xes2 = self.gridEy
xes3 = self.gridEz
ax.hold(True)
if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2])
if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2])
if nodes:
ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bs', zs=self.gridN[:, 2])
if centers:
ax.plot(self.gridCC[:, 0], self.gridCC[:, 1], 'ro', zs=self.gridCC[:, 2])
if faces:
ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2])
ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2])
ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2])
ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'g>', zs=self.gridFx[:, 2])
ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'g<', zs=self.gridFy[:, 2])
ax.plot(self.gridFz[:, 0], self.gridFz[:, 1], 'g^', zs=self.gridFz[:, 2])
if edges:
ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2])
ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2])
ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2])
ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'k>', zs=self.gridEx[:, 2])
ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'k<', zs=self.gridEy[:, 2])
ax.plot(self.gridEz[:, 0], self.gridEz[:, 1], 'k^', zs=self.gridEz[:, 2])
# Plot the grid lines
if lines:
@@ -489,14 +465,14 @@ class TensorView(object):
X = np.r_[X1, X2, X3]
Y = np.r_[Y1, Y2, Y3]
Z = np.r_[Z1, Z2, Z3]
plt.plot(X, Y, 'b-', zs=Z)
ax.grid(True)
ax.hold(False)
ax.plot(X, Y, 'b-', zs=Z)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('x3')
if showIt: plt.show()
ax.grid(True)
ax.hold(False)
if showIt: plt.show()
def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None):
assert normal in 'xyz', 'normal must be x, y, or z'
+27 -10
View File
@@ -322,7 +322,7 @@ class TreeFace(object):
if not self.isleaf: return
if self.dim == 2:
line = np.c_[self.node0.x0, self.node1.x0].T
ax.plot(line[:,0], line[:,1],'r-')
ax.plot(line[:,0], line[:,1],'b-')
if text: ax.text(self.center[0], self.center[1],self.num)
elif self.dim == 3:
if text: ax.text(self.center[0], self.center[1], self.center[2], self.num)
@@ -665,10 +665,10 @@ class TreeCell(object):
def plotGrid(self, ax, text=False):
if not self.isleaf: return
if self.dim == 2:
ax.plot(self.center[0],self.center[1],'b.')
ax.plot(self.center[0],self.center[1],'ro')
if text: ax.text(self.center[0],self.center[1],self.num)
elif self.dim == 3:
ax.plot([self.center[0]],[self.center[1]],'b.', zs=[self.center[2]])
ax.plot([self.center[0]],[self.center[1]],'ro', zs=[self.center[2]])
if text: ax.text(self.center[0], self.center[1], self.center[2], self.num)
@@ -1048,21 +1048,38 @@ class TreeMesh(InnerProducts, BaseMesh):
zP = self._getEdgeP(zEdge, xEdge, yEdge)
return sp.vstack((xP, yP, zP))
def plotGrid(self, ax=None, text=True, plotC=True, plotF=True, plotE=False, plotEx=False, plotEy=False, plotEz=False, showIt=False):
def plotGrid(self, ax=None, text=False, centers=False, faces=False, edges=False, lines=True, nodes=False, showIt=False):
self.number()
axOpts = {'projection':'3d'} if self.dim == 3 else {}
if ax is None: ax = plt.subplot(111, **axOpts)
if plotC: [c.plotGrid(ax, text=text) for c in self.cells]
if plotF: [f.plotGrid(ax, text=text) for f in self.faces]
if plotE and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edges]
if plotEx and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesX]
if plotEy and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesY]
if plotEz and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesZ]
if lines:
[f.plotGrid(ax, text=text) for f in self.faces]
if centers:
[c.plotGrid(ax, text=text) for c in self.cells]
if faces:
fX = np.array([f.center for f in self.sortedFaceX])
ax.plot(fX[:,0],fX[:,1],'g>')
fY = np.array([f.center for f in self.sortedFaceY])
ax.plot(fY[:,0],fY[:,1],'g^')
if edges:
eX = np.array([e.center for e in self.sortedFaceY])
ax.plot(eX[:,0],eX[:,1],'c>')
eY = np.array([e.center for e in self.sortedFaceX])
ax.plot(eY[:,0],eY[:,1],'c^')
if nodes:
ns = np.array([n.x0 for n in self.sortedNodes])
ax.plot(ns[:,0],ns[:,1],'bs')
ax.set_xlim((self.x0[0], self.h[0].sum()))
ax.set_ylim((self.x0[1], self.h[1].sum()))
if self.dim == 3:
ax.set_zlim((self.x0[2], self.h[2].sum()))
ax.grid(True)
ax.hold(False)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
if showIt: plt.show()
def plotImage(self, I, ax=None, showIt=True):
+3 -3
View File
@@ -11,18 +11,18 @@ def exampleLrmGrid(nC, exType):
assert exType in possibleTypes, "Not a possible example type."
if exType == 'rect':
return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
return list(ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False))
elif exType == 'rotate':
if len(nC) == 2:
X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2)
amt[amt < 0] = 0
return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt
return [X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt]
elif len(nC) == 3:
X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2)
amt[amt < 0] = 0
return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt
return [X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt]
def meshTensors(*args):
+17
View File
@@ -21,7 +21,24 @@ Each mesh code follows the guiding principles that are present in this
tutorial, but the details, advantages and disadvantages differ between
the implementations.
.. plot::
from SimPEG import Mesh, Utils, np
import matplotlib.pyplot as plt
sz = [10,10]
tM = Mesh.TensorMesh(sz)
qM = Mesh.TreeMesh(sz)
qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.3 else 0)
rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
fig, axes = plt.subplots(1,3,figsize=(14,5))
opts = {}
tM.plotGrid(ax=axes[0], **opts)
axes[0].set_title('TensorMesh')
qM.plotGrid(ax=axes[1], **opts)
axes[1].set_title('TreeMesh')
rM.plotGrid(ax=axes[2], **opts)
axes[2].set_title('LogicallyRectMesh')
Variable Locations and Terminology
+9 -9
View File
@@ -9,15 +9,6 @@ Tensor Mesh
:undoc-members:
Cylindrical 1D Mesh
===================
.. automodule:: SimPEG.Mesh.Cyl1DMesh
:show-inheritance:
:members:
:undoc-members:
Logically Rectangular Mesh
==========================
@@ -27,6 +18,15 @@ Logically Rectangular Mesh
:undoc-members:
Cylindrical 1D Mesh
===================
.. automodule:: SimPEG.Mesh.Cyl1DMesh
:show-inheritance:
:members:
:undoc-members:
Base Mesh
=========