From 48e506f4bb12cae57f980a5f89e1ad68212f96e7 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 3 Mar 2014 15:50:06 -0800 Subject: [PATCH] Clean Plot grid code and add to documentation --- SimPEG/Mesh/LogicallyRectMesh.py | 76 +++++++++++++------------- SimPEG/Mesh/TensorView.py | 94 ++++++++++++-------------------- SimPEG/Mesh/TreeMesh.py | 37 +++++++++---- SimPEG/Utils/meshutils.py | 6 +- docs/api_Mesh.rst | 17 ++++++ docs/api_MeshCode.rst | 18 +++--- 6 files changed, 128 insertions(+), 120 deletions(-) diff --git a/SimPEG/Mesh/LogicallyRectMesh.py b/SimPEG/Mesh/LogicallyRectMesh.py index d606fe78..41dd2641 100644 --- a/SimPEG/Mesh/LogicallyRectMesh.py +++ b/SimPEG/Mesh/LogicallyRectMesh.py @@ -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') diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index 4c04cda9..3cc98b61 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -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' diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index ae83eae8..8ac841d4 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -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): diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 48dcb8cf..67d747c3 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -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): diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index 46bd877b..e5a72491 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -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 diff --git a/docs/api_MeshCode.rst b/docs/api_MeshCode.rst index 090bb6c2..65820fc8 100644 --- a/docs/api_MeshCode.rst +++ b/docs/api_MeshCode.rst @@ -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 =========