diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py deleted file mode 100644 index 110fe9bd..00000000 --- a/SimPEG/DiffOperators.py +++ /dev/null @@ -1,336 +0,0 @@ -import numpy as np -from scipy import sparse as sp -from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros - - -def ddx(n): - """Define 1D derivatives, inner, this means we go from n+1 to n+1""" - return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") - - -def checkBC(bc): - """ Checks if boundary condition 'bc' is valid. """ - if(type(bc) is str): - bc = [bc, bc] - assert type(bc) is list, 'bc must be a list' - assert len(bc) == 2, 'bc must have two elements' - - for bc_i in bc: - assert type(bc_i) is str, "each bc must be a string" - assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'" - return bc - - -def ddxCellGrad(n, bc): - """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1""" - bc = checkBC(bc) - - D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") - # Set the first side - if(bc[0] == 'dirichlet'): - D[0, 0] = 2 - elif(bc[0] == 'neumann'): - D[0, 0] = 0 - # Set the second side - if(bc[1] == 'dirichlet'): - D[-1, -1] = -2 - elif(bc[1] == 'neumann'): - D[-1, -1] = 0 - return D - - -def av(n): - """Define 1D averaging operator from cell-centres to nodes.""" - return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") - - -class DiffOperators(object): - """ - Class creates the differential operators that you need! - """ - def __init__(self): - raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') - - def faceDiv(): - doc = "Construct divergence operator (face-stg to cell-centres)." - - def fget(self): - if(self._faceDiv is None): - # The number of cell centers in each direction - n = self.n - # Compute faceDivergence operator on faces - if(self.dim == 1): - D = ddx(n[0]) - elif(self.dim == 2): - D1 = sp.kron(speye(n[1]), ddx(n[0])) - D2 = sp.kron(ddx(n[1]), speye(n[0])) - D = sp.hstack((D1, D2), format="csr") - elif(self.dim == 3): - D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) - D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) - D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) - D = sp.hstack((D1, D2, D3), format="csr") - # Compute areas of cell faces & volumes - S = self.area - V = self.vol - self._faceDiv = sdiag(1/V)*D*sdiag(S) - - return self._faceDiv - return locals() - _faceDiv = None - faceDiv = property(**faceDiv()) - - def nodalGrad(): - doc = "Construct gradient operator (nodes to edges)." - - def fget(self): - if(self._nodalGrad is None): - # The number of cell centers in each direction - n = self.n - # Compute divergence operator on faces - if(self.dim == 1): - G = ddx(n[0]) - elif(self.dim == 2): - D1 = sp.kron(speye(n[1]+1), ddx(n[0])) - D2 = sp.kron(ddx(n[1]), speye(n[0]+1)) - G = sp.vstack((D1, D2), format="csr") - elif(self.dim == 3): - D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0])) - D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1)) - D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1)) - G = sp.vstack((D1, D2, D3), format="csr") - # Compute lengths of cell edges - L = self.edge - self._nodalGrad = sdiag(1/L)*G - return self._nodalGrad - return locals() - _nodalGrad = None - nodalGrad = property(**nodalGrad()) - - def setCellGradBC(self, BC): - """ - Function that sets the boundary conditions for cell-centred derivative operators. - - Examples:: - - BC = 'neumann' # Neumann in all directions - BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else - BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain, - # Dirichlet else - - """ - if(type(BC) is str): - BC = [BC for _ in self.n] # Repeat the str self.dim times - elif(type(BC) is list): - assert len(BC) == self.dim, 'BC list must be the size of your mesh' - else: - raise Exception("BC must be a str or a list.") - - for i, bc_i in enumerate(BC): - BC[i] = checkBC(bc_i) - - self._cellGrad = None # ensure we create a new gradient next time we call it - self._cellGradBC = BC - return BC - _cellGradBC = 'neumann' - - def cellGrad(): - doc = "The cell centered Gradient, takes you to cell faces." - - def fget(self): - if(self._cellGrad is None): - BC = self.setCellGradBC(self._cellGradBC) - n = self.n - if(self.dim == 1): - G = ddxCellGrad(n[0], BC[0]) - elif(self.dim == 2): - G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0])) - G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0])) - G = sp.vstack((G1, G2), format="csr") - elif(self.dim == 3): - G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0])) - G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0])) - G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0])) - G = sp.vstack((G1, G2, G3), format="csr") - # Compute areas of cell faces & volumes - S = self.area - V = self.vol - self._cellGrad = sdiag(S)*G*sdiag(1/V) - return self._cellGrad - return locals() - _cellGrad = None - cellGrad = property(**cellGrad()) - - def edgeCurl(): - doc = "Construct the 3D curl operator." - - def fget(self): - if(self._edgeCurl is None): - # The number of cell centers in each direction - n1 = self.nCx - n2 = self.nCy - n3 = self.nCz - - # Compute lengths of cell edges - L = self.edge - - # Compute areas of cell faces - S = self.area - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - - D32 = kron3(d3, speye(n2), speye(n1+1)) - D23 = kron3(speye(n3), d2, speye(n1+1)) - D31 = kron3(d3, speye(n2+1), speye(n1)) - D13 = kron3(speye(n3), speye(n2+1), d1) - D21 = kron3(speye(n3+1), d2, speye(n1)) - D12 = kron3(speye(n3+1), speye(n2), d1) - - O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) - O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) - O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) - - C = sp.vstack((sp.hstack((O1, -D32, D23)), - sp.hstack((D31, O2, -D13)), - sp.hstack((-D21, D12, O3))), format="csr") - - self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) - return self._edgeCurl - return locals() - _edgeCurl = None - edgeCurl = property(**edgeCurl()) - - def aveF2CC(): - doc = "Construct the averaging operator on cell faces to cell centers." - - def fget(self): - if(self._aveF2CC is None): - n = self.n - if(self.dim == 1): - self._aveF2CC = av(n[0]) - elif(self.dim == 2): - self._aveF2CC = sp.hstack((sp.kron(speye(n[1]), av(n[0])), - sp.kron(av(n[1]), speye(n[0]))), format="csr") - elif(self.dim == 3): - self._aveF2CC = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") - return self._aveF2CC - return locals() - _aveF2CC = None - aveF2CC = property(**aveF2CC()) - - def aveE2CC(): - doc = "Construct the averaging operator on cell edges to cell centers." - - def fget(self): - if(self._aveE2CC is None): - # The number of cell centers in each direction - n = self.n - if(self.dim == 1): - raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') - elif(self.dim == 2): - self._aveE2CC = sp.hstack((sp.kron(av(n[1]), speye(n[0])), - sp.kron(speye(n[1]), av(n[0]))), format="csr") - elif(self.dim == 3): - self._aveE2CC = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveE2CC - return locals() - _aveE2CC = None - aveE2CC = property(**aveE2CC()) - - def aveN2CC(): - doc = "Construct the averaging operator on cell nodes to cell centers." - - def fget(self): - if(self._aveN2CC is None): - # The number of cell centers in each direction - n = self.n - if(self.dim == 1): - self._aveN2CC = av(n[0]) - elif(self.dim == 2): - self._aveN2CC = sp.hstack((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") - elif(self.dim == 3): - self._aveN2CC = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveN2CC - return locals() - _aveN2CC = None - aveN2CC = property(**aveN2CC()) - - def aveN2CCv(): - doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate." - - def fget(self): - if(self._aveN2CCv is None): - # The number of cell centers in each direction - n = self.n - if(self.dim == 1): - self._aveN2CCv = av(n[0]) - elif(self.dim == 2): - self._aveN2CCv = sp.block_diag((sp.kron(av(n[1]), av(n[0])), - sp.kron(av(n[1]), av(n[0]))), format="csr") - elif(self.dim == 3): - self._aveN2CCv = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0])), - kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr") - return self._aveN2CCv - return locals() - _aveN2CCv = None - aveN2CCv = property(**aveN2CCv()) - - def getMass(self, materialProp=None, loc='e'): - """ Produces mass matricies. - - :param str loc: Average to location: 'e'-edges, 'f'-faces - :param None,float,numpy.ndarray materialProp: property to be averaged (see below) - :rtype: scipy.sparse.csr.csr_matrix - :return: M, the mass matrix - - materialProp can be:: - - None -> takes materialProp = 1 (default) - float -> a constant value for entire domain - numpy.ndarray -> if materialProp.size == self.nC - 3D property model - if materialProp.size = self.nCz - 1D (layered eath) property model - """ - if materialProp is None: - materialProp = np.ones(self.nC) - elif type(materialProp) is float: - materialProp = np.ones(self.nC)*materialProp - elif materialProp.shape == (self.nCz,): - materialProp = materialProp.repeat(self.nCx*self.nCy) - materialProp = mkvc(materialProp) - assert materialProp.shape == (self.nC,), "materialProp incorrect shape" - - if loc=='e': - Av = self.aveE2CC - elif loc=='f': - Av = self.aveF2CC - else: - raise ValueError('Invalid loc') - - diag = Av.T * (self.vol * mkvc(materialProp)) - - return sdiag(diag) - - def getEdgeMass(self, materialProp=None): - """mass matrix for products of edge functions w'*M(materialProp)*e""" - return self.getMass(loc='e', materialProp=materialProp) - - def getFaceMass(self, materialProp=None): - """mass matrix for products of face functions w'*M(materialProp)*f""" - return self.getMass(loc='f', materialProp=materialProp) - - def getFaceMassDeriv(self): - Av = self.aveF2CC - return Av.T * sdiag(self.vol) diff --git a/SimPEG/TensorView.py b/SimPEG/TensorView.py deleted file mode 100644 index 687a520d..00000000 --- a/SimPEG/TensorView.py +++ /dev/null @@ -1,335 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import matplotlib -from mpl_toolkits.mplot3d import Axes3D -from SimPEG.utils import mkvc - - -class TensorView(object): - """ - Provides viewing functions for TensorMesh - - This class is inherited by TensorMesh - """ - def __init__(self): - pass - - def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False): - """ - Mesh.plotImage(I) - - Plots scalar fields on the given mesh. - - Input: - - :param numpy.array I: scalar field - - Optional Input: - - :param str imageType: type of image ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez') or combinations, e.g. ExEy or FxFz - :param int figNum: number of figure to plot to - :param matplotlib.axes.Axes ax: axis to plot to - :param str direction: slice dimensions, 3D only ('x', 'y', 'z') - :param bool numbering: show numbering of slices, 3D only - :param str annotationColor: color of annotation, e.g. 'w', 'k', 'b' - :param bool showIt: call plt.show() - - .. plot:: examples/mesh/plot_image_2D.py - :include-source: - - .. plot:: examples/mesh/plot_image_3D.py - :include-source: - """ - assert type(I) == np.ndarray, "I must be a numpy array" - assert type(numbering) == bool, "numbering must be a bool" - assert direction in ["x", "y","z"], "direction must be either x,y, or z" - - - if imageType == 'CC': - assert I.size == self.nC, "Incorrect dimensions for CC." - elif imageType == 'N': - assert I.size == self.nN, "Incorrect dimensions for N." - elif imageType == 'Fx': - if I.size != np.prod(self.nFx): I, fy, fz = self.r(I,'F','F','M') - elif imageType == 'Fy': - if I.size != np.prod(self.nFy): fx, I, fz = self.r(I,'F','F','M') - elif imageType == 'Fz': - if I.size != np.prod(self.nFz): fx, fy, I = self.r(I,'F','F','M') - elif imageType == 'Ex': - if I.size != np.prod(self.nEx): I, ey, ez = self.r(I,'E','E','M') - elif imageType == 'Ey': - if I.size != np.prod(self.nEy): ex, I, ez = self.r(I,'E','E','M') - elif imageType == 'Ez': - if I.size != np.prod(self.nEz): ex, ey, I = self.r(I,'E','E','M') - elif imageType[0] == 'E': - plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} - fig = plt.figure(figNum) - # Determine the subplot number: 131, 121 - numPlots = 130 if plotAll else len(imageType)/2*10+100 - pltNum = 1 - ex, ey, ez = self.r(I,'E','E','M') - if plotAll or 'Ex' in imageType: - ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(ex, imageType='Ex', ax=ax_x, **options) - pltNum +=1 - if plotAll or 'Ey' in imageType: - ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(ey, imageType='Ey', ax=ax_y, **options) - pltNum +=1 - if plotAll or 'Ez' in imageType: - ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(ez, imageType='Ez', ax=ax_z, **options) - pltNum +=1 - return - elif imageType[0] == 'F': - plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} - fig = plt.figure(figNum) - # Determine the subplot number: 131, 121 - numPlots = 130 if plotAll else len(imageType)/2*10+100 - pltNum = 1 - fx, fy, fz = self.r(I,'F','F','M') - if plotAll or 'Fx' in imageType: - ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(fx, imageType='Fx', ax=ax_x, **options) - pltNum +=1 - if plotAll or 'Fy' in imageType: - ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(fy, imageType='Fy', ax=ax_y, **options) - pltNum +=1 - if plotAll or 'Fz' in imageType: - ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(fz, imageType='Fz', ax=ax_z, **options) - pltNum +=1 - return - else: - raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'") - - - if ax is None: - fig = plt.figure(figNum) - fig.clf() - ax = plt.subplot(111) - else: - assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!" - fig = ax.figure - - if self.dim == 1: - if imageType == 'CC': - ph = ax.plot(self.vectorCCx, I, '-ro') - elif imageType == 'N': - ph = ax.plot(self.vectorNx, I, '-bs') - ax.set_xlabel("x") - ax.axis('tight') - elif self.dim == 2: - if imageType == 'CC': - C = I[:].reshape(self.n, order='F') - elif imageType == 'N': - C = I[:].reshape(self.n+1, order='F') - C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:]) - elif imageType == 'Fx': - C = I[:].reshape(self.nFx, order='F') - C = 0.5*(C[:-1, :] + C[1:, :] ) - elif imageType == 'Fy': - C = I[:].reshape(self.nFy, order='F') - C = 0.5*(C[:, :-1] + C[:, 1:] ) - elif imageType == 'Ex': - C = I[:].reshape(self.nEx, order='F') - C = 0.5*(C[:,:-1] + C[:,1:] ) - elif imageType == 'Ey': - C = I[:].reshape(self.nEy, order='F') - C = 0.5*(C[:-1,:] + C[1:,:] ) - - ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T) - ax.axis('tight') - ax.set_xlabel("x") - ax.set_ylabel("y") - - elif self.dim == 3: - if direction == 'z': - - # get copy of image and average to cell-centres is necessary - if imageType == 'CC': - Ic = I[:].reshape(self.n, order='F') - elif imageType == 'N': - Ic = I[:].reshape(self.n+1, order='F') - Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] ) - elif imageType == 'Fx': - Ic = I[:].reshape(self.nFx, order='F') - Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:]) - elif imageType == 'Fy': - Ic = I[:].reshape(self.nFy, order='F') - Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:]) - elif imageType == 'Fz': - Ic = I[:].reshape(self.nFz, order='F') - Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:]) - elif imageType == 'Ex': - Ic = I[:].reshape(self.nEx, order='F') - Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1]) - elif imageType == 'Ey': - Ic = I[:].reshape(self.nEy, order='F') - Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1]) - elif imageType == 'Ez': - Ic = I[:].reshape(self.nEz, order='F') - Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:]) - - # determine number oE slices in x and y dimension - nX = np.ceil(np.sqrt(self.nCz)) - nY = np.ceil(self.nCz/nX) - - # allocate space for montage - nCx = self.nCx - nCy = self.nCy - - C = np.zeros((nX*nCx,nY*nCy)) - - for iy in range(int(nY)): - for ix in range(int(nX)): - iz = ix + iy*nX - if iz < self.nCz: - C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = Ic[:, :, iz] - else: - C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan - - C = np.ma.masked_where(np.isnan(C), C) - xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())] - yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())] - # Plot the mesh - ph = ax.pcolormesh(xx, yy, C.T) - # Plot the lines - gx = np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0]) - gy = np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1]) - # Repeat and seperate with NaN - gxX = np.c_[gx, gx, gx+np.nan].ravel() - gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel() - gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel() - gyY = np.c_[gy, gy, gy+np.nan].ravel() - ax.plot(gxX, gxY, annotationColor+'-', linewidth=2) - ax.plot(gyX, gyY, annotationColor+'-', linewidth=2) - ax.axis('tight') - - if numbering: - pad = np.sum(self.hx)*0.04 - for iy in range(int(nY)): - for ix in range(int(nX)): - iz = ix + iy*nX - if iz < self.nCz: - ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad, - '#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large') - - ax.set_title(imageType) - if showIt: plt.show() - return ph - - def plotGrid(self, 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 - :param bool faces: plot faces - :param bool centers: plot centers - :param bool edges: plot edges - :param bool lines: plot lines connecting nodes - :param bool showIt: call plt.show() - - .. plot:: examples/mesh/plot_grid_2D.py - :include-source: - - .. plot:: examples/mesh/plot_grid_3D.py - :include-source: - """ - 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') - if faces: - ax.plot(xs1[:, 0], xs1[:, 1], 'g>') - ax.plot(xs2[:, 0], xs2[:, 1], 'g^') - - # Plot the grid lines - if lines: - NN = self.r(self.gridN, 'N', 'N', 'M') - X1 = np.c_[mkvc(NN[0][0, :]), mkvc(NN[0][self.nCx, :]), mkvc(NN[0][0, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][0, :]), mkvc(NN[1][self.nCx, :]), mkvc(NN[1][0, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, 0]), mkvc(NN[0][:, self.nCy]), mkvc(NN[0][:, 0])*np.nan].flatten() - 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.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 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]) - 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]) - - # Plot the grid lines - if lines: - NN = self.r(self.gridN, 'N', 'N', 'M') - X1 = np.c_[mkvc(NN[0][0, :, :]), mkvc(NN[0][self.nCx, :, :]), mkvc(NN[0][0, :, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][0, :, :]), mkvc(NN[1][self.nCx, :, :]), mkvc(NN[1][0, :, :])*np.nan].flatten() - Z1 = np.c_[mkvc(NN[2][0, :, :]), mkvc(NN[2][self.nCx, :, :]), mkvc(NN[2][0, :, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, 0, :]), mkvc(NN[0][:, self.nCy, :]), mkvc(NN[0][:, 0, :])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, 0, :]), mkvc(NN[1][:, self.nCy, :]), mkvc(NN[1][:, 0, :])*np.nan].flatten() - Z2 = np.c_[mkvc(NN[2][:, 0, :]), mkvc(NN[2][:, self.nCy, :]), mkvc(NN[2][:, 0, :])*np.nan].flatten() - X3 = np.c_[mkvc(NN[0][:, :, 0]), mkvc(NN[0][:, :, self.nCz]), mkvc(NN[0][:, :, 0])*np.nan].flatten() - Y3 = np.c_[mkvc(NN[1][:, :, 0]), mkvc(NN[1][:, :, self.nCz]), mkvc(NN[1][:, :, 0])*np.nan].flatten() - Z3 = np.c_[mkvc(NN[2][:, :, 0]), mkvc(NN[2][:, :, self.nCz]), mkvc(NN[2][:, :, 0])*np.nan].flatten() - 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.set_xlabel('x1') - ax.set_ylabel('x2') - ax.set_zlabel('x3') - if showIt: plt.show() diff --git a/docs/examples/mesh/plot_TensorMesh.py b/docs/examples/mesh/plot_TensorMesh.py index 94c627f8..e773ed13 100644 --- a/docs/examples/mesh/plot_TensorMesh.py +++ b/docs/examples/mesh/plot_TensorMesh.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh pad = 7 padfactor = 1.4 diff --git a/docs/examples/mesh/plot_grid_2D.py b/docs/examples/mesh/plot_grid_2D.py index 258fbdfc..4c6d5da5 100644 --- a/docs/examples/mesh/plot_grid_2D.py +++ b/docs/examples/mesh/plot_grid_2D.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh h1 = np.linspace(.1,.5,3) h2 = np.linspace(.1,.5,5) diff --git a/docs/examples/mesh/plot_grid_3D.py b/docs/examples/mesh/plot_grid_3D.py index caf8dbd2..6278802c 100644 --- a/docs/examples/mesh/plot_grid_3D.py +++ b/docs/examples/mesh/plot_grid_3D.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh h1 = np.linspace(.1,.5,3) h2 = np.linspace(.1,.5,5) diff --git a/docs/examples/mesh/plot_image_2D.py b/docs/examples/mesh/plot_image_2D.py index c7bcb0d5..9fef3f1b 100644 --- a/docs/examples/mesh/plot_image_2D.py +++ b/docs/examples/mesh/plot_image_2D.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh n = 20 h = np.ones(n)/n diff --git a/docs/examples/mesh/plot_image_3D.py b/docs/examples/mesh/plot_image_3D.py index d00756a1..c67b0d84 100644 --- a/docs/examples/mesh/plot_image_3D.py +++ b/docs/examples/mesh/plot_image_3D.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh n = 20 h = np.ones(n)/n diff --git a/docs/examples/mesh/plot_nodes.py b/docs/examples/mesh/plot_nodes.py index 5d2c912b..ef11b351 100644 --- a/docs/examples/mesh/plot_nodes.py +++ b/docs/examples/mesh/plot_nodes.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from SimPEG import TensorMesh +from SimPEG.mesh import TensorMesh x0 = np.zeros(2) h1 = np.linspace(.1,.5,3)